JuliaAttic / QuBase.jl

A foundational library for quantum mechanics in Julia
Other
43 stars 6 forks source link

WIP: Consistently use AbstractQuArray and define its interface. #29

Closed acroy closed 9 years ago

acroy commented 9 years ago

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

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 commented 9 years ago

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 commented 9 years ago

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 commented 9 years ago

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 commented 9 years ago

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 commented 9 years ago

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

Does that make sense?

acroy commented 9 years ago

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 commented 9 years ago
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 commented 9 years ago

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

acroy commented 9 years ago

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 commented 9 years ago

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

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

jrevels commented 9 years ago

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 commented 9 years ago

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 commented 9 years ago

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 commented 9 years ago

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

acroy commented 9 years ago

Let's move on ...