Closed amitjamadagni closed 9 years ago
I think for sparse matrices one mostly wants to avoid calculating expm
directly, since the latter is typically dense and large. This is the reason why empv
et al are so important. For those one rather computes expm*v
where v
is some vector. That said, I think that expm(full(S))
is fine for small problems. If one is really serious about it, one needs something like expmv
for QuArrays
and functions that displace/squeeze/etc a vector. If we want to include such functions in QuBase, we can pull in Expokit.jl or ExpmV.jl for instance.
@acroy At some point, we should probably include a definition like the following to dispatch on for sparse QuArrays:
typealias SparseQuArray{B,T,N,S<:AbstractSparseArray} QuArray{B,T,N,S}
If we want to include such functions in QuBase, we can pull in Expokit.jl or ExpmV.jl for instance.
+1, and I think it makes sense to include those functions (in the form of using
the packages, of course).
I would suggest to wait with using Expokit.jl et al. We have to look into it for the propagators anyways.
+1 to SparseQuArray
.
One build passes and the other build fails. The error says "does not pass tests in QuBase" for the second build. Any insights would be helpful.
I guess the Travis failure is not related to this PR. Might be a consequence of JuliaLang/julia#10380. @jrevels : Any ideas?
@acroy done with the edits.
Looks Ok, but I want to wait for @jrevels to comment on the Travis failure. Could you add a test for displaceop
? If one applies the displacement op to the ground state, one gets a coherent state and you can test if the elements in the state vector are correct.
@acroy I get an
ERROR: InexactError()
in setindex! at sparse/sparsematrix.jl:1253
in generic_scale! at linalg/generic.jl:18
when I use the combination of complex
type for alpha
parameter in displaceop
. So I have replaced it with scale. I will work on this and see if this can be made to work if possible with scale.
Edit : And also the build is failing, when I run the tests on my machine it passes
sv1 = [0.5403023058681398,0.8414709848078965]
sv2 = [0.5403023058681398 + 0.0im,0.0 + 0.8414709848078965im]
println(displaceop(2,1)*statevec(1,FiniteBasis(2)))
println(QuArray(sv1))
println(@assert displaceop(2,1)*statevec(1,FiniteBasis(2)) == QuArray(sv1))
println(@assert displaceop(2,1*im)*statevec(1,FiniteBasis(2)) == QuArray(sv2))
println(displaceop(2,1) * statevec(1,FiniteBasis(2)) == QuArray(sv1))
# output
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.5403023058681398,0.8414709848078965]
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.5403023058681398,0.8414709848078965]
nothing
nothing
true
which is as expected. But there is an error in the first build saying the assertion fails.
I guess you have to use @test_approx_eq
or similar. Why do you "hard code" the expected result instead of calculating exp(-1/2)
or whatever for the coefficients?
@acroy we could have done something like
exp(-0.5)*(statevec(1,FiniteBasis(2))+statevec(2,FiniteBasis(2)))
2-element QuVector in FiniteBasis{Orthonormal,1}:
...coefficients: Array{Float64,1}
[0.6065306597126334,0.6065306597126334]
which is not equal to the required output. (Truncated version)
We could define coherent state as
coherent(n::Int,alpha::Number) = displaceop(n,alpha)*statevec(1,FiniteBasis(n))
# then test as
coherent(2,1) == displaceop(2,1)*statevec(1,FiniteBasis(2))
We could also have the expansion of coherent states i.e., |alpha> = exp(-mod(alpha)^2)/2)*(sigma[(alpha)^n/sqrt(n!)|n>)
A function coherentstatevec
should probably either be implemented using the displace
function (which will use expmv
) or directly using the expansion in Fock states. But your version is probably Ok for now.
Anyways, the test probably fails because your basis is too small. Using QuDOS.coherentstatevec
I get
julia> cs = QuDOS.coherentstatevec(10,1.)
QuStateVec{Array{Complex{Float64},1}}(Complex{Float64}[0.606531+0.0im,0.606531+0.0im,0.428882+0.0im,0.247615+0.0im,0.123808+0.0im,0.0553686+0.0im,0.022603+0.0im,0.00854887+0.0im,0.00299672+0.0im,0.00110007+0.0im],10,[10])
julia> cs[1]
0.6065306597114603 + 0.0im
julia> exp(-0.5)
0.6065306597126334
julia> Base.Test.@test_approx_eq_eps real(cs[1]) exp(-0.5) 1e-6
You might be able to use a similar test if you increase the basis size.
@acroy I have made an attempt at implementing coherentstatevec
using expm as well the expanding the summation. I have added the tests for displaceop
using the same. A review would be helpful.
I solved the issue with the Travis failure (see #26). Please rebase and squash after addressing the final comments.
@acroy I have tried these things out. Locally they pass for me, but the build-bot shows an error
println(coeffs(displaceop(2,1.0)'))
println(coeffs(squeezingop(lowerop(2), QuArray(eye(2)),2.0)))
println(@assert coeffs(squeezingop(lowerop(2), QuArray(eye(2)), 2.0)') == coeffs(displaceop(2,1)))
println(@assert squeezingop(lowerop(2),QuArray(eye(2)), 2.0)' == displaceop(2,1))
[0.5403023058681398 0.8414709848078965
-0.8414709848078965 0.5403023058681397]
[0.5403023058681398 0.8414709848078965
-0.8414709848078965 0.5403023058681397]
nothing
nothing
Does it make a difference if you take coeffs(displaceop(2,1.0))
in the last assert (note the 1.0)? Otherwise you have to resort to @test_approx_eq
or @test_approx_eq_eps
.
But locally for me it is passing, I guess it should pass by the build-bot as well, right ? I am missing something here ??
It's strange but it can happen. In the worst case you have to put in println statements like in your local test.
The ==
in your @test_approx_eq
commit is superfluous ...
Hah, nice! Could you please squash your commits. I will have another look tomorrow, but I think everything looks good now.
An attempt at implementing expm, full (sparse to dense) as expm is not implemented for sparse matrix. With respect to the implementation of expm with respect to sparse matrix, I have been referring to the https://github.com/JuliaLang/julia/issues/10035 which had hints at scipy implementation. The scipy implementation is based on the following algorithm (http://core.ac.uk/download/pdf/17664.pdf) as cited in the documentation. I will try a hand at implementing this if already not implemented, but that said as exponentiation is only available for dense matrix, I have implemented full function and added the expm on the top of that.
Also the
*
method using DualMatrix was giving some inconsistent results, an attempt has been made to get it work.Any comments would be helpful. Thanks