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.walter@gmail.com> wrote:
On Wed, Oct 27, 2010 at 10:50 PM, Nicolai Heitz <nicolaiheitz@gmx.de> wrote:
m 27.10.2010 02:02, schrieb Sebastian Walter:
On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen<pav@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.htm... 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-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
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion