juliaattic / qubase.jl Goto Github PK

As briefly discussed in #5 we should decide on a policy for exports. Since exports in some sense also define the package API, this might be tricky to decide now. AFAICT there are several options:

1. Don't export anything.
2. Only export QuArray and FiniteBasis and the related functions.
3. In addition to 2, also export abstract types.
4. In addition to 3, also export "helper functions/types". For example, functions creating special matrices (like raisematrix etc).

We should have some way to indicate that a QuVector or a QuMatrix is transposed/dual. I am thinking of something like JuliaLang/julia#6837, i.e. introducing a dual/transposed type. This would make it easier to define inner products and such.

references: JuliaLang/julia#6837, JuliaLang/julia#4774

As discussed in #7, we may want to define types/labels to denote quantum pictures of the workspace. Transformation of propagators and states between different pictures may be also helpful.

It seems to me like partial traces would be something required for many uses of JuliaQuantum, yet I couldn't find an implementation in QuBase.

Hacking something together really quick is not hard, but does anybody have a piece of code lying around that computes it the JuliaQuantum way?

The following error occurs when accessor functions are called on ladderops.

julia> r = raiseop(4)
4x4 QuMatrix in QuBase.FiniteBasis{QuBase.Orthonormal,1}:
...coefficients: Base.SparseMatrix.SparseMatrixCSC{Float64,Int64}

    [2, 1]  =  1.0
    [3, 2]  =  1.41421
    [4, 3]  =  1.73205

julia> coeffs(r)
ERROR: MethodError: `coeffs` has no method matching coeffs(::QuBase.QuArray{QuBase.FiniteBasis{QuBase.Orthonormal,1},Float64,2,Base.SparseMatrix.SparseMatrixCSC{Float64,Int64}})

# but the following thing works.

julia> m = [1+1im 2+2im 3+3im; 
            4.9+4im 5+5im 6+6im; 
            7+7im 8+8im 9+9im]
3x3 Array{Complex{Float64},2}:
 1.0+1.0im  2.0+2.0im  3.0+3.0im
 4.9+4.0im  5.0+5.0im  6.0+6.0im
 7.0+7.0im  8.0+8.0im  9.0+9.0im

julia> spm = sparse(m)
3x3 sparse matrix with 9 Complex{Float64} entries:
    [1, 1]  =  1.0+1.0im
    [2, 1]  =  4.9+4.0im
    [3, 1]  =  7.0+7.0im
    [1, 2]  =  2.0+2.0im
    [2, 2]  =  5.0+5.0im
    [3, 2]  =  8.0+8.0im
    [1, 3]  =  3.0+3.0im
    [2, 3]  =  6.0+6.0im
    [3, 3]  =  9.0+9.0im

julia> qspm = QuArray(spm)
3x3 QuMatrix in QuBase.FiniteBasis{QuBase.Orthonormal,1}:
...coefficients: Base.SparseMatrix.SparseMatrixCSC{Complex{Float64},Int64}

    [1, 1]  =  1.0+1.0im
    [2, 1]  =  4.9+4.0im
    [3, 1]  =  7.0+7.0im
    [1, 2]  =  2.0+2.0im
    [2, 2]  =  5.0+5.0im
    [3, 2]  =  8.0+8.0im
    [1, 3]  =  3.0+3.0im
    [2, 3]  =  6.0+6.0im
    [3, 3]  =  9.0+9.0im

julia> coeffs(qspm)
3x3 sparse matrix with 9 Complex{Float64} entries:
    [1, 1]  =  1.0+1.0im
    [2, 1]  =  4.9+4.0im
    [3, 1]  =  7.0+7.0im
    [1, 2]  =  2.0+2.0im
    [2, 2]  =  5.0+5.0im
    [3, 2]  =  8.0+8.0im
    [1, 3]  =  3.0+3.0im
    [2, 3]  =  6.0+6.0im
    [3, 3]  =  9.0+9.0im

