hmm, I have just realized that I forgot to upload the new version to pypi: it is now available on http://pypi.python.org/pypi/algopy
On Thu, Oct 28, 2010 at 10:47 AM, Sebastian Walter <sebastian.wal...@gmail.com> wrote: > 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
