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

[amitjamadagni]

WIP :

• [x] Spin operators to be implemented
• [x] Tests

Implemented

• [x] +,- on QuArray, CTranspose
• [x] momentm_op, position_op added.

[amitjamadagni]

@jrevels I was thinking a bit off on it. I realize it is costlier. Anyways thanks for this :).

[acroy]

For spin operators we might want to keep an eye on FixedSizeArrays.jl (also JuliaLang/julia#7568)?

[amitjamadagni]

@acroy I was doing something like this based on QuTiP :

function spin_j(j)
          m = reverse([-j:j-1])
          N = length(m)+1
          m_coeffs = [sqrt(j*(j+1.0) - (x+1.0)*x) for x in m]
          return QuArray(spdiagm(m_coeffs,1,N,N))
end

As far as my understanding of FixedSizeArrays.jl goes we are trying to create the array in a more efficient way with focus on the size. I am still trying to further understand but please do let me know if this is going in the right direction.

[acroy]

@amitjamadagni : for now we probably want (have) to use sparse matrices as containers as you do. I just stumbled across FixedSizeArrays and thought those could be useful for stuff like Pauli spin matrices, which are quite common. AFAICT the FixedSizeArrays need Julia 0.4 anyways, so this is something for QuBase 2.0 :-)

[amitjamadagni]

@acroy @jrevels a review would be helpful. I have added the spin operators. Thanks

[acroy]

You have to rebase against the latest master and maybe we should put the spin operators into a separate file?

[amitjamadagni]

@acroy sure okay. But could we move the entire machinery into a single file with the name as operators.jl ??

[acroy]

Let's see what Jarrett says, he has a better overview of the repo structure. I thought we could move the spin operators into a file spinops.jl or such ...

[amitjamadagni]

@acroy I am trying to rebase against the latest master. In the attempt I created a bit of mess with the commits. Any pointers relating to rebase would be helpful. Thanks.

[jrevels]

I agree with using Rationals for the spin operators and moving them to their own file ("spinops.jl" sounds good to me).

@amitjamadagni You might've already seen this before, but I found this tutorial very helpful when learning to rebase. Is there a specific thing you're having trouble with, error message, etc.?

[amitjamadagni]

@jrevels Yeah I have been trying and playing around but the problem I end up is I am unable to see the previous commits when I do a rebase -i master, to squash. I will give it a try once again and let you know.

[amitjamadagni]

@acroy and also it would helpful to hear your comments on the unicode method names. I have been able to generate sigma_x but not sigma_y and sigma_z (I guess there is no support for \sigma_y and \sigma_z).

[acroy]

Regarding the Rationals I actually thought to use them in spin_h1/2 instead of the isinteger test. It somehow feels strange to allow arbitrary floats for the spin. The matrices would of course (?) have floating point elements (convert m to float before constructing the sparse matrix).

I haven't looked too much into the Unicode stuff yet. If the other sigmas are not working, maybe we should look at this at some other point (it's not crucial, but just nice to have).

[amitjamadagni]

@acroy I guess we would still require the isinteger test as a check on the input, just having a rational type would not help as we can still have inputs like 4/3 which are not allowed and also we do not have a isrational method to check. Please do let me know if I am missing something here. Also this is a better check on the floats that we are using, I guess there are no edge cases which I can think of so. Please do let me know your thoughts on it. Thanks.

[acroy]

I thought we could test the numerator for being positive and the denominator for being 1 or 2. Might be that isinteger for Rationals actually does something like that. I don't insist on using Rationals though, it just feels odd to allow arbitrary floats and but then restrict the input (what happens if you put in 0.500001?)

[amitjamadagni]

@acroy I just ran some benchmarks. Here are the results : Using Rational we get : 0.002287712 seconds (17960 bytes allocated) Using Float we get : 0.002204122 seconds (17936 bytes allocated)

Any thoughts on this would be helpful.

[acroy]

I wouldn't expect a big performance difference for that piece of code. If we would use Rationals as elements of the matrices, it would certainly have an impact, since (almost) all linear algebra routines are optimized for floats.

One might get a (negligible?) performance improvement if one passes numerator and denominator as separate integers to those routines. If both values are knows at compile time, Julia might be able to optimize the functions (eliminate branches for instance). But I don't think that it is crucial, unless we expect those routines to be called millions of times (although then garbage collection will kill us first).

