On Wed, Oct 27, 2010 at 10:50 PM, Nicolai Heitz <nicolaihe...@gmx.de> wrote: > m 27.10.2010 02:02, schrieb Sebastian Walter: > >> On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen<p...@iki.fi> wrote: >>> Tue, 26 Oct 2010 14:24:39 -0700, Nicolai Heitz wrote: >>>>> http://mail.scipy.org/mailman/listinfo/scipy-user >>>> I contacted them already but they didn't responded so far and I was >>>> forwarded to that list which was supposed to be more appropriated. >>> I think you are thinking here about some other list -- scipy-user >>> is the correct place for this discussion (and I don't remember seeing >>> your mail there). > I was pretty sure that I put it there. Unfortunately it had a different > name there: Errors in calculation of the Hessian using numdifftools. But > I can't find the post by myself at the moment so maybe something went > wrong. >>> [clip] >>>> 1) Can I make it run/fix it, so that it is also going to work for the SI >>>> scaling? >>> Based on a brief look, it seems that uniform scaling will not help you, >>> as you have two very different length scales in the problem, >>> >>> 1/sqrt(m w^2)>> C >>> >>> If you go to CM+relative coordinates you might be able to scale them >>> separately, but that's fiddly and might not work for larger N. >>> >>> In your problem, things go wrong when the ratio between the >>> length scales approaches 1e-15 which happens to be the machine epsilon. >>> This implies that the algorithm runs into some problems caused by the >>> finite precision of floating-point numbers. >>> >>> What exactly goes wrong and how to fix it, no idea --- I didn't look into >>> how Numdifftools is implemented. > > Probably you are right. I converted it to natural units and it worked > out in the 2 ion case. Increasing the number of ions leads to problems > again and I have no idea where those problems come from. > >>>> 2) How can I be sure that increasing the number of ions or adding a >>>> somehow more complicated term to the potential energy is not causing the >>>> same problems even in natural units? >>>> >>>> 3) In which range is numdifftools working properly. >>> That depends on the algorithm and the problem. Personally, I wouldn't >>> trust numerical differentiation if the problem has significantly >>> different length scales, it is important to capture all of them >>> accurately, and it is not clear how to scale them to the same size. >>> Writing ND software that works as expected all the time is probably >>> not possible even in theory. >>> >>> Numerical differentiation is not the only game in the town. I'd look >>> into automatic differentiation (AD) -- there are libraries available >>> for Python also for that, and it is numerically stable. >>> >>> E.g. >>> >>> http://en.wikipedia.org/wiki/Automatic_differentiation#Software >>> >>> has a list of Python libraries. I don't know which of them would be >>> the best ones, though. >>> >> they all have their pro's and con's. >> Being (co-)author of some of these tools, my personal and very biased >> advice is: >> if you are on Linux, I would go for PYADOLC. it provides bindings to a >> feature-rich and well-tested C++ library. >> However, the installation is a little tricker than a "setup.py build" >> since you will need to compile ADOL-C and get Boost::Python to work. >> PYADOLC can also differentiate much more complicated code than your >> example in a relatively efficient manner. > > Is there by chance any possibility to make PYADOLC run on a (lame) > windows xp engine. If not what else would u recommend (besides switching > to Linux, what I am going to do soon).
1) PYADOLC A windows version has been requested several times now. But until recently ADOL-C wasn't available as windows version. So yes, in principle it should be possible to get it to work on windows: You will need 1) boost:python http://www.boost.org/doc/libs/1_44_0/libs/python/doc/index.html 2) ADOL-C sources http://www.coin-or.org/projects/ADOL-C.xml 3) scons http://www.scons.org/ on windows. If you want to give it a try I could help to get it to work. 2) Alternatives: You can also try the ALGOPY which is pure Python and is known to work on Linux and Windows. The installation is also very easy (setup.py build or setup.py install) I have added your problem to the ALGOPY documentation: http://packages.python.org/algopy/examples/hessian_of_potential_function.html The catch is that ALGOPY is not as mature as PYADOLC. However, if you are careful to write clean code it should work reliably. Sebastian >> For your example the code looks like: >> >> --------------------- code --------------------------- >> >> .... >> >> c=classicalHamiltonian() >> xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10) >> >> import adolc; import numpy >> >> # trace the computation >> adolc.trace_on(0) >> x = adolc.adouble(c.initialposition()) >> adolc.independent(x) >> y = c.potential(x) >> adolc.dependent(y) >> adolc.trace_off() >> >> hessian = adolc.hessian(0, xopt) >> eigenvalues = numpy.linalg.eigh(hessian)[0] >> normal_modes = c.normal_modes(eigenvalues) >> print 'hessian=\n',hessian >> print 'eigenvalues=\n',eigenvalues >> print 'normal_modes=\n',normal_modes >> --------------------- code --------------------------- >> >> and you get as output >> Optimization terminated successfully. >> Current function value: 0.000000 >> Iterations: 81 >> Function evaluations: 153 >> hessian= >> [[ 5.23748399e-12 -2.61873843e-12] >> [ -2.61873843e-12 5.23748399e-12]] >> eigenvalues= >> [ 2.61874556e-12 7.85622242e-12] >> normal_modes= >> [ 6283185.30717959 10882786.30440101] >> >> >> Also, you should use an eigenvalue solver for symmetric matrices, e.g. >> numpy.linalg.eigh. >> > Your code example looks awesome and leads to the correct results. Thank > you very much. I try to make it work on my pc as well. >> regards, >> Sebastian >> >>> -- >>> Pauli Virtanen >>> >>> _______________________________________________ >>> NumPy-Discussion mailing list >>> numpy-discuss...@scipy.org >>> http://mail.scipy.org/mailman/listinfo/numpy-discussion >>> >> _______________________________________________ >> NumPy-Discussion mailing list >> numpy-discuss...@scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
