Incorporating Wavefunction and DifferentialOperator into representations

This post is a bit late because I’ve been experiencing some computer troubles recently.

Since my last post, the represent branch has been merged into sympy master! There are still some issues with representing continuous bases that need to be ironed out, so I have started a represent2 branch to deal with some of these changes.

The first commits in the branch are some simple fixes to docstrings and cleanup of the interfaces of Wavefunction (Wf) and DifferentialOperator (DO). Of note are the addition of an _eval_expand to Wf and DO, so that we can do simplifications like expand(Wavefunction((x+y)**2, x, y)) == Wavefunction(x**2, x, y) + Wavefunction(y**2, x, y) + Wavefunction(2*x*y, x, y).

The biggest tasks for represent2 are to incorporate the new Wf and DO classes into the internal represent logic of classes in quantum. Another big task is to simplify the logic of represent. In particular, rep_innerproduct and rep_expectation will become helper functions rather than being called from the main represent logic. This simplifies the main logic, but retains these methods to be called from the internal _represent methods. This is important because these functions know about proper indexing to form representations in these very standard forms, and this means we don’t have to reproduce that logic internally.

Both of these tasks are proving to be quite tricky and time consuming. I hope to have them done in the next couple of days. If that is the case, I expect to be able to rapidly move on to implementing coordinate systems and finally getting some example textbook quantum systems up and running. With only a few weeks of GSoC left, I’ll have to work extra hard to make that happen. Modifying the represent logic is taking much longer than expected, but we’re dealing with very complicated logic here and we want to make sure its robust. I’m confident, though, that I can have a nice finished product by the end of GSoC.

Summarizing the changes in represent

As the pending pull request gets closer to being merged, I thought it might be helpful to myself and others interested in the code if I summarized the changes that I have made to the code in represent.py. Hopefully people can use this post as a reference as they start digging into that portion of the code in reviewing the PR.

Represent has behavior to handle almost any kind of expression, from Add, TensorProduct, or Mul to the basic QExpr. Most of the more complicated types of expressions, like Add, recursively call represent on smaller atoms of the expression, hoping to ultimately get to a single QExpr, which is a sort of termination point in the recursion. The logic for handling a single QExpr is what I mainly modified in my work. Previously, the QExpr handling was a simple attempt to call expr._represent. Now, however, the following sequence occurs:

• First, call get_basis. This function exists to unify different ways of calling represent. Whether represent is called with an operator or state class/instance as the basis in options, get_basis will examine this and try to return a state instance. At this point, the replace_none option is set to False. This is in place because represent_default_basis will only be called from expr._represent if there is no basis specified. Since the basis that get_basis returns may not be the one called by represent_default_basis, we want to make sure that if the basis in options is None, represent_default_basis gets called.
• After setting the basis appropriately, call expr._represent. If this is successful, return the result
• If a NotImplementedError (NIE) is raised, we set the replace_none option to True so that any future calls to get_basis will actually fill a basis into the options. At this point, any attempted call to represent_default_basis has failed, so we should now try representing in a basis state that we know.
• If the passed QExpr is a Ket or Bra, we then call rep_innerproduct, which will return an innerproduct (<x’|x>) type representation for the given Ket
• If the passed QExpr is an Operator, we call rep_expectation, which will return an expectation value like representation (<x’ | A | x>) for the given Operator

In addition to QExpr handling, there were some additions to the processing of Muls. Originally, represent is called on each arg of the Mul individually. While this still happens, there is code added to keep track of the unities that have been inserted and the current index that should be being appended to dummy kets. Details of where the unities are inserted are in this previous blog post.

Finally, before the result is returned, integrate_result is now called as well. If the expression passed to it is an Expr (presumably a continuous result), then it integrates over any unities (e.g. |x_1><x_1|) which were inserted into the quantum expression. This collapses some of the Delta functions that were originally in the expression.

To see the latest represent, check out the version in my branch here

A week of cleanup

This week was cut short by storms in the Midwest again. I lost power and internet from Monday until Wednesday, so I was stuck trying to find local cafes that still had power and also had internet. I did manage to have a very productive conversation with Brian about things that should still be cleaned up. I will post a more detailed description of the changes once I’ve finished the todo list we came up with, but for right now, some general themes that I’ve been working on include:

• Docstring cleanup throughout the PR: as this is my first time working with Sympy, I’m still figuring out how to write good docs that conform to the standards throughout the code base. As such, alot of work will go into writing better descriptions and examples, as well as cleaning up the formatting.
• More test coverage: There are many tests that have been added, but its time to start thinking more generally and come up with more complicated tests. Examples of tests that were added include support for Bras in the operator <-> state mapping as well as testing multi-dimensional Lapalacians in DifferentialOperator.
• API streamlining: the API for DifferentialOperator described in my last post was a bit cludgy, so we’ve now reduced it to be just the general expression form. This might be more verbose, but it also simplifies the internal code quite a bit. In Wavefunction, I’ve cleaned up the normalization a bit, as well as adding a normalize() function which will return a normalized Wavefunction. Things like _eval_dagger and _eval_conjugate are defined now as well.

These changes were mainly to operatorset.py, Wavefunction, and DifferentialOperator. You can see some commits here, here, and here. There’s still a bit left to do, mostly with docstrings and test coverage, and I will talk more about this later. This week was also the GSoC midterm evaluation, so I will try to post a summary of the first half and the changes that I have made as well.

Revamped DifferentialOperator and pull request cleanup

The first task that I handled this week was to generalize DifferentialOperator so that it can handle more complicated expressions like a Laplacian or something similar. Its constructor now operates in three different modes: (see the commit here)

