Planet SymPy

A few days back, I was on my way to write up a blog post detailing how impressive the new system is (and it /does/ have some great speedups in there). To be thorough, I wanted to show results for different parameters being disabled, but to my surprise I discovered clause learning was hurting. This flies in the face of everything we know about SAT solving, so the hunt for reasoning commenced.

Well long story short, the hunt continues but a likely reason is the lack of impact that the simple clause learning scheme actually has. Monitoring the clauses that cause unit propagation (a measure of how useful a clause is) reveals that the learned clauses aren't actually ever being used -- big-time bummer. It's because of this, I've disabled the conflict analysis in the latest version.

But results are still be be had. First off, how much better are we than the original DPLL solver, and how bad does the clause learning actually hurt? Well the answers are (roughly) pretty good, and not much ;). Here's a graph plotting the run-time for the three approaches (-CA means clause learning was disabled):


Note that the y-axis is on a log scale, so the differences are actually more drastic than they appear. Take-away: the new sat-solver is doing great, and the clause learning doesn't hurt that much...just doesn't help like we'd expect.

Next up is the impact of watch literals. I'll have a blog post in a day or two for the final installment of "what goes into a SAT solver", and the focus will be on watch literals. The introduction of this data structure hack heralded in a wave of killer SAT solvers, and they helped out here too. The following is a comparison (with clause learning disabled) between using watch literals and not using them:


Oh, and it should be noted that when points don't exist, it's because the solver just couldn't handle it in a reasonable amount of time.

So in the interest of actually letting SymPy users see these enhancements, I'm going to wrap up this week with that final post, and post-pone the clause learning investigation until more time is available. Next up is investigating how a truth maintenance system can be used to smartly cache queries to the assumption framework. If we can avoid repeated SAT queries (no matter how beefy the solver is), that's a win.

One final note, there will likely be an incorporation of an industrial SAT-solver interface for those who really want the power. This will be through the Numberjack Library, and be optional if you happen to have it installed. The start of it is in this branch, but we'll need to wait until the next version of Numberjack is released due to limitations it currently has.

That's all for now. Cheers.

Now that I've overcome the immediate effects of SDWS (Sage Days Withdrawal Syndrome), I should write a brief report about Sage Days 24 which took place last week at the Research Institute for Symbolic Computation (RISC), located in a small town called Hagenberg just outside Linz in Austria. Many thanks for the organizers for inviting me!

Since I'm starting as a Ph.D. student at RISC in October, it was nice to get an early view of the site and meet some current students and staff. The location is quite beautiful, my only complaint being the extremely high temperature last week (they say the weather is more normal most of the year!).

SD24 was titled "Symbolic Computation in Differential Algebra and Special Functions", which pretty accurately describes the topics covered. Some points of interest:

William Stein talked about his vision for Sage (short and long term), including goals for special functions support.

Frédéric Chyzak gave a talk about the Dynamic Dictionary of Mathematical Functions. I really like the approach of starting from a minimal definition of a special function (a differential equation with initial conditions) and generating series expansions, Chebyshev approximations (etc.) algorithmically.

Nico Temme talked about numerical evaluation of special functions, which was familiar grounds to me, but with much food for thought.

I met Manuel Kauers, my future supervisor, who gave a neat presentation of some recent research.

There were many other interesting talks and discussions as well, and I gave a talk about special function evaluation in mpmath (based on my SD23 talk). And of course, I met various other Sage developers and got some coding done as well.

I had hoped to get generalized hypergeometric functions into Sage during SD24. There is now a patch, but it still needs some more work.

I also implemented parabolic cylinder functions (PCFs) in mpmath, the last remaining chapter of hypergeometric functions listed in the DLMF (mpmath now covers chapters 1-20, 22, 24-25, partially 26-27, and 33). Code commit here.

Technically, parabolic cylinder functions are just scaled (linear combinations of) confluent hypergeometric functions or Whittaker functions, but they are a bit tricky to compute due to cancellation and branch cuts.

An example from Nico Temme's talk was to evaluate D_ν(z) with ν = -300.14 and z = 300.15 (Maple 13 takes 5 seconds to give an answer that has the wrong sign and is off by 692 orders of magnitude). The new pcfd function in mpmath instantaneously (in about a millisecond) gives the correct result:

>>> pcfd(-300.14, 300.15)
6.83322814925713e-10526
>>> mp.dps = 30
>>> pcfd('-300.14', '300.15')
6.83322814923312844480669009487e-10526

