Sean Vig's 2011 GSOC Project

Developing Wigner-3nj Symbols in SymPy

Continuing GSoC Work

This last week, I have made progress on my project working on laying the base work for the spin states and in reimplementing logic in the cg_simp method for Clebsch-Gordan coefficients.

First, I have started the work on the implementation of coupled/uncoupled spin states. Currently, this is implemented by adding a coupled property to the spin states. This can be set to True for coupled, False for uncoupled or left as None for other states. As this evolves, I will move to having uncoupled product states be represented by a TensorProduct of two spin states. The next key will be establishing represent and rewrite logic for these spin states. Part of this will be figuring out how exactly these methods will work and what they will return. Namely, the represent method, noting that when representing an uncoupled state as a coupled state, it returns states with multiple j values, which under the current logic, would return matrices of different dimension. Also, we will have to determine what represent will do to uncoupled tensor product spin states. This next week, I will likely rebase this branch against the CG branch so I can start using the Clebsch-Gordan coefficients to implement these functions as the CG pull is finalized.

With the Clebsch-Gordan coefficients, this last week I was able to get the simplification of symbolic Sum objects working. I did this using the pattern matching built into sympy with Wild and .match. The final step with this should be to rework the logic of _cg_simp_add to make it easier to add in additional symmetries.

Transitioning to Spin States

This last week, the first thing that was taken care of was finishing the x/y/z spin basis representation. Having fixed the Wigner small d-function in the Rotation class, the tests for this were put into the pull request and the pull was merged into the sympy master, making this my first pull request since starting GSoC. There is still some changes that will come for the Rotation class, namely the creation of a symbolic WignerD class which is returned by the current Rotation.D and Rotation.d functions, but that will be dealt with in a later pull request.

With the x/y/z basis stuff finally out of the way, I moved back to getting the Glebsch-Gordan coefficient/Wigner-3j symbols to a state where they can be pulled. Having fallen behind in getting the CG coefficient simplification to a suitable state and with the work on the x/y/z spin basis pushing the timeline back even further, the current goal is to merge what I have so far and move on to the coupled spin states. What I have so far is the classes for the Wigner-3j symbols and the Clebsch-Gordan coefficients which can be manipulated symbolically and evaluated, and a very rough version of the cg_simp method. Currently, this method can handle 3 numerical simplification, however the code is still messy and having more cases would be ideal. That said, in an effort to make sure the key parts of this GSoC project are covered, I’ll be moving into writing the coupled spin states.

For the spin states portion of this project, I will develop a means of writing coupled and uncoupled spin states. The uncoupled product basis states will be written using the TensorProduct, which is in the current quantum module; each of the states in the tensor product will be states as they are currently implemented. To represent coupled basis spin states, the current spin states will be modified to included a coupled parameter. This value stores the J_i’s of the spin spaces which are being coupled. In addition to the spin states being implemented, methods will be written to utilize the CG coefficients mentioned earlier to go between coupled and uncoupled basis representation. Look for more next week as this code is fleshed out.

More CG Simplification and Wrapping Up X/y/z Spin Bases

For the first part of this last week, I continued on my work to get sums of Clebsch-Gordan coefficients to simplify. Using the same general logic that I outlined in the last blog entry, besides general cleaning up of the code, most of the work at the beginning of the week was spent on trying to develop a function that could check an expression for CG coefficients matching a set of conditions.

The rationale behind trying to write such a method is that it would make it much easier to identify the times where symmetries could be utilized. With such a method, the process of checking for CG coefficients could be done in a single function and the logic for implementing CG symmetries could be handled in this one function. The current method uses lists of tuples to specify the conditions on CG coefficients. For example, if the j1 value of a CG coefficient needed to be =0, you could pass the tuple (“j1”,0), or if the m1 and m3 values match (“m1”,“m3”). All the conditions for each CG coefficient are combined into a single tuple. The current snag with simplifications of sums of products of CG coefficients. For example: $$ \sum_{\alpha,\beta} C_{\alpha\beta}^{c\gamma} C_{\alpha\beta}^{c’\gamma’} = \delta_{cc’} \delta_{\gamma\gamma’} $$

While the current method would be able to check for specific values on the CG coefficients, I have yet to come up with a good way to check that the m1 and m2 values are the same when they can take any value, as in this example. As it stands, this code still seems like quite a hack and will need some work before it is good to go.

What is left with this part of the project is:

  • Getting simplification to work with sums of products of terms (as in the example above)
  • Applying CG symmetries to perform simplifications
  • Simplification of symbolic CG sums
  • Fixing up the printing of CG terms
  • Final testing/documentation