Similarly with respect to other functions. I have tried playing around but I have not found any solution for this. It would be great to hear something on this. Thanks.

Will QuArray support storing base's expression? or a new subtype of QuArray?

my problem is,for example:

to store a specific state like

the QuArray type can store it like this at present

a = QuArray([x->1/sqrt(2)*sin(2/a)*sin(1*pi/a*x),x->1/sqrt(2)*sin(2/a)*sin(2*pi/a*x)])

however I hope the coeffs can be isolated from the base' expressions

Hello again! I've been busy with other stuff lately, and haven't been putting time into QuBase.jl.

I believe that QuBase.jl should start targeting the upcoming v0.4 release of Julia. Here are a few reasons:

• It seems that QuBase.jl won't be ready for release until after v0.4 releases anyway (which is planned to be in a couple of weeks, AFAIK).
• Switching over to v0.4 would fix #41
• We'd be able to use Traits.jl, which I think is the right way to standardize interfaces for QuStates and QuOperators so that we can dispatch on behavior without being restricted to specific storage implementations (see my comment in #32).

At some point, we need to figure out how QuArrays promote with other array types (e.g. Julia's base Array).

For example, what should the following operations return?

+(::AbstractQuArray, ::Array)
kron(::AbstractQuArray, ::Array)
*(::AbstractQuArray, ::Array)

As you have noticed, Jarrett (jrevels) is going to present some of our works in the 17th SQuInT workshop on February 19-21, 2015—I will also be there for sure. Therefore, we have to stick to a tight timeline and workflow to efficiently achieve a basic version of the project before the workshop. Here, I tentatively propose a framework on this direction in order to integrate this project into the big picture of JuliaQuantum’s roadmap. @jrevels will take over the leadership on the coordination of this project as he has been working on the antecedent package on the same line and is going to present the work. Let’s see if everyone involved happy with the short-term plan, and I’ll update the plan accordingly. You can open new issues to discuss specific topics. Thanks.

1. Basic typing system of quantum mechanics.
   1. Timeline:
      • A tentative version before Jan 15 so that other people can play with it?
      • A prerelease version before Feb 1st?
      • Technical status frozen on Feb 10?
   2. concepts to be covered (@jrevels: most can be moved from QuDirac.jl code base. Can you start another post to monitor the planned hierarchy of those types or start documenting the plan in Docile.jl? This may help avoid repeated work and conflicts for others.):
      • quantum states (bra-ket vector space notation should be enough).
      • quantum operators (with some important operators implemented in the base).
      • density operator or density matrix (may be absorbed into the quantum operator type).
      • super-operators.
      • propagators (@acroy: can you think about this part?)
      • others.
2. Unit test -- Timeline:
   • Before Feb 1st (Sunday): Think out and implement some simple tests on each part. Pay enough comments in the code file. All work should be mainly done by who defines the types.
   • All tests done by Feb 10? A second person should be involved in this testing process.
3. Systematic test and making examples -- Timeline:
   • Before Feb 1st: making some simple quantum dynamics solvers based on the QuBase.jl for demonstration purpose (@acroy: can lead this part? I believe you have enough foundation to achieve it easily. This part could be evolved into the future quantum solver project for JuliaQuantum). Finish some simple demonstrations of how to use this the basic types to make more complicated operator-state calculations (@jrevels: I believe you have enough material to cover this).
   • Before Feb 7 (Saturday): find some interesting problems in literature with plots to solve with and without the programmed quantum dynamics solvers. Again, @acroy, should be easy for you. I will also think of some interesting examples in quantum information and control field. Feedbacks to the code base if there is anything wrong. If possible, compare with other packages and do some benchmarks. Examples and results will be used for the vivid graphs in the poster and future tutorials, documentations. Yes, poster presentation like fancy plots.
4. Documentations (@jrevels) -- Timeline:
   • Generate a brief description with Docile.jl before Feb 1, so that I can link it to the website and give a feedback. If we have time, we can consider using http://readthedocs.org but optional for now.
   • Work out the framework of the post before Feb 15 with completed examples and demons. Report progress with draft file back to public.
   • Final beautification and modification upon feedbacks before Feb 17. I can print out the poster for free and bring it to Berkeley on Feb 18.
5. Future work and asides:
   • More complete types and functions. @Jutho: can you point out the hierarchy of the type system you are working on with tensor algebra before Feb 1st? We will consider those possibilities while working on the basic version of QuBase.jl as early as possible.
   • You can leave more suggestions on future below.

Edited on Jan 6: we may not necessarily stick to the timeline I proposed at the early stage of this repo. I think @jrevels and @acroy know better than I do on organizing this project. See: #7 .

Our first task in this repo, I believe, is to decide what version of Julia we'll want to try to target with our initial release.

The latest stable build, v0.3, is advantageous because it has been very thoroughly tested and is, well, stable. We'd have a more solid platform to build from and debugging would probably end up being easier as a result.

The other obvious option is to target v0.4, which means we'd be working off of the master branch until v0.4 is released. There are many drastic changes that may be delivered in v0.4, the majority of which I find extremely desirable:

• call overloading
• call to convert fallback
• Potential development of Tim Holy's Traits trick
• The resolution of the vector transpose issue
• Cleanup of A_mul_B etc.

Developing off of master will force us to contribute back to Julia's development when language bugs crop up, or we realize we desire certain features (IMO, a good opportunity for those of us who haven't already contributed heavily to core Julia). We'd also be in a perfect position once v0.4 goes stable, since we'd be developing from it all along. Unfortunately, working to support v0.4 could mean a lot of wasted work if we put effort into supporting features that get removed.