All the curve plots from DLMF 12.3 can be reproduced as follows:

f1 = lambda x: pcfu(0.5,x)
f2 = lambda x: pcfu(2,x)
f3 = lambda x: pcfu(3.5,x)
f4 = lambda x: pcfu(5,x)
f5 = lambda x: pcfu(8,x)
plot([f1,f2,f3,f4,f5], [-3,3], [0,3])

f1 = lambda x: pcfv(0.5,x)
f2 = lambda x: pcfv(2,x)
f3 = lambda x: pcfv(3.5,x)
f4 = lambda x: pcfv(5,x)
f5 = lambda x: pcfv(8,x)
plot([f1,f2,f3,f4,f5], [-3,3], [-3,3])

f1 = lambda x: pcfu(-0.5,x)
f2 = lambda x: pcfu(-2,x)
f3 = lambda x: pcfu(-3.5,x)
f4 = lambda x: pcfu(-5,x)
plot([f1,f2,f3,f4], [-8,8], [-6,6])

f1 = lambda x: pcfv(-0.5,x)
f2 = lambda x: pcfv(-2,x)
f3 = lambda x: pcfv(-3.5,x)
f4 = lambda x: pcfv(-5,x)
plot([f1,f2,f3,f4], [-8,8], [-6,6])

f1 = lambda x: pcfu(-8,x)
f2 = lambda x: pcfv(-8,x)*gamma(0.5-(-8))
f3 = lambda x: hypot(f1(x), f2(x))
plot([f1,f2,f3], [-6,6], [-120,120])

f1 = lambda x: diff(pcfu,(-8,x),(0,1))
f2 = lambda x: diff(pcfv,(-8,x),(0,1))*gamma(0.5-(-8))
f3 = lambda x: hypot(f1(x), f2(x))
plot([f1,f2,f3], [-6,6], [-220,220])

A little more information can be found in the mpmath documentation section for PCFs.

Thanks again to the NSF grant enabling this work.

This week was quite spread out in terms of what I coded in quantum.py. I made inner products work with operators and other objects in between them. We made pretty print work with various quantum classes. And finally I made a suite of validating/combining/separating functions for inner and outer products.

First of all, let's look at how pretty makes quantum objects very readable. Let's create the iconic psi state (ket):


In [5]: psi = Ket('psi')

In [6]: psi
Out[6]: |ψ>

Next let's create an operator and dagger it (and leave it unevaluated):


In [9]: a = Operator('A')

In [10]: Dagger(a)
Out[10]:
†
A

(the dagger symbol normally appears as a superscript to A, but Blogger is temperamental)

