QuCoeffs: A more controlled approach to QuArray coefficient containers

[jrevels]

This PR presents a new wrapper type, QuCoeffs, which gives us better control over the behavior of QuArrays' coefficient containers.

Specifically, this is an attempt at providing a solution to #9.

Major changes:

1. The QuCoeffs type. This type provides a wrapper for lazy conjugation/transposition on arbitrary coefficient containers; I tried to incorporate some of the ideas from WIP:Add Transpose immutable JuliaLang/julia#6837. I specifically focused on getting the implementation to make sense for vectors and matrices (nothing N>2). Having QuCoeffs also helps separate some of the logical concerns that arise for QuArrays. Now, we can think of our type structure the following way:
   • Behavior of coefficients is defined with QuCoeffs
   • Behavior of bases is defined with AbstractBasis
   • Their combined/interacting behavior is defined with QuArray
2. QuArray coefficient containers are type-restricted to QuCoeffs (which themselves can hold whatever)
3. QuArray is no longer a subtype of AbstractArray. This allows for a much, much cleaner implementation of both the QuArray and QuCoeffs types. I think I'm seeing in action the argument that @Jutho's been making for a while now - at this point, QuArray behavior really does seem to deviate from the behavior defined on AbstractArrays. We'll have to implement our own printing methods, though :(
4. Minor cleanup/reorganization related to the introduction of QuCoeffs. Convience constructors (statevec, zero, eye) have moved to their own file.

I'm really interested in getting @Jutho's eyes on this - my secret hope is that we could use TensorToolbox.jl to extend the functionality of (or maybe replace) QuCoeffs in the future.

[acroy]

Thanks for doing this! It certainly looks promising. I will also try to have a closer look (time is a bit constrained these days).

[acroy]

Do we need new types here? Aren't bools sufficient?

Allowing for conjugation and transposition separately is certainly most general, but do we have a use case for this? I thought that indicating "duality" (whatever that means in the specific case) would be enough (and maybe easier to grok for implementing inner products and such).

[jrevels]

The bool wrapper types are there for type stability purposes. Otherwise, the only type stable way to construct an instance of QuCoeffs would be to manually declare them as parameters (which is obviously annoying since you'd have to manually declare the other parameters as well) or have the constructor always set those parameters to false (which is limiting/less flexible).

As for only indicating duality - that's how I started, but then I ran into some of the problems discussed in JuliaLang/julia#6837. Mainly, figuring out which operations we want to result in copies and which we want to result in views. If we just had a duality indicator instead of separate conjugate and transposition indicators, then the behavior gets wonky/unintuitive - a transpose would cause a copy, a conjugation would cause a copy, but the ctranspose would cause a view. Thus, doing conj(transpose(arr)) would be a double copy operation, while ctranspose(arr) would be a no-copy operation, which seems ugly/confusing to me.

I'm of the opinion that it makes the most sense for all three of the above operations to "agree" to return a view, or a copy. I obviously went with the view here, and tried to reduce the inefficiencies associated with computing conjugation on the fly by having it be a no-op for non-Complex types.

[acroy]

Thanks for the clarification. I should probably have a look into JuliaLang/julia#6837 again to be on the same page.

a transpose would cause a copy, a conjugation would cause a copy, but the ctranspose would cause a view. Thus, doing conj(transpose(arr)) would be a double copy operation, while ctranspose(arr) would be a no-copy operation, which seems ugly/confusing to me.

Maybe this would indeed be too confusing. I was just thinking that we mainly want v'=dual(v) which typically just happens to be ctranspose(v)=conj(transpose(v)) (i.e. all three should return a copy, but dual would return a view).

[jrevels]

I was just thinking that we mainly want v'=dual(v) which typically just happens to be ctranspose(v)=conj(transpose(v)) (i.e. all three should return a copy, but dual would return a view).

Done. The code is definitely way cleaner and we avoid the view/copy mess, nice!

Disambiguating the dual and ctranspose functions will probably be a point of confusion for users, but hopefully nothing that thorough, explicit docs won't be able to clarify.

[Jutho]

This certainly ooks interesting. Tensors from TensorToolbox.jl are still more general in that every single tensor index can be associated to a vector space or its dual (contravariant and covariant), or even to its conjugated or conjugated dual in the case of non-Hilbert spaces. This is not necessary to represent kets, bras and operators, but is necessary for e.g. representing states as tensor networks. Unfortunately I am having much less time than I thought and hoped to spend on this, so I have not really worked on this for a while now and it doesn't seem I'll have any time throughout January either. Maybe I should just try to get more people involved, but this also requires writing some basic documentation.

[i2000s]

@Jutho I have friends working closely in tensor networks. I am trying to get more people on board on our future directions. I would appreciate if you can provide a minimum documentation and some examples on your package before the SQuInT meeting (Feb 19-21). Some of our volunteers (most are not Julia experts) will have meetings there to think about if anyone would like to do something over the coming months. Hopefully it helps.

[jrevels]

One of the advantages of not subtyping QuCoeffs to AbstractArray is that we can more easily enforce our own indexing scheme (and also are freed from having to make a bunch of functions to resolve ambiguities). I didn't implement indexing for N>2. @Jutho: You said TensorToolbox covers a more general dual indexing scheme for arbitrary N, yeah? Do you think it would make sense to focus our efforts with QuCoeffs on N<=2 case until TensorToolbox reaches a suitable state? I wouldn't want to repeat work :)

