Debugging DiracDelta integration and next steps in representations

 

This week has been pretty eventful, though not much code itself has been written. In Illinois we had some storms early in the week which knocked out my internet for a few days…but now I’m back and up and running.

I mentioned in my last post that I found a bug in the integration of DiracDelta functions. I had my first branch of GSoC merged into sympy with a fix for this, but then found out that the tests that I had written for it were failing in master for some other users! The problem that I fixed was that the when going through a Mul with multiple DiracDeltas, deltaintegrate would only keep the last one it found and completely throw out the rest. Instead, I had it collapse the first one it found with the integrated variable in it, and keep the rest with the appropriate subsitutions. This comes from the identity that

If has other delta functions inside of it, the identity still applies!

The problem with my fix turned out to be a sorting issue! The test that I had written was tailored for the results on my computer, but the args of the Mul were being sorted differently on others’ computers so that the test ended up failing. I then opened another pull request with code that used sympy’s default_sort_key in deltaintegrate to make sure all the arguments of the Mul were sorted the same way before the function was integrated, and that fixed the issue!

In addition to that, I was able to touch base with my mentor Brian Granger and he gave me some good suggestions on changes to make to my existing code and also some new directions to take in representing continuous bases. I think these will lay good groundwork for starting to actually implement multiple coordinate systems and the canonical quantum systems (particle in a box, harmonic oscillator, etc.). These changes include:

1. Rather than having an API to get the basis kets of an operator, it makes more sense to have an API to get the basis operators of a given ket. This is because in higher dimensions, position kets are actually simultaneous eigenkets of multiple operators. To this end, there will be a global function which can map between states and sets of operators, or vice versa.
2. In addition to this, rather than just taking an operator class as the specification for a basis in represent, represent will now be able to take three things as inputs for the basis: a single operator, a set of operators, or a state. It will then use the global function to map between them.
3. The default labeling systems needs to be revamped with states with multiple labels in mind. It makes more sense to define it as something like “default_args” and have it always return a tuple. More tests with time-dependent states will also help make sure this is clear

I have also been starting to think about how to handle implementing the various coordinate systems that will need to be done. It might make sense to have a general Position type class, which can then me called to have a request coordinate set returned. This needs to be fleshed out a bit more though. First, I will work on making the changes Brian suggested and submit my first big GSoC pull request.

Representing arbitrary quantum expressions in sympy

While the code in sympy.physics.quantum is currently very good at representations in discrete Hilbert spaces, representing expressions in continuous Hilbert spaces requires alot of legwork. The work I’ve done so far does a good job on representing simple expressions, but what about some long arbitrary expression?

Brian gave a good example. Let’s say we want to represent the expression

in the position basis. The easiest way to do this is just to insert unities in the appropriate places. So, we try to evaluate

where we then integrate over and

Let’s now consider a simpler case, where it’s also easier to see directly the expected result. How might we go about representing the expression

where is an eigenket of the position operator . The easiest way to do this is to just bring in a position bra from the left, to compute

Let’s make sure this is consistent with what we get if we instead insert a unity like we would with any arbitrary expression:

If we then integrate over

Since and are just dummy variables, this is indeed consistent.

How do we then implement this in represent.py? Well, first, its useful to change the get_basis_kets in HermitianOperator a bit. Instead of just taking a number of kets, it now takes a number of kets AND a start index. That way, if we recursively keep track of the index in the options of represent, when rep_innerproduct or rep_expectation are called they know which basis kets they should be using for their calculations.

If we are keeping track of an index for kets, let’s think about when we should increment this index. There are a few key cases to consider (reading the expression from right to left):

1. Ket followed by an operator : Here, the ket and the operator should be using the same index for their kets. The representation of is and representation of will start with the same ket (). In this case, we don’t increment our index in between processing the ket and the operator.
2. Bra followed by a ket : In this case, the two representations should use different kets (there’s nowhere to insert a unity), so we should increment our index before processing the ket, thus calculating
3. Operator followed by a bra : This is a case where we are inserting a unity, but we still need to increment the index we’re keeping track of. That’s because the representation of starts with an ket on the right side, so the index when we take this representation is just 1. But, for representing the bra, we take $\langle x|x_2\rangle$, so we increment the index. Thus, we can take the representation as .
4. Operator followed by an operator : We again have to increment the index even though we’re inserting a unity. The start index of the expectation value for the first operator will be one less than the second, so we end up calculating .