Inner products can have objects (mainly operators) between them now. Outer products will not have this functionality because it is not valid for them. When someone creates an inner product, they must instantiate the class because states and operators (and such) cannot automatically combine with their __mul__ method due to a very tricky reason (the subtlety lies deep within sympy's Mul class):

In Part 1 of this series of blog posts, I gave what I believed to be the prerequisites to understanding the mathematics behind the Risch Algorithm (aside from a basic understanding of derivatives and integrals from calculus). In this post, I will elaborate on what is meant by “elementary function,” a term that is thrown around a lot when talking about Risch integration.

The usual definition of elementary function given in calculus is any function that is a constant, a polynomial, an exponential (, ), a logarithm (, ), one of the standard trig functions or their inverses (sin, cos, tan, arcsin, arccos, arctan, etc.), and any combination of these functions via addition, subtraction, multiplication, division, taking powers, and composition. Thus, even a function as crazy as is elementary, by this definition.

But for the rigorous definition of an elementary function, we must take into consideration what field we are working over. Before I get into that, I need some definitions. Suppose that is the field we are working over. You can imagine that , the field of rational functions in x with rational number coefficients. As with the previous post, imagine as a function, for example, . Let be a differential extension of . We have not defined this, but it basically means that our derivation works the same in as it does in . You can imagine here that .

We say that is a primitive over if . In other words, the derivative of is does not contain , only elements of . Obviously, by the definition of a derivation (see the last post in the series), any element of is a primitive over , because the derivative of any element of a field is again an element of that field (you can see this by the definition of a derivation, also given in the last post). But also if for some , then is a primitive over , because .

We say that is a hyperexponential over if . Written another way, for some . We know from calculus that the functions that satisfy differential equations of the type are exactly the exponential functions, i.e., .

The last class of functions that needs to be considered is algebraic functions. I will not go into depth on algebraic functions, because my work this summer is only on integrating purely transcendental functions. Therefore, the only concern we shall have with algebraic functions in relation to the integration algorithm is to make sure that whatever function we are integrating is not algebraic, because the transcendental algorithms will not be valid if they are. Hopefully in a future post I will be able to discuss the Risch Structure Theorems, which give necessary and sufficient conditions for determing if a Liouvillian function (see next paragraph) is algebraic.

Now, we say that a function is Liouvillian over if is algebraic, a primitive, or a hyperexponential over . For to be a Liouvillian monomial over , we have the additional condition that . This just means that we cannot consider something like over as a Liouvillian monomial. Otherwise (I believe) we could run into undecidability problems.

We call a logarithm over if for some , i.e., . We call an exponential over if (or ) for some , i.e., . Note the difference between an exponential monomial and a hyperexponential monomial.

We can finally give the rigorous definition of an elementary extension. is an elementary extension of if there are such that and is elementary over for all . An elementary function is any element of an elementary extension of with the derivation . A function has an elementary integral over if there exists an elementary extension of and such that , i.e., .

Usually, we start with , the field of rational functions in x with rational number coefficients. We then build up an elementary extension one function at a time, with each function either being a logarithm or exponential of what we have already built up, or algebraic over it. As I noted above, we will ignore algebraic functions here. We generally start with because it is computable (important problems such as the zero equivalence problem or the problem of determining certain field isomorphisms are decidable), but the above definition lets us start with any subfield of .

Now you may be wondering: we’ve covered algebraic functions, exponentials and logarithms, and obviously rational functions are elements of , but what about trigonometric functions? Well, from a theoretical stand point, we can make our lives easier by noticing that all the common trigonometric functions can be represented as exponentials and logarithms over . For example, . You can see here that all the common trig functions can be represented as complex exponentials or logarithms like this. However, from an algorithmic standpoint, we don’t want do convert all trig expressions into complex exponentials and logarithms in order to integrate them. For one thing, our final result will be in terms of complex exponentials and logarithms, not the original functions we started with, and converting them back may or may not be an easy thing to do. Also, aside from the fact that we have different functions than we were expecting, we also will end up with an answer containing , even if our original integrand did not.

Fortunately, the integrating tangents directly is a solved problem, just like integrating algebraic, exponential, or logarithmic functions is solved. We can’t integrate functions like or directly as monomials like we can with or , because the derivatives of sin and cos are not polynomials in their respective selves with coefficients in . However, we can use a trick or two to integrate them. One way is to rewrite and proceed to integrate it as a tangent. Another alternative is to write . This function is algebraic over , but if we do not already have in our differential extension, it is transcendental, and we can rewrite it as (this is used in Bronstein’s text, so I believe what I just said is correct, though I haven’t verified it with the structure theorems just yet). These both work using the relevant identities for sin too. Of course, there is still the problem of rewriting the final integrand back in terms of sin or cos. Otherwise, you will get something like instead of for . Bronstein doesn’t elaborate on this too much in his book, so it is something that I will have to figure out on my own.

The second option I gave above leads nicely into the main point I wanted to make here about elementary functions. Notice that everywhere in the definitions above, things depend on the field we are working in. Therefore, cannot be an elementary extension over , but it can be over . Also, the error function, defined as cannot be an elementary extension over , but it can over . In fact this is how we can integrate in terms of some special functions, including the error function: by manually adding (or whatever) to our differential extension. Therefore, the usual definition of an elementary anti-derivaitve and the above Risch Algorithm definition of an elementary integral coincide only when the extension consists only of elementary functions of the form of the usual definition (note that above, our final fields are and , respectively).

Originally, I was also going to talk about Liouville’s Theorem in this blog post, but I think it has already gotten long enough (read “I’m getting tired”), so I’ll put that off until next time.

As I said in the last blog post, measurement has been implemented as a one-shot irreversible function in my code; I made it a function because measurement should not distribute to each orthogonal state inside the quantum state of the system. This could not easily be implemented in Sympy as every object (including my gates) is assumed to be distributive. This leads to some generality problems as one cannot construct a general circuit object with a mix of measurement gates and Unitary Gates. I did a one-shot measurement rather than an ensemble measurement as the one-shot measurement was needed immediately for my implementation of Shor's algorithm.

This week, I fixed these problems. In order to ensure measurement did not distribute, I had to subclass Gate into a NonDistributiveGate class which would be given special logic inside my apply method. This required me to carefully review apply and make sure that the .expand() method was used judiciously so that integer factors combined correctly while ensuring that gates did not distribute to the individual orthogonal states. This took much longer than I thought it would. I have also begun to implement a notion of ensemble measurement which will return a tuple of measurement outcomes.

Another thing I have been putting off is the creation of a set of quantum logic which implements common arithmetic quantum mechanically; this logic would be key in doing a legitimate and realistic simulation of all processes needed to do Shor's or Grover's algorithms. Making this logic required a look through the literature on reversible logic and a review of my old notes from my assembly language class. Reversible logic requires some junk bits which store extra information; since I am simulating a quantum computer and not a purely classical reversible machine, I can cheat a little and use measurement to irreversibly destroy information which bring down the number of necessary junk bits and the logic needed to manage them. Thus, I was able to implement a quantum mechanical Add, Bitshift and Multiply. I still need to create logic which can take powers and preform the mod function. Also, I might want to consider optimizing the number of gates and junk bits needed for these operations as they require quite a few.

The addition of tests to my code has been going smoothly as well. If I want to merge this code by the end of GSoC, I will need to begin to base my code off Matt Curry's code and document-test my code like crazy. These will likely be my primary concern for the next few weeks.

---

  Everyone should learn from their mistakes, right? Of course ;). Well SAT solvers are no different. If the solver discovers that the problem has no solution when a subset of the variables are set a certain way, then that is useful knowledge -- this knowledge is captured in the form of new clauses that are added to the theory.

  Now you may think that since we're doing a systematic search, we won't ever need that information again: for example if x=True / y=True / z=False led to an inconsistency and we backtracked to try other values of z, y and x, will we ever use the fact that (!x | !y | z) must hold? When using a static variable ordering, the answer is no. This is because if we backtrack, we will never see the same setting of variables again that led to the inconsistency. However static variable orderings are brutal :p.

  Let's expand on our example to see clause learning in a little more detail. Imagine the SAT solver went as follows:

• Set x = True
• From unit propagation set v1 = True and v2 = False
• Set y = True
• From unit propagation set v3 = False, v4 = True, v5 = False
• Set z = False
• From unit propagation we get the empty clause !!conflict!!

  So what do we know? Well right off the bat x = v1 = y = v4 = True and v2 = v3 = v5 = z = False is enough to make the theory unsatisfiable (so solution exists with those settings). But remember that unit propagation is a sound procedure -- the variables v1 - v5 are all implied by setting x and y to True [0]. So really just the settings to x, y, and z are enough to cause the problem. How could we have avoided it? Well what if we had the clause (!x | !y | z)? Note that this is just say at least one of the variables was set different. If we had that clause to begin with, our search would have looked like this:

• Set x = True
• From unit propagation set v1 = True and v2 = False
• Set y = True
• From unit propagation set v3 = False, v4 = True, v5 = False, z = True
• ...

  So what we do is add the clause, (!x | !y | z), to our database of clauses and continue the search. Since we're backtracking, we'll set z=True anyways (the clause wasn't really needed), but in the future (say if we backtrack higher and go down another branch, we may face the same decisions. This is when the clause really becomes useful.

  Now the method I just described uses a very simple technique to figure out the clause it should create / add -- it picks the negation of all the decisions you've made so far. However there are far more complex schemes that are ongoing research. I may incorporate these into SymPy some day, but for now the default works great ;).

  Another common concern (which I'll avoid dealing with for now) is how to detect useless clauses -- eventually learned clauses are no longer triggering and just taking up precious cycles in the solving process (and memory too). Elaborate cache flushing schemes have been proposed (and still are being actively researched), but SymPy likely won't be getting to the point where that would be necessary.

  So that's all for clause learning. It's extremely helpful in the context of SAT solving, and is finally in our SymPy solver.

[0]: I won't go into detail here, but you get the same result if you set a number of variables and then do unit prop as if you interleaved unit prop with each variable setting.

---

  So behind any reasonable SAT solver is a great heuristic. Variable State Independent Decaying Sum (VSIDS) is one such heuristic. What does a heuristic do? Well at each node in the search tree, it should tell you what variable to set, and what to set it to.

  VSIDS works by three main steps:

1. First it assigns a score to every literal (or variable setting if you will) matching the number of clauses that the literal appears in. So if a literal is in a large number of clauses, it will be picked before one that isn't.
2. Every time a clause is added, the literals in that clause get their score bumped up (by 1). Every time a clause is deleted, the scores get reduced (by 1).
3. At regular intervals, /every/ literal has their score reduced by a constant factor (divide by two is typical).

  So what does this accomplish? Well at a high level it prefers variable settings that will satisfy the majority of clauses, and over time (through steps 2 and 3) it will prefer literals that will satisfy a large number of learned clauses. This has the effect of avoiding bad parts of the search space where a solution is unlikely to exist.

  The SymPy solver now has VSIDS driving the search, with the one step of #3 missing (coming soon ;)). The beauty of VSIDS is how nicely it plays with the clause-learning of the SAT solver to make things exceptionally fast.

---

  That's all for now. I'll have another post soon showing some recent results, and one more within the following week to describe the data structures that make all this possible.

---

As I mentioned in my previous post, I did not want to implement array arguments for the C code generator by simply taking the same route as I did with the Fortran generator. In hindsight it was clear that I could do a much better implementation, and I’ve spent the week figuring out exactly how it should be done, and putting together a working implementation.  By the way it’s been beneficial to do the C and Fortran implementations in parallel, as it provides a good opportunity to test out alternative implementations.

The principal weakness of the Fortran code generator is that the support of array arguments has only limited testability. It is currently testable only on two levels:

1) The properties and behavior of IndexedElement instances,

2) The output from the Fortran code generator.