1. A single symbol or string: this represents the variable to differentiate with respect to. This a simple type operator.
2. A tuple with two entries: the first entry of the tuple is the same as in 1), but the second is an integer representing the order of the derivative to be taken. Therefore, DifferentialOperator((“x”, 2)) would represent
3. Two arguments, an Expr and a tuple: the first entry is a general expression involving a Function. The second argument is a tuple. The first entry in the tuple is the symbol of the Function in argument 1 which is to be replaced by the Wavefunction we apply this operator to. The second entry in the tuple is the symbol which the function is evaluated with before replacement. So, another way of expressing a single derivative would be DifferentialOperator(Derivative(f(x), x), (f, x)) were f = Function(‘f’) and x = Symbol(‘x’). This allows us to put any arbitrary expression in and have it applied to the Wavefunction accordingly.

Once the new DifferentialOperator was implemented, it was time to add this functionality to cartesian.py. As a result, it is now possible to represent X operators in the momentum basis and vice versa!

In addition to the DifferentialOperator changes, this week was about coming up with more tests and trying to find bugs in corner cases of representations. One bug related to integration in sympy was discovered when testing out representations that involve DifferentialOperator. DifferentialOperator is non-commutative, and as_independent, which is called by the integration routine, would treat these operators as if they did not depend on the variable of integration even when they did. I opened issue 2549 here and then shortly after had a pull request here with a solution that eventually got merged in. This meant that DifferentialOperators were now integrated over correctly, so more tests were added that took advantage of this new functionality.

I had a brief discussion with Brian on IRC about outstanding issues in the PR, and one thing that we both agreed on was the the enumerate_states function of represent.py would not work for very general kets. The way it was previously structured was to simply take the labels of kets and append indices to the end of them. This works for continuous bases, and is useful when trying to form an inner product like , but it’s completely useless for kets of discrete bases. This commit factors the actual enumeration of the states out to an internal ket function, leaving represent.py’s enumerate_states as the public interface. I think this design is ideal, because this means that each ket can decide for itself how to handle a request for a certain number of kets from the basis.

Hopefully Brian and I can discuss any leftover outstanding issues with my PR soon and we can get it merged in. After that happens, I’m going to start implementing many different coordinate systems and lay the groundwork for being able to handle many more example systems. This point in my schedule is when I already wanted to be implementing the coordinate systems, so there’s a chance that I may fall a bit behind that schedule. Luckily, I incoporated a 3 week buffer into the end of the timeline, just in case that happened, so I think I’m still on track to get a nice working product by the end of the summer.

Mapping states to operators and other pull request changes

This week has been a very busy and very productive week. I submitted my first big pull requests for the changes to representations that this blog describes, and have had many discussions with Brian about issues and changes to that code. The pull request is here. This post will try to summarize some of the issues that were raised and point to the commits that address those issues.

Originally, I had written an API in StateBase for retrieving the operator associated with that state, and a corresponding API in operator for retrieving the ket associated with that Operator. Brian quickly realized, however, that this is not the best design for the sorts of things we’ll be working with, because many kets (particularly when you get to 2D or 3D systems) will actually be simultaneous eigenkets of multiple operators. So, instead of the previous API, I wrote a new file, operatorset.py, which contains operator_to_state and state_to_operator mapping functions. This is nice, because in making the file we define a global dictionary with all of the appropriate mappings of classes, so that in one place there’s a complete record of these correspondences. In place of the _get_basis_kets function, which was important for trying to form inner products or expectations values, I added an enumerate_states method to represent.py with essentially the same functionality.

This seemed to be working quite well, but I soon stumbled upon an issue related to passing Bra states into the state_to_operator global function. The dictionary for mapping states to operators simply has ket classes and their corresponding operators. So, to map a bra state to an operator, we need to know its dual_class. If someone passes a Bra instance, this is fine, because dual_class is a property. However, if they instead pass a Bra class, we can’t call dual_class because its a property! So, one task in this pull request was to actually change dual_class to a classmethod throughout sympy.physics.quantum, which worked perfectly.

Some other concerns were raised as well. To see the full set of commits, you can look at the PR, but some other ones included:

• Commutativity of Wavefunctions and DifferentialOperators: If we represent an expression in terms of Wavefunctions and DifferentialOperators, then we want that expression to preserve its order so that we can call qapply on it. DifferentialOperators are already set to be non-commutative, but Wavefunctions inherit from Function, which has its own rules for commutativity. So, one simple fix was to override Function’s is_commutative in Wavefunction. Related to that, DifferentialOperator was changed to return a Wavefunction rather than an expression when applied to a Wavefunction. If we have two consecutive derivatives, this is very important for that to actually work, since we have an _apply_operator_Wavefunction defined in DifferentialOperator.
• default_label changed to default_args: The default_label functionality was changed to default_args, which now returns a tuple rather than a simple string. If no args are passed, the QExpr constructor will try to call default_args for that object to create a default instance. Another notable change is that, previously, I had a default_label defined in StateBase which would try, by default, to simply lowercase the default labels of its corresponding operator. It turns out that this is not at all general, so we instead opt for having each Ket and Operator define its own default_args. This will be the most flexible framework for being able to handle all kinds of complicated states. Some examples using this API are in cartesian.py
• default_args added to many base classes: Once the default_args API was well defined, I added default_args for Operator, Ket, and TimeDependentKet, with some tests. In Bra, I put a default_args which simply returns the default_args of its dual_class. This can be overridden for special cases, but  I think it is a good default case. It is probably redundant to have to define default_args in both a Bra and Ket class.

Many changes have been made, and I think the PR is almost ready to be merged. Brian also raised some concerns about being able to handle more general DifferentialOperators, like Laplacians. This will require some more thinking on my part, and will probably be a task for next week.