Browse files

improved known_bugs_and_pitfalls according to the new information tha…

…t i got from the nump

mailing list
started finite differences module (for testing)
  • Loading branch information...
1 parent 5081a59 commit 4cd7a50f48686e9f7c3286681952f6f90e3c03d3 @b45ch1 committed Jan 19, 2010
Showing with 117 additions and 47 deletions.
  1. +8 −47 KNOWN_BUGS_AND_PITFALLS.rst
  2. +58 −0 SConstruct.old
  3. +19 −0 adolc/finitedifferences.py
  4. +32 −0 adolc/tests/test_finitedifferences.py
View
55 KNOWN_BUGS_AND_PITFALLS.rst
@@ -33,51 +33,12 @@ Pitfalls
In [9]: y
Out[9]: array([ 4., 5., 6.])
-
- This is not a bug of PYADOLC but the way `*=` is implemented in numpy. One can see similar behaviour for the dtype int and float::
-
- In [1]: import numpy
-
- In [2]: x = numpy.ones(2,dtype=int)
-
- In [3]: y = 1.3 * numpy.ones(2,dtype=float)
-
- In [4]: z = x * y
-
- In [5]: z
- Out[5]: array([ 1.3, 1.3])
-
- In [6]: x *= y
-
- In [7]: x
- Out[7]: array([1, 1])
-
- In [8]: x.dtype
- Out[8]: dtype('int32')
-
- that means that the inplace operation `x *= y ` is *not* the same as `x = x * y`.
- This also means that::
-
- x = adouble(numpy.array([1,2,3],dtype=float))
- It is inconsistent to the Python behaviour and therefore a little surprising::
-
- In [9]: a = 1
-
- In [10]: b = 1.3
-
- In [11]: c = a * b
-
- In [12]: c
- Out[12]: 1.3
-
- In [13]: a *= b
-
- In [14]: a
- Out[14]: 1.3
-
- This is intended behaviour of numpy, but it leads to incorrect computations since no exception or
- warning is raised by numpy. For more info see
- http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg03236.html
-
-
+ That means that the inplace operation `x *= y ` is *not* the same as `x = x * y`.
+
+ This is not a bug of PYADOLC but a bug in numpy's implementation of the augmented
+ assignment statemetns `*=`, etc. for arrays of objects.
+
+ Numpy tries to cast the dtype of `y` to the dtype `x`. If x has dtype `float` then on each element
+ y[i].__float__() is called.
+
View
58 SConstruct.old
@@ -0,0 +1,58 @@
+import distutils.sysconfig
+import os
+import numpy
+
+
+# -1: CUSTOMIZE THIS TO FIT YOUR SYSTEM !!!
+python_include_path= '/usr/include/python2.5'
+
+adolc_include_path = '/u/walter/workspace/ADOL-C-2.1.0/ADOL-C/src'
+adolc_library_path = '/u/walter/workspace/ADOL-C-2.1.0/ADOL-C/src/.libs'
+
+colpack_include_path = '/u/walter/workspace/colpack/build/include'
+colpack_library_path = '/u/walter/workspace/colpack/build/lib'
+
+LIBS = ['adolc',
+ 'boost_python',
+ 'colpack',
+ ]
+LIBPATH = [
+ adolc_library_path,
+ colpack_library_path,
+ ]
+INCLUDEPATH = [
+ adolc_include_path,
+ python_include_path,
+ colpack_include_path,
+ ]
+
+
+
+# 0: setup the command line parsing
+AddOption('--prefix',
+ dest='prefix',
+ nargs=1, type='string',
+ action='store',
+ metavar='DIR',
+ help='installation prefix')
+
+
+env = Environment(
+ PREFIX = GetOption('prefix'),
+ TMPBUILD = '/tmp/builddir',
+ CPPPATH=[distutils.sysconfig.get_python_inc(),numpy.get_include()] + INCLUDEPATH,
+ CXXFLAGS="-ftemplate-depth-100 -DBOOST_PYTHON_DYNAMIC_LIB -O2 -Wall",
+ LIBPATH=["/usr/lib/python2.5/config"] + LIBPATH,
+ LIBS= LIBS,
+ RPATH = LIBPATH, #include information where shared libraries can be found to avoid errors like: "ImportError: libboost_python-gcc42-mt-1_34_1.so.1.34.1: cannot open shared object file: No such file or directory"
+ SHLIBPREFIX="", #gets rid of lib prefix
+)
+
+Export('env')
+Export('adolc_include_path')
+SConscript('adolc/SConscript')
+SConscript('adolc/sparse/SConscript')
+
+env.Install( target='./build/adolc/', source = ['adolc/__init__.py','adolc/wrapped_functions.py', 'adolc/cgraph.py','adolc/_adolc.so'])
+env.Install( target='./build/adolc/sparse/', source = ['adolc/sparse/__init__.py', 'adolc/sparse/_colpack.so', 'adolc/sparse/wrapped_functions.py'])
+
View
19 adolc/finitedifferences.py
@@ -0,0 +1,19 @@
+
+import numpy
+
+
+def jacobian(ffcn, x, epsilon = 10**-8):
+ """
+ Computes the Jacobian of the function ffcn by finite differences in point x.
+ """
+
+ N = numpy.size(x)
+ xv = numpy.zeros((N,N))
+ V = epsilon * numpy.eye(N)
+
+ for n in range(N):
+ xv[:,n] = x
+
+ J = (ffcn(xv + V) - ffcn(xv))/epsilon
+ return J
+
View
32 adolc/tests/test_finitedifferences.py
@@ -0,0 +1,32 @@
+from numpy.testing import *
+import numpy
+import numpy.random
+
+from adolc.finitedifferences import *
+
+
+
+
+class FiniteDifferencesTest(TestCase):
+ def test_jacobian(self):
+ N = 10
+ A = numpy.random.rand(N)
+
+ def f(x):
+ return numpy.dot(A, x)
+
+ x = numpy.random.rand(N)
+
+ assert_array_almost_equal(jacobian(f,x),A, decimal=7)
+
+
+
+
+
+
+if __name__ == '__main__':
+ try:
+ import nose
+ except:
+ print 'Please install nose for unit testing'
+ nose.runmodule()

0 comments on commit 4cd7a50

Please sign in to comment.