The missing piece in there, is the testability of complex expressions composed of IndexedElement instances. It is important to be able to test the consistency of indices in an expression, as it is very easy to create ambigous or even meaningless expessions of IndexedElement objects. For instance, the innocent looking expression: (which is considered undefined in my development branch)

x(i) + x(j)

would be ambigous at best. It could be (1) a scalar being the sum of components i and j of array x or (2) a symmetric rank 2 object a(i, j) where every component is made by adding components i and j of x. To avoid this kind of ambiguities, it is necessary to implement checks of index consistency. And it is preferable to implement it in the indexed module so that these checks are available whenever the module is in use. This gives the bonus of another level of testability:

3) The contraction structure of indexed expressions and the algorithms used to analyze the expressions.

With this in place, the reliability of array based code from the C code generator depends directly on how robust the analysis of indexed expressions is.   The third level of testability is extremely important to make Sympy generated code a viable option for scientific applications.

In my branch codegen_C on github, these thoughts have materialized in a new file sympy/tensor/indexed_methods.py, as well as in modifications of the C code printer.   I have implemented testable high-level interfaces that can analyze expressions and provide information about the index structure. The code printers should now be able to generate loops over arrays in a much more reliable fashion.

After last week’s breakthrough, work this week has been very slow. I started working on the Parametric Risch Differential Equation Problem, which is almost identical to the Risch Differential Equation Problem in how it is solved, except there are a few extra steps. Unfortunately, because it is so similar, Bronstein breezes through the description. This is fine for the parts that are the same, but he is a little unclear on how some of the new parts fit in. Also, his pseudocode has a line more or less saying

if r1 = … = rn = 0 then
    N = -1
else
    N = max(deg(r1), …, deg(rn))

for i from 0 to N
    for j from 1 to m
        Mij = coefficient(rj, t^i)

where M is a matrix. It is not very clear what this is supposed to mean in the case where N = -1. Obviously, you can’t have a a matrix with negative dimensions. Clearly, this means that this particular function doesn’t apply somehow in this case, but I am not really even sure where it fits in to the whole algorithm at this point in reading. After reading a few more pages in, it gives a few hints here and there on how it is to be used, but never is it explicitly shown, in pseudocode or otherwise. So for now, I think my best bet is to read ahead and get a fuller understanding of the complete function before I try implementing anything (this is what I had been doing before, but I caught up to myself).

Also, on an unrelated note, I just found out today that I passed my Google Summer of Code midterm evaluation. This means that I will receive half of my stipend for the program (the other half comes after passing the final evaluation at the end of the summer), and that I can continue working on my project in the program.

EDIT:

Later in the text, it runs through an example and says “… , hence M and A are 0 by 0 matrices.” So obviously, that is what was meant.

This was a good week. I coded some stuff that wasn't nearly as abstract as Hilbert spaces, and my midterm was completed.

Overall, I would say that I'm quite satisfied with my progress. We are still working on the many products between quantum objects, but that doesn't mean I can't work on the objects themselves or functions. As we saw in my previous post, basic symbolic quantum mechanics is coming together and it is very nice to see so.