My point was only that it is somewhat unintuitive/confusing to be able to pass arbitrary floats as spins. Moreover, there might actually be edge case, like passing 0.500001 or such.

[jrevels]

I can understand passing the numerator or denominator separately.

But if we do go for allowing a single number argument for spin, I think it would be odd for the user to enter something like spin_jx(1/2) and have it throw an error.

I agree that we should test for edge cases like @acroy (and maybe start increasing/decreasing precision until we it hits a failure mode) and then use that knowledge to help make a decision. I'd do it myself but I'm in transit right now...

[amitjamadagni]

@acroy my bad I had forgot to mention that the input check was for Rational but then I had converted the internals to float. The memory used is more in case of rationals when compared to float, can this be a case for consideration ?? Also from my understanding we should be setting the precision so that we consider the edge cases, right ??

[acroy]

We should check what isinteger is doing under the hood. Maybe it already uses round or such and we don't have to worry about the edge cases...

A more sophisticated (but probably overly complicated) approach might be something like

abstract Spin{S} # S is just the numerator of the spin value
immutable HalfSpin{S} <: Spin{S} end
immutable IntegerSpin{S} <: Spin{S} end

plus constructors and some convenience functions. Then we exactly know what kind of spin we have to deal with. Of course the problem would be shifted to implementing the constructor taking a float (or conversion from float), but this could be fixed/changed later without affecting anything else. The question is if we need those types for something else? Anyways it might be a nice exercise.

@amitjamadagni Where those timings obtained by running the function once? If so you basically measure compilation time/memory. Have a look at the performance tips to read an explanation. If you ran it twice already, just ignore me ;- )

[amitjamadagni]

@acroy the instances of isinteger in floatfuncs.jl and rational.jl from base are as follows :

isinteger(x::FloatingPoint) = (trunc(x)==x)&isfinite(x)

isinteger(x::Rational) = x.den == 1

I have searched for trunc but could not find correct references to it. I will work on the other idea as well and send in a commit.

[amitjamadagni]

@acroy I have made an attempt at implementing immutable types HalfSpin and IntegerSpin. Any thoughts on this would be helpful.

[acroy]

Cool! That looks basically like I anticipated it. We will need some constructors/conversions, say from float to HalfSpin, but for this we could just take your previous implementation.

[amitjamadagni]

@acroy I will send in a commit with the implementation 1 and then I guess this would be ready for merge ??

[acroy]

You only need your precious approach for the HalfSpin constructor and the we will see. Don't forget to add tests. Maybe using commutator relations might be a good approach for testing?

[amitjamadagni]

