## Cool results in represent

This post right now is just a quick update with an example python session showing the cool things we can do in represent now, after a very hectic week of code writing. I’m currently in the process of finalizing things in the way that represent works, so I will post a much more detailed post at some point later, but right now I just want to show that we can actually do quantum mechanics in continuous bases! (In this case, with the particle in a box system).

Here is an example python session with notes on the output:

>>> from sympy.physics.quantum import *
>>> from sympy.physics.quantum.piab import *
>>> wf = represent(PIABKet())
>>> wf
Wavefunction(2**(1/2)*(1/L)**(1/2)*sin(pi*n*x_1/L), (x_1, 0, L))
>>> wf.norm
1

We get a particle-in-the-box Wavefunction when it is represented and it’s properly normalized.

>>> represent(XOp()*PIABKet())
Wavefunction(2**(1/2)*x_2*(1/L)**(1/2)*sin(pi*n*x_2/L), x_2)

We now get a wavefunction with an extra factor of x!

>>> represent(PxOp()*PIABKet(), basis=XKet)
Wavefunction(-2**(1/2)*hbar*I*pi*n*(1/L)**(1/2)*cos(pi*n*x_2/L)/L, x_2)

The PxOp actually takes the derivative of the wavefunction correctly! (Momentum operators in the position basis are differential operators).

>>> represent(PIABBra()*XOp()*PIABKet(), basis=XKet)
0

Here, we insert two unities. The first one integrated collapses a delta function, but the second one integrated actually computes the expectation value of x for the particle-in-a-box wavefunction (which is what you expect from representing <psi|X|psi>).

EDIT: As you’ll see in the comments, Raoul pointed out that this is actually incorrect, and I will be looking into it!

There are still a few kinks being worked out, but we’re very close to having a nice finished product which is why I am reserving a longer post for tomorrow. I am quite pleased with the results we’re seeing so far though!

## An altered internal API for default representations and fun with Wavefunctions

This week was a very, very busy week. While there were only a handful of commits, they were probably some of the most crucial ones yet in bringing this home and being able to start writing example systems. With pencils down only two weeks away, it is seriously crunchtime.

A very nice development was the incorporation of Wavefunctions (Wf) and Differential Operators (DO) into the main represent logic. If you now try to represent anything in cartesian.py, you will get expressions that actually include Wf and have the DO applied appropriately! (I will soon be adding support for a qapply flag, which if set to false will leave all DOs intact without applying them.) Support for this included adding calls to qapply as well a function which unwraps all Wf around expressions before integrating over unities, then rewrapping the final result in a single Wf object.

With this done, I moved on to another quite important task in the represent logic that I discussed with Brian. This involves dealing with the ambiguous case in represent where no basis is specified in the initial call to represent. The current way that the quantum modules are set up to handle this logic involves the states having an internal _represent_default_basis method, which simply calls the appropriate _represent_FOO method internally in the class. This entire process is basically a black box to the outer represent logic. This becomes a problem in continuous bases, where we actually need to know which basis was chosen for representing (in order to know if there were any unities inserted and integrate appropriately). This is also a more general problem in that if you have an arbitrary expression, you don’t know if all of the internal _represent_default_basis calls actually represent the individual QExprs in the same basis! You could, in fact, have one QExpr represented in a different basis than another and end up with some very weird things.

So, there are really two problems exposed by this API weakness, and we need to address them both. First, the main represent logic needs to know what the default basis for a given QExpr is. Second, in the case of arbitrary quantum expressions, we need to choose one basis to represent all of the individual QExprs in.

Brian and I have discussed strategies to address both of these issues. The first is a change in the internal API for default representations. Rather than using the current _represent_default_basis, classes should now contain an internal _get_default_basis, which simply returns  the class which is meant to be the default basis for that QExpr. This way, we can still replicate the previous behavior, but the main represent logic now knows which basis is the default for that class. Because this involves changing an already quite large base of classes in the spin and quantum computing modules, this was quite an arduous task. After much testing and grappling with errors, the spin and cartesian classes now follow the new represent conventions for this, with the quantum computing classes soon to follow. (I should note, that as these changes were made, the spin classes were also changed to follow additional new conventions for represent. The first is that internally, the basis option is converted to a basis state rather than left as an operator in represent. This means that all _represent_Op methods had to be changed to _represent_Ket methods. Second, rather than being able to specify a single operator, you must now specify the complete commuting set of operators for a given eigenket. This means, for the spin classes, rather than being able to pass basis=Jx, you now pass basis=set([J2, Jx]).)

The second change is to have represent choose which basis to represent in during the first step of the recursion. With the current algorithm, _represent_default_basis was simply called at the lowest level of the recursion, on the individual QExpr. Now, the basis is chosen based on the default basis of the first QExpr in the expression, and all QExprs in the arbitrary expression are represented in this basis. This means all the representations will be consistent and we aren’t left with any strange final expressions. This change in the represent logic is next on my TODO list.

It is really nice to see all of the pieces of a now very robust representation logic coming together. The schedule I have been following is not quite what I initially planned, but the changes to represent logic for continuous bases ended up being way more subtle than what I initially anticipated. With all the possibilities for different quantum systems out there, I really envision that at the end of the day this revamped represent will be able to handle what is thrown at it. By GSoC pencils down, I will have at least a few different nice example systems implemented and some very clear documentation. Documenting this is very important to me because I’ve put a ton of time into making sure it will work in very general cases (with Brian’s amazing help and discussion), so making sure users know how it works and how to use its full power is the most important thing to me.

## 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 $\frac{d}{dx}$ 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 $\frac{d^2}{dx^2}$
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 $\langle x_1 | x \rangle$, 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.