prathmeshrmadhu / numdifftools Goto Github PK

Dear Per,

I have filed a scipy ticket suggesting numerical differentiation routines,
possibly from numdifftools, should be included in scipy:
http://projects.scipy.org/scipy/ticket/1510

Maybe you can participate in the discussion on that ticket or even make it 
happen?

Cheers, Christoph

The tests fail for me, both on 0.6.0 from pypi and trunk, with both python 2.7 
and 3.4.

Output is attached from test_all.py on trunk with Python 2.7.8, numpy 1.8.1, 
scipy 0.14.0.

I have had a couple of people ask me about this bug because my package 
exercises that part of numdifftools (see last couple of messages 
http://groups.google.com/group/pymc/browse_thread/thread/3fc6c26701cd7943). The 
bug is fixed in the source. Any chance a new version could be released that 
includes this fix?

More of an open-ended question than an issue...I have an application where I 
need to calculate the Hessian of a multivariate function many, many times. It's 
a highly non-linear, non-convex function, with many local extrema, and so the 
best step size is not immediately obvious. Numdifftools seems eminently suited 
for this, but there is a comment in the code that the Hessian function:

"is NOT a tool for frequent use on an expensive to evaluate objective function, 
especially in a large number of dimensions. Its computation will use roughly  
O(6*n^2) function evaluations for n parameters."

What would you recommend if you DO have an expensive objective function with a 
large (~10) number of dimension. I'm considering modifying the code to pick up 
some additional efficiency, but I wanted to see if you, as the author, had any 
input first. 

For example, one idea is to determine the "best" step size for some arbitrary 
set of inputs, and then fix the step size for the remainder of the Hessian 
evaluations. Even if I have upwards of 1% error in the calculation, I just want 
to get the step size on the right order of magnitude. That way I don't need to 
iteratively determine the best step size for each and every Hessian evaluation. 

Thanks for your thoughts.

What steps will reproduce the problem?
1. Run this script : http://pastebin.com/Qutr3Gmf

What is the expected output? What do you see instead?

scipy.misc.derivative works well on this example, but not numdifftools 

What version of the product are you using? On what operating system?

a fresh checkout from svn trunk 

# $ svnversion
42M


Please provide any additional information below.

BTW, could you please also 

# run 2to3 -w numdifftools/core.py

and commit the result, so that numdifftools will work with python3

 trying to easy_install Numdifftools, I get the following error:


john-salvatiers-macbook:~ johnsalvatier$ easy_install Numdifftools
Searching for Numdifftools
Reading http://pypi.python.org/simple/Numdifftools/
Best match: Numdifftools 0.2.1
Downloading 
http://pypi.python.org/packages/source/N/Numdifftools/Numdifftools-
0.2.1.zip#md5=42a9f9756e4fbc56de05424f2bba0b98
Processing Numdifftools-0.2.1.zip
Running Numdifftools-0.2.1/setup.py -q bdist_egg --dist-dir 
/var/folders/B8/B85hpAz-GTGNEmYQ9zGZwE+++TI/-Tmp-/easy_install-
lXQU4f/Numdifftools-0.2.1/egg-dist-tmp-Tu9Rgs
Traceback (most recent call last):
  File 
"/Library/Frameworks/Python.framework/Versions/Current/bin/easy_install", 
line 8, in <module>
    load_entry_point('setuptools==0.6c9', 'console_scripts', 
'easy_install')()
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 1671, in main
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 1659, in with_ei_usage
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 1675, in <lambda>
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/distutils/
core.py", line 151, in setup
    dist.run_commands()
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/distutils/
dist.py", line 974, in run_commands
    self.run_command(cmd)
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/distutils/
dist.py", line 994, in run_command
    cmd_obj.run()
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 211, in run
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 446, in easy_install
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 476, in install_item
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 655, in install_eggs
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 930, in build_and_install
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/command/easy_install.py", 
line 919, in run_setup
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/sandbox.py", line 27, in 
run_setup
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/sandbox.py", line 63, in run
  File 
"/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/setuptools-0.6c9-py2.5.egg/setuptools/sandbox.py", line 29, in 
<lambda>
  File "setup.py", line 21, in <module>
ImportError: No module named numdifftools

I already discussed the problem in the numpy and scipy forum and they couldn't 
help. We assume that numdifftools is not working properly in the scaling of SI 
units as the following example shows.

What is the expected output? What do you see instead?
see describtion below

What version of the product are you using? On what operating system?
numdifftools 0.3.1 on windows XP

Please provide any additional information below.

I am using Python 2.6 and the following packages:

1) numpy
2) scipy
3) numdifftools

I want to use
Python to calculate the eigenvalues of the Hamiltonian given in the
code. Therefore I have to do the following steps:

1) define the Hamiltonian or the potential respectively (-->  function:
potential)
2) minimize the potential ( I am using scipy.optimize.fmin to calculate
the minimum of the potential)
(2b) Expand the Hamiltonian around the minimum position. This step is
not done in the code because it is only necessary for the purpose that
one understand the physical background and why one have to do step 3)
3)   Calculate the Hessian matrix of the potential, that means calculate
the second derivatives of the potential at the point of the minimum
position (-->  numdifftools.Hessian)
4) Calculate the eigenvalues of the Hessian (-->scipy.linalg.eigvals)

The problem can be solved analytically except of the calculation of the
minima in step 2:

Now I have the following problem:
The Hessian matrix evaluated at the minimum position is wrong.

What I checked so far:

1) The potential seems to be calculated correctly and the minimum
position for two ions (N=2) is fine.
2) The Hessian matrix calculation works fine for several other functions.
3) The Hessian matrix is correct when I set the Coulomb interaction to
zero (C=0)
4) The Hessian matrix is correct when I set C to some simple arbitrary
integer numbers like 1, 5 ,8 ....
5) The Hesse matrix gives the correct number of the first part of the
diagonal elements even with C=2.3e-28 ! But in that case the offdiagonal
elements and the part of the diagonal element which refers to the second
derivative in the Coulomb term is wrong! The offdiagonal term should be
C/(abs(x_1-x_2))**3 = 2.6e-12 but it is 3.4e-24!
5) In principal Python can divide those numbers (2*2.3e-28)/(5.6e-6)**3
6) I played around with the accurateness of the calculation of the
Hessian by varying the arguments numTerms, method and metOrder but that
didn't change anything in my case.

My source code is attached. 

# import a tool to use / as a symbol for normal division
from __future__ import division

#import system data
import sys, os

#Loading the required packages
import scipy as sp
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import numdifftools as nd

# and subpackages
from scipy import *
from scipy import linalg, optimize, constants

#-------------------------------------------------------------------------------
----------
#       Hamiltonian      H=sum_i(p_i2/(2m)+ 1/2 * m * w2 x_i2)+ 
sum_(i!=j)(a/|x_i-x_j|)
#-------------------------------------------------------------------------------
----------

class classicalHamiltonian:
    def __init__(self):

        self.N = 2                            #N is a scalar, it's the number of ions in the chain
        f = 1000000                            #f is a scalar, it's the trap frequency
        self.w = 2*pi*f                         #w is a scalar, it's the angular velocity corresponding to the trap frequency
        self.C = (4*pi*constants.epsilon_0)**(-1)*constants.e**2    #C is a scalar, it's the Coulomb constant times the electronic charge in SI
        self.m = 39.96*1.66e-27                        #m is the mass of a single trapped ion in the chain




    def potential(self, positionvector):                     #Defines the potential that is going to be minimized

        x= positionvector                         #x is an 1-d array (vector) of lenght N that contains the positions of the N ions
        w= self.w
        C= self.C
        m= self.m


        #First we consider the potential of the harmonic osszilator
        Vx = 0.5 * m * (w**2) * sum(x**2)               

        for i in range(len(x)):
            for j in range(i+1, len(x)):
                Vx += C * (abs(x[i]-x[j]))**(-1)    #then we add the coulomb interaction

        return Vx

    def initialposition(self):        #Defines the initial position as an estimate for the minimize process

        N= self.N
        x_0 = r_[-(N-1)/2:(N-1)/2:N*1j]
        return x_0

    def normal_modes(self, eigenvalues):    #the computed eigenvalues of the matrix Vx are of the form (normal_modes)2*m.
        m = self.m
        normal_modes = sqrt(eigenvalues/m)
        return normal_modes

#C=(4*pi*constants.epsilon_0)**(-1)*constants.e**2
c=classicalHamiltonian()
#print c.potential(array([-0.5, 0.5]))
xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10)
hessian = nd.Hessian(c.potential)
eigenvalues = linalg.eigvals(hessian(xopt))
normal_modes = c.normal_modes(eigenvalues)

What steps will reproduce the problem?
Installing numdifftools from source.

I have lates svn-version of scipy installed (scipy 0.10.0.dev). Installing 
numdifftools from source causes the setup to download a scipy version 
scipy>=0.8. The installed version obviously requires this condition (0.10>0.8), 
but I guess it is detected as 0.1 instead to 0.10.

What steps will reproduce the problem?

Using python2.7 and numdifftools 0.6.0

import numpy as np
import numdifftools as nd

def(x):
    return np.sum(x**2)

hess=nd.Hessian(f)
N=3
x0=np.zeros(3)

print hess(x0)

will raise 

Numdifftools-0.6.0-py2.7.egg/numdifftools/core.py:295: RuntimeWarning: divide 
by zero encountered in log2
  stepNom = np.maximum(np.log2(np.abs(x0)), 0.02)
[[ 2.  0.  0.]
 [ 0.  2.  0.]
 [ 0.  0.  2.]]


What is the expected output? What do you see instead?

Expected output is just the ndarray without the error message

On what operating system?

Linux version 3.2.0-4-amd64 ([email protected]) (gcc version 4.6.3 
(Debian 4.6.3-14) ) #1 SMP Debian 3.2.63-2+deb7u2

Often times other information must be passed to a function beyond the 
evaluation point. Optimization functions typically allow for an args=() option 
in the function call, which passes all variables in the tuple as additional 
arguments to the function being optimized. Can the same feature be built into 
numdifftools (I'm particularly interested in the Hessian feature).

Could you please make numdifftools Python 3 compatible?

Maybe all that is needed is fix the print and import statements?

$ python3.2 setup.py install --user
...
byte-compiling 
build/bdist.macosx-10.8-x86_64/egg/numdifftools/test/test_hessian.py to 
test_hessian.cpython-32.pyc
  File "build/bdist.macosx-10.8-x86_64/egg/numdifftools/test/test_hessian.py", line 69
    print c.potential(array([-0.5, 0.5]))
          ^
SyntaxError: invalid syntax


$ python3.2 -c 'import numdifftools as nd; nd.test(coverage=True)'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/deil/Library/Python/3.2/lib/python/site-packages/Numdifftools-0.4.0-py3.2.egg/numdifftools/__init__.py", line 1, in <module>
ImportError: No module named info

I also noticed that the numdifftools setup.py file declares that it works with 
Python 2.6 only:
http://code.google.com/p/numdifftools/source/browse/trunk/setup.py#62
But actually at the moment it works fine with Python 2.7 already.