Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

WIP : Spin Operators, Tests. WD : +,-,ops added. #22

Closed
wants to merge 9 commits into from
Closed

WIP : Spin Operators, Tests. WD : +,-,ops added. #22

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0fe8886
WIP : Operators, Docs, Tests. WD : +,-,ops added.
amitjamadagni Mar 29, 2015
98fec4b
Use of Type parameters.
amitjamadagni Mar 30, 2015
f3b5be2
Use of scale!
amitjamadagni Mar 30, 2015
d9574c2
Use of scale, AbstractQuArray usage.
amitjamadagni Mar 30, 2015
722fbc5
WIP : Spin operators.
amitjamadagni Apr 1, 2015
0dc1e40
Spin Operators.
amitjamadagni Apr 2, 2015
d9007d7
Check on input.
amitjamadagni Apr 2, 2015
62145d3
Merge branch 'master' into op_machinery
amitjamadagni Apr 3, 2015
879e9f6
Added : Operators
amitjamadagni Apr 10, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions src/arrays/arraymath.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,9 @@ end

*(dm1::DualMatrix, dm2::DualMatrix) = (dm1.qarr*dm2.qarr)'

+(qarr1::Union(QuArray,CTranspose), qarr2::Union(QuArray,CTranspose)) = bases(qarr1) == bases(qarr2) ? (return QuArray(coeffs(qarr1)+coeffs(qarr2), bases(qarr1))) : "Not compatible"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you want to use error("Bases not compatible.") here. Should we use type parameters (for the basis and N) to restrict the functions in the first place?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@acroy something like

function +{B<:AbstractBasis,N}(qarr1::Union(QuArray{B,N}, CTranspose{B,N}), qarr2::Union(QuArray{B,N},CTranspose{B,N}))

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that way we get e MethodError "for free" if someone tries to add say a QuVector and a QuMatrix. Nevertheless we have to check that the bases are the same (I am not sure if == works though).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think == is the right thing to do here, since it can be overloaded to work properly for each basis type - "value" equality is what we're looking for here, yes? Also, all we need to do to support mixed-basis type operations is to simply define some promotion rules between the various basis types, of which there aren't that many.

-(qarr1::Union(QuArray,CTranspose), qarr2::Union(QuArray,CTranspose)) = bases(qarr1) == bases(qarr2) ? (return QuArray(coeffs(qarr1)-coeffs(qarr2), bases(qarr1))) : "Not compatible"

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know that we already have this pattern in the code, but maybe we should be annotating arguments as AbstractQuArray{B,N} instead of Union(QuArray{B,N},CTranspose{B,N}) for functions like these. It's both less wordy and more generic.

Thoughts, @acroy?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess we could, as QuArray and CTranspose are subtypes of AbstractQuArray. But could we compare bases if we used such a construct ??

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. Your new method signature for -, for example, would look like this:

-{B<:OrthonormalBasis,N}(qarr1::AbstractQuArray{B,N},qarr2::AbstractQuArray{B,N})

...and the bases function is defined for both QuArray and CTranspose.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jrevels : That would be nicer. I guess we haven't really specified the interface for AbstractQuArrays yet, but this would mean all types <:AbstractQuArray should support coeffs and rawcoeffs etc, right? Then we can generalize all the functions as you suggest and it will be easy to add new types with some specialized behavior (density matrices!).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep - I think bases/rawbases should be in there as well.

Edit: Though defining the interface for AbstractQuArrays will probably need its own PR.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. That should be a separate PR and maybe we should also document the interface somewhere.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually we don't want to always return a QuArray here if we allow the arguments to be any subtype of AbstractQuArray. We probably need some promotion rules to get this right. But as we said above all of this should be done in a separate PR. For now it is Ok to leave it like that. @jrevels ?

