WIP: Consistently use AbstractQuArray and define its interface.

[acroy]

This is a first stab at using AbstractQuArray instead of Union(QuArray,CTranspose) whenever possible. The interface for an AbstractQuArray is so-far basically defined by the functions

• rawcoeffs, coeffs
• coefftype
• rawbases, bases
• copy
• similar_type (see comments for a definition)

One of the open issues is type promotion of different AbstractQuArrays for example in arithmetic operations. Right now + and - always returns a QuArray, which is probably not desirable. Somewhat related is the question of how we want to implement new subtypes of AbstractQuArray (for example density matrices). Any ideas and suggestions are welcome.

[jrevels]

This is good stuff!

What do you think of removing coeffs and bases from the required functions list, and simply have default definitions on AbstractQuArray (we could use the ones currently defined for QuArray)? Then, they could be overloaded only when needed (with CTranspose, for example).

Right now + and - always returns a QuArray, which is probably not desirable.

We should definitely overload these when both arguments are of type CTranspose to avoid taking the eager transpose since A' + B' = (A + B)'.

Somewhat related is the question of how we want to implement new subtypes of AbstractQuArray (for example density matrices)

This a good example to play with. Assuming the implementation will look somewhat like:

type QuDensity{B,T,V} <: AbstractQuMatrix{B,T}
   state:::QuVector{B,T,V}
end

...we would have to have a promote_rule definition like so:

function Base.promote_rule{B1,B2,T1,T2,M,V}(::Type{QuMatrix{B1,T1,M}}, ::Type{QuDensity{B2,T2,V}})
    return QuMatrix{promote_type(B1,B2), promote_type(T1,T2), find_promotion(M,V)}
end

So we'd have to define promotion for bases (which I believe could make sense). The tough part is coding up the solution to what I called find_promotion above. I believe that might not be easy to generally define - the implementation may have to be specific to each promote_rule method we come up with, rather than being a separate method, but you get my drift.

The nice thing about using coeffs is that it works by falling back to promotion between underlying coefficient representations, even when their direct relationships to each other are not necessarily well-typed. For example, in the "lazy" representation of QuDensity above, we can only parameterize the type of a vector. Can we deduce what the matrix version of that vector type will be without actually eagerly performing the outer product of the coefficients therein (without resorting to reflection a la Base.return_types)?

We'll definitely need conversion and promotion rules for the sake of genericness, though...otherwise, simple things like storing collections of Q<:AbstractQuArrays could get quite inefficient.

In this same vein, how do we want to handle conversion/promotion between QuArrays and base Array types?

EDIT::

This gets even trickier if we try to get more generic with the QuDensity definition, e.g. if we try to do something like:

type QuDensity{B,T,V<:AbstractQuVector} <: AbstractQuMatrix{B,T}
   state::V
   QuDensity(state::AbstractQuVector{B,T}) = new(state)
end

[acroy]

What do you think of removing coeffs and bases from the required functions list, and simply have default definitions on AbstractQuArray (we could use the ones currently defined for QuArray)?

I like that. Done in one of the last commits.

We should definitely overload these when both arguments are of type CTranspose to avoid taking the eager transpose

Also done.

I am still trying to figure out the promotion business. As you say, it would be nice to rely on the promotion rules for the coeff containers ...

Another thing I realized only now is that CTranspose wouldn't work for the density matrix, it only works with QuArrays. Maybe we should generalize it along the lines of your last suggestion for QuDensity?

[acroy]

In this same vein, how do we want to handle conversion/promotion between QuArrays and base Array types?

That shouldn't be hard. From our side to base, conversion basically means calling coeffs I would say. And in the other direction we would call QuArray(coeffs). Promotion might be trickier, but I am not sure if we really want/need that?

[jrevels]

Another thing I realized only now is that CTranspose wouldn't work for the density matrix, it only works with QuArrays. Maybe we should generalize it along the lines of your last suggestion for QuDensity?

Definitely, that's a good idea.

Promotion might be trickier, but I am not sure if we really want/need that?

