Mailman 3 March 2000 - NumPy-Discussion

Pyfort 3.1 at Source Forge
by Paul F. Dubois 01 Mar '00

01 Mar '00

Pyfort, a Python-Fortran connection tool, is now available at pyfortran.sourceforge.net, instead of xfiles.llnl.gov. Documentation is available at the site. Version 3.1 uses distutils for easy extension building. A sample package is in the test subdirectory. Support for the pgf90 compiler on Linux is now available. Note that the limitation to non-explicitly interfaced routines has not been removed yet, but it will handle using pgf90-compiled F90 files in your Fortran.

1 0

Re: [Gmath-devel] Re: [Numpy-discussion] Derivatives
by Hassan Aurag 01 Mar '00

01 Mar '00

I have to agree with you on most counts that Numerical Recipes in C=20 is not a full-blown encyclopedia on all subtleties of doing numerical=20 computations. However it does its job greatly for a big class of stuff I am=20 interested in: minimization, non-linear system of equations solving=20 (The newton routine given there is good, accurate and fast.) There are errors and problems as in most other numerical books. In=20 truth, I don't think there is anything fully correct out there. When trying to make computations you have to do a lot of testing and=20 a lot of thought even if the algorithm seems perfect. That applies to=20 all books, recipes et al. I have just discovered that tan(pi/2) =3D 1.63317787284e+16 in=20 numerical python. And we all know this is bad. It should be infinity=20 period. We should define a symbol called infinity and put the correct=20 definition of tan, sin .... for all known angles then interpolate if=20 needed for the rest, or whatever is used to actually compute those thing= s. Peace <snip> - ------- End of Forwarded Message I believe you're actually talking about the math module, not the numeric module (I'm not aware of tan or pi definitions in numeric, but I haven't bothered to double check that). Never the less, I think it has relevance here as Numeric is all about doing serious number crunching. This problem is caused by the lack of infinite precision in python. Of course, how is it even possible to perform infinite precision on non-rational numbers? The obvious solution is to allow the routine (tan() in this case) to recognize named constants that have relevance in their domain (pi in this case). This would fix the problem: math.tan(math.pi) = -1.22460635382e-16 but it still doesn't solve your problem because the named constant would have the mathematical operation performed on it before it's passed into the function, ruining whatever intimate knowledge of the given named constant that routine has. Perhaps you could get the routine to recognize rational math on named constants (the problem with that solution is how do you not burden other routines with the knowledge of how to process that expression). Assuming you had that, even for your example, should the answer be positive or negative infinity? Another obvious solution is just to assume that any floating point overflow is infinity and any underflow is zero. This obviously won't work because some asymptotic functions (say 1/x^3) will overflow or underflow at values for x for which the correct answer is not properly infinity or zero. It is interesting to note that Matlab's behaviour is the same as Python's, which would indicate to me that there's not some easy solution to this problem that Guido et. al. overlooked. I haven't really researched the problem at all (although now I'm interested), but I'd be interested if anyone has a proposal for (or reference to) how this problem can be solved in a general purpose programming language at all (as there exists the distinct possibility that it can not be done in Python without breaking backwards compatibility). Terry -----BEGIN PGP SIGNATURE----- Version: PGPfreeware 5.0i for non-commercial use Charset: noconv iQA/AwUBOL2OUqfuGVwXgOQkEQJTpQCggOuFT2ZVavzMhy+jZgoehnrK5uIAoMzO D5OOdLtBvT97ee7vkckO+0Qt =SmqL -----END PGP SIGNATURE-----

1 0

Re: [Gmath-devel] Re: [Numpy-discussion] Derivatives
by Hassan Aurag 01 Mar '00

01 Mar '00