# scaling
Base.scale!(num::Number, qarr::QuArray) = (scale!(num, rawcoeffs(qarr)); return qarr)
Base.scale!(num::Number, ct::CTranspose) = CTranspose(scale!(num', ct.qarr))
Expand Down Expand Up @@ -90,18 +93,18 @@ normalize(qarr::Union(QuArray,CTranspose)) = normalize!(copy(qarr))

# General tensor product definitions for orthonormal bases
function tensor{B<:OrthonormalBasis,T1,T2,N}(qarr1::AbstractQuArray{B,T1,N}, qarr2::AbstractQuArray{B,T2,N})
return QuArray(kron(coeffs(qarr1), coeffs(qarr2)),
return QuArray(kron(coeffs(qarr1), coeffs(qarr2)),
map(tensor, bases(qarr1), bases(qarr2)))
end

function tensor{B<:OrthonormalBasis,T1,T2,N}(qarr1::CTranspose{B,T1,N}, qarr2::CTranspose{B,T2,N})
return QuArray(kron(rawcoeffs(qarr1), rawcoeffs(qarr2)),
return QuArray(kron(rawcoeffs(qarr1), rawcoeffs(qarr2)),
map(tensor, rawbases(qarr1), rawbases(qarr2)))'
end

function tensor{B<:OrthonormalBasis}(ket::QuVector{B}, bra::DualVector{B})
return QuArray(kron(coeffs(ket), coeffs(bra)),
bases(ket,1),
return QuArray(kron(coeffs(ket), coeffs(bra)),
bases(ket,1),
bases(bra,1))
end

Expand Down
24 changes: 18 additions & 6 deletions src/arrays/ladderops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
###########################
# Ladder Operator Methods #
###########################
# n specifies which particle (a.k.a tensor product
###########################
# n specifies which particle (a.k.a tensor product
# factor) the operator acts on
raiseop(b::AbstractBasis, n=1) = QuArray(raisematrix(size(b), n), b, b)
raiseop(lens, n=1) = raiseop(FiniteBasis(lens), n)

lowerop(b::AbstractBasis, n=1) = QuArray(lowermatrix(size(b), n), b, b)
lowerop(lens, n=1) = lowerop(FiniteBasis(lens), n)

Expand All @@ -22,10 +22,22 @@
before_eye(lens, pivot) = speye(prod(lens[[1:pivot-1]]))
after_eye(lens, pivot) = speye(prod(lens[[pivot+1:length(lens)]]))
function eye_sandwich(lens, pivot, op)
return kron(before_eye(lens, pivot),
op,
return kron(before_eye(lens, pivot),
op,
after_eye(lens, pivot))
end

function positionop(n::Int)
cop = raiseop(n)
return (cop+cop')/sqrt(2.)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should use scale! here since + already makes a copy and / would make another one (same for momentumop below).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@acroy Ah! Once again, I comment only to find that you've beat me to it by the time I'm done writing. I'm starting to believe that you are, in fact, a ninja...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, I just realized we won't be able to perform an in-place scaling of the Float64 array by a complex number; see my edit on my comment below.

end

function momentumop(n::Int)
cop = raiseop(n)
return im*(cop-cop')/sqrt(2.)
end

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A couple of points here:

  1. We should support define these operator functions so that they can target specific factors of product states. See how raiseop and lowerop for an example of this; you could just reuse eye_sandwhich to make it work. Eventually, though, we should replace the usage of that method with specialized implementations that don't make temporary arrays.
  2. I believe it would be more efficient to scale the resulting QuArrays in place than make a copy, e.g.
function momentumop(n::Int)
    cop = raiseop(n)
    return scale!(im/sqrt(2.), cop-cop')
end

Edit: Actually, we can't scale the momentum operator in-place with an imaginary number, as it'll throw an InexactError since the Float64 sparse matrix won't be able to hold complex numbers. Thus, I would do:

function positionop(n::Int)
    cop = raiseop(n)
    return scale!(1/sqrt(2.), cop+cop')
end

function momentumop(n::Int)
    cop = raiseop(n)
    return scale(im/sqrt(2.), cop-cop')
end

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jrevels if I my understanding of the first comment is right it meant that we should be able to track intermediate steps in producing the final answer. In specific have functions which produce the intermediate results and reuse them as necessary ?? But in these case it is a straight forward implementation. Please do correct me if I am missing something.

export raiseop,
lowerop
lowerop,
positionop,
momentumop