It seems that when we try and implement products between these objects (inner and outer products included too now unfortunately), many difficulties arise in the design decisions. That's okay though, sometimes things just catch us off guard sometimes. I know that in a few weeks everything will be nicely glued together (in a tensor product sense).

So far this week I have been playing around with learning and implementing Quantum Algorithms. Since I have already created the framework to build a Quantum Fourier Transform (QFT), I decided to start work on implementing Shor's algorithm for prime factorization which uses the QFT to amplify certain select probabilities. This makes for a good demonstration that everything built up until now is working correctly.

Here's a quick overview of what the algorithm needs to factor a number. It starts by picking a random number, a, less than N and makes sure that it is co-prime to N using Euclid's (classical) Algorithm. A superposition of many states |0>...|k>...|2**n> is made using Hadamard Gates (where n is the number of qbits in a quantum register). It then does the controlled-mod operation (a**k % N) of this register and stores the result in a second register. The second register is measured, collapsing the state of the first register into only the values which result in a particular value for the second register (this is due to the first and second register being entangled by the controlled-mod operation). We can take the QFT of this to derive the order of a in modular N arithmetic. This can be used to derive prime factors of N.

As we can see, the algorithm requires me to add a few new things to my code. First, it requires a system of measurement. To start, I have made measurement a function that takes in a state and the qbit it wishes to measure. Eventually, I will make this a Gate object. In addition, a controlled-mod gate is needed to make the algorithm work; eventually, I will want to develop a way to build this out of elementary gates, but for now I have simply overwritten the apply method for a controlled mod gate that knows how to directly apply itself to a state.

For now, I am applying the Quantum Fourier Transform gate by gate, which makes it rather slow (5-qbit long numbers like 21 take way to much time and RAM). I plan on making this a gate in and of itself that knows how to apply itself using a FFT, but also knows how to decompose itself into the individual gates that make it up. That way, user can choose to apply gates in a way they see fit.

Here is an example of my code factoring 15:

In [14]: shor(15)
a= 11
N= 15
Hadamard Gates Applied
controlled Mod'd
measured 2
QFT'd
measured 1
|1000000000001011>
r= 2
Out[14]: (5, 3)
#comments such as "Hadamard Gates Applied" are so that one can keep track of where in the algorithm we are

All in all, this week is going quite well.

It’s been an exciting week.  My work on Fortran code generation made it into the official master branch of Sympy.  This includes printing of free-form Fortran, a brand new FCodeGen object, and a module for indexed objects.  The indexed objects allows generation of code involving arrays.

Also, Brian introduced me to the Python library called Theano.   It “…allows you to define, optimize, and evaluate mathematical expressions involving multi-dimensional arrays efficiently.”  Now, that description could very well have been a statement about my doings so far this summer.  The Theano project has an emphasis on optimization, and is also able to perform calculations on a GPU if available, so it would be very nice if Theano could be used together with Sympy.   After a little discussion with Brian and Ronan I was inspired to examine the possibility to use Theano as a back-end.

The expression tree structures in Theano and Sympy differ slightly, so I wrote a little module to convert the tree of a Sympy expression to a data structure with identical topology as the corresponding Theano expression.  The idea is that expression trees on equivalent form can be converted between Sympy and Theano by a trivial mapping of the nodes.  However, I was not able to follow this track further as I had to prioritize the review process of the Fortran code generation.

I’ve spent the rest of the time to bring the C generator up to par with the Fortran facilities.  The current result is here, but it is not ready yet as the tests of array arguments are still XFAIL-ing.   I don’t want to quick fix this by copying the approach I used for Fortran arrays, because there is need for a more robust implementation. This demands improvements in the indexed module, in particular functions to inspect and validate expressions.  The code printers should rely on these functions to determine the correct implementation of the mathematical expressions.

The indexed module has become a central piece in my project as the way to represent arrays in Sympy.  I think that the Indexed objects can be useful not only for generating loops in C and Fortran as I have done so far, but also for interfacing with other array-based systems, such as Theano, and also directly with Numpy.  So  improvements to the indexed module would benefit the interaction with all these array-based systems, and I should give it a high priority.

Yesterday I came back home from Sage Days 23 in Leiden. I don't really have time to write a detailed travel report right now, but overall it was very, very nice. Those interested can check out William Stein's photos (1, 2, 3, 4) of the event (I can be seen in a few of the photos, wearing a blue t-shirt during days 1-3 and a yellow t-shirt during day 4).