I have to agree with you on most counts that Numerical Recipes in C=20 is not a full-blown encyclopedia on all subtleties of doing numerical=20 computations. However it does its job greatly for a big class of stuff I am=20 interested in: minimization, non-linear system of equations solving=20 (The newton routine given there is good, accurate and fast.) There are errors and problems as in most other numerical books. In=20 truth, I don't think there is anything fully correct out there. When trying to make computations you have to do a lot of testing and=20 a lot of thought even if the algorithm seems perfect. That applies to=20 all books, recipes et al. I have just discovered that tan(pi/2) =3D 1.63317787284e+16 in=20 numerical python. And we all know this is bad. It should be infinity=20 period. We should define a symbol called infinity and put the correct=20 definition of tan, sin .... for all known angles then interpolate if=20 needed for the rest, or whatever is used to actually compute those thing= s. Peace <snip> - ------- End of Forwarded Message I believe you're actually talking about the math module, not the numeric module (I'm not aware of tan or pi definitions in numeric, but I haven't bothered to double check that). Never the less, I think it has relevance here as Numeric is all about doing serious number crunching. This problem is caused by the lack of infinite precision in python. Of course, how is it even possible to perform infinite precision on non-rational numbers? The obvious solution is to allow the routine (tan() in this case) to recognize named constants that have relevance in their domain (pi in this case). This would fix the problem: math.tan(math.pi) = -1.22460635382e-16 but it still doesn't solve your problem because the named constant would have the mathematical operation performed on it before it's passed into the function, ruining whatever intimate knowledge of the given named constant that routine has. Perhaps you could get the routine to recognize rational math on named constants (the problem with that solution is how do you not burden other routines with the knowledge of how to process that expression). Assuming you had that, even for your example, should the answer be positive or negative infinity? Another obvious solution is just to assume that any floating point overflow is infinity and any underflow is zero. This obviously won't work because some asymptotic functions (say 1/x^3) will overflow or underflow at values for x for which the correct answer is not properly infinity or zero. It is interesting to note that Matlab's behaviour is the same as Python's, which would indicate to me that there's not some easy solution to this problem that Guido et. al. overlooked. I haven't really researched the problem at all (although now I'm interested), but I'd be interested if anyone has a proposal for (or reference to) how this problem can be solved in a general purpose programming language at all (as there exists the distinct possibility that it can not be done in Python without breaking backwards compatibility). Terry -----BEGIN PGP SIGNATURE----- Version: PGPfreeware 5.0i for non-commercial use Charset: noconv iQA/AwUBOL2OUqfuGVwXgOQkEQJTpQCggOuFT2ZVavzMhy+jZgoehnrK5uIAoMzO D5OOdLtBvT97ee7vkckO+0Qt =SmqL -----END PGP SIGNATURE-----

1 0

RE: [Numpy-discussion] Derivatives
by Rodgers, Kevin 01 Mar '00

01 Mar '00

Whenever I see somebody implementing code based on "Numerical Recipies", I feel obligated to send them the following link: http://math.jpl.nasa.gov/nr/nr.html YMMV, as always . . . Kevin Rodgers Northrop Grumman Ryan Aeronautical Center krodgers(a)ryanaero.com > -----Original Message----- > From: Hassan Aurag [SMTP:aurag@crm.umontreal.ca] > Sent: Wednesday, March 01, 2000 6:33 AM > To: Gmath-devel(a)lists.sourceforge.net; > numpy-discussion(a)lists.sourceforge.net > Subject: [Numpy-discussion] Derivatives > > > > Hi, > > attached is a file called Derivative.py. > > It computes derivatives and is based on an algorithm found in > Numerical Recipes in C. > > What to do you think about it and has anyone started a "serious" > calculus oriented subpackage for Numerical Python in general? > > I mean: derivatives, partial derivatives, jacobian, hessian > implemented fast and precise. > > On another note, why isn't infinity defined in NumPy? > > Why is tan(pi/2) a number even if big? Shouldn't it be infinity? > > > > > > > _______________________________________________ > Numpy-discussion mailing list > Numpy-discussion(a)lists.sourceforge.net > http://lists.sourceforge.net/mailman/listinfo/numpy-discussion

