Planet SymPy

February 08, 2010

Jason Gedge

We got monkeys!

Yep, we do have monkeys. Blender monkeys, to be exact. I whipped up a simple loader for Wavefront OBJ models. Only loads the basic geometry now, so I have to work on the material stuff. One problem is my lack of a shader that does more advanced illumination, so that's something I have to work on. The two screenshots I've posted only do per-pixel lighting with Lambertian reflectance. I also want

by Jason G (noreply@blogger.com) at February 08, 2010 01:01 PM

February 05, 2010

Fredrik Johansson

Mpmath 0.14 released

I've just released version 0.14 of mpmath! Some highlights in this release have been covered in previous development updates on this blog:

  • mpmath in Sage to become 3x faster -- a Cython extension module will be added to Sage that greatly speeds up mpmath (the Sage patch is currently being reviewed; if all goes to plan, the patch will be accepted soon and the mpmath version in Sage will be updated to 0.14 at the same time).

  • Zeta evaluation with the Riemann-Siegel expansion -- Juan Arias de Reyna contributed code for very efficient and robust evaluation of the Riemann zeta function for arguments with large imaginary part.

  • Various features discussed in this update: 3D surface plotting (wrapping matplotlib, matrix calculus (transcendental functions of a matrix argument), Borel summation of divergent hypergeometric series, and options for whether to use Mathematica's conventions for hypergeometric functions have been added.

  • Improved accuracy for hypergeometric functions with large parameters, and analytic continuation implemented for 3F2, 4F3, and higher hypergeometric functions.

  • Support for using Python floats/complexes for faster low-precision math. This is very handy for plotting, multidimensional integration, and other things requiring many function evaluations.

There are many other changes as well. See the CHANGES file and the list of commits for more details. Thanks to Juan Arias de Reyna, Vinzent Steinberg, Jorn Baayen and Chris Smith who contributed patches, and various people who tested and submitted bug reports (thanks also to anyone else I forgot to mention).

I'm sorry that it took several months to finish this release! This was partly due to large internal reorganizations done in order to support floats and the Cython backend in Sage. I've also had to take some time off to focus on my studies and other things.

An example: spherical harmonics


Here I will demonstrate three new features with a single example: 3D plotting, one of the newly added special functions (spherical harmonics), and low-precision evaluation with the fp context. Of course, everything here can be done in arbitrary precision with mp as well; fp is just faster for plotting.

The following code visualizes spherical harmonics using a 3D polar plot. I've chosen to visualize the real part; the code is easily edited to show the absolute value (for example) instead.

from mpmath import fp

def Y(m,n):
def g(theta,phi):
R = fp.re(fp.spherharm(m,n,theta,phi))**2
x = R*fp.cos(phi)*fp.sin(theta)
y = R*fp.sin(phi)*fp.sin(theta)
z = R*fp.cos(theta)
return [x,y,z]
return g

fp.splot(Y(0,0), [0,fp.pi], [0,2*fp.pi])

The first few spherical harmonics are thus:
Y(0,0)
Y(1,0)
Y(1,1)
Y(2,0)
Y(2,1)
Y(2,2)
Y(3,0)
Y(3,1)
Y(3,2)
Y(3,3)


Some numerical examples



A selection of evaluations that were not implemented, failed, gave inaccurate results or were extremely slow in previous versions of mpmath:


>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True

# reflection formula for Barnes G-function
>>> barnesg(-100000+10000j)
(-7.332534002219787256775675e+22857579769 +
1.304469872717249495403882e+22857579770j)

# Riemann zeta function, large height
>>> zeta(0.5+100000000j)
(-3.362839487530727943146807 + 1.407234559646447885979583j)

# Accurate evaluation of large-degree Bernoulli polynomial
>>> bernpoly(1000,100)
4.360903799140670486890619e+1996

# Computation of Euler numbers
>>> eulernum(20)
370371188237525.0
>>> eulernum(50, exact=True)
-6053285248188621896314383785111649088103498225146815121L
>>> eulernum(2000000000000000000)
3.19651108713502662532039e+35341231273461153426

# Accurate evaluation of hypergeometric functions with large parameters
>>> hyp2f1(1000,50,25,-1)
-2.886992344761576738138864e-274
>>> legenp(0.5, 300, 0.25)
-1.717508549497252387888262e+578
>>> besseli(2, 10, derivative=100)
821.2386210064833486609255

# Evaluation of 4F3 and 4F2
>>> hyper([1,2.125,3.25,4.5], [5.25,6.5,7], 1+2j)
(1.006737556189825231039221 + 0.3058770612264986390003873j)
>>> hyper([1,2.125,3.25,4.5], [5.25,6.5], 1+2j)
(0.3220085070641791195103201 + 0.5251231752161637627314328j)

# Matrix functions
>>> A = matrix([[2,3,1+j],[1,0,-1],[2,1,5]])
>>> mnorm(A - logm(expm(A)))
9.540212722670391941827867e-25
>>> mnorm(A - expm(logm(A)))
4.375286550134565812513542e-26
>>> nprint(expm(2*A))
[ (-2855.4 + 10733.0j) (-1871.95 + 8597.9j) (-7721.42 + 12906.7j)]
[(-3910.68 - 1749.82j) (-3152.69 - 1272.11j) (-4460.18 - 3558.54j)]
[ (17039.0 + 22621.8j) (14320.3 + 17405.9j) (13587.7 + 35087.8j)]
>>> nprint(expm(A)**2)
[ (-2855.4 + 10733.0j) (-1871.95 + 8597.9j) (-7721.42 + 12906.7j)]
[(-3910.68 - 1749.82j) (-3152.69 - 1272.11j) (-4460.18 - 3558.54j)]
[ (17039.0 + 22621.8j) (14320.3 + 17405.9j) (13587.7 + 35087.8j)]
>>> nprint(chop(sqrtm(A)*sqrtm(A)))
[2.0 3.0 (1.0 + 1.0j)]
[1.0 0.0 -1.0]
[2.0 1.0 5.0]

# A convenience function for precise evaluation of exp(j*pi*x)
>>> expjpi(36); exp(j*pi*36)
(1.0 + 0.0j)
(1.0 + 4.702307256907087875509661e-25j)

# Some fp functions have specialized implementations
# for improved speed and accuracy
>>> fp.erfc(5); float(mp.erfc(5))
1.5374597944280347e-12
1.5374597944280349e-12
>>> fp.erf(5); float(mp.erf(5))
0.99999999999846256
0.99999999999846256
>>> fp.ei(-12-3j); complex(mp.ei(-12-3j))
(4.6085637890261843e-07-3.141592693919709j)
(4.6085637890261901e-07-3.141592693919709j)
>>> fp.zeta(-0.5); float(fp.zeta(-0.5)) # real-argument fp.zeta
-0.20788622497735468
-0.20788622497735468
>>> fp.ci(40); float(mp.ci(40))
0.019020007896208769
0.019020007896208765
>>> fp.gamma(4+5j); complex(mp.gamma(4+5j))
(0.14965532796078551+0.31460331072967695j)
(0.14965532796078521+0.31460331072967596j)


For more examples and images, see the posts linked at the top of this post!

by Fredrik Johansson (fredrik.johansson@gmail.com) at February 05, 2010 09:16 PM

February 03, 2010

Fredrik Johansson

mpmath in Sage to become 3x faster

I've blogged previously about the potential speedups attainable by rewriting core parts of mpmath in Cython (here and here), but all I had back then was some standalone classes and functions that couldn't easily be integrated. Well, I now finally have a fully working, fully integrated Cython implementation of the mpf and mpc types (plus related utility functions) -- and it successfully runs the entire mpmath test suite!

On my computer (Ubuntu 64 bit, AMD Athlon 64 3700+), running import mpmath; mpmath.runtests() in Sage takes this much time:

Before [mpmath (svn trunk) + Sage 4.3.1]: 62.75 seconds

After [mpmath (svn trunk with modifications) + Cython mpmath types + Sage 4.3.1 + a Sage performance patch by Robert Bradshaw]: 19.87 seconds

I have essentially only implemented arithmetic operations and conversions, so core transcendental functions (as well as square roots and powers) still fall back to the Python code, and they account for the majority of the running time. According to cProfile, about 20% of the total time is spent in exp() and gamma() alone. Hypergeometric series evaluation is also still done in Python. Replacing these functions with Cython versions should give a significant boost, and will be very easy to do in the future (in terms of code infrastructure). So realistically, the running time for the unit tests can be cut much further, probably in half.

The code isn't public yet. I will soon commit the changes to the mpmath svn repository and create a patch for Sage with the Cython extensions. The code just needs a little more cleanup, and there are no doubt a couple of subtle issues (such as corner-case memory leaks) that need fixing... but essentially, it's working, and the tests pass, so I expect it to be releasable quite soon.

Update: after adding just a little more Cythonized wrapper code, the time is now down to 17.64 seconds. I have still not touched exp, log, cos, sin, gamma, or any other core transcendental functions. Expect further improvements.

Update 2 (2010-02-03): preliminary version of the Cython code here

by Fredrik Johansson (fredrik.johansson@gmail.com) at February 03, 2010 03:04 AM

February 02, 2010

Jason Gedge

Game Design

So, one thing that I've been working on over the past week is reviving my "game engine" that I started working on a couple of years ago. It's not far, but I'm currently happy with the way things are going:Pretty simple, I know. Shows off some basic texturing and per-pixel lighting, but that's about all I got for now. Currently this is all done in OpenGL, but I've abstracted many of the concepts

by Jason G (noreply@blogger.com) at February 02, 2010 03:33 AM

February 01, 2010

Fabian Pedregosa

Scikit-learn 0.1

Today I released the first public version of Scikit-Learn (release notes).

It’s a python module implementing some machine learning algorithms, and it’s shaping quite good. For this release I did not want to do any incompatible changes, so most of them are just bug fixes and updates.

For the next release, however, some more radical changes are planned, and definitely something should be done about the (incredibly long) namespace, having to tape

from scikits.learn.machine.manifold_learning.regression.neighbors import Neighbors

each time you want to perform a nearest-neighbor algorithms is just not practical!

Here is a nice screenshot,

Scikit - learn

by fabian at February 01, 2010 01:32 PM

January 30, 2010

Jason Gedge

Current Work

I thought that since I haven't posted in awhile, I'd put something new up to let everyone know what's been up for the past few months.Well first, posts have been delayed because up until mid-December I was busy finishing the last of my courses required for my M.Sc. After that, I spent holidays both relaxing and working on my first paper. The paper focuses on improving stereo matching for

by Jason G (noreply@blogger.com) at January 30, 2010 12:10 AM

January 27, 2010

Fredrik Johansson

Using Sage numbers in mpmath

I've written a very basic mpmath context for computing "natively" with Sage's real and complex numbers. It can be tried out by upgrading mpmath in Sage to the svn trunk, applying this patch this patch to fix the helper code in Sage, and then importing this file.

Some examples:


----------------------------------------------------------------------
| Sage Version 4.3.1, Release Date: 2010-01-20 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: from mpsage import SageContext
sage: from mpmath import mp, timing
sage: mp.prec = 150
sage: sc = SageContext(150)
sage:
sage: print sc.quad(lambda t: sc.exp(-t**2), [0,sc.inf])
0.88622692545275801364908374167057259139877471
sage: print mp.quad(lambda t: mp.exp(-t**2), [0,mp.inf])
0.88622692545275801364908374167057259139877473
sage: timing(sc.quad, lambda t: sc.exp(-t**2), [0,sc.inf])
0.40903496742248535
sage: timing(mp.quad, lambda t: mp.exp(-t**2), [0,mp.inf])
0.43610286712646484
sage:
sage: print sc.zeta(4+5*I)
0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897*I
sage: print mp.zeta(4+5*I, method='euler-maclaurin')
(0.95121830700949569861514024576822624312072806 + 0.024932824507357722259855155244740571939365897j)
sage: timing(sc.zeta, 4+5*I)
0.015111517906188966
sage: timing(mp.zeta, 4+5*I, method='euler-maclaurin')
0.028736209869384764
sage:
sage: print sc.zeta(1+100000000*I)
1.8349308953278466065175598876062914152382527 + 1.0691847904236128969174476736978194200591565*I
sage: print mp.zeta(1+100000000*I)
(1.8349308953278466065175598876062914152563156 + 1.069184790423612896917447673697819420203942j)
sage: timing(sc.zeta, 1+100000000*I)
1.0952680110931396
sage: timing(mp.zeta, 1+100000000*I)
2.0537140369415283


Unfortunately, there are some issues (besides the fact that some methods are missing, so not everything works yet).

This context doesn't provide variable precision, so the user has to manually allocate sufficient extra precision to compensate for rounding errors.

I first tried to do it, but variable precision is very inconvenient to implement using Sage's way of managing precision. There is no direct way to perform operations with a given target precision (independent of the inputs), and the second best option (which is to perform coercions to the target precision everywhere) is very slow, besides losing accuracy. The only way to implement variable precision in a reasonable way is to bypass Sage's RealNumber and ComplexNumber (at least their public interfaces) and wrap MPFR directly using a precision model similar to what MPFR and mpmath uses, where the precision of the result is always specified and independent of the inputs.

Secondly, there is not that much of a speedup (in the examples above, the speedup is at most about 2x). This is mainly due to the fact that the context uses wrapper classes for RealNumber and ComplexNumber, with all interface and conversion code written in Python. So the overhead is about the same as for the corresponding code in vanilla mpmath (where it accounts for about 50% of the total overhead). The reason RealNumber and ComplexNumber can't be used directly, even in a fixed-precision setting, is that mpmath in many places multiplies numbers by floats (usually exact float literals like 0.5), and Sage always coerces to the number with less precision. This could presumably be fixed by replacing all float and complex constants in mpmath, but I'm not in the mood to do that right now.

There should be a significant performance benefit if direct-from-Cython classes and conversion methods were used. It's also very important to optimize certain core functions; for example, quadrature would benefit greatly from a fast fdot implementation.

All things considered, I'm probably not going to continue much more down this particular path. It's better to write a fully compatible Cython context with classes designed directly for mpmath. Trying to wrap Sage's numbers did however help identify a few problems with the existing interfaces, so I might extend the work on this context a little more just to find more such issues.

by Fredrik Johansson (fredrik.johansson@gmail.com) at January 27, 2010 08:27 PM

January 19, 2010

Fredrik Johansson

Zeta evaluation with the Riemann-Siegel expansion

I'm very grateful to Juan Arias de Reyna who has contributed a module implementing zeta function evaluation using the Riemann-Siegel formula in mpmath. This finally permits rapid evaluation of the Riemann zeta function for arguments with large imaginary part.

To follow tradition on this blog, pictorial proof shall be given. Here is a plot of a segment of the critical line, ζ(1/2+ti) with t between 1010+50 and 1010+55:

A complex plot of showing the critical strip, t ranging between 108+40 and 108+45 (note: the y scale is reversed):



Juan Arias de Reyna, who is a professor of mathematics at the University of Seville, has done a thorough job with this code. He has even proved rigorous error bounds for his algorithm (subject to assumptions that the underlying functions in mpmath being well-implemented). The news is that the bounds -- documented in an as-yet unpublished paper -- are valid off the critical line. The code also computes derivatives (up to 4th derivatives), although not as rigorously but still very robustly.

I integrated the code (and added a few optimizations) during the last few days. The zeta function in mpmath now automatically switches between Borwein's algorithm (close to the real line), Euler-Maclaurin summation, and the Riemann-Siegel expansion.

Some example values and timings on my laptop, computing ζ(1/2+10ni) for n up to 12:

>>> from timeit import default_timer as clock
>>> mp.dps = 25
>>>
>>> for n in range(13):
... t1 = clock(); y = zeta(0.5 + 10**n*j); t2 = clock()
... print n, y, round(t2-t1,3)
...
0 (0.1439364270771890603243897 - 0.7220997435316730891261751j) 0.02
1 (1.544895220296752766921496 - 0.1153364652712733754365914j) 0.003
2 (2.692619885681324090476096 - 0.02038602960259816177072685j) 0.042
3 (0.3563343671943960550744025 + 0.9319978312329936651150604j) 0.073
4 (-0.3393738026388344575674711 - 0.03709150597320603147434421j) 0.434
5 (1.073032014857753132114076 + 5.780848544363503984261041j) 0.167
6 (0.07608906973822710000556456 + 2.805102101019298955393837j) 0.146
7 (11.45804061057709254500227 - 8.643437226836021723818215j) 0.181
8 (-3.362839487530727943146807 + 1.407234559646447885979583j) 0.336
9 (-2.761748029838060942376813 - 1.677512240989459839213205j) 0.877
10 (0.3568002308560733825395879 + 0.286505849095836103292093j) 2.695
11 (0.6436639255801185727194357 + 0.1168615914616853418448829j) 8.583
12 (2.877961809278403355251079 - 3.206771071318398923493412j) 26.934


Similar results off the critical line:

>>> for n in range(13):
... t1 = clock(); y = zeta(1.0 + 10**n*j); t2 = clock()
... print n, y, round(t2-t1,3)
...
0 (0.5821580597520036481994632 - 0.9268485643308070765364243j) 0.021
1 (1.390287313237401426796005 - 0.109785153066302056909746j) 0.004
2 (1.632833506686711866610705 - 0.06813120384181249010120548j) 0.043
3 (0.9409368682927533108010138 + 0.04522665207209509908865644j) 0.083
4 (0.4973279229716308441790286 - 0.5878238243194009766923214j) 0.598
5 (1.618122122846936796567759 + 1.070441041470623686626035j) 0.233
6 (0.9473872625104789104802422 + 0.5942199931209183283333071j) 0.195
7 (2.859846483332530337008882 + 0.491808047480981808903986j) 0.409
8 (1.83493089532784660651756 + 1.069184790423612896917448j) 0.455
9 (0.9038018561650776977609945 - 1.189857828822373901473908j) 1.393
10 (0.5418173564211820524034624 + 0.635303581895880322679247j) 3.824
11 (0.5365466615361310937110304 - 0.1234443975100346650640542j) 12.031
12 (0.8225630719679733497367277 - 0.4484762282040223492401122j) 38.061


The implementation also supports use of the fp context. Double precision unavoidably becomes insufficient as the imaginary part approaches 1015, but it has the advantage of speed in the range where it works:

>>> for n in range(13):
... t1 = clock(); y = fp.zeta(0.5 + 10**n*1j); t2 = clock()
... print n, y, round(t2-t1,3)
...
0 (0.143936427077-0.722099743532j) 0.007
1 (1.5448952203-0.115336465271j) 0.001
2 (2.69261988568-0.0203860296026j) 0.003
3 (0.356334367195+0.931997831233j) 0.004
4 (-0.339373802616-0.0370915059691j) 0.123
5 (1.07303201485+5.78084854433j) 0.005
6 (0.076089072636+2.80510210471j) 0.006
7 (11.4580404601-8.64343725721j) 0.006
8 (-3.36283920435+1.40723433071j) 0.011
9 (-2.76174796643-1.67750591108j) 0.028
10 (0.356829034142+0.286525774475j) 0.083
11 (0.64314751322+0.116713738713j) 0.256
12 (2.8689206645-3.21135962379j) 0.808


We have done some comparison with Mathematica, and the mpmath version appears to be about as fast (a bit faster or a bit slower, sometimes substantially faster, depending on circumstance). The most expensive part of the computation occurs in a simple internal function that adds lots of ns terms. I think for Sage, it will be very easy to switch to a Cython version of this function which should improve speed by a large factor.

But most importantly, Mathematica's Zeta[] is notoriously buggy for large imaginary arguments. As a first example, here Mathematica 7.0 computes two entirely different values at different precisions:

In[3]:= N[Zeta[1/4+10^12 I], 15]
Out[3]= -0.0125397 + 0.0139723 I
In[4]:= N[Zeta[1/4+10^12 I], 30]
Out[4]= 358.066828240154490750947835567 - 580.623633400912069466146291259 I


With mpmath:

>>> mp.dps = 15; zeta(0.25+1e12j)
(358.066828240154 - 580.623633400912j)
>>> mp.dps = 30; zeta(0.25+1e12j)
(358.066828240154490750947835567 - 580.623633400912069466146291259j)


As a second example, if Mathematica is asked for derivatives, it's more likely than not to return complete nonsense, and even increasing the precision doesn't help:

In[2]:= N[Zeta'[1/2+10^6 I], 15]

883 883
Out[2]= 2.48728166483172 10 - 7.66644043045624 10 I

In[3]:= N[Zeta'[1/2+10^6 I], 25]

940 940
Out[3]= 1.586685034587255948191759 10 + 2.158475809806136995106119 10 I

In[4]:= N[Zeta'[1/2+10^6 I], 30]

1022
Out[4]= -1.071044014417407205715473623855 10 -

1021
> 2.73478073192015960107298455871 10 I


For contrast, mpmath computes derivatives perfectly:

>>> mp.dps = 15; zeta(0.5+1e6j, derivative=1)
(11.6368040660025 - 17.127254072213j)
>>> mp.dps = 25; zeta(0.5+1e6j, derivative=1)
(11.63680406600252145919591 - 17.12725407221299600357895j)
>>> mp.dps = 30; zeta(0.5+1e6j, derivative=1)
(11.6368040660025214591959071246 - 17.1272540722129960035789468265j)
>>> diff(zeta, 0.5+1e6j) # Numerical verification
(11.6368040660025214591959071246 - 17.1272540722129960035789468265j)


That's all for now. I'm back in school again, so maybe I won't have as much time for programming in the near future. But my schedule is flexible, so we'll see. A new release of mpmath shouldn't be far away.

Update: I should point out that the bugs in Mathematica's Zeta[] only seem to occur in recent versions. Mathematica 5 and older versions do not seem to have these problems.

by Fredrik Johansson (fredrik.johansson@gmail.com) at January 19, 2010 10:38 PM

January 13, 2010

Fredrik Johansson

YAMDU (yet another mpmath development update)

At the moment, I'm able to get quite a bit of work done on mpmath. This includes general code cleanup, gardening the issue tracker, and finishing features that were works-in-progress for a long while.

Today's topics are:

3D plotting
Matrix calculus
Borel summation of divergent hypergeometric series
Optional consistency with Mathematica
Internal reorganization

3D plotting


Jorn Baayen submitted a patch that adds 3D surface plotting (along with some other minor changes to the visualization module). Thanks a lot! I previously blogged about 3D plotting using matplotlib.

You can now easily produce plots of functions of the form z = f(x,y):

splot(lambda x, y: 10*sinc(hypot(x,y)), [-10,10], [-10,10])



You can just as easily plot more general parametric functions x,y,z = f(u,v). The following plots a Möbius strip:

def f(u,v):
d = (1+0.5*v*cos(0.5*u))
x = d*cos(u)
y = d*sin(u)
z = 0.5*v*sin(0.5*u)
return x,y,z

splot(f, [0,2*pi], [-1,1])



(I'm not sure why there are black spots in the image. This seems to be a matplotlib bug.)

Matrix calculus


It's now possible to compute exponentials, sines, cosines, square roots, logarithms, and complex powers (Az) of square matrices:


>>> mp.dps = 25; mp.pretty = True
>>> A = matrix([[2,3,1+j],[1,0,-1],[2,1,5]])
>>> nprint(expm(A))
[ (27.6034 + 41.2728j) (24.9924 + 31.1189j) (8.67256 + 68.271j)]
[(-15.0962 + 1.51713j) (-10.1713 + 1.30035j) (-28.7003 - 0.401789j)]
[ (125.368 + 37.7417j) (98.9812 + 26.9693j) (158.616 + 85.1593j)]
>>> nprint(logm(expm(A)))
[(2.0 - 5.35398e-26j) (3.0 + 6.9172e-26j) (1.0 + 1.0j)]
[(1.0 + 4.94422e-26j) (4.42102e-25 - 1.14411e-25j) (-1.0 - 1.95948e-27j)]
[(2.0 + 1.58966e-26j) (1.0 - 8.57028e-27j) (5.0 + 3.23653e-27j)]
>>> nprint(sqrtm(A)**2)
[(2.0 - 8.25839e-30j) (3.0 - 8.40399e-30j) (1.0 + 1.0j)]
[(1.0 + 4.63764e-30j) (4.84564e-30 - 3.27721e-31j) (-1.0 + 7.47261e-31j)]
[(2.0 - 4.31871e-30j) (1.0 + 1.78726e-31j) (5.0 + 2.54582e-29j)]
>>> nprint(powm(powm(A,1+j),1/(1+j)))
[(2.0 - 1.12995e-26j) (3.0 - 8.93158e-27j) (1.0 + 1.0j)]
[(1.0 - 1.54352e-26j) (9.23906e-27 - 1.67262e-26j) (-1.0 - 2.62243e-27j)]
[(2.0 + 9.97431e-27j) (1.0 + 1.56341e-26j) (5.0 - 1.79194e-26j)]


The code also works with the double precision fp context, which obviously is much faster for large matrices. When computing square roots and logarithms, most of the time is spent on matrix inversion, which can be accelerated by substituting in numpy:


>>> from mpmath import fp
>>> A = fp.randmatrix(20)
>>> B = fp.sqrtm(A)
>>> timing(fp.sqrtm, A)
3.6470590535948006
>>>
>>> import numpy
>>> fp.inverse = lambda A: fp.matrix(numpy.linalg.inv(A.tolist()))
>>>
>>> timing(fp.sqrtm, A)
1.0815329881311735
>>> C = fp.sqrtm(A)
>>>
>>> fp.mnorm(A-B**2)
3.5001514923930131e-14
>>> fp.mnorm(A-C**2)
2.6462503380426079e-14


It could be much faster still by doing everything in numpy. Probably in the future fp should be improved to seamlessly use numpy internally for all available matrix operations. Patches are welcome!

Borel summation of divergent hypergeometric series


The hypergeometric series of degree p > q+1 are divergent for |z| > 0. Nevertheless, they define asymptotic expansions for analytic functions which exist in the sense of Borel summation. Previously, mpmath knew only how to evaluate 2F0 (whose Borel sum has a closed form), but now it can evaluate the functions of higher degree as well.

To illustrate, the cosine integral Ci(z) has an asymptotic series representation in terms of 3F0:



This representation is efficient for very large values of |z|, where the argument of the hypergeometric series is small. But with z = 0.5, say, the series does not yield a single digit. With a value of z = 10 or so, the series only gives about 3 digits with an optimal truncation:

>>> z = 10.0
>>> for n in range(1,16):
... print sum(rf(0.5,k)*rf(1,k)*rf(1,k)*(-4/(z**2))**k/fac(k) for k in range(n))
...
1.0
0.98
0.9824
0.98168
0.9820832
0.98172032
0.9821993216
0.981327538688
0.9834198176768
0.977017443971072
1.00134646405284
0.888946391275078
1.50939479300832
-2.52351981825774
27.9653146429137


Of course, there are ways to compute the cosine integral using convergent series. But if we pretend that we only know about the asymptotic series, then the Borel regularization means that we can evaluate the function to high precision even for small z:

>>> mp.dps = 25; mp.pretty = True
>>> z = mpf(0.5)
>>> u = -4/z**2
>>> H1 = hyper([1,1,1.5],[],u)
>>> H2 = hyper([0.5,1,1],[],u)
>>>
>>> print log(z)-log(z**2)/2-cos(z)/z**2*H1+sin(z)/z*H2
-0.1777840788066129013358103
>>> print ci(z)
-0.1777840788066129013358103


The z = 10 series above, to high precision:

>>> hyper([0.5,1,1],[],-4/(10.0**2))
0.9819103501017016869905255
>>> mp.dps = 40
>>> hyper([0.5,1,1],[],-4/(10.0**2))
0.9819103501017016869905255431829554224704


It's instructive to visualize the optimal truncations of an asymptotic series compared to the exact solution. Notice how, for an alternating series like this, the truncations alternate between over- and undershooting:

def optimal_asymp(z):
u = -4/(z**2)
term_prev = inf
s = 0.0
for k in range(30):
term = rf(0.5,k)*rf(1,k)*rf(1,k)*u**k/fac(k)
if abs(term) abs(term_prev):
s += term
term_prev = term
else:
break
return s

def exact(z):
u = -4/(z**2)
return hyper([0.5,1,1],[],u)

mp.dps = 10
plot([optimal_asymp, exact], [0,10])



The catch with this feature? Computing the Borel regularization involves evaluating (nested) numerical integrals with a hypergeometric function in the integrand. Except for special parameters (where degree reductions happen internally -- the Ci series above happens to be such a case), it's not very fast.

Optional consistency with Mathematica


I often try to follow Mathematica's conventions regarding limits, special values and branch cuts of special functions. This simplifies testing (I can directly compare values with Mathematica), and usually Mathematica's conventions seem well-reasoned.

There are exceptions, however. One such exception concerns Mathematica's interpretation of 2F1(a,b,c,z) for b = c a negative integer. Mathematica interprets this as a polynomial, which doesn't make much sense since the denominator of the zero term is zero. It's not even consistent with how Mathematica evaluates this function for a symbolic variable:

In[3]:= Hypergeometric2F1[3,-2,-2,2]

Out[3]= 31

In[4]:= Hypergeometric2F1[3,b,b,2]

Out[4]= -1


Maybe there is a good reason for doing what Mathematica does, but it's at the very least not documented anywhere. I've now changed mpmath to instead interpret the two parameters as eliminating each other and giving a 1F0 function (which is what Mathematica does for a generic value). The Mathematica-compatible value can be recovered by explicitly disabling parameter elimination:

>>> mp.pretty = True
>>>
>>> hyp2f1(3,-2,-2,2)
-1.0
>>> hyp2f1(3,-2,-2,2,eliminate=False)
31.0


On a related note, I've fixed the Meijer G-function to switch between its z and 1/z forms automatically to follow Mathematica's definition of the function. This introduces discontinuities in the function for certain orders. The user can explicitly choose which form to use (so as to obtain a continuous function) with series=1 or series=2.

Internal reorganization


In a long overdue change, I've moved many of the modules in mpmath into subdirectories. The multiprecision library code is now located in mpmath/libmp; special functions are in mpmath/functions, calculus routines are in mpmath/calculus, and routines for matrices and linear algebra are in mpmath/matrices.

This shouldn't affect any documented interfaces, but it will break external code that directly uses mpmath internal functions. The mpmath interfaces in SymPy and Sage will need some editing for the next version update. The breakage should be straightforward to fix (mostly just a matter of changing imports).

Since the next version of mpmath is definitely going to break some code, I might use the opportunity to do some other cosmetic interface changes as well. Ideally, after the next version, the interface will be stable until and beyond mpmath 1.0 (whenever that happens). But the version number is still at 0.x for a reason -- no compatibility guarantees.

by Fredrik Johansson (fredrik.johansson@gmail.com) at January 13, 2010 01:54 PM

January 07, 2010

Fredrik Johansson

Accurate hypergeometric functions for large parameters

Holiday season means more free time for coding (yay!). In this commit, I've fixed the remaining major accuracy issues with hypergeometric functions in mpmath. The effect, generally speaking, is that mpmath will now deal much more robustly with large parameters of hypergeometric-type functions.

Previously, you would obtain an accurate value with parameters of magnitude ≈ 10−1 - 101 (say), and often for large parameters as well, but for some functions the accuracy would quickly deteriorate if parameters were increased (or moved closer to singularities). You could still obtain any accuracy by simply increasing the working precision, but you had to manually figure out the amount of extra precision required; the news is that mpmath now automatically gives a value with full accuracy even for large parameters (or screams out loud if it fails to compute an accurate value).

This doesn't mean that you can't trick mpmath into computing a wrong value by choosing sufficiently evil parameters, but it's much harder now than simply plugging in large values. Abnormally large values will now mainly accomplish abnormal slowness (while mpmath implements asymptotic expansions with respect to the argument of hypergeometric functions, evaluation for asymptotically large parameters is a much harder problem as far as I'm aware -- so mpmath in effect accomplishes it by brute force.)

The most trivial change was to change the series summation to aim for control of the relative error instead of the absolute error. This affects alternating series, where large parameters lead to very small function values. Previously, something like the following would happen:


>>> mp.dps=5; hyp1f1(-5000, 1500, 100)
1.6353e-14
>>> mp.dps=15; hyp1f1(-5000, 1500, 100)
8.09813050863682e-25
>>> mp.dps=30; hyp1f1(-5000, 1500, 100)
-1.38318247777802583583082760724e-39


This isn't disastrously bad since the absolute error is controlled (summing this series naively in floating-point arithmetic might give something like 1e+234 due to catastrophic cancellation). But now, you get the relatively accurate answer right away, which is much nicer:


>>> mp.dps = 5; hyp1f1(-5000, 1500, 100)
9.1501e-185
>>> mp.dps = 15; hyp1f1(-5000, 1500, 100)
9.15012134245639e-185
>>> mp.dps = 30; hyp1f1(-5000, 1500, 100)
9.15012134245639443114220499541e-185


Of course, if the value is too small, the evaluation loop must abort eventually. The default behavior now is to raise an exception if relative accuracy cannot be obtained. The user can force a return value by either permitting a higher internal precision before aborting, or specifying a size 2−zeroprec below which the value is small enough to be considered zero:


>>> hyp1f1(-8000, 4000, 1500)
Traceback (most recent call last):
...
ValueError: hypsum() failed to converge to the requested 53 bits of accuracy
using a working precision of 3568 bits. Try with a higher maxprec,
maxterms, or set zeroprec.
>>>
>>>
>>> hyp1f1(-8000, 4000, 1500, maxprec=10000)
4.36754212717293e-1350
>>> hyp1f1(-8000, 4000, 1500, accurate_small=False)
-1.99286380611911e-25
>>> hyp1f1(-8000, 4000, 1500, zeroprec=1000)
0.0


Since exceptions are inconvenient, it will be necessary to add some more symbolic tests for exact zeros for certain functions (particularly orthogonal polynomials) where exact zeros arise as important special cases.

Also, cancellation detection in hypercomb is now the default for all functions. This fixes, among other things, a bug in the Meijer G-function reported by a user via email. Before:


>>> meijerg([[],[]],[[0],[]],0.1)
0.90483741803596
>>> meijerg([[],[]],[[0,0],[]],0.1)
1.47347419562219
>>> meijerg([[],[]],[[0,0,0],[]],0.1)
1.61746025085449
>>> meijerg([[],[]],[[0,0,0,0],[]],0.1) # suspicious
13194139533312.0


After:

>>> meijerg([[],[]],[[0],[]],0.1)
0.90483741803596
>>> meijerg([[],[]],[[0,0],[]],0.1)
1.47347419562219
>>> meijerg([[],[]],[[0,0,0],[]],0.1)
1.61745972140487
>>> meijerg([[],[]],[[0,0,0,0],[]],0.1)
1.56808223438324


Another important change is more correct handling parameters very close to negative integers, particularly those appearing in denominators of series. Previously, unless the integer n was small enough or the precision was set high enough, this situation would yield a bogus value (that of a prematurely truncated series):


>>> n = -20
>>> nh = fadd(n, 1e-100, exact=True)
>>> hyp0f1(n, 0.25) # nonsense
0.987581857511351
>>> hyp0f1(nh, 0.25) # same nonsense
0.987581857511351


Now, correctly:


>>> n = -20
>>> nh = fadd(n, 1e-100, exact=True)
>>> hyp0f1(n, 0.25)
Traceback (most recent call last):
...
ZeroDivisionError: pole in hypergeometric series
>>> hyp0f1(nh, 0.25)
1.85014429040103e+49


Finally, and probably most importantly, the limit evaluation code has been changed to adaptively decrease the size of perturbation until convergence. Under fairly general assumptions, the maximum accurate perturbation at precision p can easily be shown to be 2−(p+C); the problem is that the parameter-dependent constant C isn't generally known. Previously C was just set to a small value (10), and naturally this would break down for some functions when parameters were increased beyond roughly that magnitude.

I briefly considered the idea of estimating C analytically in terms of the parameters, and while I think this can be done rigorously, it seems difficult -- especially to do it tightly (grossly overestimating C would murder performance). The algorithm implemented now is quasi-rigorous, and although there is some slowdown (sometimes by a fair amount), the improved reliability is definitely worth it. A user can also manually supply a perturbation size of their own taste, thereby overriding adaptivity. If the user supplies an intelligent value, this gives both the speed of the old code and full rigor. Probably some intelligent choices for particular functions could be made automatically by mpmath too, to recover the old speed in common cases.

An example of a situation where this becomes necessary is for Legendre functions with certain large parameter combinations. With the old code:

>>> mp.dps = 15; legenp(0.5, 300, 0.25) # bad
4.19317029788723e+626
>>> mp.dps = 30; legenp(0.5, 300, 0.25) # bad
3.72428336871098039162972441784e+611
>>> mp.dps = 60; legenp(0.5, 300, 0.25) # bad
2.93794154954090326636196697693611381787845107728298382876544e+581
>>> mp.dps = 100; legenp(0.5, 300, 0.25) # only accurate to a few digits
-1.71750854949725238788826203712778687036438365374945625996246145924802366061559
6579520831362887006032e+578


With the new code:

>>> mp.dps = 15; legenp(0.5, 300, 0.25)
-1.71750854949725e+578
>>> mp.dps = 30; legenp(0.5, 300, 0.25)
-1.71750854949725238788826203713e+578
>>> mp.dps = 60; legenp(0.5, 300, 0.25)
-1.71750854949725238788826203712778687063419097363472262860567e+578
>>> mp.dps = 15; legenp(0.5, 300, 0.25, hmag=300) # faster, with manual perturbation
-1.71750854949725e+578


Another example is for Bessel functions differentiated to high order. With the old code:

>>> mp.dps = 15; besseli(2, 10, derivative=100) # bad
2.63560662943646e+34
>>> mp.dps = 30; besseli(2, 10, derivative=100) # bad
23408889310840499424.9813614712
>>> mp.dps = 60; besseli(2, 10, derivative=100) # only accurate to a few digits
821.238621006501815018537509753672810563116338269219460709828


With the new code:

>>> mp.dps = 15; besseli(2, 10, derivative=100)
821.238621006483
>>> mp.dps = 30; besseli(2, 10, derivative=100)
821.238621006483348660925541651
>>> mp.dps = 60; besseli(2, 10, derivative=100)
821.238621006483348660925541650744338655411830854860813048862


In related news, I've given all hypergeometric-type functions the ability to accept kwargs and pass them on to hypercomb and hypsum. This means that tuning and error control options will be available for a large number of functions (Bessel functions, etc), which could be useful for some users who need to do advanced calculations. I have yet to document these features thoroughly (the interface will also perhaps be tweaked before the next release).

by Fredrik Johansson (fredrik.johansson@gmail.com) at January 07, 2010 01:54 PM

Fabian Pedregosa

scikit-learn project on sourceforge

This week we created a sourceforge project to host our development of scikit-learn. Although the project already had a directory in scipy’s repo, we needed more flexibility in the user management and in the mailing list creation, so we opted for SourceForge.

To be honest, after using git and Google Code for bug tracking, I was not very excited about using subversion/sourceforge again. On the other hand, we needed some sort of compromise that would allow a very heterogeneous range of developers to work together, and after some (surprisingly civilized) emails and some chatting with Gael, we agreed that SourceForge was indeed the best choice.

In case you are interested, there’s a (preliminary) web page with more info. You might also want to have a look at the previous project’s web page.

by fabian at January 07, 2010 01:17 PM

December 29, 2009

Aaron Meurer

Automatically Remove Trailing Whitespace in XCode


I like XCode, and I use it to edit all of my source for SymPy. But, like many editors, it likes to auto-indent new lines to the level of indentation of the previous line. This is a useful feature, but it makes for training whitespace out the wazoo, since blank lines will be indented in. I am constantly finding myself using SymPy’s strip_whitespace script to clean up my files.

This bugged me enough that I Googled a solution, and found this. It is a simple XCode plugin that, among other things, adds an option to strip trailing whitespace on save. Just install in the PlugIns folder in the XCode package and enable the option in the new Google pane of the XCode preferences.

by asmeurer at December 29, 2009 11:56 PM

December 28, 2009

Official SymPy blog

SymPy 0.6.6 released

SymPy 0.6.6 has been released on December 20, 2009. It is available at

http://sympy.org


The source distribution can be downloaded from:
http://sympy.googlecode.com/files/sympy-0.6.6.tar.gz

You can get the Windows installer here:
http://sympy.googlecode.com/files/sympy-0.6.6.win32.exe

And the html documentation here:
http://sympy.googlecode.com/files/sympy-0.6.6-docs-html.zip

About SymPy

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python.

Changes since last stable release


  • many documentation improvements, including docstrings and doctests
  • new assumptions system (GSoC) (See assumptions documentation for more information or have a look at Fabian's blog.)
    • Note: This is going to replace the old assumption system. It is encouraged to use it for new code, however it is not completely finished and parts of sympy have yet to be rewritten to use it; this scheduled for the 0.7 release.
  • improvements to test runner
  • printing improvements (especially LaTeX, but also mathml and pretty printing)
  • discriminant of polys
  • block diagonal methods for matrices
  • vast improvements to solving of ODEs (GSoC) (See ODE documentation for full details or Aaron's blog).
  • logcombine function
  • improvements to sets
  • better trigonometric simplification
  • improvements to piecewise functions
  • improvements to solve() and nsolve()
  • improvements to as_numer_denom()
  • much better quartic and cubic polynomial rootfinding
  • code refactoring and cleanup
  • physics: coupled clusters and wick expansion
  • matrices: symbolic QR solving
  • mpmath updated
  • pyglet updated
  • many, many bug fixes and small improvements

The following 22 people have contributed patches to this release (sorted alphabetically):

  • Aaron Meurer
  • Alan Bromborsky
  • Andy R. Terrel
  • Bill Flynn
  • Chris Smith
  • Eh Tan
  • Fabian Pedregosa
  • Fredrik Johansson
  • Jorn Baayen
  • Julio Idichekop Filho
  • Kevin Goodsell
  • Łukasz Pankowski
  • Luke Peterson
  • Øyvind Jensen
  • Ondrej Certik
  • Oscar Benjamin
  • Priit Laes
  • Renato Coutinho
  • Ronan Lamy
  • Ryan Krauss
  • Ted Horst
  • Toon Verstraelen
  • Vinzent Steinberg

The following 10 people helped reviewing patches:

  • Aaron Meurer
  • Andy R. Terrel
  • Chris Smith
  • Fabian Pedregosa
  • Fredrik Johansson
  • Luke Peterson
  • Mateusz Paprocki
  • Ondrej Certik
  • Ronan Lamy
  • Vinzent Steinberg

by Vinzent Steinberg (noreply@blogger.com) at December 28, 2009 01:57 PM

December 22, 2009

Fredrik Johansson

Analytic continuation of 3F2, 4F3 and higher functions

As of a recent commit, mpmath can evaluate the analytic continuation of the generalized hypergeometric function p+1Fp for any p. Previously 2F1 (the Gaussian hypergeometric function) was supported -- see earlier blog posts -- but not 3F2 and higher. This addition means that the generalized hypergeometric function is finally supported essentially everywhere where it is "well-posed" (in the sense that the series has nonzero radius of convergence), so it is a rather significant improvement. Unfortunately, the implementation is still not perfect, but I decided to commit the existing code since it is quite useful already (and long overdue).

As proof of operation, I deliver plots of 3F2 and 4F3, requiring both |z| < 1 and |z| > 1:

from mpmath import *
f1 = lambda z: hyp3f2(1,2,3,4,5,z)
f2 = lambda z: hyper([1,2,3,4],[5,6,7],z)
plot([f1,f2], [-2,2])




A portrait of 3F2 restricted to the unit circle:

plot(lambda x: hyp3f2(1,2,3,4,5,exp(2*pi*j*x)), [-1,1])



A numerical value of 5F4 with z on the unit circle, with general complex parameters:

>>> mp.dps = 50
>>> print hyper([1+j,0.5,-2j,1,0.5+0.5j],[0.5j,0.25,-j,-1-j],-1)
(-1.8419705729324526212110109087877199070037836117341 -
0.57799141426978758652682670969879368618437786161612j)


High precision values of 3F2 at z = 1 and z = 1 + ε:

>>> print hyp3f2(1,1,2,3.5,1.5,1)
2.2603241503205748814116375624327119385797950825522
>>> print hyp3f2(1,1,2,3.5,1.5,'1.0001')
(2.2626812356987790952469291649495098300894035980837 -
0.00092505569713084902188662960467216683326695892509251j)


A complex plot of 3F2:

>>> mp.dps = 5
>>> cplot(lambda z: hyp3f2(2.5,3,4,1,2.25,z), [-2,2], [-2,2],
... points=50000)



A bit of theoretical background: the hypergeometric series pFq has an infinite radius of convergence when p ≤ q, so it can in principle be evaluated then by adding sufficiently many terms at sufficiently high precision (although in practice asymptotic expansions must be used for large arguments, as mpmath does). In the balanced case when p = q+1, i.e. for 1F0(a,z), 2F1(a,b,c,z), 3F2(...), 4F3(...), ..., the series converges only when |z| 1 (plus in some special instances). This is due to the fact that the hypergeometric function has a singularity (a pole or branch point, depending on parameters) at z = 1, in turn owing to the fact that the hypergeometric differential equation is singular at z = 1. (The reason that the p = q+1 and not p = q case is "balanced" is that there is an extra, implicit factorial in the hypergeometric series.)

The main ingredient of the analytic continuation of p+1Fp is the inversion formula which replaces z with 1/z and thus handles |z| > 1. This was easy to implement -- the only complication is that integer parameters result in singular gamma factors, but the mechanisms to handle those automatically were already in place. There was no particular reason why I hadn't added that code already.

The tricky part is the unit circle, and the close vicinity thereof, where neither series converges (quickly). I've been looking for good ways handle this, with mixed results.

When I posed the problem on this blog several months back, a reader suggested this paper by Wolfgang Bühring which gives a series for the analytic continuation around z = 1. My finding from trying to implement it is that the rate of convergence of the series generally is poor, and therefore it is not immediately effective for computation. However, convergence acceleration with nsum improves the situation considerably in many cases. Some parameter combinations render the convergence acceleration useless, but even then, it can give a few correct digits, so it is better than nothing (although the implementation should probably warn the user when the result probably is inaccurate). I'm unfortunately not aware of any parameter transformations that would substantially improve convergence. The current implementation uses this method for 3F2 when |z-1| is small; it should work for 4F3 and higher too, but the series coefficients are much more complicated (involving multiply nested sums), so that's yet to be done.

For the rest of the unit circle, I've settled for simply using convergence acceleration directly. This essentially just amounts to passing the hypergeometric series to nsum, which applies Richardson extrapolation and iterated Shanks transformations. The Shanks transformation is actually perfect for this -- it's almost tailor made for convergence acceleration of the balanced hypergeometric series -- and is able to sum the series even outside the circle of convergence, using just a few terms. This covers most of the unit circle -- the catch is just that the acceleration asymptotically deteriorates and ultimately becomes useless close to z = 1, so complementary method (such as the Bühring's is still required there).

Mathematica seems to support unit-circle evaluation of hypergeometric functions quite well. Unfortunately, I don't know how it does it. According to Wolfram's Some Notes On Internal Implementation page,


The hypergeometric functions use functional equations, stable recurrence relations, series expansions, asymptotic series and Padé approximants. Methods from NSum and NIntegrate are also sometimes used.


This looks similar to what I'm doing -- high order Padé approximants and Mathematica's NSum should be equivalent to the nsum-based series acceleration in mpmath. But probably Mathematica employs some more tricks as well.

I've also tested direct integration of the hypergeometric differential equation and of Mellin-Barnes integral representations, but these approaches don't seem to work much better (at least not without many improvments), and at best seem to be relatively slow.

by Fredrik Johansson (fredrik.johansson@gmail.com) at December 22, 2009 06:51 PM

Fabian Pedregosa

Winter in Paris is not funny

This week I arrived to the place where I will be working the following two years: Neurospin.

Neurospin

It’s a research center located 20 km from Paris, and so far things are going smoothly: the place is beautiful, work is great and food is excellent. Well OK, I do miss some things from Spain and weather is horrible, but from now on it can only get better, I suppose.

by fabian at December 22, 2009 05:36 PM

December 15, 2009

Fabian Pedregosa

Learning, Machine Learning

My new job is about managing an open source package for machine learning in Python. I’ve had some experience with Python now, but I am a total newbie in the field of machine learning, so my first task will be to find a good reference book in the subject and start reading.

The books I’ve come across so far have been:

The classic “Artificial Intelligence, A Modern Approach” by Rusell & Norvig, also known by its initials AIMA has been a very valuable introduction. It has four chapters devoted to Machine Learning, and while it does not go very deeply into the maths nor provide detailed analysis of the algorithm, it can be read fairly easily and has very interesting historical notes at the end of each chapter.

“Data Mining, Practical Machine Learning Tools and Techniques” by Ian H. Witten and Eibe Frank. It is easy going, a bit more in-depth but does not go very deep into the maths (and when it does, the section is marked as ‘optional’). The second part of the book, called “The Weka Machine Learning Workbench” describes Weka, a machine learning framework written in Java. For me, this is an invaluable resource, as Weka seems to be fairly well designed and will surely provide some design inspiration.

“Pattern Classification”, by Richard O Duda, Peter E Hart, David G Stork. A bit harsh at the beginning, this book is slowly becoming my favorite. Doesn’t hide the maths behind the algorithms, has complete sections on computational complexity, lots of exercises at the end of each chapter … I just wish it came with answers to (selected) exercises.

Interestingly, this books starts describing the most specific methods and finishes with the most general ones. In AIMA, you’ll find the reverse scheme.

by fabian at December 15, 2009 09:34 PM

December 13, 2009

Ondřej Čertík

ESCO 2010 conference

An interesting conference 2nd European Seminar on Coupled Problems (ESCO) will be held on June 28 — July 2, 2010 in Pilsen, Czech Republic.
Among the topics are solving PDEs and applications and using Python for scientific computing. In particular, Gaël Varoquaux is the keynote speaker.

Unfortunately, it was later announced that the SciPy 2010 conference is going to be at the same time, which is really unfortunate. But here are some reasons why you should consider going to ESCO 2010 instead:
  • If you like numerical calculations (finite elements, differences, volumes, ...) and solving partial differential equations and other problems and also programming in Python, together with C/C++ or Fortran, you will have a chance to meet some of the top people in the field. SciPy conference usually has people who solve PDEs (e.g. SciPy 09 had about 6), but ESCO 2010 will have about 60. So ESCO wins.

  • Robert Cimrman, who you probably know from the scipy and numpy mailinglists, also the author of the sfepy FEM package in Python, lives in Pilsen, so he'll gladly show you some good Pilsen pubs. SciPy 2010 is going to be in Austin and while Austin has cool pubs too, I must be fair and I liked that (I was there at the Sage 08 days), but it's just not comparable, the beer is better in Pilsen, it's a historic city and there are more pubs.

  • Pilsen is close to Prague, so you will have the chance to visit it. You should walk in the old town, have couple beers etc. Here you can see some photos that Gaël took when we met in Prauge. Again, this is incomparable with Austin.

  • It is held in the Pilsner Urquell Brewery. When Pavel Šolín announced that at the SciPy 09 conference, Jarrod asked "Ah, in a beer pub?". So let me be clear. The word pilsner (type of the beer) is coming from the Czech city Pilsen (Plzeň in Czech). Pilsner Urquell is not some beer pub (e.g. even Reno where I live now has a beer pub), it's The Brewery. Austin is a cool place (and Texas steaks are really good), but as you can see now, it absolutely cannot compete with Pilsen.


If you have time, I can fully recommend to go to ESCO 2010.

by Ondřej Čertík (noreply@blogger.com) at December 13, 2009 08:49 PM

November 13, 2009

Aaron Meurer

How to get both 32-bit and 64-bit Python in Snow Leopard


We had some discussion on one of the Python issues about whether my Python in Snow Leopard should be 32-bit or 64-bit. I originally thought that it was tied to what the kernel was, but I turned out to be wrong.

From what I discovered, the important thing is what the Python was compiled as. You can tell what your Python has been compiled as by running:

>>> import sys
>>> from math import log
>>> log(sys.maxsize, 2)

If this is just under 31, then it is 32 bit. If it returns 63, then it is 64. An easier way to tell it to run:

>>> 2**40

If you get 1099511627776L, then you have 32-bit Python, if you get 1099511627776, you have 64-bit Python (notice that the number is long in 32-bit Python, because it is larger than maxint).

This test won’t work in Python 3 because all integers are “long” by default, but the first part will still work.

So why does this matter, you ask? Well, aside from the fact that much longer numbers are not long (anything less than 2**63 – 1 = 9223372036854775807), there is the issue of hashing.

In 64-bit Python:

>>> hash('a')
12416037344

but in 32-bit Python

>>> hash('a')
-468864544

SymPy uses hash values to order arguments, so often it happens that behavior in one architecture will not show up in the other. These problems are often hard to track and fix, but the worst is when things work fine on the machine you are working on. This actually happened to me with my GSoC project. I was renumbering the arbitrary constants in the printing order in an expression, but it turned out that the printing order of an expression can be dependent on .args order, so I had to modify the tests to canonize the numbering first.

So here comes the crux of the post. It turns out that on Mac OS X, if you install the binary from python.org (Mac Installer Disk Image), this installs a 32-bit Python (for compatibility reasons) in /Library/Frameworks/Python.framework/Versions/2.6/bin/python2.6
. However, if you install Python using 64-bit fink in Snow Leopard, it will compile it from source into 64-bit, and install it into /sw/bin/python2.6.

So now I have an easy way to test both architectures without having to ssh into some other machine, which was what I was doing before.

by asmeurer at November 13, 2009 01:16 AM

October 15, 2009

Dale Peterson

New rolling torus code

I changed the code to the rolling torus example pretty significantly. Little fixes in how the arrows representing the body fixed unit vectors are now fixed. Additionally, I added some code that calculates the kinetic, potential and total energy of the system, and plots it. The scipy solver seems to do a pretty good (not perfect) job of maintaining constant energy. I think tightening up the error tolerances would improve this, but the simulation is well behaved at this point so I didn’t bother.

by admin at October 15, 2009 09:32 PM

October 10, 2009

Dale Peterson

Gyrostat

I just added an example of use for the derivation of the equations of a gyrostat using PyDy.  It is located in the examples/gyrostat/ subdirectory.

In case you aren’t familiar, here is the definition of a gyrostat:

A gyrostat is a mechanical system which is comprised of more than one body and yet has the rigid body property that its inertia components are time independent constants.

The reason I became interested in gyrostats is because in modeling the bicycle, the rear frame and the rear wheel can be treated as a gyrostat, and in doing so, the number of parameters that appear in the equations of motion will be reduced by two.  The same can be done for the front fork and handlebar assembly and front wheel.

In addition to the Python script that generates the equations, I took the time to do a complete write up of the model in LaTeX and generate a nice pdf of it.

by admin at October 10, 2009 04:20 AM

September 10, 2009

Fredrik Johansson

Python floats and other unusual things spotted in mpmath

I've just put up a branch (named "mp4") containing changes to mpmath that I've been working on for several days now. This branch includes some rather large changes, including support for fixed-precision (machine-precision) arithmetic, and a new implementation of the multiprecision type. The code is a bit rough at the moment, and not all changes are final, so don't expect this in a release in the immediate future.

New mpf implementation


The mp4 branch uses a new mpf type that handles both real and complex values. It makes complex arithmetic quite a bit faster (up to twice the speed), although real arithmetic is a little slower. Some of the slowdown is due to wrapping old-format functions instead of providing new ones, so there is a bit of unnecessary overhead. The biggest advantage is not that the implementation becomes faster, but that it becomes simpler. For example, a nice feature is that complex calculations that yield real results don't need to be converted back to a real type manually. Unfortunately, writing code that also supports Python float and complex (see below) somewhat reduce the usefulness of such features.

Imaginary infinities are a little more well-behaved (i.e. j*(j*inf) gives -inf and not a complex NaN), and the representation is flexible enough that it might be possible to support unsigned infinities, infinities with an arbitrary complex sign, and possibly even zero with an arbitrary complex sign -- hey, perhaps infinities and zeros with order and residue information attached to them so you can evaluate gamma(-3)**2 / gamma(-5) / gamma(-2) directly? (Well, I'm not sure if this is actually possible, and if it is, it's probably way overkill :-)

This new implementation will not necessarily make it -- it depends on whether the pros outweigh the cons. I would very much appreciate review by others. Regardless of the outcome, writing it has been very useful for identifying and fixing problems in the interface between "low level" and "high level" code in mpmath.

This should particularly simplify the future support for a C or Cython-based backend. In fact, the current fixed-precision backend (see below) only uses about 200 lines of code -- and although it's not complete yet, it handles a large set of calculations. Adding a fast multiprecision backend based on Sage or GMPY can probably be done in an evening now. I'm not going to do such a thing right now, as I want to fix up the existing code before adding anything new; should someone else be interested though, then by all means go ahead.

I've also inserted a few other optimizations. As noted before on this blog, many of the special functions in mpmath are evaluated as linear combinations of hypergeometric series. All those functions basically fall back to a single routine which sums hypergeometric series. The mp4 branch contains a new implementation of this routine that is up to about twice as fast as before. The speed is achieved by code generation: for every different "type" (degree, parameter types) of hypergeometric series, a specialized version is generated with various bookkeeping done in advance so it doesn't have to be repeated in every iteration of the inner loop.

As an example, the following runs at 1.7x the speed:

>>> 1/timing(mpmath.hyp1f1, 1.5, 2.2, 1.35)
5620.8844813722862
>>> 1/timing(mp4.hyp1f1, 1.5, 2.2, 1.35)
9653.1737629459149


Nearly all tests pass with the new mpf type used instead of the old one -- the only significant missing piece is interval arithmetic.

Fixed-precision arithmetic



Several people have requested support for working with regular Python floats and complexes in mpmath. Often you only need a few digits of accuracy, and mpmath's multiprecision arithmetic is unnecessarily slow. Many (though not all) of the algorithms in mpmath work well in fixed precision; the main problem with supporting this feature has been that of providing an appropriate interface that avoids cluttering the code.

In the mp4 branch, this is solved by adding a fixed-precision context. This context uses the same "high-level" methods for special functions, calculus, linear algebra, etc, as the multiprecision context, while low-level functions are replaced with fixed-precision versions. For example exp is just a wrapper around math.exp and cmath.exp.

While the default multiprecision context instance is called mp, the default fixed-precision context instance is called fp. So:


>>> from mp4 import mp, fp
>>> mp.pretty = True
>>> mp.sqrt(-5); type(_)
2.23606797749979j
<class 'mp4.ctx_mp.mpf'>
>>> fp.sqrt(-5); type(_)
2.2360679774997898j
<type 'complex'>



The fixed-precision context is still missing a lot of low-level functions, so many things don't work yet. Let's try a couple of calculations that do work and see how they compare.

Lambert W function


The Lambert W function was the first nontrivial function I tried to get working, since it's very useful in fixed precision and its calculation only requires simple functions. The lambertw method is the same for mp as for fp; the latter just uses float or complex for the arithmetic.


>>> mp.lambertw(7.5)
1.56623095378239
>>> mp.lambertw(3+4j)
(1.28156180612378 + 0.533095222020971j)
>>> fp.lambertw(7.5)
1.5662309537823875
>>> fp.lambertw(3+4j)
(1.2815618061237759+0.53309522202097104j)

The fixed-precision results are very accurate, which is not surprising since the Lambert W function is implemented using a "self-correcting" convergent iteration. In fact, the multiprecision implementation could be sped up by using the fixed-precision version to generate the initial value. The speed difference is quite striking:

>>> 1/timing(mp.lambertw, 7.5)
1249.5319808144905
>>> 1/timing(mp.lambertw, 3+4j)
603.49697841726618
>>> 1/timing(fp.lambertw, 7.5)
36095.559380378654
>>> 1/timing(fp.lambertw, 3+4j)
20281.934235976787

Both the real and complex versions are about 30x faster in fixed precision.

Hurwitz zeta function


The Hurwitz zeta function is implemented mainly using Euler-Maclaurin summation. The main ingredients are Bernoulli numbers and the powers in the truncated L-series. Bernoulli numbers are cached, and can be computed relatively quickly to begin with, so they're not much to worry about. In fact, I implemented fixed-precision Bernoulli numbers by wrapping the arbitrary-precision routine for them, so they are available to full 53-bit accuracy. As it turns out, the fixed-precision evaluation achieves nearly full accuracy:


>>> mp.hurwitz(3.5, 2.25); fp.hurwitz(3.5, 2.25)
0.0890122424923889
0.089012242492388816
>>> t1=timing(mp.hurwitz,3.5,2.25); t2=timing(fp.hurwitz,3.5,2.25); 1/t1; 1/t2; t1/t2
473.66504799548278
7008.0267335004173
14.7953216374269


There is a nice 15x speedup, and it gets even better if we try complex values. Let's evaluate the the zeta function on the critical line 0.5+ti for increasing values of t:


>>> s=0.5+10j
>>> mp.hurwitz(s); fp.hurwitz(s)
(1.54489522029675 - 0.115336465271273j)
(1.5448952202967554-0.11533646527127067j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
213.37891598750548
4391.4815202596583
20.580672180923461

>>> s=0.5+1000j
>>> mp.hurwitz(s); fp.hurwitz(s)
(0.356334367194396 + 0.931997831232994j)
(0.35633436719476846+0.93199783123336344j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
8.1703453931669383
352.54843617352128
43.149759185011469

>>> s = 0.5+100000j
>>> mp.hurwitz(s); fp.hurwitz(s)
(1.07303201485775 + 5.7808485443635j)
(1.0730320148426455+5.7808485443942352j)
>>> t1=timing(mp.hurwitz,s); t2=timing(fp.hurwitz,s); 1/t1; 1/t2; t1/t2
0.17125208455924443
9.2754935956407891
54.162806949260492


It's definitely possible to go up to slightly larger heights still. The Euler-Maclaurin truncation in the Hurwitz zeta implementation is not really tuned, and certainly has not been tuned for fixed precision, so the speed can probably be improved.

Hypergeometric functions


Implementing hypergeometric series in fixed precision was trivial. That the multi-precision implementation is actually quite fast can be seen by comparing direct evaluations:


>>> mp.hyp1f1(2.5, 1.2, -4.5)
-0.0164674858506064
>>> fp.hyp1f1(2.5, 1.2, -4.5)
-0.016467485850591584
>>> 1/timing(mp.hyp1f1, 2.5, 1.2, -4.5)
6707.6667199744124
>>> 1/timing(fp.hyp1f1, 2.5, 1.2, -4.5)
12049.135305946565


The float version is only about twice as fast. Unfortunately, hypergeometric series suffer from catastrophic cancellation in fixed precision, as can be seen by trying a larger argument:


>>> mp.hyp1f1(2.5, 1.2, -30.5)
6.62762709628679e-5
>>> fp.hyp1f1(2.5, 1.2, -30.5)
-0.012819333651375751


Potentially, checks could be added so that the fixed-precision series raises an exception or falls back to arbitrary-precision arithmetic internally when catastrophic cancellation occurs. However, it turns out that this evaluation works for even larger arguments when the numerically stable asymptotic series kicks in:


>>> mp.hyp1f1(2.5, 1.2, -300.5)
1.79669541078302e-7
>>> fp.hyp1f1(2.5, 1.2, -300.5+0j)
1.7966954107830195e-07


The reason I used a complex argument is that the asymptotic series uses complex arithmetic internally, and Python has an annoying habit of raising exceptions when performing complex-valued operations involving floats (a proper workaround will have to be added). In this case the speedup is close to an order of magnitude:


>>> 1/timing(mp.hyp1f1, 2.5, 1.2, -300.5)
532.36666412814452
>>> 1/timing(fp.hyp1f1, 2.5, 1.2, -300.5+0j)
4629.4746136865342


Another example, a Bessel function. This function is calculated using a hypergeometric series and also a couple of additional factors, so the speedup is quite large:


>>> mp.besselk(2.5, 7.5)
0.000367862846522012
>>> fp.besselk(2.5, 7.5)
0.00036786284652201188
>>> 1/timing(mp.besselk, 2.5, 7.5)
3410.5578142787444
>>> 1/timing(fp.besselk, 2.5, 7.5)
21879.520083463747


The fixed-precision context does not yet handle cancellation of singularities in hypergeometric functions. This can be implemented as in the multiprecision case by perturbing values, although accuracy will be quite low (at best 5-6 digits; sometimes an accurate result will be impossible to obtain).

Numerical calculus


Root-finding works, as does numerical integration, with nice speedups:


>>> fp.quad(lambda x: fp.cos(x), [2, 5])
-1.8682217014888205
>>> mp.quad(lambda x: mp.cos(x), [2, 5])
-1.86822170148882
>>> 1/timing(fp.quad, lambda x: fp.cos(x), [2, 5])
1368.3622602114056
>>> 1/timing(mp.quad, lambda x: mp.cos(x), [2, 5])
55.480651962846274

>>> fp.findroot(lambda x: fp.cos(x), 2)
1.5707963267948966
>>> mp.findroot(lambda x: mp.cos(x), 2)
1.5707963267949
>>> 1/timing(fp.findroot, lambda x: fp.cos(x), 2)
19160.822293284604
>>> 1/timing(mp.findroot, lambda x: mp.cos(x), 2)
1039.4805452292442


The tests above use well-behaved object functions; some corner cases are likely fragile at this point. I also know, without having tried, that many other calculus functions utterly don't work in fixed precision (not by algorithm, nor by implementation). Some work will be needed to support them even partially. At minimum, several functions will have to be changed to use an epsilon of 10-5 or so since full 15-16-digit accuracy requires extra working precision which just isn't available.

Plotting


The plot methods work, and are of course way faster in fixed-precision mode. For your enjoyment, here is the gamma function in high resolution; the following image took only a few seconds to generate:

>>> fp.cplot(fp.gamma, points=400000)



With the default number of points (fp.cplot(fp.gamma)), it's of course instantaneous.

by Fredrik Johansson (fredrik.johansson@gmail.com) at September 10, 2009 11:14 PM

September 07, 2009

Aaron Meurer

Google Summer of Code 2009 Wrap Up


Sorry about the extreme delay with this. I of course have been busy with classes.

Note that this will just be a summary of the summer, with my comments looking back on it. If you want more details on each individual thing that I implemented, look back on my previous blog posts.

Let me start from the beginning. Around late February to early March of this year, I discovered the existence of Google Summer of Code. I knew that I wanted to do some kind of work this summer, preferably an internship, so it piqued my interest. At that time, the mentoring organizations were still applying for GSoC 2009, so I could only look at the ones from 2008. Most of them were either Linux things or Web things, neither of which I had any experience in or am I much interested in. I took a free course in Python at my University the previous semester, and it was the programming language that I knew best at the time. I had learned some Java in my first semester CS class (did I mention that this was my first year at college?), and I hated it, and I was still learning C for my second semester CS class. So I looked at what the Python Foundation had to offer. I am a double major in math and computer science, so I looked under the math/science heading. That’s when I saw SymPy.

I should not that I have been ahead in Math. It was my second semester, and I was taking Discrete Mathematics, Ordinary Differential Equations, Basic Concepts of Math, and Vector Analysis. So I looked for project ideas on the SymPy page that related to what I knew. The only one that I saw, other than core improvements, was to improve the ODE solving capabilities. I got into contact with the community and looked at the source, finding that it was only capable of solving 1st order linear equations and some special cases of 2nd order linear homogeneous equations with constant coefficients. I already at that point knew several methods from my ODE course, and I knew much of what I would learn.

Application Period
The most difficult part of the Google Summer of Code, in my opinion, is the application period. For starters, you have to do it while you are still in classes, so you pretty much have to do it in your free time. Also, if you have never applied for Google Summer of Code before, you do not really know what a good application should look like. I have long had my application available on the SymPy Wiki, and I will reference it here a few times. First off, it was recommended to me by some of the SymPy developers that I put as many potential things that I could do in the summer in my application as I though I could do. I was still only about half way through my ODEs course when I wrote the application, but I had the syllabus so I knew the methods I would be learning at least by name. So that is exactly what I did: I packed my application with every possible thing that I knew we would be learning about in ODEs.

After I felt that I had a strong application, and Ondrej had proofread it for me, I submitted it. There were actually two identical applications, one for the Python Software Foundation, and one for Portland State University. This is because SymPy was not accepted as a mentoring organization directly, so they had to use those two foundations as proxies. A requirement of acceptance is to submit a patch that passes review. I decided to add a Bernoulli solver, because Bernoulli can be solved in the general case much like the 1st order linear ODE, which was already implemented.

After I applied, there was an acceptance period. I used that period to become aquatinted with the SymPy community and code base. A good way to do this is to try to fix EasyToFix issues. I found issue 694, which is to implement a bunch of tests from a paper by Michael Wester for testing computer algebra systems. The tests cover every possible thing that a full featured CAS could do, so it was a great way to learn SymPy. The issue is still unfinished, so working on it is still a good way to learn how to use SymPy.

Also, it was important to learn git, SymPy’s version control system. The learning curve it pretty steep if you have never used version control system before, but once you can use it, it becomes a great tool at your disposal.

Acceptance
After being accepted, I toned down my work with SymPy to work on finishing up my classes. My classes finished a few weeks before the official start, so I used that period to get a jump start on my project.

The GSoC Period
For the start of the period, I followed my timeline. I implemented 1st order ODEs with homogeneous coefficients and 1st order exact ODEs. These were both pretty simple to do, as I expected.

The next thing I wanted to do was separable. My goal at this point was to get every relevant exercise from my textbook to work with my solvers. One of the exercises from my book (Pg. 56, No. 21) was dy=e^{x + y}dx. I soon discovered that it was impossible for SymPy to separate e^{x + y} \rightarrow e^{x}e^{y}, because the second would be automatically combined in the core. I also discovered that expand(), which should have been the function to split that, expanded using all possible methods indiscriminately. Part of my separatevars() function that I was writing to separate variables in expressions would be to split things like x + yx \rightarrow x(1 + y) and 2 x y + x^{2} + y^{2} \rightarrow (x + y)^{2}, but expand()
as it was currently written would expand those.

So I spent a few weeks hacking on the core to make it not auto-combine exponents. I came up with a rule that exponents would only auto-combine if they had the same term minus the coefficient, the same rule that Add uses to determine what terms should auto combined by addition. So it would combine e^{x}e^{x} \rightarrow e^{2x}, but e^{x}e^{y} would be left alone. It turns out that some of our algorithms, namely the Gruntz limit algorithm, relies on auto-combining. We already had a function that could combine exponents, powsimp(), but it also combined bases, as in x^zy^z \rightarrow (xy)^z, so I had to split the behavior so that it could act only as auto-combining once did (by the way, use powsimp(expr, combine='exp', deep=True) to do this). Then, after some help from Ondrej on pinpointing the exact location of the bugs, I just applied the function there. The last thing I did here was to split the behavior of expand, so that you could do expand(x*(y + 1), mul=False) and it would leave it alone, but expand(exp(x + y), mul=False) would return exp(x)*exp(y). This split behavior turned out to be useful in more than one place in my code.

This was the first non bug fix patch of mine that was pushed in to SymPy, and at the time of this writing, it is the last major one in the latest stable version. It took some major rebasing to get my convoluted commit history ready for submission, and it was during this phase that I git finally clicked for me, especial the git rebase command. This work took a few weeks from my ODEs time, and it became clear that I would not be doing every possible thing from my application. The reason that I included so much in my application was that my project was non-atomic. I could implement a little or a lot and still have a working useful module.

If you look at my timeline on my application, you can see that the first half is symbolic methods, and the second half is other methods, things like series. It turns out that we didn’t really learn much about systems of ODEs in my course and we learned very little about numerical methods (and it would take much more to know how to implement them). We did learn series methods, but they were so annoying to do that I came to hate them with a passion. So I decided to just focus on symbolic methods, which were my favorite anyway. My goal was to implement as many as I could.

After I finished up separable equations, I came up with an idea that I did not have during the application period. dsolve() was becoming cluttered fast with all of my solution methods. The way that it worked was that it took an ODE and it tried to match methods one by one until it found one that worked, which it then used. This had some drawbacks. First, as I mentioned, the code was very cluttered. Second, the ODEs methods would have to be applied in a predetermined order. There are several ODEs that match more than one method. For example, 2xy + (x^2 + y^2)\frac{dy}{dx}=0 has coefficients that are both homogeneous of order 2, and is also exact, so it can be solved by either method. The two solvers return differently formatted solutions for each one. A simpler example is that 1st order ODEs with homogeneous coefficients can be solved in two different ways. My working solution was to try them both and then apply some heuristics to return the simplest one. But sometimes, one way would create an impossible integral that would hand the integration engine. And it made debugging the two solvers more difficult because I had to override my heuristic. This also touches on the third point. Sometimes the solution to an ODE can only be represented in the form of an unevaluatable integral. SymPy’s integrate() function is supposed to return an unevaluated Integral class if it cannot do it, but all too often it will just hang forever.

The solution I came up with was to rewrite dsolve using a hints method. I would create a new function called classify_ode() that would do all of the ODE classification, removing it from the solving code. By default, dsolve would still use a predetermined order of matching methods. But you could override it by passing a “hint” to dsolve for any matching method, and it would apply that method. There would also be options to only return unevaluated integrals when applicable.

I ended up doing this and more (see the docstrings for classify_ode() and dsolve() in the current git master branch), but before I could I needed to clean up some things. I needed to rewrite all of dsolve() and related functions. Before I started the program, there were some special cases in dsolve for second order linear homogeneous ODEs with constant coefficients and one very special case ODE for the expanded form of \frac{d^2}{dx^2}(xe^{-y}) = 0.

So the first thing I did was implement a solver for the general homogeneous linear with constant coefficients case. These are rather simple to do: you just find the roots of the characteristic polynomial built off of the coefficients, and then put the real parts of the roots in front of the argument of an exponential and the imaginary parts in front of the arguments of a sine and cosine (for example, 3 \pm 2i would give C1e^{3x}\sin{2x} + C2e^{3x}\cos{2x}. The thing was, that if the imaginary part is 0, then you only have 1 arbitrary constant on the exponential, but if it is non-zero, you get 2, one for each trig function. The rest falls out nicely if you plug 0 in for b into $e^{ax}(C1\sin{bx} + C2\cos{box})$ because the sine goes to 0 and the cosine becomes 1. But you would end up with C1 + C2 instead of just C1 in that case. I had already planned on doing arbitrary constant simplification as part of my project, so I figured I would put this on hold and do that first. Then, once that was done, the homogeneous case would be reduced to 1 case instead of the usual 2 or 3.

My original plan was to make an arbitrary constant type that automatically simplified itself. So, for example, if you entered C1 + 2 + x with C1 an arbitrary constant, it would reduce to just C1 + x. I worked with Ondrej, including visiting him in Los Alamos, and we build up a class that worked. The problem was that, in order to have auto-simplification, I had to write the simplification directly into the core. Neither of us liked this, so we worked a little bit on a basic core that would allow auto-simplification to be written directly in the classes instead of in the Mul.flatten() and Add.flatten() methods. It turns out that my constant class isn’t the only thing that would benefit from this. Things like the order class (O(x)) and the infinity class (oo) are auto-simplified in the core, and things could be much cleaner if they happened in the classes themselves. Unfortunately, modifying the core like this is not something that can happen overnight or even in a few weeks. For one thing, it needed to wait until we had the new assumptions system, which was another Google Summer of Code project running parallel to my own. So we decided to shelf the idea.

I still wanted constant simplification, so I settled with writing a function that could do it instead. There were some downsides to this. Making the function as general as the classes might have been would have been far too much work, so I settled on making it an internal-only function that only worked on symbols named C1, C2, etc. Also, unlike writing the simplification straight into Mul.flatten() which was as simple as removing any terms that were not dependent on x, writing a function that parsed an expression and simplified it was considerably harder to write. I managed to churn out something that worked, and so I was ready to finish up the solver I had started a few paragraphs ago.

After I finished that, I still needed to maintain the ability to solve that special case ODE. Apparently, it is an ODE that you would get somewhere in deriving something about relativity, because it was in the relativity.py example file. I used Maple’s excellent odeanalyser() function (this is where I go the idea for my classify_ode())to find a simple general case ODE that it fit (Liouville ODE). After I finished this, I was ready to start working on the hints engine.

It took me about a week to move all classification code into classify_ode(), move all solvers into individual functions, separate simplification code into yet other functions, and tie it all together in dsolve(). In the end, the model worked very well. The modularization allowed me to do some other things that I had not considered, such as creating a special “best” hint that used some heuristics that I originally developed for first order homogeneous which always has two possible solutions to try to give the best formatted solution for any ODE that has more than one possible solution method. It also made debugging individual methods much easier, because I could just use the built in hint calls in dsolve() instead of commenting out lines of code in the source.

This was good, because there was one more method that I wanted to implement. I wanted to be able to solve the inhomogeneous case of a nth order linear ode with constant coefficients. This can be done in the general case using the method of variation of parameters. It was quite simple to set up variation of parameters up in the code. You only have to set up a system of integrals using the Wronskian of the general solutions. It would usually be a very poor choice of a method if you were trying to solve an ODE by hand because taking the Wronskian and computing n integrals is a lot of work. But for a CAS, the work is already there. I just have to set up the integrals.

It soon became clear that even though, in theory, the method of variation of parameters can solve any ODE of this type, in practice, it does not always work so well in SymPy. This is because SymPy have very poor simplification, especially trigonometric simplification, so sometimes there would be a trigonometric Wronskian that would be identically equal to some constant, but it could only simplify it to some very large rational function of sines and cosines. When these were passed to the integral engine, it would cause it to fail, because it could not find the integral for such a seemingly complex expression.

In addition, taking Wronskians, simplifying them, and then taking n integrals is a lot of work as I said, and even when SymPy could do it, it took a long time. There is another method for solving these types of equations called undetermined coefficients that does not require integration. It only works on a class of ODEs where the right hand side of the ODE is a simple combination of sines, cosines, exponentials, and polynomials in x. It turns out that these kinds of functions are common anyway, so most ODEs of this type that you would encounter could be solved with this method. Unlike variation of parameters, undetermined coefficients requires considerable setup, including checking for different cases. This would be the method that you would want to use if you had to solve the ODE by hand because, even with all the setup, it only requires solving a system of linear equations vs. solving n integrals with variation of parameters, but for a CAS, it is the setup that matters, so this was a difficult prospect.

I spent the last couple of weeks writing up the necessary algorithms to setup the required system of linear equations and handling the different cases. After I finally worked out all of the bugs, I ran some profiling against my variation of parameters solver. It turned out that for ODEs that had trigonometric solutions (which take longer to simplify), my undetermined coefficients solver was an order of magnitude faster than the variation of parameters solver (and that is just for the ODEs that the variation of parameters engine could even solve at all). For ODEs that only had exponentials, it was still 2-4 times faster.

I finished off the summer by writing extensive documentation for all of my solvers and functions. Hopefully someone who uses SymPy to solve ODEs can learn something about ODE solving methods as well as how to use the function I wrote when they read my documentation.

Post-GSoC
I plan on continuing development with SymPy now that the Google Summer of Code period is over. SymPy is an amazing project, mixing Python and Computer Algebra, and I want to help it grow. I may even apply again in a future year to implement some other thing in SymPy, or maybe apply as a mentor for SymPy to help someone else improve it.

Advice
What follows is some general advice for someone who wants to apply for Google Summer of Code. Some of the advice pertains specifically to SymPy, and some of it is general advice that I think would apply to any project.

- Get involved early. As soon as you decide that you want to participate in Google Summer of Code, start getting involved in the project. Get into contact with them and discuss possible projects. If you are looking before the participating organizations are announced, look at the organizations from previous years. For some organizations, it will vary; for others (like Python), it is almost given that they will be accepted every year.

- Some projects (including SymPy) require you to send in a patch that passes review to be accepted. This will give you a change to start familiarizing yourself with the code base. If you are applying to SymPy, the Wester example I mentioned above is a really good way to learn what SymPy can do and how it works.

- Subscribe to the mailing list, and once you are comfortable with it, participate. Also, it is a good idea to idle in IRC (SymPy is on freenode at #sympy). This will help you get to know the main contributors for the project.

- For you application, see if the people in the project you are applying for will review it. If they like your project idea, they will try to help you write a good application so you can be accepted and you can implement it. If they don’t like your idea, then they will tell you and you should change it, otherwise you will not be accepted, no matter how well written your proposal is. I have my proposal on the wiki (see link above). I am not saying that it is necessarily a very good proposal, but it did get accepted. If you are applying to SymPy, Ondrej will proofread your applications for you.

- If you are an IRC fan, there is also #gsoc on freenode, where you can ask all your GSoC related questions. Be warned that it does get pretty noisy in the application period, especially right before the applications are due and right before proposals are accepted.

- I cannot stress this one enough. If you have never worked with a version control system before, it is perhaps more important to spend your time learning it than it is to learn the code base for your project. These things have a steep learning curve if you have never used them before. Once you master them though, they can make your life much easier. Also, the sooner you learn to use them well, the easier your life will be later on down the road. I spent a good part of the last week of GSoC cleaning up my commit history from the first half of the summer when I bad very poor committing/log habits. If your project uses git, such as SymPy does, you might look at this tutorial. If it uses something else, good luck. Seriously, git is the only good version control system. See this video.

- Expect to spend only about half of the summer actually implementing stuff. You may think that you are a good programmer and that your code will not be so buggy that you will need to spend that much time fixing bugs, and you may be right. But the fact is, you will be working on code bases written by may programmers that are not so good. You will need to fix several already existing bugs to make your code work, which means that you will need to learn the code base well, learn how to read other people’s code, and how to fix bugs that you had no part in creating. You will be glad if a bug is in your code because you will usually know immediately what causes it and how to fix it. But if a bug is somewhere else, you will need to find it, figure out why it happens, what is supposed to happen, and how to fix it without breaking anything else. This is also why it is important to be active in the developer community.

- Good luck.

by asmeurer at September 07, 2009 06:49 AM

September 06, 2009

Jason Gedge

Components, systems, subsystems, entities...*collapses*

So I've recently been reading into various articles and forum topics on Component-Based design, or Entity Systems (ES). This concept was probably first used in games when the first Dungeon Siege came out. Apparently it's closest neighbour is possibly that of Aspect-Oriented programming, which I know zero about. At first it seems like a complicated idea, but after I looked through a bunch of info

by Jason G (noreply@blogger.com) at September 06, 2009 12:31 AM

September 05, 2009

Fabian Pedregosa

Summer of Code is over

Google Summer of Code program is officially over. It has been four months of intense work, exciting benchmarks and patch reviewing. It was a huge pleasure working with you guys!

As for the project, I implemented a complete logic module and then an assumption system for sympy (sympy.logic, sympy.assumptions, sympy.queries). I even had time to make the logic module fast. On top of this, there’s the refine module. It is there where you can see some nice examples and where all the power of sympy.queries and sympy.logic is exposed.

Although this sounds good, there are some things that I did not complete on time. I could not remove the old assumption system. There are simply too many things that depend on this to remove it on one move. However, I agreed with Ondrej that we both would be working on this the days 15-30 September. This has to be done because we definitely do not want to make a sympy release with two different assumption systems!

PD: a more detailed report lives here

by fabian at September 05, 2009 10:21 AM

August 27, 2009

Jason Gedge

Pointgrey Cameras

For anyone who doesn't know my research, I'm looking into stereo vision algorithms in an underwater camera array. What I'm hoping to do is a rough scene reconstruction that has improved results over blindly using an existing stereo algorithm.Now on to what this post is about. I found a little subtle and, as far as I can tell, undocumented feature of the Pointgrey cameras. If you use the

by Jason G (noreply@blogger.com) at August 27, 2009 12:57 PM

August 25, 2009

Ondřej Čertík

SciPy 2009 Conference

I attended the SciPy 2009 conference last week at Caltech.

When I first presented about SymPy at the scipy 07 conference exactly 2 years ago, it was just something that we started, so noone really used that. Last week, there were already 6 other people at the conference who contributed one or more patches to SymPy: Robert Kern, Andrew Straw, Pauli Virtanen, Brian Granger, Bill Flynn, Luke Peterson.

I gave a SymPy tutorial, main presentation and Luke gave a PyDy + SymPy lightning talk (seek to 16:04).

I also gave my experience with designing a traits GUI for FEM lightning talk (seek to 5:35).

My advisor Pavel Solin gave a talk about Hermes and FEMhub and other things that we do in our group in Reno.

Besides that, it was awesome to meet all the guys from the scientific Python community, meet old friends and also get to know other people that I only knew from the lists. I was pleased to meet there people who solve PDE using Python, we had many fruitful discussions together, and as a result, I already created FEMhub spkg packages for FiPy, other will follow. Our aim is to create a nice interface to all of them, so that we can easily test a problem in any PDE package and see how it performs.

Overall, my impression is very positive. I am very glad I chose Python as my main language many years ago, all the essential parts are now getting there, e.g. numerics (numpy, scipy, ...), symbolic (sympy, Sage, ...), 2d plotting (matplotlib, ...), 3d plotting (mayavi), GUI (traits UI, which supports both GTK, QT on linux and native widgets on Mac and Windows, and Sage notebook for web), excellent documentation tool with equations support (Sphinx), lots of supporting libraries, like sparse solvers and then very easy way to wrap C/C++ code using Cython and to speed up critical parts of the Python code using Cython. It's not that each of those libraries is the best in the world --- in fact, not a single one is --- but together as an ecosystem, plus the high level of (free) support on the lists for all of those libraries, this in my opinion makes Python a number one choice for scientific computing, together with C, C++ and sometimes Fortran for CPU intensive tasks and/or legacy libraries.

by Ondřej Čertík (noreply@blogger.com) at August 25, 2009 12:04 AM

August 24, 2009

Dale Peterson

New examples of use with PyDy

In the last several days I’ve made a lot of updates to tools for using PyDy and have come very close to settling on the way that PyDy is used to derive the equations of motion.  I have added and updated several examples, including the double pendulum and a rigid body with two reaction wheels used for attitude control.  Animations for each example have been added using visual python, so visualization of the dynamics is very tangible.

There are two main things that I am still trying to iron out with PyDy.  The first is how to handle ignorable coordinates and how to allow for the user to control whether or not they are included in the output equations of motion.  For example, for the symmetric rolling disc, there are 4 ignorable coordinates:  two for the location in the ground plane, one for the heading and one for the spin.  For purposes of stability analysis, the kinematic differential equations associated with these coordinates are not needed.  However, for purposes of animation, these equations are needed.  My goal is to make it easy to clearly specify whether or not these equations are desired in the output equations.

The second major thing that still needs work is handling dependent generalized speeds in nonholonomic systems.  When the motion equations are generated, they will involve these dependent generalized speeds, and their time derivatives, but these quantities can be computed from the constraint equations and therefore can be left implicit in the final equations of motion.  The derivatives of the dependent generalized can be intelligently computed by some careful formations of the gradients of the components of the matrix that relates the dependent speeds to the independent ones.

I will be working on both of these tasks this week and will be writing examples for the Whipple bicycle model in addition to a spinning top and the rattleback.

by admin at August 24, 2009 05:11 AM

August 23, 2009

Ondřej Čertík

SymPy on the google phone

Here is a proof that sympy works on the google phone:



Thanks to Rae S. Yip from Caltech and Nicolas Pinto for taking the picture.

by Ondřej Čertík (noreply@blogger.com) at August 23, 2009 12:51 AM

August 19, 2009

Fabian Pedregosa

Speed improvements for ask() (sympy.queries.ask)

I managed to overcome the overhead in ask() that arises when converting between symbol and integer representation of sentences in conjunctive normal.

The result went beyond what I expected. The test suite for the query module got 10x times faster in my laptop. From 26 seconds, it descended to an impressive 2.03 secs. There is still room for improvement, but it is no longer “so desperately slow”.

I’ll submit those patches soon to sympy’s trunk, but in the meantime they are in my logic branch:

git pull http://fseoane.net/git/sympy.git logic

by fabian at August 19, 2009 10:36 PM

Jason Gedge

First Look

So what's this score editing application I've been talking about in my posts? We've decided to give a pre-beta teaser. Minus a couple of rendering issues, here's what she looked like a couple of days ago (click image for larger preview):For our first release we focused on guitar tabs, and making sure the application as a whole works together properly. Unfortunately this release will be a private

by Jason G (noreply@blogger.com) at August 19, 2009 11:55 AM

August 18, 2009

Fabian Pedregosa

Logic module (sympy.logic): improving speed

Today I’ve been doing some speed improvements for the logic module. More precisely, I implemented an efficient internal representation for clauses in conjunctive normal form.

In practice this means a huge performance boost for all problems that make use the function satisfiable() or dpll_satisfiable(). For example, test_dimacs.py has moved from 2.7 seconds to an impressive 0.3 sec, and ask() runs on average 3x times faster, although both problems still have an overhead because of the conversion to this new representation that can be avoided in most times.

Now, the details. Traditionally, dpll (the algorithm that we use for deciding satisfiability) used to store clauses as arrays of symbols, and this worked fine, but sadly comparing symbols in sympy is slow, and this algorithm does a lot of comparisons … but we can map each sympy symbol to a unique integer, and with minor modifications to the algorithm we get these performance gains.

Now, the code. You can pull from my branch logic:

git pull http://fseoane.net/git/sympy.git logic

There are now some obvious performance tweaks we can do:

- in ask(), we can skip the conversion to integer representation by ‘precompiling’ known_facts_dict into this representation. This should be easy and will probably give performance boosts of several orders of magnitude.

- this integer representation is very similar to the one used in dimacs CNF files, so a parser that directly converts CNF files to this integer representation should make solving CNF files much faster.

I would like to give some credit to Ronan Lamy, who sent a patch some time ago, and although I did not include it (yet) into main sympy branch, it inspired me for these modifications.

by fabian at August 18, 2009 09:35 PM

August 17, 2009

Aaron Meurer

Undetermined Coefficients


[Sorry for the delay in this post. I was having some difficulties coming up with some of the rationales below. Also, classes have started, which has made me very busy.]

If there was one ODE solving method that I did not want to implement this summer, it was undetermined coefficients. I didn’t really like the method too much when we did it my my ODE class (though it was not as unenjoyable as series methods). The thing that I never really understood very well is to what extent you have to multiply terms in the trial function by powers of x to make them linearly independent of the solution to the general equation. We did our ODEs homework in Maple, so I would usually just keep trying higher powers of x until I got a solution. But to implement it in SymPy, I had to have a much better understanding of the exact rules for it.

From a user’s point of view, the method of undetermined coefficients is much better than the method of variation of parameters. While it is true that variation of parameters is a general method and undetermined coefficients only works on a special class of functions, undetermined coefficients requires no integration or advanced simplification, so it is fast (very fast, as well shall see below). All that the CAS has to do is figure out what a trial function looks like, plug it into the ODE, and solve for the coefficients, which is a system of linear equations.

On the other hand, from the programmer’s point of view, variation of parameters is much better. All you have to do is take the Wronskian of the general solution set and use it to set up some integrals. But the Wronskian has to be simplified, and if the general solution contains sin’s and cos’s, this requires trigonometric simplification not currently available in SymPy (although it looks like the new Polys module will be making a big leap forward in this area). Also, integration is slow, and in SymPy, it often fails (hangs forever).

Figuring out what the trial function should be for undetermined coefficients is way more difficult to program, but having finnally finished it, I can say that it is definitely worth having in the module. Problems that it can solve can run orders of magnitude faster than the variation of parameters, and often variation of parameters can’t do the integral or returns a less simplified result.

So what is this undetermined coefficients? Well, the idea is this: if you knew what each linearly independent term of the particular solution was, minus the coefficients, then you could just set each coefficient as an unknown, plug it into the ODE, and solve for them. It turns out that resulting system of equations is linear, so if you do the first part right, you can always get a solution.

The key thing here is that you know what form the particular solution will take. However, you don’t really know this ahead of time. All you have is the linear ode a_ny^{(n)}(x) + \dots + a_1y'(x) + a_0y(x) = F(x) (as far as I can tell, this only works in the case where the coefficients a_i are constant with respect to x. I’d be interested to learn that it works for other linear ODEs. At any rate, that is the only one that works in my branch right now.). The solution to the ode is y(x) = y_g(x) + y_p(x), where y_g(x) is the solution to the homogeneous equation f(x) \equiv 0, and y_p(x) is the particular solution that produces the F(x) term on the right hand side. The key here is just that. If you plug y_p(x) into the left hand side of the ode, you get F(x).

It turns out that this method only works if the function F(x) only has a finite number of linearly independent derivatives (I am unsure, but this might be able to work in other cases, but it would involve much more advanced mathematics). So what kind of functions have a finite number of linearly independent solutions? Obviously, polynomials do. So does e^x, \cos{x}, and \sin{x}. Also, if we multiply two or more of these types together, then we will get a finite number of linearly independent solutions after applying the product rule. But is that all? Well, if we take the definition of linear independence from linear algebra, we know that a set of n vectors \{\boldsymbol{v_1}, \boldsymbol{v_2}, \boldsymbol{v_3}, \dots, \boldsymbol{v_n}\}, not all zero, are linearly independent only if a_1\boldsymbol{v_1} + a_2\boldsymbol{v_2} + a_3\boldsymbol{v_3} + \dots + a_n\boldsymbol{v_n}=0 holds only when a_1 \equiv 0, a_2 \equiv 0, a_3 \equiv 0, \dots, a_n \equiv 0, that is, the only solution is the trivial one (remember, this is the definition of linear independence). They are linearly dependent if there exist weights a_1, a_2, a_3, \dots, a_n, not all 0, such that the equation a_1\boldsymbol{v_1} + a_2\boldsymbol{v_2} + a_3\boldsymbol{v_3} + \dots + a_n\boldsymbol{v_n}=0 is satisfied. Using this definition, we can see that a function f(x) will have a finite number of linearly independent derivatives if it satisfies a_nf^{(n)}(x) + a_{n - 1}f^{(n - 1)}(x) + \dots + a_1f'(x) + a_0f(x) = 0 for some n and with a_i\neq 0 for some i. But this is just a homogeneous linear ODE with constant coefficients, which we know how to solve. The solutions are all of the form ax^ne^{b x}\cos{cx} or ax^ne^{b x}\sin{cx}, where a, b, and c are real numbers and n is a non-negative integer. We can set the various constants to 0 to get the type we want. For example, for a polynomial term, b will be 0 and c will be 0 (use the cos term).

So this gives us the exact form of functions that we need to look for to apply undetermined coefficients, based on the assumption that it only works on functions with a finite number of linearly independent derivatives.

Well, implementing it was quite difficult. For every ODE, the first step in implementation is matching the ODE, so the solver can know what methods it can apply to a given ODE. To match in this case, I had to write a function that determined if the function matched the form given above, which was not too difficult, though not as trivial as just grabbing the right hand side in variation of parameters. The next step is to use the matching to format the ODE for the solver. In this case, it means finding all of the finite linearly independent derivatives of the ODE, so that the solver can just create a linear combination of them solve for the coefficients. This was a little more difficult, and it took some lateral thinking.

At this point, there is one more thing that needs to be noted. Since the trial functions, that is, the linearly independent derivative terms of the right hand side of the ODE, are of the same form as the solutions to the homogeneous equation, it is possible that one of the trial function terms will be a solution to the homogeneous equation. If this happens, plugging it into the ODE will cause it to go to zero, which means that we will not be able to solve for a coefficient for that term. Indeed, that term will be of the form C1*\textrm{term} in the final solution, so even if we had a coefficient for it, it would be absorbed into this term from the solution to the homogeneous equation. For example, variation of parameters will give a coefficient for such terms, even though it is unnecessary. This is a clue that Maple uses variation of parameters for all linear constant coefficient ODE solving, because it gives the unnecessary terms with the coefficients that would be given by variation of parameters, instead of absorbing them into the arbitrary constants.

We can safely ignore these terms for undetermined coefficients, because their coefficients will not even appear in the system of linear equations of the coefficients anyway. But, without these coefficients, we will run into trouble. It turns out that if a term x^ne^{ax}\sin{bx} or x^ne^{ax}\cos{bx} is repeated solution to the homogeneous equation, and x^{n + 1}e^{ax}\sin{bx} or x^{n + 1}e^{ax}\cos{bx} is not, so that n is the highest x power that makes it a solution to the homogeneous equation, and if the trial solution has x^me^{ax}\sin{bx} or x^me^{ax}\cos{bx} terms, but not x^{m + 1}e^{ax}\sin{bx} or x^{m + 1}e^{ax}\cos{bx} terms, so that m is the highest power of x in the the trial function terms, then we need to multiply these trial function terms by x^{n + m} to make them linearly independent with the solutions of the homogeneous equation.

Most references simply say that you need to multiply the trial function terms by “sufficient powers of x” to make them linearly independent with the homogeneous solution. Well, this is just fine if you are doing it by hand or you are creating the trial function manually in Maple and plugging it in and solving for the coefficients. You can just keep upping the powers of x until you get a solution for the coefficients. Creating those trial functions in Maple, plugging them into the ODE, and solving for the coefficients is exactly what I had to do for my homework when I took ODEs last spring, and this “upping powers” trial and error method is exactly the method I used. But when you are doing it in SymPy, you need to know exactly what power to multiply it by. If it is too low, you will not get solution to the coefficients. If it is too high, you can actually end up with too many terms in the final solution, giving a wrong answer.

Fortunately, my excellent ODEs textbook gives the exact cases to follow, and so I was able to implement it correctly. The textbook also gives a whole slew of exercises, all for which the solutions are given. As usual, this helped me to find the bugs in my very complex and difficult to write routine. It also helped me to find a match bug that would have prevented dsolve() from being able to match certain types of ODEs. The bug turned out to be fundamental to the way match() is written, so I had to write my own custom matching function for linear ODEs.

The final step in solving the undetermined coefficients is of course just creating a linear combination of the trial function terms, plugging it into the original ODE, and setting the coefficients of each term on each side equal to each other, which gives a linear system. SymPy can solve these easily, and once you have the values of the coefficients, you can use them to build your particular solution, at which point, you are done.

The results were astounding. Variation of parameters would hang on many simple inhomogeneous ODEs because of poor trig simplification of the Wronsikan, but my undetermined coefficients method handles them perfectly. Also, there is no need to worry about absorbing superfluous terms into the arbitrary constants as with variation of parameters, because they are removed from within the undetermined coefficients algorithm.

But the biggest thing was speed. Here are some benchmarks on some random ODEs from the test suite. WordPress code blocks are impervious to whitespace, as I have mentioned before, so no pretty printing here. Also, it truncates the hints. The hints used are 'nth_linear_constant_coeff_undetermined_coefficients' and 'nth_linear_constant_coeff_variation_of_parameters':

In [1]: time dsolve(f(x).diff(x, 2) - 3*f(x).diff(x) - 2*exp(2*x)*sin(x), f(x), hint='nth_linear_constant_coeff_undetermined_coefficients')
CPU times: user 0.07 s, sys: 0.00 s, total: 0.08 s
Wall time: 0.08 s
Out[2]:
f(x) == C1 + (-3*sin(x)/5 - cos(x)/5)*exp(2*x) + C2*exp(3*x)

In [3]: time dsolve(f(x).diff(x, 2) - 3*f(x).diff(x) - 2*exp(2*x)*sin(x), f(x), hint='nth_linear_constant_coeff_variation_of_parameters')
CPU times: user 0.92 s, sys: 0.01 s, total: 0.93 s
Wall time: 0.94 s
Out[4]:
f(x) == C1 + (-3*sin(x)/5 - cos(x)/5)*exp(2*x) + C2*exp(3*x)

In [5]: time dsolve(f(x).diff(x, 4) - 2*f(x).diff(x, 2) + f(x) - x + sin(x), f(x), hint='nth_linear_constant_coeff_undetermined_coefficients')
CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s
Wall time: 0.06 s
Out[6]:
f(x) == x - sin(x)/4 + (C1 + C2*x)*exp(x) + (C3 + C4*x)*exp(-x)

In [7]: time dsolve(f(x).diff(x, 4) - 2*f(x).diff(x, 2) + f(x) - x + sin(x), f(x), hint='nth_linear_constant_coeff_variation_of_parameters')
CPU times: user 5.43 s, sys: 0.03 s, total: 5.46 s
Wall time: 5.52 s
Out[8]:
f(x) == x - sin(x)/4 + (C1 + C2*x)*exp(x) + (C3 + C4*x)*exp(-x)

In [9]: time dsolve(f(x).diff(x, 5) + 2*f(x).diff(x, 3) + f(x).diff(x) - 2*x - sin(x) - cos(x), f(x), 'nth_linear_constant_coeff_undetermined_coefficients')
CPU times: user 0.10 s, sys: 0.00 s, total: 0.10 s
Wall time: 0.11 s
Out[10]:
f(x) == C1 + (C2 + C3*x - x**2/8)*sin(x) + (C4 + C5*x + x**2/8)*cos(x) + x**2

In [11]: time dsolve(f(x).diff(x, 5) + 2*f(x).diff(x, 3) + f(x).diff(x) - 2*x - sin(x) - cos(x), f(x), 'nth_linear_constant_coeff_variation_of_parameters')


The last one involves a particularly difficult Wronskian for SymPy (run it with hint=’nth_linear_constant_coeff_variation_of_parameters_Integral’, simplify=False).

Wall time comparisons reveal amazing speed differences. We’re talking orders of magnitude.

In [13]: 0.94/0.08
Out[13]: 11.75

In [14]: 5.52/0.06
Out[14]: 92.0

In [15]: oo/0.11
Out[15]: +inf


Of course, variation of parameters has the most difficult time when there are sin and cos terms involved, because of the poor trig simplification in SymPy. So let’s see what happens with an ODE that just has exponentials and polynomial terms involved.

In [16]: time dsolve(f(x).diff(x, 2) + f(x).diff(x) - x**2 - 2*x, f(x), hint='nth_linear_constant_coeff_undetermined_coefficients')
CPU times: user 0.10 s, sys: 0.00 s, total: 0.10 s
Wall time: 0.10 s
Out[17]:
f(x) == C1 + x**3/3 + C2*exp(-x)

In [18]: time dsolve(f(x).diff(x, 2) + f(x).diff(x) - x**2 - 2*x, f(x), hint='nth_linear_constant_coeff_variation_of_parameters')
CPU times: user 0.19 s, sys: 0.00 s, total: 0.19 s
Wall time: 0.20 s
Out[19]:
f(x) == C1 + x**3/3 + C2*exp(-x)

In [20]: time dsolve(f(x).diff(x, 3) + 3*f(x).diff(x, 2) + 3*f(x).diff(x) + f(x) - 2*exp(-x) + x**2*exp(-x), f(x), hint='nth_linear_constant_coeff_undetermined_coefficients')
CPU times: user 0.09 s, sys: 0.00 s, total: 0.09 s
Wall time: 0.09 s
Out[21]:
f(x) == (C1 + C2*x + C3*x**2 + x**3/3 - x**5/60)*exp(-x)

In [22]: time dsolve(f(x).diff(x, 3) + 3*f(x).diff(x, 2) + 3*f(x).diff(x) + f(x) - 2*exp(-x) + x**2*exp(-x), f(x), hint='nth_linear_constant_coeff_variation_of_parameters')
CPU times: user 0.29 s, sys: 0.00 s, total: 0.29 s
Wall time: 0.29 s
Out[23]:
f(x) == (C1 + C2*x + C3*x**2 + x**3/3 - x**5/60)*exp(-x)


The wall time comparisons here are:

In [24]: 0.20/0.10
Out[24]: 2.0

In [25]: 0.29/0.09
Out[25]: 3.22222222222

So we don’t have orders of magnitude anymore, but it is still 2 to 3 times faster. Of course, most ODEs of this form will have sin or cos terms in them, so the order of magnitude improvement over variation of parameters can probably be attributed to undetermined coefficients in general.

Of course, we know that variation of parameters will still be useful, because functions like \ln{x}, \sec{x} and \frac{1}{x} do not have a finite number of linearly independent derivatives, and so you cannot apply the method of undetermined coefficients to them.

There is one last thing I want to mention. You can indeed multiply any polynomial, exponential, sin, or cos functions together and still get a function that has a finite number of linearly independent solutions, but if you multiply two or more of the trig functions, you have to apply the power reduction rules to the resulting function to get it in terms of sin and cos alone. Unfortunately, SymPy does not yet have a function that can do this, so to solve such a differential equation with undetermined coefficients (recommended, see above), you will have to apply them manually yourself. Also, just for the record, it doesn’t play well with exponentials in the form of sin’s and cos’s or the other way around (complex coefficients on the arguments), so you should back convert those first too.

Well, this concludes the first of two blog posts that I promised. I also promised that I would write about my summer of code experiences. Not only is this important to me, but it is a requirement. I really hope to get this done soon, but with classes, who knows.

by asmeurer at August 17, 2009 10:33 PM

Los Alamos “Sprint”


Last weekend, Luke came to visit Ondrej in Los Alamos, so I decided to drive him up from Albuquerque and visit him again. It was nice meeting Luke and seeing Ondrej again.

Aside from coding (the main thing that I did was fix an ugly match bug that was preventing dsolve() from recognizing certain ODEs), we visited the atomic museum in Los Alamos, the Valles Caldera, and some of the surrounding hot springs.

Here are some pictures that Luke took with his iPhone. Stupid WordPress seems to insist on flipping some of them (I can’t fix it):
Luke Me and Luke Me Me and Luke Me and Ondrej Me and Ondrej The Valles Caldera Ondrej and Me The road Ondrej and Me

This is one of three posts that I plan on doing this week. I just finished my GSoC project today/last night, so I will be blogging about that. I plan on doing a post on the method of Undetermined Coefficients, as well as some other things that I managed to do. The other post will be my general musings/advice for GSoC. That will probably be my last post here in a while. I plan on continuing work with SymPy, but I get very busy with classes, so I most likely won’t be doing much until next summer.

by asmeurer at August 17, 2009 06:42 PM

Fabian Pedregosa

Refine module

This commit introduced a new module in sympy: the refine module.

The purpose of this module is to simplify expressions when they are bound to assumptions. For example, if you know that x>0, then you can simplify abs(x) to x. This code was traditionally embedded into the core, but now this will be part of an external module (sympy.refine) upon which the core has no dependencies.

In a not very original move, I named the main function in this module refine(). It’s syntax is very straightforward: first argument is an expression and second argument are assumptions. Some examples (from isympy):

In [1]: refine(1+abs(x), Assume(x, Q.positive))
Out[1]: 1 + x

In [2]: refine(exp(I*x*pi), Assume(x, Q.odd))
Out[2]: -1

In [3]: refine(exp(I*x*pi), Assume(x, Q.even))
Out[3]: 1

Right now the module lacks some rules, but the design (very similar to the query module) will make adding these rules an easy task.

by fabian at August 17, 2009 05:20 PM

August 16, 2009

Ondřej Čertík

Los Alamos Sprint

Last weekend, Luke and Aaron came to visit me here in Los Alamos and we accomplished the following:

* documentation doctests fixed
* pexpect wrappers to maxima
* couple match bugs fixed
* lots of patches reviewed and pushed in
* made pydy simplify trig expressions
* pexpect wrappers to autolev
* work on the odes module
* visited hot springs

Pictures of the last item were requested, so Aaron (left), Luke:


Aaron, Ondrej:


Ondrej, Aaron:

by Ondřej Čertík (noreply@blogger.com) at August 16, 2009 05:42 PM

August 15, 2009

Jason Gedge

Java Annotations: A [Somewhat] Brief Introduction

So in my last post where I described a messaging system we implemented, I also mentioned our use of annotations. I thought it would be appropriate to write a follow-up post with a brief introduction to them, so here it is. I'm going to talk about annotations in the Java sense, but a lot of this propagates to other [reflective] languages too.Annotations are often defined as "notes of explanation

by Jason G (noreply@blogger.com) at August 15, 2009 03:35 PM

August 13, 2009

Fredrik Johansson

August 12, 2009

Dale Peterson

Developments in PyDy

I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. One of the more useful methods I implemented is the output_eoms() method in the NewtonianReferenceFrame class. Once the equations of motion have been derived, this method creates python functions for them (in the form that scipy.odeint wants them in) and writes them to a file that you specify. Once you have this, you can import the function and use it for numerical integration. Additionally, this method also outputs a function that calculates the generalized speeds given the coordinates and their time derivatives, which is useful if you know the initial conditions of the coordinates time derivatives but need to determine the initial conditions of the generalized speeds when you want to perform a numerical integration. Finally, this method also optionally outputs an animate function, which takes a variable number of Points and (Axis, Angle) tuples and outputs a function which given the coordinates and parameters, will compute the inertial position and orientation of Points and/or Reference Frames. This is useful when you want to animate your simulation.
Both the rolling_disc.py and rigidbody.py files in the pydy/examples folder demonstrate the usage of these new methods.

One issue that has been repeatedly rearing its ugly head is the use of Function instances instead of Symbol instances. For everything except time differentiation, PyDy expressions could work fine with Symbols instead of Functions. I have become to realize that everything just works faster and more efficiently in Sympy with Symbols, so I will be working in the next week or so to replace the use of Functions in PyDy with Symbol instances instead. This should allow lots of things to work much better without having to call a bunch of helper methods to get expressions in their simplest form, or having to repeatedly substitute and back substitute with Symbols and Functions. For time differentiation, I will keep track of a list of symbols that represent the generalized coordinates, and then any expression that needs time differentiation will be examined to see if it has any occurrences of those coordinates. For expressions involving those coordinates, the partial derivative will be taken with respect to each coordinate, and then multiplied by the symbol that will be used to represent the time derivative of that coordinate, in essence a manual (as opposed to letting Sympy take care of it by using sympy Function) application of the chain rule. The user will be required to eclare these coordinate symbols in a special way, specifically using the .setcoords() method. Finally, I will implement method that will replace the symbols with Function instances (depending on time), so that if one wants to play with an expression interactively, all that they will have to do is call this function so that things like q1 = Symbol(’q1′) get replaced with q1 = Function(’q1′)(t), and then time differentiation will work just fine.

by admin at August 12, 2009 10:26 PM

Aaron Meurer

Testing implicit solutions to ODEs


So, the hard deadline for GSoC it Monday, so this will probably be my last post until then (I am very busy trying to finish up the ode module by then). But this is one of those things that you just have to blog about.

So I have this checksol function in test_ode.py that attempts to check it the solutions to odes are valid or not. It was a relic of the old ode module. For that, it would just substitute the solution into the ode and see if it simplified to 0. That is what it still does, if the solution is solved for f(x) (the function for all of my ode tests). But if the solution is implicit in f, either because solve() is not good enough to solve it or because it cannot be solved, then that method obviously does not work. So what I was trying to do is what my textbook suggested. Take the derivative of the solution implicitly n times, where n is the order of the ode, and see if that is equal to the ode. Basically, I was subtracting the ode from it and seeing if it reduced to 0.

However, it wasn’t really working at all for most of my implicit solutions, even the really simple ones. I ended up XFAILing most of my implicit checksol tests. I think every single homogeneous coefficients had an implicit solution, and none of them were working with checksol().

So I started to ask around on IRC to see if anyone had any better ideas for testing these. Ondrej couldn’t think of anything. Luke and Chris worked on an example that I gave them, and it seemed to be that it wasn’t correct (which I didn’t believe for a second, because the solution was straight out of my text, and both homogeneous coefficients integrals produced that same solution). It turns out that we were mixing up \log{\frac{y}{x}} and \log{\frac{x}{y}} terms. One of those appeared in the ode and the other appeared in the solution (the ode was y dx  + x\log{\frac{y}{x}}dy - 2x dy = 0 and the solution is \frac{y}{1 + \log{\frac{x}{y}}}=C, number 9 from my odes text, pg. 61.

So Chris had a novel idea. For 1st order odes, you can take the derivative of the solution and solve for \frac{dy}{dx}, which will always be possible, because differentiation is a linear operator. Then substitute that into the original ode, and it will reduce.

So we were talking about this on IRC later, and I had an epiphany as to why my original method wasn’t working. After trying it manually on an ode, I found that I had to multiply through the solution’s derivative by \frac{x}{f(x)} to make it equal to the ode. Then, that reminded me of an important solution method that I didn’t have time to implement this summer: integrating factors. I remember that my textbook had mentioned that there is a theorem that states that every 1st order ODE that is linear in the derivative has a unique integrating factor that makes it exact. And I realized, the derivative of the solution will be equal to the ODE if and only if the ODE is exact. I checked my exact tests and verified my hunch. I had to XFAIL all of my implicit homogeneous coefficients solutions, but all of my exact checksols were working just fine.

So I refactored my checksol function to do this, and it now can check almost every one of my failing checksols. The exceptions are some where trigsimp() cannot simplify the solution to 0 (we have a poor trigsimp), a second order solution (the above trick only works on 1st order odes, I believe), and some other simplification problems.

The only down side to this new routine is that it is kind of slow (because of the simplification). I am going to have to skip a test of only 6 solutions because it takes 24 seconds to complete.

by asmeurer at August 12, 2009 09:27 PM

August 11, 2009

Fredrik Johansson

Torture testing special functions

Today I implemented a set of torture test for special functions in mpmath.

The function torturer takes a function, say f = airyai or f = lambda z: hyp1f1(1.5, -2.25, z), and attempts to evaluate it with z = a · 10n as n ranges between -20 and 20, with a ≈ 1, 1+i, i, -1+i, -1. (The "≈" indicates multiplication by an additional quasi-random factor near unity.)

Thus it tests a wide range of magnitudes (both tiny and huge), with pure positive, negative, imaginary, and complex arguments; i.e. covering most distinct regions of the complex plane. (It doesn't cover the lower half-plane, but that's generally just a waste of time since most function algorithms are agnostic about the sign of the imaginary part -- this should be done for some functions in the future, however.)

Each evaluation is also performed at a range of precisions, between 5 and 150 digits by default (for fast functions the max precision is set much higher). The results at two successive precisions are compared to verify that they agree relatively to the lesser of the precisions (with a single-digit error allowed for roundoff error).

This doesn't guarantee correctness, of course -- I might have implemented the wrong formula for the function -- but it's a strong check that whatever formula has been implemented is being evaluated accurately. It does provide some amount of correctness testing for those functions that use more than one algorithm depending on precision and/or the magnitude of the argument (and this includes the hypergeometric functions). If one algorithm is used at low precision and another at high precision, then an error in either will be revealed when comparing results.

My original intention with was just to check robustness of the asymptotic expansions of hypergeometric functions. Testing with a continuous range of magnitudes ensures that there aren't any regions where convergence gets impossibly slow, a loop fails to terminate, an error estimate is insufficient, etc. (And this needs to be done at several different levels of precision.) I then ended up including many other functions in the list of torture tests as well.

The results? I discovered (and fixed!) some 10-20 bugs in mpmath. Most were of the kind "only 43 accurate digits were returned where 50 were expected" and due to inadequate error estimates and/or insufficient extra precision being allocated (mostly in argument transformations). I actually knew about many of the problems, and had just neglected to fix them because they wouldn't show up for "normal" arguments (and in particular, there weren't any existing tests that failed).

Lesson to be learned: if it isn't tested, it probably doesn't work (fully). Semi-automatic exhaustive or randomized testing (the more "evil" inputs, the better), is a necessary complement to testing hand-picked values.

There is more to be done:


  • Some functions don't yet pass the torture tests

  • Rather than just large/small, the relevant stress measure for some functions is closeness to a negative integer, closeness to the unit circle, etc

  • Functions of more than one argument should be tested with randomly selected values for all parameters



But already now, the torture tests have done a good job. Tellingly, functions I implemented recently are much less buggy than functions I implemented long ago (even just three months ago). Periodically rewriting old code as one learns more is bound to improve quality :-)

Speaking of bugs, Mathematica doesn't pass the mpmath torture test suite (far from it). An example of something that will fail:


In[47]:= N[HypergeometricPFQ[{1},{-1/4,1/3},10^10],15]

General::ovfl: Overflow occurred in computation.

General::ovfl: Overflow occurred in computation.

General::ovfl: Overflow occurred in computation.

General::stop: Further output of General::ovfl
will be suppressed during this calculation.

(... nothing for 5 minutes)

Interrupt> a

Out[47]= $Aborted


And mpmath (instantaneously):


>>> mp.dps = 15; mp.pretty = True
>>> hyp1f2(1,(-1,4),(1,3),10**10)
-3.53521222026601e+86866


By the way, the current set of torture tests takes about 11 minutes to run on my computer with psyco and gmpy (much longer without either). It's possibly a better benchmark than the standard tests.

Update: running the torture tests with multiprocessing on sage.math takes only 28 seconds (which is the time for the slowest individual test). Go parallel!

by Fredrik Johansson (fredrik.johansson@gmail.com) at August 11, 2009 05:28 PM

August 10, 2009

Freddie Witherden

The Status of 0.2

Last Wednesday I pushed out the first release candidate of 0.2. The main motivation behind the release was to get the new caching code (needed by matplotlib) out into the wild as soon as possible.

However, soon after my mentor, Michael Droettboom, sent me a patch adding support for the MathML torture tests to mathtex. These tests resulted in several issues being filed against mathtex. Further, Michael was kind enough to provide patches for many of then — which I hastily added to mathtex.

The result? The mathtex of today has much better rendering than the mathtex of a week ago. However, it does also mean that the differences between 0.2 RC 1 and 0.2 will be substantial — enough to warrant a version-bump. It is for this reason that I am planning to skip the 0.2 release series and go straight to 0.3 RC 1.

I expect to tag the release in the next day or so, but there is at least one issue which I would like to squish before doing so.

by Freddie Witherden (noreply@blogger.com) at August 10, 2009 10:59 PM

Fabian Pedregosa

Query module - finally in trunk

The query module is finally in the main SymPy repository. I made substantial changes since last post, most of them at the user interface level (thanks to Vinzent and Mateusz for many insightful comments).

Main function is ask(), which replaces the old expression.is_* syntax. You can ask many things. For example, you can ask if a given expression is integer, prime or real:

>>> ask(2, Q.integer)
True
>>> ask(x**2, Q.integer)
None
>>> ask(x**2, Q.integer, Assume(x, Q.integer))
True
>>> ask(sqrt(2)*x, Q.integer, Assume(x, Q.integer))
>>> ask(I*x, Q.real, Assume(x, Q.imaginary))
True
>>> ask(x*y, Q.prime, Assume(x, Q.prime) & Assume(y, Q.prime))
False

as you see, it returns True when it is sure that the expression is an integer, None if it does not know, and False if it is certainly not an integer.

The second argument, which we will call ‘key’ specifies what we want to ask about. For example, Q.integer ask wether expression is an integer, Q.negative wether it’s a negative, etc.
For a complete list of all the keys available, see doc/src/modules/queries.txt in sympy codebase.

It also accepts an optional third argument where you can specify assumptions that symbols in expr satisfy. That can be any kind of boolean expression involving assumptions, for example:

Assume(x, Q.positive) & Assume(x, Q.integer)

or

Assume(x, Q.positive) | Assume(x, Q.negative)

, or even

NAnd(Assume(x, Q.positive), Assume(x, Q.integer)

by fabian at August 10, 2009 07:51 PM

Aaron Meurer

Homogeneous coefficients corner case


Before I started the program, I implemented Bernoulli equations. But the general solution to Bernoulli equations involves raising something to the power of \frac{1}{1-n}, where n is the power of the dependent term (see the Wikipedia page for more info). This works great, as I soon discovered, unless n == 1. Then you get something to the power of \infty. So I had to go in and remove the corner case.

So you think that after that I would have been more careful after that about checking that if general solution that divides by something I would test to see if that something is not zero before returning it as a solution.

Well, as I was just trying to implement some separable equation tests, I was going through the exercises of my ode text as I usually do for tests, and I came across xy' - y = 0. If you recall, this equation also has coefficients that homogeneous of the same order (1). From the general solution to homogeneous coefficients, you would plug it into \int{\frac{dx}{x}}=\int{\frac{-Q(1,u)du}{P(1,u)+uQ(1,u)}}+C where u = \frac{y}{x} or \int{\frac{dy}{y}}=\int{\frac{-P(u,1)du}{uP(u,1)+Q(u,1)}}+C where u = \frac{x}{y} (here, P and Q are from the general form P(x,y)dx+Q(x,y)dy=0). Well, it turns out that if you plug the coefficients from my example into those equations, the denominator will become 0 for each one. So I (obviously) need to check for that P(1,u)+uQ(1,u) and uP(u,1)+Q(u,1) are not 0 before running the homogeneous coefficients solver on a differential equation.

by asmeurer at August 10, 2009 05:30 PM

August 09, 2009

Fredrik Johansson

3D visualization of complex functions with matplotlib

Hooray! matplotlib 0.99 is out and it has 3D plotting, finally!

I've shown a lot of color plots of complex functions on this blog to demonstrate complex functions in mpmath. These plots are informative, but sometimes a 3D plot (typically of the function's absolute value) gives a much better view. A big advantage of 3D plots over 2D color plots is that far fewer evaluation points are required for a good high-resolution image, and this helps when visualizing the slower functions in mpmath.

There will probably be a 3D plot function in a future version of mpmath (or two functions; for two-variable real, and complex functions), similar in style to the existing matplotlib wrappers plot and cplot. Until I've figured out the details, I'll share a couple of test plots.

Coulomb wave function of a complex argument, F(6,4,z):



Principal-branch logarithmic gamma function:



Gamma function:



Imaginary part of Lambert W function, 0th branch:



Riemann zeta function in the critical strip:



Those are all wireframe plots. It's also possible to do surface plots. Using the Riemann zeta function again, a surface plot looks as follows:



Surface + wireframe plot:



None of these is perfect. I'd like to be able to do a surface plot with a widely spaced wireframe mesh to pronounce the geometry, and smooth colored surface between the meshes. This doesn't seem possible with matplotlib because the surface plot doesn't do smoothing (even with a stride selected); overlaying a wireframe works decently, but the wireframe doesn't render with occlusion and this looks bad for some functions. Since the pure wireframe plot is much faster, I think I prefer it for now.

For complex functions, it'd also be nice with a color function separate from the geometry function -- then you could plot the phase as color in the same picture.

Color helps for visualizing complicated structure, especially for phase plots. Below, I've plotted the phase (actually the sine of the phase, to make it continuous) of a Jacobi theta function θ3(1+4j/3,q) (restricted to |q| 0.8, because it gets extremely messy for larger |q|):



The phase of the Barnes G-function:



I could do many more, but that will probably enough for this blog :-)

To reproduce any of these, use the following script with trivial modifications. It's a straightforward extension of the example scripts in the matplotlib documentation.


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import pylab
import numpy as np
import mpmath
mpmath.dps = 5

# Use instead of arg for a continuous phase
def arg2(x):
return mpmath.sin(mpmath.arg(x))

#f = lambda z: abs(mpmath.loggamma(z))
#f = lambda z: arg2(mpmath.exp(z))
#f = lambda z: abs(mpmath.besselj(3,z))
f = lambda z: arg2(mpmath.cos(z))

fig = pylab.figure()
ax = Axes3D(fig)
X = np.arange(-5, 5, 0.125)
Y = np.arange(-5, 5, 0.125)
X, Y = np.meshgrid(X, Y)
xn, yn = X.shape
W = X*0
for xk in range(xn):
for yk in range(yn):
try:
z = complex(X[xk,yk],Y[xk,yk])
w = float(f(z))
if w != w:
raise ValueError
W[xk,yk] = w
except (ValueError, TypeError, ZeroDivisionError):
# can handle special values here
pass
print xk, xn

# can comment out one of these
ax.plot_surface(X, Y, W, rstride=1, cstride=1, cmap=cm.jet)
ax.plot_wireframe(X, Y, W, rstride=5, cstride=5)

pylab.show()

by Fredrik Johansson (fredrik.johansson@gmail.com) at August 09, 2009 11:23 PM

August 07, 2009

Fabian Pedregosa

django, change language settings dynamically

After some failed attempts, I just found how to change the language settings dynamically in django, and I thought it could be useful to someone. Just use function activate() from django.utils.translation. For example:

from django.utils.translation import activate

activate(‘es-ES’)

will change global language settings to ‘es-ES’ (Spain spanish). I use it because I have two forms and I want that errors appear in different languages, but I do not want to get into gettext

by fabian at August 07, 2009 02:12 PM

August 05, 2009

Fredrik Johansson

Coulomb wave functions and orthogonal polynomials

Work has been a bit slow over the last two weeks, partially due to the fact that I was away during one of them. Nevertheless, I've been able to add a couple of more functions to mpmath, as described in the following post below. With these functions committed, I'm probably going to stop adding features and just do a new release as soon as possible (some documentation and general cleanup will be required first).

Coulomb wave functions

I have, due to a request, implemented the Coulomb wave functions (commit). The Coulomb wave functions are used in quantum physics; they solve the radial part of the Schrödinger equation for a particle in a 1/r potential.

The versions in mpmath are fully general. They work not only for an arbitrarily large or small radius, but for complex values of all parameters. Some example evaluations:


>>> mp.dps = 25; mp.pretty = True
>>> coulombf(4, 2.5, 1000000)
-0.9713831935409625197404938
>>> coulombf(2+3j, 1+j, 100j)
(2.161557009297068729111076e+37 + 1.217455850269101085622437e+39j)
>>> coulombg(-3.5, 2, '1e-100')
2.746460651456730226435998e+252


The definitions of the Coulomb wave functions for complex parameters are based mostly on this paper by N. Michel, which describes a special-purpose C++ implementation. I'm not aware of any other software that implements Coulomb wave functions, except for GSL which only supports limited domains in fixed precision.

The Coulomb wave functions are numerically difficult to calculate: the canonical representation requires complex arithmetic, and the terms vary by many orders of magnitude even for small parameter values. Arbitrary-precision arithmetic, therefore, is almost necessary even to obtain low-precision values, unless one uses much more specialized evaluation methods. The current mpmath versions are simple and robust, although they can be fairly slow in worst cases (i.e. when they need to use a high internal precision). Timings for average/good cases are in the millisecond range,


>>> timing(coulombf, 3, 2, 10)
0.0012244852348553437
>>> timing(coulombg, 3, 2, 10)
0.0068572681723514609


but they can climb up towards a second or so when large cancellations occur. Someone who needs faster Coulomb wave functions, say for a limited range of parameters, could perhaps find mpmath useful for testing a more specialized implementation against.

The speed is plenty for basic visualization purposes. Below I've recreated graphs 14.1, 14.2 and 14.5 from Abramowitz & Stegun.











mp.dps = 5

# Figure 14.1
F1 = lambda x: coulombf(x,1,10)
F2 = lambda x: coulombg(x,1,10)
plot([F1,F2], [0,15], [-1.2, 1.7])

# Figure 14.3
F1 = lambda x: coulombf(0,0,x)
F2 = lambda x: coulombf(0,1,x)
F3 = lambda x: coulombf(0,5,x)
F4 = lambda x: coulombf(0,10,x)
F5 = lambda x: coulombf(0,x/2,x)
plot([F1,F2,F3,F4,F5], [0,25], [-1.2,1.6])

# Figure 14.5
F1 = lambda x: coulombg(0,0,x)
F2 = lambda x: coulombg(0,1,x)
F3 = lambda x: coulombg(0,5,x)
F4 = lambda x: coulombg(0,10,x)
F5 = lambda x: coulombg(0,x/2,x)
plot([F1,F2,F3,F4,F5], [0,30], [-2,2])


Also, a plot for complex values, cplot(lambda z: coulombg(1+j, -0.5, z)):



Orthogonal polynomials

I've added the (associated) Legendre functions of the first and second kind, Gegenbauer, (associated) Laguerre and Hermite functions (commits 1, 2, 3, 4). This means that mpmath finally supports all the classical orthogonal polynomials. As usual, they work in generalized form, for complex values of all parameters.

A fun thing to do is to verify orthogonality of the orthogonal polynomials using numerical quadrature:

>>> mp.dps = 15
>>> chop(quad(lambda t: exp(-t)*t**2*laguerre(3,2,t)*laguerre(4,2,t), [0,inf]))
0.0
>>> chop(quad(lambda t: exp(-t**2)*hermite(3,t)*hermite(4,t), [-inf,inf]))
0.0


As another basic demonstration, here are some Legendre functions of the second kind, visualized:
F1 = lambda x: legenq(0,0,x)
F2 = lambda x: legenq(1,0,x)
F3 = lambda x: legenq(2,0,x)
F4 = lambda x: legenq(3,0,x)
F5 = lambda x: legenq(4,0,x)
plot([F1,F2,F3,F4,F5], [-1,1], [-1.3,1.3])




Implementing these functions is a lot of work, I've found. The general case is simple (just fall back to generic hypergeometric code), but the special cases (singularities, limits, asymptotically large arguments) are easy to get wrong and there are lots of them to test and document.

My general methodology is to implement the special functions as the Wolfram Functions site defines them (if they are listed there), test my implementation against exact formulas, and then test numerical values against Mathematica. Unfortunately, Wolfram Functions often leaves limits undefined, and is not always consistent with Mathematica. In fact, Mathematica is not even consistent with itself. Consider the following:


In[277]:= LaguerreL[-1, b, z] // FunctionExpand

Out[277]= 0

In[278]:= LaguerreL[-1, -2, z] // FunctionExpand

z 2
E z
Out[278]= -----
2



It's tempting to just leave such a case undefined in mpmath and move on (and perhaps fix it later if it turns out to be used). Ideally mpmath should handle all singular cases correctly and consistently, and it should state in the documentation precisely what it will calculate for any given input so that the user doesn't have to guess. Testing and documenting special cases is very time-consuming, however, and although I'm making progress, this work is far from complete.

by Fredrik Johansson (fredrik.johansson@gmail.com) at August 05, 2009 01:58 AM

August 03, 2009

Jason Gedge

Java: not always so cross-platform

There's no doubting the fact that Java really makes one's life far easier in general when creating a piece of software that works on various platforms. But even with such power, there's always little things that really can make one's life a pain when creating a piece of software that is intended to be high quality.One of the major sources of pain in Java is cross-platform GUIs. It's nice that

by Jason G (noreply@blogger.com) at August 03, 2009 03:53 PM

Becoming a Jedi, err...JNI Master!

I myself am far from being a JNI master, since I only started doing my first JNI a few days ago. The problem we were having is that there is no way in Java to get some native behavior on OS X. For example, the closest to native you can get with Open/Save file dialogs is what you get from java.awt.FileDialog, which isn't very much.Since we're developing in Swing it is key that if we do anything,

by Jason G (noreply@blogger.com) at August 03, 2009 03:45 PM

August 01, 2009

Freddie Witherden

Rendering Improvements

Over the past couple of days I have been working to improve the rendering quality of mathtex. Rather than describe the changes with several paragraphs of prose I thought it more apt to have a before-and-after picture instead.

Before:

After:As can be seen the result is much closer to that produced by LaTeX.

by Freddie Witherden (noreply@blogger.com) at August 01, 2009 11:07 PM

Aaron Meurer

Variation of Parameters and More


Well, the last time I posted a project update, I had resigned myself to writing a constant simplifying function and putting the Constant class on the shelf. Well, just as I suspected, it was hell writing it, but I eventually got it working. Already, what I have in dsolve() benefits from it. I had many solutions with things like \frac{x}{C_1} or -C_1x in them, and they are now automatically reduced to just C_1x. Of course, the disadvantage to this, as I mentioned in the other post, is that it will only simplify once. Also, I wrote the function very specifically for expressions returned by dsolve. It only works, for example, with constants named sequentially like C1, C2, C3 and so on. Even with making it specialized, it was still hell to write. I was also able to get it to renumber the constants, so something like C2*sin(x) + C1*cos(x) would get transfered to C1*sin(x) + C2*cos(x). It uses Basic._compare_pretty() (thanks to Andy for that tip), so it will always number the constants in the order they are printed.

Once I got that working, it was just little work to finish up what I had already started with solving general linear homogeneous odes (a_ny^{(n)} + a_{n-1}y^{(n-1)} + \dots + a_2y'' + a_1y' + a_0y = 0 with a_i constant for all i). Solving these equations is easy. You just set up a polynomial of the form a_nm^n + a_{n-1}m^{n-1} + \cdots + a_2m^2 + a_1m + a_0 = 0 and find the roots of it. Then you plug the roots into an exponential times x^i for i from 1 to the multiplicity of the root (as in Cx^ie^{root \cdot x}). You usually expand the real and complex parts of the root using Euler’s Formula, and, once you simplify the constants, you get something like x^ie^{realpart \cdot x}(C_1\sin{(impart \cdot x)} + C_2\cos{(impart \cdot x)}) for each i from 1 to the multiplicity of the root. Anyway, with the new constantsimp() routine, I was able to set this whole thing up as one step, because if the imaginary part is 0, then the two constants will be simplified into each other. Also, SymPy has some good polynomial solving, so I didn’t have any problems there. I even made good use of the collect() function to factor out common terms, so you get something like (C_1 + C_2x)e^{x} instead of C_1e^{x} + C_2xe^{x}, which for larger order solutions, can make the solution much easier to read (compare for example, ((C_1 + C_2x)\sin{x} + (C_3 + C_4x)\cos{x})e^{x} with the expanded form, C_1e^{x}\sin{x} + C_2xe^{x}\sin{x} + C_3\cos{x}e^{x} + C_4x\cos{x}{e^x} as the solution to {\frac {d^{4}}{d{x}^{4}}}f \left( x \right) -4\,{\frac {d^{3}}{d{x}^{3}}}f \left( x \right) +8\,{\frac {d^{2}}{d{x}^{2}}}f \left( x \right) -8\,{\frac {d}{dx}}f \left( x \right) +4\,f \left( x \right) =0).

I entered all 30 examples from the relevant chapter of my text (Ordinary Differential Equations by Morris Tenenbaum and Harry Pollard), and the whole thing runs in under 2 seconds on my machine. So it is fast, though that is mostly due to fast polynomial solving in SymPy.

So once I got that working well, I started implementing variation of parameters, which is a general method for solving all equations of form a_ny^{(n)} + a_{n-1}y^{(n-1)} + \dots + a_2y'' + a_1y' + a_0y = F(x). The method will set up an integral to represent the particular solution to any equation of this form, assuming that you have all n linearly independent solutions to the homogeneous equation a_ny^{(n)} + a_{n-1}y^{(n-1)} + \dots + a_2y'' + a_1y' + a_0y = 0. The coefficients a_i do not even have to be constant for this method to work, although they do have to be in my implantation because otherwise it will not be able to find general solution to the homogeneous equation.

So, aside from doing my GSoC project this summer, I am also learning Linear Algebra, because I could not fit it in to my schedule next semester and I need to know it for my Knot Theory class. It turns out that it was very useful in learning the method of variation of parameters. I will explain how the method works below, but first I have a little rant.

Why is the Wikipedia article on variation of parameters the only website anywhere that covers variation of parameters in the general case? Every other site that I could find only covers 2nd order equations, which I understand is what is taught in most courses because applying it to anything higher can be tedious and deriving the nth order case requires knowledge of Cramer’s Rule, which many students may not know. But you would think that there would at least be sites that discuss what I am about to discuss below, namely, applying it to the general case of an nth order inhomogeneous linear ode. Even the Wolphram MathWorld article only explains the derivation for a second order linear ODE, mentioning at the bottom that it can be applied to nth order linear ODEs. I did find a website called Planet Math that covers the general case, but it wasn’t on the top of the Google results list and took some digging to find. It also has problems of its own, like being on a very slow server and some of the LaTeX on the page not rendering among them.

This partially annoys me because the Wikipedia article is not very well written. You have to read through it several times to understand the derivation (I will try to be better below). The Planet Math site is a little better, but like I said, it took some digging to find, and I actually found it after I had written up half of this post already.

But it is also part of a larger attitude that I am finding more and more of where anything that is not likely to be directly applied is not worth knowing and thus not worth teaching. Sure, it is not likely that any person doing hand calculations will ever attempt variation of parameters on an ode of order higher than 2 or 3, but that is what computer algebra systems like SymPy are for. Unfortunately, it seems that they are also in a large part for allowing you to not know how or why something mathematically is true. What difference does it make if variation of parameters can be applied to a 5th order ODE if I have to use Maple to do actually do it anyway. As long as the makers of Maple know how to apply variation of parameters to a nth order ODE, I can get along just fine. At least with SymPy, the source is freely available, so anyone who does desire to know how things are working can easily see. Anyway, I am done ranting now, so if you were skipping that part, this would be the point to start reading again.

So you have your linear inhomogeneous ODE: a_ny^{(n)} + a_{n-1}y^{(n-1)} + \dots + a_2y'' + a_1y' + a_0y = F(x). a_n cannot be zero (otherwise it would be a n-1 order ODE), so we can and should divide through by it. Lets pretend that we already did that, and just use the same letters. Also, I will rewrite a_n as a_n(x) to emphasize that the coefficients do not have to be constants for this to work. So you have your linear inhomogeneous ODE: y^{(n)} + a_{n-1}(x)y^{(n-1)} + \dots + a_2(x)y'' + a_1(x)y' + a_0(x)y = F(x). So, as I mentioned above, we need n linearly independent solutions to the homogeneous equation y^{(n)} + a_{n-1}(x)y^{(n-1)} + \dots + a_2(x)y'' + a_1(x)y' + a_0(x)y = 0 to use this method. Let us call those solutions y_1(x), y_2(x), \dots, y_n(x). Now let us write our particular solution as y_p(x) = c_1(x)y_1(x), c_2(x)y_2(x), \dots, c_n(x)y_n(x). Now, if we substitute our particular solution in to the left hand side of our ODE, we should get F(x) back. So we have (y_p)^{(n)} + a_{n-1}(x)(y_p)^{(n-1)} + \dots + a_2(x)y_p'' + a_1(x)y_p' + a_0(x)y_p = F(x). Now, let me rewrite y_p as a summation to help keep things from getting too messy. I am also going to write c_i instead of c_i(x) on terms for additional sanity. Every variable is a function of x. y_p(x) = \sum_{i=1}^{n} c_i y_i. The particular solution should satisfy the condition of the ODE, so
y_p^{(n)} + a_{n-1}y_p^{(n-1)} + \dots + a_2y_p'' + a_1y_p' + a_0y_p = F(x).

(\sum_{i=1}^{n} c_i y_i)^{(n)} + a_{n-1}(\sum_{i=1}^{n} c_i y_i)^{(n-1)} + \dots + a_2(\sum_{i=1}^{n} c_i y_i)^{(2)} +
a_1(\sum_{i=1}^{n} c_i y_i)^{(1)} + a_0\sum_{i=1}^{n} c_i y_i = F(x).

Now, if we apply the product rule to this, things will get ugly really fast, because we have to apply the product rule on each term as many times as the order of that term (the first term would have to be applied n times, the second, n-1 times, and so on). But there is a trick that we can use. In the homogeneous case, there is no particular solution, so in that case the c_i terms must all vanish identically because the solutions are linearly independent of one another. Thus, if we plug the particular solution into the homogeneous case, we get

(\sum_{i=1}^{n} c_i y_i)^{(n)} + a_{n-1}(\sum_{i=1}^{n} c_i y_i)^{(n-1)} + \dots + a_2(\sum_{i=1}^{n} c_i y_i)^{(2)} +
a_1(\sum_{i=1}^{n} c_i y_i)^{(1)} + a_0\sum_{i=1}^{n} c_i y_i = 0.

We already know that if we plug the y_i terms in individually of the c_i terms, that the expression will vanish identically because the y_i terms are solutions to the homogeneous equation. The product rule on each term will be evaluated according to the Leibniz Rule, which is that (c_i \cdot f_i)^{(n)}=\sum_{k=0}^n {n \choose k} c_i^{(k)} y_i(x)^{(n-k)}. Now the c_i y_i^{(n)} terms will vanish because we can factor out a c_i and they will be exactly the homogeneous solution. Because the expression is identically equal to zero, the remaining terms must vanish as well. If we assume that each \sum_{i=1}^n c_i' y_i^{(j)}=0 for each j from 0 to n-2, then this will take care of this; the terms with higher derivatives on c_i will also be 0, if this is true, though we do not need them for our derivation. In other words,
 c_1' y_1  + c_2' y_2 + \cdots + c_n' y_n = 0
c_n' y_1' + c_n' y_2' + \cdots + c_n' y_n' = 0
\vdots
c_n' y_1^{(n-2)} + c_n' y_2^{(n-2)} + \cdots + c_n' y_n^{(n-2)} = 0.

So, turning back to our original ODE with the particular solution substituted in, we have
(\sum_{i=1}^{n} c_i y_i)^{(n)} + a_{n-1}(\sum_{i=1}^{n} c_i y_i)^{(n-1)} + \dots + a_2(\sum_{i=1}^{n} c_i y_i)^{(2)} +
a_1(\sum_{i=1}^{n} c_i y_i)^{(1)} + a_0\sum_{i=1}^{n} c_i y_i = F(x).
But we know that most of the terms of this will vanish, from our assumption above. If we remove those terms, what remains is \sum_{i=1}^{n} c_i' y_i^{(n-1)} = F(x). So this is where it is nice that I learned Cramer’s Rule literally days before learning how to do Variation of Parameters in the general case. We have a system of n equations (the n-1 from above, plus the one we just derived), of n unknowns (the c_i terms). The determinant that we use here is used often enough to warrant a name: the Wronskian. We have that c_i' = \frac{W_i(x)}{W(x)}, or c_i = \int \frac{W_i(x)}{W(x)}, where W_i(x) is the Wronskian of the fundamental system with the ith column replaced with \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ F(x) \end{bmatrix}. So we finally have y_p = \sum_{i=1}^n \int \frac{W_i(x)}{W(x)} y_i.

Well, that’s the theory, but as always here, that is only half of the story. A Wronskian function is already implemented in SymPy, and finding W_i(x) simply amounts to F(x) times the Wronskian of the system without the ith equation, all times (-1)^i. So implementing it was easy enough. But it soon became clear that there would be some problems with this method. Sometimes, the SymPy would return a really simple Wronskian, something like -4e^{2x}, but other times, it would return something crazy. For example, consider the expression that I reported in SymPy issue 1562. The expression is (thanks to SymPy’s latex() command, no thanks to WordPress’s stupid auto line breaks that have forced me to upload my own image. If it wasn’t such a pain, I would do it for every equation, because it looks much nicer.):
Crazy Trig Wronskian (SymPy).
This is the Wronskian, as calculated by SymPy’s wronskian() function, of
\begin{bmatrix}x \sin{x}, & \sin{x}, & 1, & x \cos{x}, & \cos{x}\end{bmatrix}, which is the set of linearly independent solutions to the ODE {\frac {d^{5}}{d{x}^{5}}}f \left( x \right) +2\,{\frac {d^{3}}{d{x}^{3}}}f \left( x \right) +{\frac {d}{dx}}f \left( x \right) -1. Well, the problem here is that, as verified by Maple, that complex Wronskian above is identically equal to -4. SymPy’s simplify() and trigsimp() functions are not advanced enough to handle it. It turns out that in this case, the problem is that SymPy’s cancel() and factor() routines do not work unless the expression has only symbols in it, and that expression requires you to cancel and factor to find the \cos^2{x} + \sin^2{x} (see the issue page for more information on this). Unfortunately, SymPy’s integrate() cannot handle that unsimplified expression in the denominator of something, as you could imagine, and it seems like almost every time that sin’s and cos’s are part of the solution to the homogeneous equation, the Wronskian becomes too difficult for SymPy to simplify. So, while I was hoping to slip along with only variation of parameters, which technically solves every linear inhomogeneous ODE, it looks like I am going to have to implement the method of undetermined coefficients. Variation of parameters will still be useful, as undetermined coefficients only works if the expression on the right hand side of the equation, F(x) has a finite set of linearly independent derivatives (such as sin, cos, exp, polynomial terms, and combinations of them (I’ll talk more about this whenever I implement it).

The good news here is that I discovered that I was wrong. I had previously believed that among the second order special cases were cases that could only be handled by variation of parameters or undetermined coefficients, but it turns out I was wrong. All that was implemented were the homogeneous cases for second order linear with constant coefficients. In addition to this, there was one very special case ODE that Ondrej had implemented for an example (examples/advanced/relativity.py). The ODE is
-2({\frac{d}{dx}}f(x)){e^{-f(x)}}+x({\frac{d}{dx}}f(x))^{2}{e^{-f(x)}}-x({\frac{d^{2}}{d{x}^{2}}}f(x)){e^{-f(x)}}, which is the second derivative of xe^{-f(x)} with respect to x. According to the example file, it is know as Einstein’s equations. Maple has a nice odeadvisor() function similar to the classify_ode() function I am writing for SymPy that tells you all of the different ways that an ODE can be solved. So, I plugged that ODE into it and got a few possible methods out that I could potentially implement in SymPy to maintain compatibility with the example equation. The chief one is that the lowest order of f in the ODE is 1 (assuming you divide out the e^{-f(x)} term, which is perfectly reasonable as that term will never be 0. You can then make the substitution u = f'(x), and you will reduce the order of the ODE to first order, which in this case would be a Bernoulli equation, the first thing that I ever implemented in SymPy.

But I didn’t do that. Reduction of order methods would be great to have for dsolve(), but that is a project for another summer. Aside from that method, Maple’s odeadvisor() also told me that it was a Liouville ODE. I had never heard of that method, and neither it seems has Wikipedia or “Uncle Google” (as Ondrej calls it). Fortunately, Maple’s Documentation has a nice page for each type of ODE returned by odeadvisor(), so I was able to learn the method. The method relies on Lie Symmetries and exact second order equations, neither of which I am actually familiar with, so I will not attempt to prove anything here. Suffice it to say that if an ODE has the form
{\frac{d^{2}}{d{x}^{2}}}y(x)+g(y(x))({\frac{d}{dx}}y(x))^{2}+f(x){\frac{d}{dx}}y(x)=0, then the solution to the ODE is
\int^{y(x)}{e^{\int g(a){da}}}{da}+C1\int{e^{-\int f(x){dx}}}{dx}+C2=0
You could probably verify this by substituting the solution into the original ODE. See the Maple Documentation page on Liouville ODEs, as well as the paper they reference (Goldstein and Braun, “Advanced Methods for the Solution of Differential Equations”, see pg. 98).

The solution is very straight forward–as much so as first order linear or Bernoulli equations, so it was a cinch to implement it. It looks like quite a few differential equations generated by doing F''(y(x), x) for some function or x and y F(y(x), x) generates equations of that type, so it could be actually useful for solving other things.

Before I sign off, I just want to mention one other thing that I implemented. I wanted my linear homogeneous constant coefficient ODE solver to be able to handle ODEs for which SymPy can’t solve the characteristic equation, for whatever reason. SymPy has RootOf() objects similar to Maple that let you represent the roots of a polynomial without actually solving it, or even being able to solve it, but a you can only use RootOf’s if you know that none of the roots are repeated. Otherwise, you would have to know which terms require an additional x^i to preserve linear independence. Well, it turns out that there is a way to tell if a polynomial has repeated roots without solving for them. There is a number associated with every polynomial of one variable called the discriminant. For example, the discriminant of the common quadratic polynomial ax^2 + bx + c is the term under the square root of the famous solution b^2 - 4ac. It is clear that a quadratic has repeated roots if and only if the discriminant is 0. Well, the same is true for the discriminant of any polynomial. I am not highly familiar with this (ask me again after I have taken my abstract algebra class next semester), but apparently there is something called the resultant, which is the product of the differences of roots between two polynomials and which can also be calculated without explicitly finding the roots of the polynomials. Clearly, this will be 0 if and only if the two polynomials share a root. So the discriminant is built from the fact that a polynomial has a repeated root iff it shares a root with its resultant. So it is basically the resultant of a polynomial and its derativave, times an extra factor. It is 0 if and only if the polynomial has a repeated root.

Fortunately, SymPy’s excelent Polys module already had resultants implemented (quite efficiently too, I might add), so it was easy to implement the discriminant. I added it as issue 1555. If you are a SymPy developer and you have somehow managed to make yourself read this far (bless your heart), please review that patch.

Well, this has turned out to be one hella long blog post. But what can I say. You don’t have to read this thing (except for possibly my mentor. Sorry Andy). And I haven’t been quite updating weekly like I am supposed to be, so this compensates. If you happened upon this blog post because, like me, you were looking for a general treatment of variation of parameters, I hope you found my little write up helpful. And if you did, and you now understand it, could you go ahead and improve the Wikipedia article. I’m not up to it.

by asmeurer at August 01, 2009 05:39 PM

July 29, 2009

Freddie Witherden

YAPP (Yet Another Progress Report)

It only just hit that it has been well over a week since my last blog post. This is unfortunate, as I like to blog at least once a week. Without further ado here is what has been happening in mathtex this week.

Firstly, the matplotlib integration was tidied up. Mathtex is currently an optional component of matplotlib and has completely replaced the internal TeX rendering engine. I am hopeful that with a few tweaks it can be merged in the near future into the mainline matplotlib repository.

Secondly, the first public version of mathtex was released — 0.1 RC 1. However, due to an (embarrassing) slip-up on my part it was quickly superseded by 0.1 RC 2. For those who are interested it can be downloaded from http://code.google.com/p/mathtex/downloads/list Assuming all is well I hope to have a final release out by the end of the week.

Finally, the mailing list has been more active than ever! Over the past few days there have been several interesting discussions on equation breaking, build errors and how to use mathtex from a C/C++ application. Although not broad enough to constitute a blog post the topics make interesting reads nether the less.

by Freddie Witherden (noreply@blogger.com) at July 29, 2009 12:52 PM

July 28, 2009

Dale Peterson

Recursive methods for computing equations of motion

In the last week, I added methods to the NewtonianReferenceFrame class of PyDy which start at the top of the reference frame tree and the point tree and do several things. The first method, recursive_subs(expr), is designed to take kinematic differential equations (which arise through a combination of linear nonholonomic constraints and/or holonomic constriants, and definitions of generalized speeds) and substitute them into all expressions for relative angular velocity between ReferenceFrames and relative velocity between Points. In doing so, the angular velocity and velocity expressions will involve only the independent time derivatives of coordinates. The second method, recursive_acc(), parses the Point/ReferenceFrame trees and forms the relative angular acceleration and acceleration. The third method, apply_gravitational_forces(), allows one to apply a gravitational forces to every particle and rigid body in the system with one command.

Currently, the NewtonianReferenceFrame class is acting like a container class for information about the system which is somewhat global in nature. For example, definitions of generalized coordinates, their time derivatives, and generalized speeds must all be assigned as attributes of this class in order for it to work correctly.

On the todo list is to write fr() and frstar(), which will also be recursive functions which will compute the equations of motion in one command. More work on the functions dealing with constraints is also needed, and there are at least two approaches that I could potentially implement for this. In particular, the approach taken with nonholonomic systems versus holonomic systems may necessitate different approaches. Currently, the work flow is as follows:
1) Form angular velocity and velocity expressions (linear in time derivatives of coordinates)
2) Form constraint equations arising from differentiated holonomic constraints and/or simple motion (rolling) constraints (m equations, linear in time derivatives of coordinates).
3) Define n-m generalized speeds as linear combinations of time derivatives of coordinates.
4) Solve generalized speed equations (item 3) for n-m ‘independent’ time derivatives of coordinates, in terms of the n-m generalized speeds and the m ‘dependent’ time derivatives of coordinates.
5) Solve constraint equations for m ‘dependent’ time derivatives of coordinates in terms of the ‘independent’ time derivatives of coordinates.
6) The complete set of kinematic differential equations will be the equations arising from item 5 and item 4.
7) Substitute these expressions into the expressions for angular velocity and linear velocity so that the constraints and definitions for generalized speeds are imposed, and so that the constrained partial velocities and partial angular velocities can be formed.

In taking this approach, the minimal set of n-m dynamic equations of motion will be formed, and they will obey all the constraints. These n-m equations, combined with the n kinematic differential equations comprise the complete set of motion equations that describe the system.

PyDy’s examples are being updated to reflect this new functionality, especially the rolling torus example. Additionally I have added a bicycle model which should soon have the nonlinear equations of motion for a very standard bicycle model.

by admin at July 28, 2009 04:07 PM

July 24, 2009

Jason Gedge

Listening In

So it came to our attention recently that our application was making abundant use of the Observer/Listener pattern. For those not familiar with this pattern, you'd use this guy when you want the outside world to know about state changes in an object. This pattern is used often when developing with various architectural patterns, such as Model-View-Controller (MVC). Other examples, in Java,

by Jason G (noreply@blogger.com) at July 24, 2009 07:17 PM

July 21, 2009

Fabian Pedregosa

can we merge now, pleeease ?

Three months after I began to write sympy.queries, I feel it’s about time to include it in sympy’s trunk, so today I sent for review 4 patches that implement the complete query module.

It’s been a lot of fun, but it has also caused me some headaches … specially last month trying to have the code as much bug-free as possible.

Next steps is to improve performace of query() and write the refine module.

PD: This weekend I’ll be at Leipzig for EuroSciPy, hope to meet there some SymPy developers!

by fabian at July 21, 2009 08:36 PM

July 20, 2009

Dale Peterson

New trig functions

Last week I spent most of my time writing and testing new trigonometric functions for Sympy. Sympy already had sin, cos, tan, cot, asin, acos, atan, acot, but did not have sec, csc, asec, or acsc. So I added those four new functions and am in the process of writing very extensive tests for them.

The reason for writing these new functions is so that the trig functions have the same behavior as the ones in Mathematica. There are a few things that Mathematica does that is very nice, one of them being that any arguments that have a pi shift (as in sin(a + b*pi)), will always return another trig function whose pi shift is within the interval between 0 and pi/4. For example,

>>> sin(11*pi/8)
-cos(pi/8)
>>> x = Symbol(‘x’)
>>> sin(x + 11*pi/8)
-cos(-x + pi/8)
>>> sin(acot(x))
1/(x*(1 + x**(-2))**(1/2))
>>> sin(pi/10)
-1/4 + 5**(1/2)/4

Getting all the behavior to work with the inverse trig functions is also implemented, as well as the special values of the trig functions for arguments which are integer multiples of pi/10. The next task is going to be to get all of this code working with the rest of Sympy so it can be put into place.

This behavior is essential to making the trigonometric simplification routine in Sympy to be more capable, and will greatly benefit both Sympy and PyDy.

by admin at July 20, 2009 06:09 PM

Aaron Meurer

Modifying a list while looping through it in Python


Here is an interesting thing that I found. x is a SymPy Symbol:
Code block 1

I would have expected to get a = [], but it only removes the first item. And yes, x + 1 passes the condition:

>>> (x + 1).has(x)
True

Clearly, it is a bad idea to modify a list while I am looping through it. I should instead be doing something like this:
Code block 2
But I am intrigued as to why exactly this fails. If any of the readers of the blog thinks that he know, please post in the comments. Or, if I figure it out, I will post an update. Also, here is a similar example, with a strange result:
Code block 3

(Sorry for the images by the way. WordPress’s so called “code” blocks are impervious to indentation.)

UPDATE (a few minutes later):
Well, I figured it would come to me as to why this was happening, and it didn’t take long. While I haven’t read the actual Python Language Reference, this is what I am assuming is happening. This is all just my guessing on how Python is implemented. Please correct me if I am wrong.

So, obviously, a Python list is just a C array. It is probably an array of pointers, which is the only way I can see that would let it be mutable with different objects (this is how compiled languages with dynamic typing like Objective-C more or less pull it off). Now C does not have the for i in list syntax that Python has (nor does any other language that I know of. That is one of the reasons that Python is so awesome!), so if you want to recurse a list (C array), you have to do the usual for (i=0; i<=len(list); i++) { from C (or it probably uses a while loop, which would allow for things like iterators, but a for loop in C is literally just a wrapper around a while loop anyway). Then of course, inside of the loop, you just have list[i] blocks. So when I was going through my list, for example, the list of numbers in the last example, it would hit item 0 (the first item), remove it, which would amount to rebuilding the list as [2, 3, 4, 5], then it would hit item 1, which is now 3, remove it, rebuilding the list, and so on. So the even numbered elements remain because it skips every element after one that it removes. CPython must have good error handling, because eventually this would cause the indices to go beyond the length of the list. It seems to me that this behavior is not very well defined. Personally, I think that whatever you are looping through in a for loop should become immutable within the loop block. I checked Python 3.1, and the behavior is exactly the same.

Based on this, .remove() rebuilds the list each time. I would have thought it would just set the value in the array to Null, but I guess that would make it more difficult to test equality with an equivalent list that doesn’t have Null values. It is good to know that .remove() does that, because it means that can be an expensive operation.

by asmeurer at July 20, 2009 06:35 AM

July 19, 2009

Freddie Witherden

Branching Off

As the more astute (or stalker) readers may have noticed there have not been any commits to the mathtex subversion repository over the past couple of days. The reason for this is because over the past couple of days (and probably for a few more) I have been working on re-integrating mathtex back into matplotlib.

The work — for those that are interested — is being done in the mathtex branch of the matplotlib subversion repository. This reintegration is especially important as although the success rate for summer of code projects is high (around 80% I believe) a surprising number of the projects never 'go gold' and make it into actual releases. Instead they live their days out in branches and separate repositories.

To ensure that mathtex does not succumb to this fate I have started integrating it back into matplotlib as soon as possible. This has several advantages: firstly it gives the API a good work-out to ensure that there are no glaring issues — reducing the chances that the API will change in the future; secondly it allows for mathtex to have a much wider userbase than it would have otherwise.

While it should not take long to port all of the existing backends over to mathtex the precise manner in which mathtex will be bundled with matplotlib is still up for debate. (But, it is almost certain that it will be bundled, resulting in no extra dependencies for matplotlib.) Once this is done — and tested — mathtex should be on much firmer ground.

by Freddie Witherden (noreply@blogger.com) at July 19, 2009 12:09 AM

July 16, 2009

Official SymPy blog

SymPy 0.6.5 released

SymPy 0.6.5 has been released on July 17, 2009. It is available at

http://sympy.org

Source distribution can be downloaded from
http://sympy.googlecode.com/files/sympy-0.6.5.tar.gz

Windows binaries can be downloaded from
http://sympy.googlecode.com/files/sympy-0.6.5.win32.exe


About SymPy

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python and does not require any external libraries.


Changes since last stable release

This release has been marked by improved documentation,
C code generation, solve and dsolve improvements, mpath update, a new
logic module and the start of Google's Summer of Code program.


Major changes include:

- Geometric Algebra Improvements
- Upgrade GA module with left and right contraction operations
- Add intersection test for the vertical segment,
reimplementation of convex_hull (Florian Mickler)

- Implement series() as function (Barry Wardell)

- Core improvements
- Refactor Number._eval_power (Fabian Pedregosa)
- fix bugs in Number._eval_power (Chris Smith)

- Matrix improvements:
- Improve jacobian function, introduce vec and vech (Ben
Goodrich)


- Solver improvements:
- solutions past linear factor found in tsolve (Chris Smith)
- refactor sympy.solvers.guess_solve_strategy (Fabian Pedregosa)
- Small cleanups to the ODE solver and tests. (Priit Laes)
- Fix corner case for Bernoulli equation. (Priit Laes)


- Improvements on partial differential equations solvers (Priit
Laes)
- Added separation of variables for PDEs (Priit Laes)


- Expand improvements (Aaron Meurer)
- Refactoring
- exp(x)*exp(y) is no longer automatically combined into
exp(x+y), use powsimp for that


- Documentation improvements:
- Test also documentation under doc/ (Fabian Pedregosa)
- Added many docstrings
- Fix Sphinx complaints/warnings/errors (Priit Laes)
- Doctest coverage (Ondrej Certik)


- New logic module (Fabian Pedregosa)
- Efficient DPLL algorithm (Fabian Pedregosa)


- LaTeX printer improvements:
- Handle standard form in the LaTeX printer correctly (Freddie
Witherde)
- Latex: _print_Mul fix (issue 1282) (Nicolas Pourcelot)
- Robust printing of latex sub and superscripts (Toon Verstraelen)
- sorting _print_Add output using a main variable (Ryan Krauss)
- Matrix printing improvements (Ryan Krauss)


- MathML printing improvements:
- MathML's printer extended (Thomas Sidoti)


- Testing framework improvements
- Make tests pass without the "py" module (Ondrej Certik)


- Polynomial module improvements:
- Fixed subresultant PRS computation and ratint() (Mateusz
Paprocki)
- Removed old module sympy.polynomials (Priit Laes)


- limit fixes:
Compute the finite parts of the limit of a sum by direct
substitution (Ronan Lamy)


- Test coverage script (Ronan Lamy)


- Code quality improvements (remove string exceptions, code quality
test improvements) (Tomasz Buchert)


- C code generation (Toon Verstraelen)


- Update mpmath (Fredrik Johansson, Vinzent Steinberg)


The following people submitted code for this release (sorted by number
of patches):

* Fabian Pedregosa
* Ondrej Certik
* Priit Laes
* Vinzent Steinberg
* Toon Verstraelen
* Tomasz Buchert
* Nicolas Pourcelot
* Aaron Meurer
* Florian Mickler
* Jochen Voss
* Alan Bromborsky
* Barry Wardell
* Riccardo Gori
* Chris Smith
* Freddie Witherden
* Ben Goodrich
* Ronan Lamy
* Akshay Srinivasan
* Johann Cohen-Tanugi
* Luke Peterson
* Mateusz Paprocki
* Ted Horst
* Thomas Sidoti
* Vinay Kumar

by Fabian Pedregosa (noreply@blogger.com) at July 16, 2009 04:09 PM

Aaron Meurer

Constant stuff


So I was able to get a working version of the Constant class, but because the code cluttered up the internal Add and Mul classes too much, Ondrej convinced me that to make a function that does the simplification instead, and I reluctantly agreed. After begining work on it, I realized that it will be much easier to make it just an internal function that handles the special cases presented by dsolve(). That means that it will only handle arbitrary constants that are independent of one variable, and it will only work with constants that are named as “C1″, “C2″, and so on.

If we ever get the sympyx core that Ondrej and I worked on when I was in Los Alamos in, it will be easy for my to use a Constant class, because it will have handler logic that will allow for the Constant class to exist independent of Add and Mul. It already can exist independent of Pow with a minor code addition, but simplifying powers is easier than simplifying addition and multiplcation because exponentiation is neither commutative nor associative, meaning that you don’t have to worry about absorbing stuff on the other side of something, like 2 + x + C.

See my constant-Mul branch for my working version of a Constant class the implements in Mul and Add. See my constant-function branch for my work on the internal function.

Because I have decided to make thing simple and make the function internal only, I should have things up and running soon. Then, it will be simple to fix up my nth order homogeneous stuff that I already have so that it works, and then to implement variation of parameters!

by asmeurer at July 16, 2009 05:06 AM

July 15, 2009

Fredrik Johansson

Hurwitz zeta function, Dirichlet L-series, Appell F1

I've added three more functions to mpmath since the last blog update: the Hurwitz zeta function, Dirichlet L-series and the Appell hypergeometric function.

Hurwitz zeta function


The Hurwitz zeta function is available with the syntax hurwitz(s, a=1, derivative=0). It's a separate function from zeta for various reasons (one being to make the plain zeta function as simple as possible). With the optional third argument, it not only computes values, but arbitrarily high derivatives with respect to s (more on which later).

It's quite fast at high precision; e.g. the following evaluation takes just 0.25 seconds (0.6 seconds from cold cache) on my laptop:


>>> mp.dps = 1000
>>> hurwitz(3, (1,3)) # (1,3) specifies the exact parameter 1/3
27.56106119970080377622787797740750928454209531301488108286446031855695818017193
31971447814492522557242833928839896347960163079006102651117997856881812399406784
05731893349475372409101251211346682545314192272823745384665266049022078198043979
99746337871597509019667252833751410875882452970087673701166818534638337151352084
77353765613742227179887954515352387981214398169985388076469409432688912355506092
00433349370785407936365607196038965500743780377965599471761097326271922616743812
30503180321801486524844073873393780084553765789276561564723000726671702263215399
29351207878018547946361756911396554507918933344824138599697851666075594487957278
34823450478752821675163024463415647205608959739100071026336912826930683261155867
06265136825191025873785920661492516925137079914816949085911457084632783280741306
98897458344746469272828267484989756595259917290800848768074029766417786815732919
07917349747600320626141004448429987450854391175868138793015691343081417209625034
88644344900903488585429308918210445984048


With s = π instead, it takes 3.5 seconds. For comparison, Mathematica 6 took 0.5 and 8.3 seconds respectively (on an unloaded remote system which is faster than my laptop, but I'm not going to guess by how much).

The function is a bit slow at low precision, relatively speaking, but fast enough (in the 1-10 millisecond range) for plotting and such. Of course, it works for complex s, and in particular on the critical line 1/2+it. Below, I've plotted |ζ(1/2+it,1)|2, |ζ(1/2+it,1/2)|2, |ζ(1/2+it,1/3)|2; the first image shows 0 ≤ t 30 and the second is with 1000 ≤ t 1030.





To show some more complex values, here are plots of ζ(s, 1), ζ(s, 1/3) and ζ(s, 24/25) for 100 ≤ Im(s) ≤ 110. I picked the range [-2, 3] for the real part to show that the reflection formula for Re(s) 0 works:



In fact, the Hurwitz zeta implementation works for complex values of both s and a. Here is a plot of hurwitz(3+4j, a) with respect to a:



The evaluation can be slow and/or inaccurate for nonrational values of a, though (in nonconvergent cases where the functional equation isn't used).

Zeta function performance


Right now the implementations of the Riemann zeta function and the Hurwitz zeta function in mpmath are entirely separate. In fact, they even use entirely different algorithms (hurwitz uses Euler-Maclaurin summation, while zeta uses the Borwein approximation). This is useful for verifying correctness of either function.

Regarding performance, it appears that the Borwein approximation is slightly faster for small |Im(s)| while Euler-Maclaurin is massively faster for large |s|.

The following benchmark might give an idea:


>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> timing(zeta, 0.5+10j); timing(hurwitz, 0.5+10j)
0.0027386162207766627
0.0082901877193889192
>>> timing(zeta, 0.5+100j); timing(hurwitz, 0.5+100j)
0.0075319349246901089
0.041037198861866478
>>> timing(zeta, 0.5+1000j); timing(hurwitz, 0.5+1000j)
0.13468499488063479
0.18062463246027072
>>> timing(zeta, 0.5+10000j); timing(hurwitz, 0.5+10000j)
7.2977593520964241
1.0488757649366072
>>> timing(zeta, 0.5+10000j); timing(hurwitz, 0.5+10000j)
0.84306636737350971
1.0601629536714867


As the last evaluation shows, the Borwein algorithm is still not doing too poorly at 0.5+10000j from a warm cache, but at this point the cache is 50 MB large. At 0.5+100000j my laptop ran out of RAM and Firefox crashed (fortunately Blogger autosaves!) so it clearly doesn't scale. The Euler-Maclaurin algorithm caches a lot of Bernoulli numbers, but the memory consumption of this appears to be negligible. With hurwitz (and the function value, for those interested),


>>> timing(hurwitz, 0.5+100000j)
8.0361576680835149
>>> hurwitz(0.5+100000j)
(1.07303201485775 + 5.7808485443635j)


and the memory usage of the Python interpreter process only climbs to 12 MB.

Just by optimizing the pure Python code, I think the Hurwitz zeta function (at low precision) could be improved by about a factor 2 generally and perhaps by a factor 3-4 on the critical line (where square roots can be used instead of logarithms -- zeta uses this optimization but hurwitz doesn't yet); with Cython, perhaps by as much as a factor 10.

The real way to performance is not to use arbitrary-precision arithmetic, though. Euler-Maclaurin summation for zeta functions is remarkably stable in fixed-precision arithmetic, so there is no problem using doubles for most applications. As I wrote a while back on sage-devel, a preliminary version of my Hurwitz zeta code for Python complex was 5x faster than Sage's CDF zeta (in a single test, mind you). If there is interest, I could add such a version, perhaps writing it in Cython for Sage (for even greater speed).

Derivatives


The Hurwitz zeta function can be differentiated easily like so:


>>> mp.dps = 25
>>> hurwitz(2+3j, 1.25, 3)
(-0.01266157985314340398338543 - 0.06298953579777517606036907j)
>>> diff(lambda s: hurwitz(s, 1.25), 2+3j, 3)
(-0.01266157985314340398338543 - 0.06298953579777517606036907j)


For simple applications one can just as well use the numerical diff function as above with reasonable performance, but this isn't really feasible at high precision and/or high orders (numerically computing an nth derivative at precision p requires n+1 function evaluations, each at precision (n+1)×p).

The specialized code in hurwitz requires no extra precision (although it still needs to perform extra work) and therefore scales up to very high precision and high orders. Here for example is a derivative of order 100 (taking 0.5 seconds) and one of order 200 (taking 4 seconds):


>>> hurwitz(2+3j, 1.25, 100)
(2.604086240825183469410664e+107 - 1.388710675207202633247271e+107j)
>>> hurwitz(2+3j, 1.25, 200)
(2.404124816484789309285653e+274 + 6.633220818921777955999455e+273j)


It should be possible to calculate a sequence of derivatives much more quickly than with separate calls, although this isn't implemented yet. One use for this is to produce high-degree Taylor series. This could potentially be used in a future version of mpmath, for example to provide very fast moderate-precision (up to 50 digits, say) function for multi-evaluation of the Riemann zeta function on the real line or not too far away from the real line. This in turn would speed up other functions, such as the prime zeta function and perhaps polylogarithms in some cases, which are computed using series over many Riemann zeta values.

Dirichlet L-series


With the Hurwitz zeta function in place, it was a piece of cake to also supply an evaluation routine for arbitrary Dirichlet L-series.

The way it works it that you provide the character explicitly as a list (so it also works for other periodic sequences than Dirichlet characters), and it evaluates it as a sum of Hurwitz zeta values.

For example (copypasting from the docstring), the following defines and verifies some values of the Dirichlet beta function, which is defined by a Dirichlet character modulo 4:


>>> B = lambda s, d=0: dirichlet(s, [0, 1, 0, -1], d)
>>> B(0); 1./2
0.5
0.5
>>> B(1); pi/4
0.7853981633974483096156609
0.7853981633974483096156609
>>> B(2); +catalan
0.9159655941772190150546035
0.9159655941772190150546035
>>> B(2,1); diff(B, 2)
0.08158073611659279510291217
0.08158073611659279510291217
>>> B(-1,1); 2*catalan/pi
0.5831218080616375602767689
0.5831218080616375602767689
>>> B(0,1); log(gamma(0.25)**2/(2*pi*sqrt(2)))
0.3915943927068367764719453
0.3915943927068367764719454
>>> B(1,1); 0.25*pi*(euler+2*ln2+3*ln(pi)-4*ln(gamma(0.25)))
0.1929013167969124293631898
0.1929013167969124293631898


Appell hypergeometric function


The function appellf1 computes the Appell F1 function which is a hypergeometric series in two variables. I made the implementation reasonably fast by rewriting it as a single series in hyp2f1 (with the side effect of supporting the analytic continuation with that for 2F1).

A useful feature of the Appell F1 function is that it provides a closed form for some integrals, including those of the form



with arbitrary parameter values, for arbitrary endpoints (and that's a quite general integral indeed). Comparing with numerical quadrature:


>>> def integral(a,b,p,q,r,x1,x2):
... a,b,p,q,r,x1,x2 = map(mpmathify, [a,b,p,q,r,x1,x2])
... f = lambda x: x**r * (x+a)**p * (x+b)**q
... def F(x):
... v = x**(r+1)/(r+1) * (a+x)**p * (b+x)**q
... v *= (1+x/a)**(-p)
... v *= (1+x/b)**(-q)
... v *= appellf1(r+1,-p,-q,2+r,-x/a,-x/b)
... return v
... print "Num. quad:", quad(f, [x1,x2])
... print "Appell F1:", F(x2)-F(x1)
...
>>> integral('1/5','4/3','-2','3','1/2',0,1)
Num. quad: 9.073335358785776206576981
Appell F1: 9.073335358785776206576981
>>> integral('3/2','4/3','-2','3','1/2',0,1)
Num. quad: 1.092829171999626454344678
Appell F1: 1.092829171999626454344678
>>> integral('3/2','4/3','-2','3','1/2',12,25)
Num. quad: 1106.323225040235116498927
Appell F1: 1106.323225040235116498927


Also incomplete elliptic integrals are covered, so you can for example define one like this:


>>> def E(z, m):
... if (pi/2).ae(z):
... return ellipe(m)
... return 2*round(re(z)/pi)*ellipe(m) + mpf(-1)**round(re(z)/pi)*\
... sin(z)*appellf1(0.5,0.5,-0.5,1.5,sin(z)**2,m*sin(z)**2)
...
>>> z, m = 1, 0.5
>>> E(z,m); quad(lambda t: sqrt(1-m*sin(t)**2), [0,pi/4,3*pi/4,z])
0.9273298836244400669659042
0.9273298836244400669659042
>>> z, m = 3, 2
>>> E(z,m); quad(lambda t: sqrt(1-m*sin(t)**2), [0,pi/4,3*pi/4,z])
(1.057495752337234229715836 + 1.198140234735592207439922j)
(1.057495752337234229715836 + 1.198140234735592207439922j)


The Appell series isn't really faster than numerical quadrature, but it's possibly more robust (the singular points need to be given to quad as above to obtain any accuracy, for example). Mpmath doesn't yet have incomplete elliptic integrals; the version above could be a start, although I'll probably want to try a more canonical approach.

by Fredrik Johansson (fredrik.johansson@gmail.com) at July 15, 2009 05:43 AM