Now that we know how to insert unities, we also need to keep track of which kets we inserted so that we can integrate over the appropriate coordinates. In the case of the position basis, there will be many delta functions that come out of that, but most of them collapse out when you actually do the integral. I first tried using integrate() on the resulting expression directly, but found a bug where the delta functions were incorrectly collapsed. An example:

>>> x, x_1, x_2 = symbols(‘x,x_1,x_2’)

>>> f = x_2 * DiracDelta(x-x_2) * DiracDelta(x_2-x_1)

>>> integrate(f, (x_2, -oo, oo))

x_1

We expect, instead, either x*DiracDelta(x-x_1) or x_1*DiracDelta(x-x_1). As a first approach to this, and collapsing dirac deltas in general, I wrote a collapse_deltas function, which, given a list of unities, would collapse the delta functions. I later discovered that a function in the integrals module, deltaintegrate(), already accomplished much of this, so I instead created a new branch and pushed a bug fix for the bug I show above.

Phew. That was alot. The code for the updated represent logic can be seen in my represent branch here. The pull request for my DiracDelta bug fix is here.

Further changes to representations

This past week brought further changes to how representations are done in sympy.physics.quantum, as well as some basic tests and doctests.

First, I added some basic tests for the default label functionality that I talked about last week. Those are here.

Beyond that, I added code to represent which would allow for calculation of operator representations, in the form of <x’ | A | x>. There is still a question, which was discussed over the list, of how to deal with representations of arbitrary expressions. I think the idea is simply to try to apply a bra from the basis you are attempting to represent in on the left of the expression. One small subtlety is that this would mean that the calling represent with the position basis on A | x> would give you the same result as calling represent on just A (that is, <x’ | A | x>). The code for that, and tests, is in this commit.

Next, I set about making some basic Wavefunction and DifferentialOperator classes. There was some debate as to whether Wavefunction is needed at all, and I think there are some justifications for it. First of all, it would mean that all representations in continuous bases could have the same class Rather than returning simple expressions, a function, which could be evaluated at any point would be returned. Additionally, this class can have extra convenience functions for calculating normalization constants or calculating quantum probabilities. It could also internally handle things like indexing of states. Finally, it makes writing DifferentialOperator a bit easier. In fact, it is as simple as implementing an _apply_operator_Wavefunction method, to calculate the derivative. There has still been debate over the implementation of Wavefunction, whether it should inherit from Lambda or Function. I have been trying to sort out implementation issues like this, because Function itself is not callable so I need to define my own __call__ method, which essentially duplicates functionality that already exists in Lambda. Initial code for this is here, with the first tests here. Some changes that I made today after reading discussions on the list, are here.

In the next week I’d really like to nail down the implementation of Wavefunction, as well as an example using the already written code in piab.py and changing represent to use Wavefunction. Hopefully by the end of the week I can have a nice example and some pretty plots!

Default labels for states in sympy.physics.quantum

This week was my first week working on sympy for GSoC, and my first step was to try to change represent.py, the piece of the sympy.physics.quantum code, to handle representations of operators and kets that may be derivatives or arbitrary functions. This is important for my work because many operators represented in the position or momentum bases are actually differential operators (for example, momentum in the position basis).

I first needed to discuss with Brian and Ondrej, my mentors, about what exactly represent.py was doing, in the quantum mechanical sense. We concluded that represent should try, if possible, to form an inner product, something like <x | psi>. For operators, it should return the result of <x | A | x’>.

As a first step to this, though, it was necessary to have some default labels for states. For example, if we want to represent something in the basis of XOp, the x operator of the position basis, then what labels for those basis kets should we have (|x_1>, |x_2>, etc.? or |x>, |x’> …). These kets need some default label because when we call represent, all we specify is the operator whose basis we wish to represent the expression in. So, to this end, I added default_label functions to Operator and QExpr, as well as some functions to states and operators so that eigenkets can identify the operator they are associated with and operators know which kets are their eigenkets.

Finally, after this, I changed represent. It currently just looks for a _represent_BasisClass style function in the expression passed to it. Now, however, if this fails, it also tries to form an inner product between an eigenket of the basis and the provided expression, if the expression is a ket or a bra.

The code for all of this is located here: https://github.com/lazovich/sympy/commits/represent

Still to be implemented is the <x | A | x’> functionality in represent. Also, I will be creating a DifferentialOperator class to be able to represent operators like d/dx which can be applied to wavefunctions. Feedback on the code is greatly appreciated from those sympy people who are interested!