Yeah, this is basically what I was getting at - I didn't know whether or not we should have it, and if we did, how it would work (e.g. would Array*QuArray --> QuArray, or Array?). I guess if we're unsure right now then we could hold off for a bit and come back to it later (though we should have an issue tracking it if that's what we do).

[acroy]

Ok. Let's do the generalized CTranspose in a separate PR and open an issue for Array/QuArray promotion.

Regarding this PR: I think the following behavior would be nice

• any operation on something with a subtype of AbstractQuArray should return the same subtype if possible (by subtype I mean QuArray, CTranspose, QuDensity, etc)
• any operation on two things with the same subtype should also return something of this type if applicable (+,-, etc)
• operations on two things with different subtypes use promote_type to figure out the
  type of the result, which is by default QuArray

Does that make sense?

[acroy]

I think we basically have to rely on the container promotion, otherwise we will never get to do our actual project. But you are right, we should be careful with the promotions. I think your example actually reveals a bug (or an omission of defining the respective promote_rule)?!

Regarding our promotions: Let's say we have a function foo, which operates on containers and foos them. Now we want to foo QA<:AbstractQuArray. Currently we do

foo_me(qa::AbstractQuArray) = QuArray(foo(coeffs(qa)), bases(qa))

which takes advantage of our QuArray constructor to determine the required container and element types. But this always returns a QuArray. I think what we need is something like this

similar_type{B,T,N,A<:AbstractArray}(qa::QuArray{B,T,N}, ::Type{A}) = QuArray{B,eltype(A),N,A}

which needs to be defined for each new subtype of AbstractQuArray. If we require in addition that each subtype CA has a constructor CA{B,T,N,A}(coeffs, bases) we can write in general

function foo_me(qa::AbstractQuArray)
    fc = foo(coeffs(qa))
    QAT = similar_type(qa, typeof(fc))
    return QAT(fc, bases(qa))
end

(for subtypes like CTranspose or QuDensity, which really hold a QuArray instead of a container, the required constructor can just construct the QuArray first and then store that).

For operations acting on two instances of the same subtype we can also make use of this idea

function bar_me{QA1<:AbstractQuArray,QA2<:AbstractQuArray}(qa1::QA1,qa2::QA2)
    @assert bases(qa1) == bases(qa2)
    bc = bar(coeffs(qa1), coeffs(qa2)) # operation on two containers, which returns a container
    QAT = similar_type(qa1, typeof(bc)) # shouldn't matter if qa1 or qa2 is used
    return QAT(bc, bases(qa1))
end

Mhm, actually this should also work if only the basis type is the same, but I don't know if one can enforce that ...

That leaves the case with two different subtypes. For this we need the promotion. But if the stuff above works, we only need to get the primary type and the basis promotion right. The container type will be determined by similar_type. For instance

Base.promote_rule{B1,B2,T1,T2,N,A1,A2}(::Type{CTranspose{B1,T1,N,A1}}, ::Type{QuArray{B2,T2,N,A2}}) = QuArray{promote_type(B1,B2),promote_type(T1,T2),N,promote_type(A1,A2)}

should do it. To make use of it we have to modify similar_type such that it takes only types as arguments, but this shouldn't be a problem.

(I haven't tried your "fudging" approach, but it looks (hacky but) reasonable. The problem is that the promoted type might depend on the operation you do. I would rather like to rely on the inventor of foo and bar to return a useful type. Otherwise it could be hard to add new operations ...)

[jrevels]

similar_type{B,T,N,A<:AbstractArray}(qa::QuArray{B,T,N}, ::Type{A}) = QuArray{B,eltype(A),N,A}

I quite like this. Like similar, but provides a constructor instead of an instance, nice.

For operations acting on two instances of the same subtype we can also make use of this idea

function bar_me{QA1<:AbstractQuArray,QA2<:AbstractQuArray}(qa1::QA1,qa2::QA2)
    @assert bases(qa1) == bases(qa2)
    bc = bar(coeffs(qa1), coeffs(qa2)) # operation on two containers, which returns a container
    QAT = similar_type(qa1, typeof(bc)) # shouldn't matter if qa1 or qa2 is used
    return QAT(bc, bases(qa1))
end

Mhm, actually this should also work if only the basis type is the same, but I don't know if one can enforce that ...

If you wanted to just enforce that the basis type is the same, the signature would be bar_me{B}(qa1::AbstractQuArray{B} ,qa2::AbstractQuArray{B}). However, I think whether or not it would need to be a value equality vs. type equality would depend on bar...

I think this idea works great, though, if I understand you correctly. Slightly modifying it, we could do:

similar_type{Q<:QuArray}{::Type{Q}} = QuArray
similar_type{QC<:CTranspose}{::Type{QC}} = CTranspose
⋮
# a similar_type definition for each Q<:AbstractQuArray
⋮

function bar_me{B}(qa1::AbstractQuArray{B} ,qa2::AbstractQuArray{B})
     @assert bases(qa1) == bases(qa2)
     bc = bar(coeffs(qa1), coeffs(qa2)) # operation on two containers, which returns a container
     QAT = similar_type(promote_type(typeof(qa1), typeof(qa2)))
     return QAT(bc, bases(qa1))
end

...where promote_type is defined as you've said.

I think the above modified version of similar_type should be a little easier to define for generalized wrapper types like QuDensity or CTranspose. Then we can just say that each relevant type should have the appropriate outer constructor to infer the types of the coefficient/basis arguments.

This is a cool approach!

[acroy]

Even better! I will try your suggestion and see how it works ...

[acroy]

Seems to work nicely. Once #31 has landed we will have to revisit all those * methods involving CTranspose. To figure out if it really works as intended we should add a new type (QuDensityMatrix or QuStateVec) and see how it goes?

I also fixed the deprecation warnings on 0.4-dev. For this I had to mess with labelbasis.jl; it would be good if you could check that.

[acroy]

Dunno. The warning said "replace int(...) by round(Int, ...)".

EDIT: But Int(...) should also work. I am fine with both.

[jrevels]

That's because int rounded off floating point numbers. Int attempts to truncate, but throws an InexactError if information will be lost. The latter behavior is actually what we want here, in this case - the numbers fed into for these functions Int should always be whole, and an error should be thrown if somehow they aren't.

[acroy]

Alright, I will change that. What else should be done in this PR? Shall I add a new type to test the promotion etc?

[acroy]

Could/should we have a macro for this construct (like @delegate you mentioned once)? This code pattern is probably quite common and having a macro would be nice for users to turn their favorite matrix function into a QuMatrix function.

[jrevels]

So like @applycoeffs f(qa) --> QuArray(f(coeffs(qa)), bases(qa))? That could work. OTOH, making users write out QuArray(f(coeffs(qa), bases(qa)) might result in more "deliberate" code, as the intent about the return type of f is clearer. I could go either way, honestly.

[acroy]

Ok. Then we leave it like that for now. It's more of a convenience feature anyways, which we can also do later if we want.

[acroy]

I did some cleanup and added a few missing exports. There is currently a ambiguity warning for *, which I don't quite understand. Moreover it seems to work nevertheless. Any ideas?

Apart from that it seems to work. I defined

type QuStateVec{B,T,A} <: QuBase.AbstractQuVector{B,T}
    state::QuVector{B,T,A}

    QuStateVec(state::QuVector{B,T,A}) = new(state)
end

QuStateVec{B,T}(coeffs::AbstractVector{T}, basis::B) = QuStateVec(QuArray(coeffs, basis))
QuStateVec{B,T,A}(qa::QuVector{B,T,A}) = QuStateVec{B,T,A}(qa)

QuBase.coefftype{B,T,A}(::QuStateVec{B,T,A}) = A

QuBase.rawcoeffs(qs::QuStateVec) = rawcoeffs(qs.state)
QuBase.rawbases(qarr::QuStateVec) = rawbases(qs.state)
QuBase.rawbases(qarr::QuStateVec, i) = rawbases(qs.state, i)

QuBase.similar_type{Q<:QuStateVec}(::Type{Q}) = QuStateVec

QuBase.typerepr(::QuStateVec) = "QuStateVec"

Base.copy(qs::QuStateVec) = QuStateVec(copy(qa.state))
Base.promote_rule{B1,B2,T1,T2,A1,A2}(::Type{QuStateVec{B1,T1,A1}}, ::Type{QuVector{B2,T2,A2}}) = QuArray{promote_type(B1,B2),promote_type(T1,T2),1,promote_type(A1,A2)}

and so far it behaved as expected (e.g. returns a QuStateVec if I add two of them or multiply one by a QuMatrix). Note that I had to provide typerepr, which we should either add to the interface or define a fallback.

I would say if everyone is happy (and we figured out the ambiguity warning), we merge this and add the new types in a separate PR ...

[acroy]

@jrevels : Thanks! Looks like your suggestion does the job. Is there anything else that should be done?

[jrevels]

Why the addition of the semicolons here?

[acroy]

This is needed for 0.4. I think it has to do with concatenating vs non-concatenating behavior.

EDIT: For instance

julia> [1:3]
WARNING: [a] concatenation is deprecated; use [a;] instead

[jrevels]

Ah, okay. I think this PR is ready to merge then :)

[acroy]

Cool :-) Then let's do it and see how it works in the real world ... Just one more question: what shall we do about typerepr? Should every type provide one or do we want a default typerepr(::AbstractQuArray)?

[jrevels]

I'd be cool with a default that looked like:

typerepr{Q<:AbstractQuArray}(::Q) = "$Q"

...and then overload it for specific types when appropriate.

[acroy]

Ok, I will add this and then merge ...

[acroy]

In the end I used typerepr(qa::AbstractQuArray) = "$(typeof(qa).name)", which was easier to add and gives nice output (could be that we don't need the other typerepr at all).

[jrevels]

Oh, nice! I didn't know about the name field. In that case, yeah, typerepr can be removed. If you want to remove it now before merging you can, otherwise I'll do so later.

[acroy]

I learned about it only by accident. Seems that the more individual names are nicer (e.g. QuVector is only printed as QuArray). So let's leave it like that then.

[acroy]

Let's move on ...