Just to get discussion going, I vote we focus the master branch at the latest stable v0.3 and make a v0.4 branch to try and keep up with the more recent developments. Then we can switch over whenever the time comes. Thoughts?

It would be great if we could use Docile package for documentation. As we are basing our initial attempt on Julia release version 0.3, we could use Docile as a package, this dependency would not be needed in the future as the next version of Julia has documentation in built, which is a merge of Docile. I have looked through this package which is on similar lines in terms of development depending on Julia 0.3 and using Docile HttpServer.jl. It would be great to hear some thoughts on this. (I have created a branch with initial commits which has edits in the REQUIRE as well as added a doc folder with README.md file).

It would be great if QuBase could support creation/annihilation operators correctly when there are multiple sites and the creation operators on different sites must anti-commute.

See discussion here.

There is a problem with our type promotion for QuArrays:

julia> sigmax
2x2 QuMatrix in FiniteBasis{Orthonormal}:
...coefficients: Array{Float64,2}
[0.0 1.0
 1.0 0.0]

julia> sigmay
2x2 QuMatrix in FiniteBasis{Orthonormal}:
...coefficients: Array{Complex{Float64},2}
Complex{Float64}[0.0 + 0.0im -0.0 - 1.0im
                 0.0 + 1.0im 0.0 + 0.0im]

julia> promote(sigmax,sigmay)
(2x2 QuMatrix in FiniteBasis{Orthonormal}:
...coefficients: Array{Float64,2}
[0.0 1.0
 1.0 0.0],2x2 QuMatrix in FiniteBasis{Orthonormal}:
...coefficients: Array{Complex{Float64},2}
Complex{Float64}[0.0 + 0.0im -0.0 - 1.0im
                 0.0 + 1.0im 0.0 + 0.0im])

The issue seems to be that (with Julia 0.3)

julia> promote_type(Array{Float64,2},Array{Complex128,2})
Array{T,N} (constructor with 9 methods)

So maybe this is actually an issue of Julia Base (0.3)? Any ideas?

Although there are already some ideas in #3, we should have more fine-grained discussion on the content of QuBase. Maybe we can use this issue to collect some ideas and refine the roadmap in #3.