This part of the project has unfortunately fallen behind the preliminary schedule by a bit, as it was due to be finished up last week. I’ll outline what I’m currently working on finish up next, but hopefully I can finish the CG stuff ASAP so I can move on to working on the spin stuff which is the true meat of the project and try to get back on schedule.

After meeting with my project mentor, Ondrej, on Wednesday, it was decided that the focus would shift to finishing up the work I’d started on x/y/z spin bases and representation of spin states that I’d started before GSoC had officially started.

The first order of business was identifying an error in the Wigner small-d function, which is used extensively in the changing of spin bases. With Ondrej noting that the small-d function was defined only on a small interval and then me discovering the bug in the Rotation.d method, we were able to address this. However, no sooner had this been done than Ondrej is able to work out a better equation for the small-d function, which will likely replace the current implementation.

Other than this, most of the work this week on the x/y/z basis representation was in documentation, testing and generally cleaning up the code to be pulled. The current pull request (my first work to be submitted since the start of GSoC) is still open here. While this pull integrates the current work on basis representation, after this pull there is still some work that will need to be done testing both the Wigner small-d and the D functions, for both symbolic and numerical values, and ensuring they return the correct results. Because the representation code relies so heavily on these functions, it is imperative that these functions evaluate properly. Once these are fully tested, there will also likely need to be more tests to ensure all the representation code returns the right values for as many odd cases as would be necessary to test. Hopefully I can finish this up soon and move on to other work that still needs to be done.

More Improvements to the Simplification Method

I was out of town for the beginning of this week, so I don’t have as much to report, nevertheless in the last few days, I have made some good improvements to the cg_simp method, allowing it to simplify cases other than just the simple sums, tho it still only handles the same case as before involving the sum of . While it is entirely possible I’m doing something stupid in performing all the checks, as far as I can tell, it works for the given case. Because the code itself is not yet very clear and it is not always straightforward what is happening, I’ll explain what I have implemented.

First note that for the simplifications that will be made, it is required to have a sum of terms. The first for loop in the method constructs two lists cg_part and other_part, the former consisting of all terms of the sum containing a CG coefficient and the latter the other terms.

Next, we iterate over the list cg_part. The number of CG coefficients is computed, which determines which simplifications can be made (currently, the only implemented simplification uses terms with 1 CG coefficient in each sum term). Those terms with the proper number of CG coefficients are then considered for the simplification. The way this will work is: based on the properties of the CG coefficient in the term, it will search the rest of the list for other terms that can be used in the simplification, and if all the terms exist, will simplify the terms.

Turning to the implementation, when iterating through the list, the first thing to do is determine the properties of the CG coefficient terms, that is to extract from the term in the sum the CG coefficient itself and the numerical leading terms. Here it is also noted the sign of these numerical terms.

Next, the rest of the list is checked to see if a simplification can be made using the determined term. To keep track of this, a list, cg_index, is initialized with False as the element of the list. In checking the later terms, we preform a similar decomposition as with the first term, that is splitting up the CG coefficient from the other components of the term, determining the CG coefficient, the leading terms and the sign of the terms. If the properties of these are correct, then the corresponding element of cg_index is updated with a tuple (term, cg, coeff) where term is the term in cg_part (so this element can be removed later), the CG coefficient and the leading numerical coefficient of the coefficient.

Now, if all the elements of cg_index are changed, the simplification is preformed. When this happens, first we find the minimum coefficient of the chosen CG coefficients, which determines the number of times we can apply the simplification. Then the replacing happens; for each element in cg_index (which is a tuple) the first element of the tuple is popped off cg_part, then, if the term is not eliminated by the simplification, a new term is created and added to cg_part, and finally a constant is added to other_part, completing the simplification.

Looking at the code, this method is very straightforward, but should be robust and scalable for treating cases of sums with numerical leading coefficients, and now that the i’s have been dotted and the t’s have been crossed on testing this method, implementing new cases should come rapidly in the next couple days. However, one place where this will still need some work is in implementing symbolic simplification, both in dealing with symbolic leading terms on the CG coefficients and symbolic CG coefficients themselves. This will take a bit of thought and likely a bit of help to complete, but this is one thing I hope to work on in the next week. In addition, as the simplification comes into place, I’ll work on polishing out the last of the details to get the classes for the Wigner3j/CG coefficients working properly.

Implementing Clebsch-Gordan Symmetries and Sum Properties

In this first week of the GSoC project, I focused on implementing methods that would simplify terms with Clebsch-Gordan coefficients. This still has a long was to go, but I will outline what I have done so far.

The first step was implementing means of dealing with sums of single coefficients. This would hopefully look something like:

>>> Sum(GC(a, alpha, 0, 0, a, alpha), (alpha, -a, a+1)

The first implementation of this used an indexing system that was able to index single coefficients, which could then be processed. This allowed the simplification function to act properly in simple numerical cases, so it could do things like:

>>> cg_simp(CG(1,-1,0,0,1,-1)+CG(1,0,0,0,1,0)+CG(1,1,0,0,1,1)+a)

The problem with this implementation is doing something as simple as having one of the terms have a constant coefficient would break it. In addition, there would be no clear way to extend this to sums involving products of multiple Clebsch-Gordan coefficients.

To deal with this, I started working a solution that could deal with having constant coefficients and products of coefficients. Currently implemented is method which creates list of tuples containing information about the Clebsch-Gordan coefficients and the leading coefficients of the Clebsch-Gordan coefficients. Currently, the only implemented logic is only able to deal with the case in that could be dealt with in the previous implementation, however, this should be able to expand to encompass more exotic cases.

Another thing that was touched on this last week was treating symmetries. These are quite simple to implement, as they need only return new Clebsch-Gordan coefficients in place of old ones, just with the parameters changed in correspondence with the symmetry operation. The key will be using these symmetries to help in simplifying terms. This will be based on the development of better logic in the simplification method and the implementation of some means of determining if these symmetries can be used to apply some property of the Clebsch-Gordan coefficients that can simplify the expression.

I will be out of town this next week on a vacation, and will not be able to get work in, but I will continue working on this when I return, with the intention of getting it to a state that can be pushed within the next couple weeks.

Official GSoC Start

This week marks the official start of the Google Summer of Code. While I started getting my feet wet last week after finishing the last of my finals and grading, the bulk of the work has just started turning out. I’ll quick cover what I have from this last week and what I’m looking to get working this week.

Before the start, I worked out expanding the functionality of the currently implemented x/y/z bases, which work I have here. The previous implementation only allowed for evaluating inner products between states in the same basis and representing the states in the Jz basis and then only with j=½ states. Using the Wigner D-function, implemented with the Rotation class, I implemented represent to go between the x/y/z bases for any j values. Both _eval_innerproduct and _rewrite_as were then created to take advantage of the represent function to extend the functionality of the inner product and to implement rewrite between any bases and any arbitrary j values.

This seems like this is some documentation and tests away from being pushed, but there is something buggy with the Rotation.d function, implementing the Wigner small d-matrix. I noticed when trying to do

>>> qapply(JzBra(1,1)*JzKet(1,1).rewrite('Jx')

and I wasn’t getting the right answer. As it turns out, the Rotation.d function, which uses Varshalovich 4.3.2 Eq 7, does not give the right answer for Rotation.d(1,1,0,-pi/2) or Rotation.d(1,0,1,pi/2). Namely, there is something wrong with the equation that doesn’t change the sign of the matrix element when reversing the sign of the beta Euler angle. Running all four differential representations given by Varshalovich for the small d matrix, Eq 7-10, give the wrong result, so the derivations of these will need to be checked to fix this. I have a bug report up here.

As for what I will be implementing this week, I already have the basics of the Wigner3j and CG class implemented, the work for this going up here. This includes creating the objects, some of the printing functionality and numerical evaluation of the elements using the functions. The meat of the class that I’m currently working on is the cg_simp function, which will simplify expressions of Clebsch-Gordan coefficients. I currently have one case working, that is Sum(CG(a,alpha,0,0,a,alpha),(alpha,-a,a)) == 2a+1 which is Varshalovich 8.7.1 Eq 1. There are still some things to smooth out with the implementation, but I should have that worked out a bit better, in addition to some more simplifications by the end of the week.

That’s all I have for now, watch for updates within the week as to what I’ve gotten done and what I have yet to do.

The Best Things in Life Are Free (and Open Source)

Hi all,

This adventure into blogging is to document the work I will be doing this summer with SymPy, an open source symbolic mathematics library written in Python. The project will be done as a part of the Google Summer of Code program. This summer, I will be developing a symbolic class for creating Clebsch-Gordan coefficients and will develop the spin algebra in the existing quantum physics module to utilize these coefficients. The gory details can be read in my application. My mentor for this project will be Ondřej Čertík. The project officially starts May 24, but I’ll be diving in once I finish up the last of my finals May 14. I will be documenting my progress on the project throughout the summer through this blog, and anyone interested in this is free to watch for updates once the project is underway. All my work will be pushed to my SymPy fork on github.

That’s enough for now, I’ll be checking back in once I get started in the summer. Now, back to trying to graduate.