I didn't get quite as much coding done as I probably had time for (in part because I had a talk to prepare, and in part because I spent a bunch of time playing with MPIR and FLINT 2 and gaining some insights but otherwise not accomplishing much productive).

The code I did write includes a complete Cython implementation of real exp (giving up to a 7x speedup) and complex sqrt for mpmath, as well as code for fast evaluation of the gamma function at rational arguments. I will submit the Cython code to Sage soon; possibly after rewriting the remaining power code in Cython as well.

Anyway, prompted by a recent question on sage-support (and because I finally got inspiration to do it while at the airport), I've now implemented Bessel function zeros in mpmath.

Bessel function zeros are very important in applied mathematics because they are needed for Fourier–Bessel series, used to solve partial differential equations in cylindrical coordinates. I'd like to write down a concrete computational example, but it'd be slightly involved and I don't really have time right now. Perhaps in a later blog post.

The code permits computing the m-th nonnegative simple zero of J_ν(z), J'_ν(z) Y_ν(z), or Y'_ν(z) for any positive integer index m and real order ν ≥ 0. One can use it for large m, large ν, or both:

>>> besseljzero(1,100)
314.9434728377671624580656
>>> besseljzero(100,1)
108.836165898409774363098
>>> besseljzero(100,100)
459.5295465754674695983379

A word of caution: I've chosen the convention Abramowitz & Stegun of including z = 0 for J'₀, which is incompatible with SciPy and the Mathematica BesselZeros package. The A&S convention makes much more sense to me because it makes the zeros a continuous function of the order. This case (ν=0, derivative=1) should be the only incompatibility.

It's a bit tricky to compute the zeros of Bessel functions, because (as far as I know) no globally valid approximations are known; not any reasonably simple ones anyway.

Asymptotic expansions are known with respect to either the order or the index. In my implementation, I've chosen to simply use the asymptotic expansion with respect to the index n. If the n is too small, the code simply chooses a larger index m for which the asymptotic expansion does work, and then isolates all the zeros between 0 and the m-th by counting sign changes. By caching this result, all the consecutive zeros for the same function order can then be computed quickly as well.

This approach is simple and quite robust (up to any silly implementation errors on my part), although the code could certainly be optimized more for large orders.

The following graphic shows the derivative of Y₅₀(z) with its zeros marked:



A neat property of Bessel functions is that they can be expanded as infinite products over their zeros:



These products converge slowly, but work with the aid of convergence acceleration:

>>> v, z = mpf(3), mpf(0.75)
>>> (z/2)**v/gamma(v+1)*
...     nprod(lambda k: 1-(z/besseljzero(v,k))**2, [1,inf])
...
0.008484383423274108843927552
>>> besselj(v,z)
0.008484383423274108843927552

>>> (z/2)**(v-1)/2/gamma(v)*
...     nprod(lambda k: 1-(z/besseljzero(v,k,1))**2, [1,inf])
...
0.0331364636065541211306414
>>> besselj(v,z,1)
0.0331364636065541211306414

I will do more work on Bessel functions in mpmath and Sage, so watch for future updates.

<object height="385" width="480"><param name="movie" value="http://www.youtube.com/v/6xxNEjwMbKw&amp;hl=en_US&amp;fs=1&amp;hd=1"/><param name="allowFullScreen" value="true"/><param name="allowscriptaccess" value="always"/><embed allowfullscreen="true" allowscriptaccess="always" height="385" src="http://www.youtube.com/v/6xxNEjwMbKw&amp;hl=en_US&amp;fs=1&amp;hd=1" type="application/x-shockwave-flash" width="480"></embed></object>

My VPS host just upgraded their servers to Ubuntu 10.04, so I just spent today saving all my old stuff, wiping my VPS, and reinstalling everything.  I was using nginx as my webserver previously, but now I’m using lighttpd, and it was pretty easy to setup.  I’m still ironing out the bugs, but I was able to get most everything up in running in several hours, so I hope to have the rest dialed in very soon.

So, although several other SymPy projects were accepted by Google's Summer of Code program, my own proposal unfortunately wasn't. My plans for summer 2010 have been updated as follows: unless my application to the University of Washington's Inverse Problems program (a summer research opportunity for undergraduates) is accepted, I'll be back home in California doing just what I said I would in the project proposal. I put my chances of being selected by the Inverse Problems people at roughly 20%, so I'll most likely be working on SymPy.