2 1

Re: [Numpy-discussion] Derivatives
by Hassan Aurag 01 Mar '00

01 Mar '00

The answer is yes! It is an 1e10 and not 1e-10. At 0.0 you got to pick the interval yourself. You can't just use the starting point x > Date: Wed, 1 Mar 2000 16:06:28 +0100 (MET) > From: Fredrik Stenberg <su96-fst(a)nada.kth.se> > To: Hassan Aurag <aurag(a)CRM.UMontreal.CA> > Subject: Re: [Numpy-discussion] Derivatives > MIME-Version: 1.0 > > > Hi, > > > > attached is a file called Derivative.py. > > > > It computes derivatives and is based on an algorithm found in > > Numerical Recipes in C. > > > > What to do you think about it and has anyone started a "serious" > > calculus oriented subpackage for Numerical Python in general? > > > > I mean: derivatives, partial derivatives, jacobian, hessian > > implemented fast and precise. > > > > On another note, why isn't infinity defined in NumPy? > > > > Why is tan(pi/2) a number even if big? Shouldn't it be infinity? > > > > > > > > I tried your algoritm on sin(x) and i got some rather interesting > results. > > ######### EXAMPLE############# > from math import sin > > def f(x): > return sin(x) > > import Derivative > > print Derivative.Diff(f,0.0) > > ########RESULT################ > -2.03844228853e-10 > > It should be approx 1.0 > > > > I found the error i think.. > check row 28 in Derivative > h = random()/1e-10 > should that be h = random()/1e+10?? > > > Fredrik >

1 0

Derivative II
by Hassan Aurag 01 Mar '00

01 Mar '00

Forgot to attach file. Here it is #This is a file created by GMatH. ##! #I have also added this automatic import as a convenience. #Soon, you will be able to use a Druid to guide you through the import process! import sys from Numeric import zeros from RandomArray import random CON = 1.4 CON2 = CON*CON BIG = 1.0e30 SAFE = 2.0 MyError = "The choosen x is too small!" MyRealError = "The choosen x is not real!" def Diff(f, x): global CON, CON2, BIG, SAFE def func(g, t): return apply(g, (t,)) if type(x) == type(1.j): raise MyRealError, x if float(x) == 0.0: h = random()/1e-10 if float(x) != 0.0: h = float(x)/10 a = zeros((10,10),'d') if float(h) == 0.0: raise MyError, x try: test = func(f, x) test = func(f, x+h) test = func(f, x-h) except: print "Is function", str(f), "defined and continous around x =", x,"?" print sys.exc_type, sys.exc_value return None else : hh = h a[0][0] = (func(f, x + hh) - func(f, x - hh))/(2.0*hh) err = BIG ans = "CANNOT COMPUTE DERIVATIVE" for i in range(10): hh = hh/CON a[0,i] = (func(f, x + hh) - func(f, x - hh))/(2.0*hh) fac = CON2 for j in range(1, i+1): a[j][i]=(a[j-1][i]*fac - a[j-1][i-1])/(fac - 1.0) fac = CON2*fac errt = max(abs(a[j][i] - a[j-1][i]), abs(a[j][i] - a[j-1][i-1])) if errt <= err: err = errt ans = a[j][i] if abs(a[i][i] - a[i-1][i-1]) >= SAFE*err: return ans return ans def PDiff(f, var_index, *args): return 0

1 0

Derivatives
by Hassan Aurag 01 Mar '00

01 Mar '00

Hi, attached is a file called Derivative.py. It computes derivatives and is based on an algorithm found in Numerical Recipes in C. What to do you think about it and has anyone started a "serious" calculus oriented subpackage for Numerical Python in general? I mean: derivatives, partial derivatives, jacobian, hessian implemented fast and precise. On another note, why isn't infinity defined in NumPy? Why is tan(pi/2) a number even if big? Shouldn't it be infinity?

1 0

NumPy-Discussion March 2000