@acroy sorry I could not get that comment :(

[acroy]

Sorry, I meant the previous implementation (isinteger) only applies to the HalfSpin constructor, which takes a float/Rational as input. If you have that and take the other comments into account we will see where we are. Moreover, you will need some tests. One way might be to use the commutation relations (simplest case are Pauli matrices) as tests.

[amitjamadagni]

@acroy I have made an attempt at the following :

• Type independence of Spin
  • Type conversion to spintypes
• Commutator method
• Tests using the commutator method

I have tried to work around on the making the HalfSpin type without den but I am unable to crack the problem. Any leads on this would be helpful.

[acroy]

Regarding the Spin types I thought something like this might work

abstract Spin{S}

immutable HalfSpin{S} <: Spin{S} end
immutable IntegerSpin{S} <: Spin{S} end

# constructors
HalfSpin(j::Integer) = j > 0 ? HalfSpin{j}() : error("Spin must be positive")

function HalfSpin(j::Float64)
    if (j > 0) && !isinteger(j) && isinteger(j*2)
        HalfSpin(int(j*2))
    else
        error("Spin must be positive and a multiple of 1/2.")
    end
end

IntegerSpin(j::Integer) = j > 0 ? IntegerSpin{j}() : error("Spin must be positive")

function IntegerSpin(j::Float64)
    if (j > 0) && isinteger(j)
        IntegerSpin(int(j))
    else
        error("Spin must be positive and integer.")
    end
end

# conversion
Base.convert{S}(::Type{IntegerSpin{S}}, j::Integer) = IntegerSpin(j)
Base.convert{S}(::Type{HalfSpin{S}}, j::Float64) = HalfSpin(j)

# extracting the spin
spin{S}(is::IntegerSpin{S}) = S
spin{S}(his::HalfSpin{S}) = S/2

[acroy]

@jrevels : Any thoughts on the Spin types. Are they overkill? We could also go back to Amit's initial approach and introduce the spins separately...

[amitjamadagni]

@acroy I have tried to play around with the code, here are few doubts that I have :

1. Why is it that we need convert, if we have already defined the methods for integer and float, we have something like

hs = HalfSpin(1/2)
#giving out 
HalfSpin{1}()

which is what is expected out of convert.

1. I have tried this

x = IntegerSpin(2.0)
println(x)
println(Spin(x))

This gives an error like

ERROR: cannot define function Spin; it already has a value

I guess I am missing something here. Please do let me know what I am missing. Thanks.

[acroy]

@amitjamadagni : I think it is just a typo, it should be spin(x).

[amitjamadagni]

@acroy in sense it fails to compile the file, giving out the error. Is it like we are extracting the value from the type ?? And also any thoughts on convert ??

[acroy]

@amitjamadagni : Not sure what you mean, the code I posted works for me as expected:

julia> x = HalfSpin(1/2); println(x); println(spin(x))
HalfSpin{1}()
0.5

julia> x = IntegerSpin(1); println(x); println(spin(x))
IntegerSpin{1}()
1

julia> x = IntegerSpin(2.); println(x); println(spin(x))
IntegerSpin{2}()
2

[acroy]

Maybe convert is not necessary. The problem is that we want to avoid functions, whose return type depends on the value(s) of the argument(s). This would cause type instability. I guess we should then rename the spin functions to something like spin_value and instead introduce

spin{S}(is::Spin{S}) = is 
spin(j::Integer) = IntegerSpin(j)
spin(j::Float64) = HalfSpin(j)

Those would be like your convert_to_spintypes. Then spin(1.0) wouldn't work of course, but OTOH there is no reason why a float should be accepted for this.

Now that I am writing this, I realize that the IntegerSpin type is actually superfluous. We can do everything in terms of HalfSpin:

abstract Spin{S}

immutable HalfSpin{S} <: Spin{S} end

# constructors
HalfSpin(j::Integer) = j > 0 ? HalfSpin{2*j}() : error("Spin must be positive")

function HalfSpin(j::Float64)
    if (j > 0) && isinteger(j)
        HalfSpin(int(j))
    elseif (j > 0) && isinteger(j*2)
        HalfSpin{int(2*j)}()
    else
        error("Spin must be positive and a multiple of 1/2.")
    end
end

# extracting the spin
spin_value{S}(his::HalfSpin{S}) = S/2

spin{S}(is::Spin{S}) = is 
spin(j::Integer) = HalfSpin(j)
spin(j::Float64) = HalfSpin(j)

(Maybe we don't even need the abstract Spin type). Now the type of all functions is independent of the actual value and everything converts the same way.

[amitjamadagni]

@acroy I have added the latest implementation and modified other convenience functions. Yeah the doubts on convert method stand clarified. Also regarding the error I made a very grave mistake. Sorry for that.

[acroy]

Don't worry, nobody is perfect :-) Let's see what Jarrett thinks of this, but I think this looks pretty much Ok (in the end we might want to remove Implementation 1 and squash the commits).

[amitjamadagni]

@jrevels a review would be helpful. Thanks :)

[amitjamadagni]

@acroy @jrevels I have made an attempt. A review would be helpful.

[acroy]

I think it is basically Ok. If you still have energy, you could add some comments to that spin function, in particular reminding us that it is not type stable. It would also be great if you could squash your commits.

[amitjamadagni]

@acroy I am unable to squash some of the previous commits.

[acroy]

Why?

[amitjamadagni]

I dont get the previous 8 commits before the last when i do a git rebase -i master.

[acroy]

It should probably be git rebase -i HEAD~9.

[amitjamadagni]

@acroy Thanks for the help. Done squashing.

[acroy]

Mhm, I can't merge it. Maybe you have to rebase again?

[amitjamadagni]

@acroy seems like I am stuck in a loop. Can I close this and open an another pull request ??