The following four points are probably uncontroversial and partly contained in #5.

1. Abstract types for bases and quantum arrays + interface.
2. QuArray and FiniteBasis related:
   • array-like functions/operations: +,-,* etc (it would be nice to have mutating functions as well, like A_mul_B!)
   • some matrix functions (expm)
   • basis-related functions/operations
3. Operator related:
   • convenience constructors: creation/annihilation (raiseop/lowerop), momentum and position,
     displacement and squeezing operators
   • more operators: angular momentum operators, Pauli spin matrices, etc
   • special operations: partial trace and partial transpose
4. State related:
   • convenience constructors: state vectors (statevec) and density matrices, coherent and squeezed states
   • spectral decomposition of density matrices
   • state properties: purity, fidelity
   • expectation value and variance of a given operator and state

(Some of the functions in 3 & 4 could probably be taken from QuDOS.)
Additionally we could include the following aspects, but might put them into separate packages eventually:

1. Quantum information related:
   • special states: W, GHZ, Bell states etc
   • entanglement measures: negativity, concurrence etc
   • more functions for quantum correlations and state distinguishability
2. State representation/visualization:
   • Husimi function, Wigner distribution, Glauber P function
   • Gaussian states, Yurke states, stretch states
   • More visualization/plot functions to explore the geometry of quantum states and operators (Such as states on the Bloch sphere).

Basically anything else should IMO first be developed in a separate packages. In particular with propagators we need to experiment a bit more to find a reasonable (and generic) interface. (Note: If we do it right, we can use QuArrays with ODE solvers from ODE.jl and exmpv from Expokit.jl without much trouble. So the propagation itself is "easy", but finding a good way to switch between methods and pass parameters is tricky).

[EDIT: I added some points from @Qi's comment below.]

Will it be better to add some short cuts for quantum computing in QuBase.jl or just make it a single package? Cause I currently use them as a single dependency in QuComputStates.jl. But most of the types uses super-types in QuBase.jl

If it is, I would happy to create a PR to merge them into QuBase.jl. :-)

Currently the functionality is not supported and it would be great if we could have a work around as it would improve the efficiency in few cases. In sense for a many body system the sparse support would increase the efficiency and performance when used in conjunction to the solvers in QuDynamics. There have been some discussions in the past but any pointers on this would be great. Consider the following :

julia> v = statevec(1, FiniteBasis(2))
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0]

julia> @time tensor(v,v)
  0.112194 seconds (104.74 k allocations: 4.530 MB)
4-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0,0.0,0.0]

julia> v_s = SparseVector([1;0])
Sparse vector of length 2 with 1 Int64 nonzero entries:
  [1]  =  1


julia> function Base.kron(v1::SparseVector, v2::SparseVector)
           return kron(SparseMatrixCSC(v1), SparseMatrixCSC(v2))
       end

julia> @time kron(v_s, v_s)
  0.082897 seconds (42.55 k allocations: 1.709 MB)
4×1 sparse matrix with 1 Int64 nonzero entries:
	[1, 1]  =  1

julia> @time tensor(v,v)
  0.000153 seconds (26 allocations: 960 bytes)
4-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0,0.0,0.0]

julia> @time kron(v_s, v_s)
  0.000010 seconds (20 allocations: 1.203 KB)
4×1 sparse matrix with 1 Int64 nonzero entries:
	[1, 1]  =  1

Though kron for SparseVectors is still not in base, the above has been constructed using the following issues :
JuliaLang/julia#19426
https://github.com/JuliaInv/jInv.jl/pull/26/files

It would be great if we could make QuArray entirely return a sparse construction instead of a dense version, or create an option if the coeffs construct should be sparse or dense with the default being set to a sparse construct. Something like :

julia> statevec(1, FiniteBasis(2))
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: SparseVector
  [1]  =  1