I would also be down to contribute to TensorToolbox if there was a little bit of documentation to help new contributors get going; no pressure, though. I, myself, won't have time to focus on anything besides QuBase/QuDirac until March (after SQuINT) anyway.

[jrevels]

I did have another thought:

Despite the obvious performance tradeoffs, one of the main reasons lazy transposition is so heavily discussed in Base right now is that it removes the conflation between row vectors and matrices, by allowing 1-D array transposition to return a 1-D array. We actually already have a way of disambiguating row vectors and matrices for QuArrays, though - the basis information (we set N via the basis field).

With this in mind, we could actually have eager duality for QuCoeffs, keeping track of conjugation/transposition in the type parameters (like the original incarnation of this PR). This would obviate the need for a special dual function, and more importantly, allow us to easily extend functions from Base without getting caught up in the indexing issues that come with lazy transposition.

Once v0.4 releases, we could move to whatever lazy behaviors get decided on by that time.

The counter-argument to the above is that we really desire the performance tradeoffs that come with lazy transposition, which I would understand, but would require quite a bit of work that might be better spent on the transposition issues in Base.

Thoughts? I'm going to draft up a competing PR which implements what I'm talking about above, and we can decide between the two approaches.

[i2000s]

@jrevels I think we should let performance and intuition talk. Better to have something to compare in comprehensive tests before making the final decision, which should become an important criterion for judging. Maybe you can branch out those different ideas on the repo. All sounds interesting!

BTW, regarding naming those types, I think it would be good to stick to the nomenclature that has been widely adopted in quantum mechanics textbooks. Certainly, we possibly will have to invent some new inner types for translating data inside Julia whenever it is necessary. I start wondering: are those names (QuCoeff and QuArray, etc) necessary and intuitive to introduce to our end users? Our users may not care about how arrays, vectors or matrices are employed and manipulated in side of the package. state (bra and ket), operator, superoperator, tensor_product_state and such should be enough for intuitive naming. It will also be easy and quick for future users to fall in love with the Julia library for quantum siblings. For a general policy, maybe we really want to make the typing system as simple as possible, and what we need to put more work on is to solve the compatibility, data manipulating, visibility and performance problems using dispatch and other wonderful features built in Julia rather than to keep inventing complicated types into the library. But we possibly can work on to make some alias and promotion rules in the end.

[jrevels]

@i2000s That's why the QuCoeffs type is not exported. It is not intended to be manipulated by an end user.

For a general policy, maybe we really want to make the typing system as simple as possible, and what we need to put more work on is to solve the compatibility, data manipulating, visibility and performance problems using dispatch and other wonderful features built in Julia rather than to keep inventing complicated types into the library.

If you have concrete suggestions for simplifying our type hierarchy, I'm all ears.

Edit: And the most helpful form for these suggestions would be in a PR, if you have the time.

[i2000s]

@jrevels Hah, you did not export those things. Obviously, I have lagged off for a while. Sorry for my ignorance.