The Google Summer of Code proposal for this project (actually, proposals - I submitted a copy each to Portland State University and the Python Software Foundation) has been submitted.

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

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.

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.

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

This project is officially winner of the "Concurso universitario de software libre de Granada" (something like college open source software contest of Granada)

To celebrate it, we've made a little video showing the design process

<object height="300" width="400"><param name="allowfullscreen" value="true"/><param name="allowscriptaccess" value="always"/><param name="movie" value="http://vimeo.com/moogaloop.swf?clip_id=3326096&amp;server=vimeo.com&amp;show_title=1&amp;show_byline=1&amp;show_portrait=0&amp;color=&amp;fullscreen=1"/><embed allowfullscreen="true" allowscriptaccess="always" height="300" src="http://vimeo.com/moogaloop.swf?clip_id=3326096&amp;server=vimeo.com&amp;show_title=1&amp;show_byline=1&amp;show_portrait=0&amp;color=&amp;fullscreen=1" type="application/x-shockwave-flash" width="400"></embed></object>
Empathy project from fabian seoane on Vimeo.

And yes, I promise that I'll put this app online soon.

Today has been a day for doing some design tweaks with my partner Angel Soler, just to give empathy a more dynamic look:



It's been also very interesting discussing with Ondrej and other members of the SymPy community the best way to output formulas. Proposed solutions so far are:

1. jsMath: this uses sympy's latex output and transforms it using javascript into something the browser can read (mathml or images, depending on the browser). Disadvantages: you have to load an entire javascript library, which takes quite a lot of time, and CPU goes up to 90% during 5 seconds.
2. MathML: sympy supports mathml output natively and mathml is a universal language for displaying formulas. This would be the ideal solution, if it wasn't because of ... 1. Sympy outputs "Content" MathML, whereas firefox (and other browsers) knows only how to display another variant of MathML, "Presentation" MathML. You can transform Content MathML to Presentation MathML (and sympy knows how to do it), but it is an expensive opperation (maybe too expensive to be done on every request), and 2: MathML support in browsers is partial. To view MathML in IE you have to install a plugin (MathPlayer) which, by the way is not open source.
   
3. Image. This would be a nice solution, and is the one chosen by wolfram's integrator as the main output. Sympy's png printer is a bit slow, since it does the following: expression -> latex -> dvi -> png.
   
4. ASCII (pretty). This is the easiest solution. It renders ok (although i'm having some issues with wrapping), but is not as pretty as the image output.
   

This is the first release to include SymPy Plot.

0.4.1 source tarball
0.4.1 windows installer
SymPy homepage
Changelog

PyOpenGL is required for SymPy plotting support. While *nix users should have no problem (apt-get install python-opengl), Windows users will find it difficult, if not impossible, to use PyOpenGL with Python 2.5. I recommend setting up your Python environment as follows:

1. Install Python 2.4.4 to C:\Python24 (the default location)
   I normally like to install all of my programs in C:\Program Files, but unfortunately PyOpenGL depends on the default install path.
2. Install Numeric-23.7.win32-py2.4.exe
3. Install numarray-1.1.1.win32-py2.4.exe
4. Install PyOpenGL-2.0.2.01.py2.4-numpy23.exe
5. Install sympy-0.4.1.win32.exe

6. >>> from sympy import Symbol
   >>> from sympy.modules.graphing import Plot
   >>> x = Symbol('x')
   >>> Plot(x**2, [x, -1, 1, 10])

   

If you get a dll not found error relating to GLUT, obtain glut32.dll and put it in C:\Python24\Lib\site-packages\OpenGL.

I've written a final SoC report and posted it in the SymPy wiki. Click here to read it. You will also find similar documents from other SymPy SoC participants here.

Well, it seems we're in the homestretch of the SoC. I've been working on the FFLU decomp. Strangely I ran into some problems with modules not being installed properly but that finally regulated itself. Now I'm in the debugging phase and hopefully we're have a nice new fraction-free LU function.

I'm thinking that maybe I can squeeze in another Fraction-free inversion algorithm before next week...

*edit* Did I say inversion? I meant QR. So LU seems to be working but my QR one is still buggy. Oh well.

So for whatever reason I kept saying block LU, when I really meant fraction-free LU. Those are, of course, significantly different. I've been looking at this paper suggested by someone on the mailing list. It seems reasonably straight-forward (even the proof, which is surprising). I hope to finish the LU one soon and hopefully the QR one before the project finishes up.