julia> statevec(1, FiniteBasis(2), dense=true)
2-element QuVector in QuBase.FiniteBasis{QuBase.Orthonormal}:
...coefficients: Array{Float64,1}
[1.0,0.0]

One other option would be to just add sparse functionality on top of the current functionality, which I am not sure would be a good option. Any pointers on this would be useful.

Hi!

Is the package still under development?

Currently, S<:AbstractStructure is pretty vaguely defined (my own fault), and I think in light of upcoming development goals (e.g. GSoC proposal) we should give it a review. To close this issue, I would like to either:

• clearly define its purpose, and give it a better name if we can think of one, or
• get rid of it, and use a type hierarchal structure for bases instead of parameterization to define the properties we want.

To begin, I'll describe what I was initially thinking when I proposed the S<:AbstractStructure type family.

My idea was that we could use type parameterization + dispatch to "easily" define the behavior of the wave functions that form the basis for a given QuArray, without having to implement new concrete subtypes of AbstractBasis/QuArray. What do I mean by this? I think an example might work best. Using the hydrogen atom, a potential use case would be defining types like this:

abstract Orthonormal <: AbstractStructure
abstract HydrogenSphPos <: Orthonormal

wavefunc(::Type{HydrogenSphPos}, i) = # ψ_nlm(r, θ, ϕ) 

Here, wavefunc is overloaded to return the basis wave function for hydrogen atom in spherical coordinates corresponding to the ith position in a QuVector. In other words, each index i maps to a unique quantum number tuple (n,l,m) which specifies the basis wave function for the ith coefficient in the QuVector.

What does this buy us?

First, since HydrogenSphPos <: Orthonormal, we know we can skip actually processing the basis wave functions when evaluating products of objects like QuArray{FiniteBasis{HydrogenSphPos}}.

Second, we could provide a method like the below, which would allow evaluation of the superposed state given by the QuArray at the given args:

waveeval{S}(qv::QuVector{FiniteBasis{S}}, args...) = sum([qv[i] * wavefunc(S, i)(args...) for i=1:length(qv)]  

Obviously, more optimization could be done so as to not use a temporary array, but you get my drift. When we move forward to v.04, we could overload call and have beautiful syntax for evaluating states:

Base.call(qv::QuVector, args...) = waveeval(qv, args...)

#now, we can say the following:

Ψ = # some QuVector{FiniteBasis{HydrogenSphPos}}
r, θ, ϕ = # positional parameters
Ψ(r, θ, ϕ) # evaluate entire state at positional parameters

Since the above example doesn't cover it, let's say a user defines a structural type that is not S<:Orthonormal. Then, we could a provide a function like the following that could be overloaded for arbitrary type S<:AbstractStructure:

inner_rule{S}(::Type{S}, i, j) = # provides an evaluation of the inner product of wavefunc(::Type{S}, i) and wavefunc(::Type{S}, j)

...which could then be used in the evaluation of products of non-orthonormal objects.

In conclusion, the above proposal exposes a couple of functions to the user (wavefunc and inner_rule) that allows them to define/extend the behavior of QuArrays without having to dig into the internals.

Key points for this issue:

1. The main problem I see with the above approach is that the amount of closures/functionals involved would probably result in poor performance for certain operations, specifically what I called waveeval above.
2. Does anybody here know how other libraries handle this kind of functionality? I'm not familiar with the methods used for simulation of this sort, and always kind of assumed that the basis functions were either hard-coded or were able to be abstracted away for the systems people are interested in simulating. If anybody can point me to documentation I could learn from, I'd be grateful.
3. If we do continue to use S<:AbstractStructure, what would be a better name for it? "AbstractStructure" is terribly vague, but I haven't been able to think of any better. "BasisFunction", maybe? I'd be worried it might get confused for a subtype of Function, then...

Currently, the CTranpose type only wraps QuArrays. To close this issue, we should generalize the CTranspose type to work for any subtype Q<:AbstractQuArray.

See the discussion in this PR.