qutip / quantumtoolbox.jl Goto Github PK
View Code? Open in Web Editor NEWQuantum Toolbox in Julia
Home Page: https://qutip.org/QuantumToolbox.jl/dev
License: BSD 3-Clause "New" or "Revised" License
Quantum Toolbox in Julia
Home Page: https://qutip.org/QuantumToolbox.jl/dev
License: BSD 3-Clause "New" or "Revised" License
The ProgressBar
was just implemented. However, it currently supports only sesolve
and mesolve
. The support for mcsolve is not straightforward since it has to support Distributed EnsembleProblem is a type stable way. For the moment, the use of RemoteChannel
to perform this, seems to make the ODE type unstable.
There are still some functions (compare to qutip
) unsupported yet.
qobj/operators.jl
commutator
fdestroy
fcreate
jmat
spin_Jx
spin_Jy
spin_Jz
spin_Jm
spin_Jp
spin_J_set
displace
num
position
momentum
phase
squeeze
tunneling
qft
identity
: this will cause conflict with Base.identity
qdiags
: this can be implement easily using Qobj
together with LinearAlgebra.diagm
or SparseArrays.spdiagm
in Julia
squeezing
: I'm not sure whether the current definition for generalized squeeze operator in QuTiP
is common in the literature. So I decide not to implement it for now.charge
: I didn't find a common definition of this, maybe this can be implemented if someone needs it in the future.To make the syntax more similar to QuTiP
, I suggest to add the following function support:
Qobj(A; kwargs...) = QuantumObject(A; kwargs...)
: This is just for us to create QuantumObject
in a lazy way (and more similar to QuTiP
). (f3c95b5)tensor(A...) = kron(A...)
: does the same thing as kron
(9949e38)⊗(A, B) = kron(A, B)
: take advantage of QuantumOptics
(9949e38)The last support allows us to express tensor(A, B, C)
as A ⊗ B ⊗ C
.
If @albertomercurio agrees with this, I can make a PR for this issue.
Consider a new eigen solve interface that can return not just eigenvalues and (unitary) tranformation matrix, but also the Vector{QuantumObject}
.
Also, support the similar syntax like QuTiP
: eigenenergies()
and eigenstates()
, and also make a keyword argument [like sparse=true
(or false
) in QuTiP
] as a option for using sparse eigensolvers or not.
Finally, also should support to solve the eigen problems for Liouvillian, which the return eigenstates should be operator-ket
after the issue #37 is solved
There are still some functions (compare to qutip
) unsupported yet.
TODO:
qobj/states.jl
pos
in basis
and fock
as 0
zero_ket
fock_dm
coherent_dm
thermal_dm
maximally_mixed_dm
spin_state
spin_coherent
bell_state
singlet_state
triplet_states
w_state
ghz_state
Implement symmetry-decomposed evolution for U(1) strong Liouvillian symmetry. Based on this paper.
A matrix with size=(n, 1)
should also be considered as a Ket
But currently it shows an error, for example:
julia > v = Matrix([2 3])' # this is a vector but in `Matrix`-type
julia > size(v)
(2, 1)
julia > Qobj(v)
DimensionMismatch: matrix is not square: dimensions are (2, 1)
Also, I suggest to remove the following support:
julia > v = [1, 2]
julia > isbra(Qobj(v, type=BraQuantumObject))
true
User should do this instead:
julia > isbra(Qobj(v'))
true
In QuTiP, the superoperators have different representations:
Maybe we can make a new abstract type AbstractSuperRepresentation
, and re-define the superoperator:
abstract type AbstractSuperRepresentation end
struct GeneralSuperRepresentation <: AbstractSuperRepresentation end
struct ChoiSuperRepresentation <: AbstractSuperRepresentation end
struct KrausSuperRepresentation <: AbstractSuperRepresentation end
struct StinespringSuperRepresentation <: AbstractSuperRepresentation end
struct ChiSuperRepresentation <: AbstractSuperRepresentation end
struct SuperOperatorQuantumObject <: QuantumObjectType
superrep::AbstractSuperRepresentation
end
# default representation
SuperOperatorQuantumObject() = SuperOperatorQuantumObject(GeneralSuperRepresentation())
const SuperOperator = SuperOperatorQuantumObject()
const ChoiSuperOperator = SuperOperatorQuantumObject(ChoiSuperRepresentation())
const KrausSuperOperator = SuperOperatorQuantumObject(KrausSuperRepresentation())
const StinespringSuperOperator = SuperOperatorQuantumObject(StinespringSuperRepresentation())
const ChiSuperOperator = SuperOperatorQuantumObject(ChiSuperRepresentation())
When the users apply the basic operations (:+
, :-
, :*
) between GeneralSuperRepresentation
and other representation, it should work. However, if the users apply them on different representations, it should print a warning message just like qutip.
function Base.:(==)(S1::SuperOperatorQuantumObject, S2::SuperOperatorQuantumObject)
if (S1 != S2)
if !(S1 == SuperOperator || S2 == SuperOperator)
# print warning
end
end
return true
end
Also need to define some new functions same as qutip:
superrep(Q)
(return superrep
of Qobj
)to_super(Q)
to_choi(Q)
to_kraus(Q)
to_stinespring(Q)
to_chi(Q)
iscp(Q)
istp(Q)
ishp(Q)
iscptp(Q)
There is only one steadystate
solver: the direct solver.
It could be very useful to add other solvers, such as all the iterative solvers that are available through the LinearSolve.jl package and the internal eigensolver.
The eigensolver uses the internal sparse solver, that uses the Krylov subspace, taking only the eigenstate with zero eigenvalue.
Perhaps, one can also add the support for a cache (it could be very easy using the LinearSolve.jl package).
Make a logo for the package, possibly using the QuantumToolbox.jl
package.
What I have in mind is for example creating the Wigner function a triangular cat state
and possibly mapping the colors following the Julia logo.
Then, we can add a documentation page demonstrating how this logo was made, showing step by step the code.
Support operator-bra and operator-ket like the ones in QuTiP
.
Note that we wish to remain the dims::Vector{Int}
to store only the Hilbert space dimension.
That is, if we have H = tensor(destroy(2), destroy(3))
. Part of the output information of it's Liouvillian should be
julia > liouvillian(H)
Quantum Object: type=SuperOperator dims=[2, 3] size=(36, 36)
and so does OperatorBraQuantumObject
and OperatorKetQuantumObject
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
We should try to make struct QuantumObject
(immutable).
This allows us to have backward compatibility for Julia v1.6
, which includes the Long-term support (LTS) release.
Implement SciMLOperators.MatrixOperator
in all ODE-based solvers
Add @async
computations in GPU implementation of the time evolution. Benchmark GPU code with NVIDIA Nsight.
These are the missing Qobj.attribute
functions in qutip
.
TODO:
addressed by #150
qobj/arithmetic_and_attributes.jl
qutip
, move to a new file qobj.functions.jl
tr
return real value if the input Qobj
is Hermitianaddressed by #155
dot
, same as matrix_element
isherm
diag
expm
logm
purity
proj
trans
unit
The current solver for finding the Floquet steady-state uses the matrix continued fractions (see, e.g., Supplementary Materials of 1 to obtain an effective Liouvillian, and then compute the steady-state to it.
Another method would solve instead the full linear system at once:
This is a tridiagonal linear system of the form
where
and
This will allow to simultaneously obtain all the
The steadystate_floquet
function has already the argument lf_solver
that can be used to specify which method to use to get the floquet steadystate.
Macrì, Vincenzo, et al. "Spontaneous scattering of Raman photons from cavity-QED systems in the ultrastrong coupling regime." Physical Review Letters 129.27 (2022): 273602. ↩
The current existing steadystate
solving methods doesn't depend on the initial condition (state), but there are some cases which leads to different stationary states because of the initial states.
By the multiple-dispatch
feature of Julia
, we can have a new steadystate
method with a given initial state.
The stationary state can be obtained based on DiffEqCallbacks.jl
and the solvers in OrdinaryDiffEq.jl
.
For the terminate condition of finding the stationary state, it would be better to use:
either the following condition is true
or
Hi, folks! How receptive is the team of the following two suggestions:
Support the function permute
in QuTiP
, something like:
function permute(A::QuantumObject, order::Vector{Int}) ... end
a = Qobj(rand(2, 2))
b = Qobj(rand(3, 3))
c = Qobj(rand(4, 4))
d = Qobj(rand(5, 5))
abcd = tensor(a, b, c, d)
permute(abcd, [4, 1, 3, 2]) == tensor(d, a, c, b) # true
Should support for the following types of quantum object:
We should gather all the import
and using
keywords to the file src/QuantumToolbox.jl
.
We should separate all the export
keyword to individual files where the structures and functions are defined.
@albertomercurio
Here is a summary for the things we need to address after moving the repository into qutip.org. Let me know if there are other things we need to do.
TODO list:
Project.toml
as qutip develop team or something elsedocs/make.jl
docs/src/benchmarks/benchmark_history.md
LISENCE
(should we change to the similar LISENCE as qutip ?)README.md
secrets.CODECOV_TOKEN
still workssecrets.GITHUB_TOKEN
still workssecrets.DOCUMENTER_KEY
still worksGeneral
Support :
julia> QuantumToolbox.versioninfo()
But don't export it, otherwise it will have conflict with Base.versioninfo()
Here is an example output from my package:
__
/ \
__ __ __ \__/ __
| | | | / \ / \
| | | | ______ ______ \__/_ _\__/
| |___| |/ __ \ / __ \ / ' \/ \
| ___ | |__) | / \ | _ _ |
| | | | ____/| ( ) | / \ / \ |
| | | | |____ | \__/ | | | | | |
|__| |__|\______) \______/|__| |_| |_|
Julia framework for Hierarchical Equations of Motion
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © NCKU-QFORT 2022 and later.
Lead developer : Yi-Te Huang
Other developers:
Simon Cross, Neill Lambert, Po-Chen Kuo and Shen-Liang Yang
Package information:
====================================
HierarchicalEOM Ver. 1.4.0
LinearSolve Ver. 2.22.1
OrdinaryDiffEq Ver. 6.69.0
SteadyStateDiffEq Ver. 2.0.1
SciMLOperators Ver. 0.3.7
FastExpm Ver. 1.1.0
JLD2 Ver. 0.4.41
ProgressMeter Ver. 1.9.0
PrecompileTools Ver. 1.2.0
System information:
====================================
Julia Version: 1.10.0
OS : Linux (x86_64-linux-gnu)
CPU : 144 × Intel(R) Xeon(R) Gold 6240L CPU @ 2.60GHz
Memory : 256 GB
WORD_SIZE: 64
LIBM : libopenlibm
LLVM : libLLVM-15.0.7 (ORCJIT, cascadelake)
BLAS : libopenblas64_.so (ilp64)
source code link: https://github.com/NCKU-QFort/HierarchicalEOM.jl/blob/main/src/HeomBase.jl#L167
The ASCII art (logo) is unnecessary lol.
But the package compat and system (platform) info might be important (maybe can also add the thread number stuff same as standard julia versioninfo
output)
I noticed that the current sesolve
and mesolve
will only return the final state corresponds to the given time list parameter (tl
). For example,
> using QuantumToolbox
> H = sigmaz()
> ψ0 = sqrt(0.5) * (basis(2, 0) + basis(2, 1))
> tlist = 0:1:10
> sol = sesolve(H, ψ0, tlist)
> sol.times
# this returns a vector which contains single value 10 (last element in tlist)
> sol.states
# this returns a vector which contains only the final state
I expect to obtain:
> sol.times
# all the elements in tlist
> sol.states
# all the states corresponding to each time step
I'm not familiar with the DiffEqCallbacks.jl
API, but I think it must be related to the parameter saveat
, as said in the docs from DifferentialEquations.jl
:
saveat
: Denotes specific times to save the solution at, during the solving phase. The solver will save at each of the timepoints in this array in the most efficient manner available to the solver. If onlysaveat
is given, then the argumentssave_everystep
anddense
arefalse
by default. Ifsaveat
is given a number, then it will automatically expand totspan[1]:saveat:tspan[2]
. For methods where interpolation is not possible,saveat
may be equivalent totstops
. The default value is[]
.
So currently, we need to do it like this
> sol = sesolve(H, ψ0, tlist; saveat=tlist)
which needs to provide the tlist
twice and becomes a bit redundant.
cu()
Arrays
We can define the function as
cu(A::QuantumObject{<:AbstractArray{T}}; word_size::Int=32) where T<:Number
The keyword argument word_size
can decide whether to change the element type for the output.
For example,
julia> H = sigmax()
Quantum Object: type=Operator dims=[2] size=(2, 2) ishermitian=true
2×2 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
⋅ 1.0+0.0im
1.0+0.0im ⋅
julia> cu(H)
Quantum Object: type=Operator dims=[2] size=(2, 2) ishermitian=true
2×2 CuSparseMatrixCSC{ComplexF32, Int32} with 2 stored entries:
⋅ 1.0+0.0im
1.0+0.0im ⋅
julia> cu(H, word_size=64)
Quantum Object: type=Operator dims=[2] size=(2, 2) ishermitian=true
2×2 CuSparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
⋅ 1.0+0.0im
1.0+0.0im ⋅
Our package still got many things to do, but once we are close to a standard release. We can start to think about writing change logs by utilizing Changelog.jl
.
We can follow the same structure in Documenter.jl
.
TODO:
CHANGELOG.md
in root folderdocs/make.jl
sol.states
In qutip
, the results is stored depends on e_ops
:
e_ops
is specified: don't return states at each time point (even the final state, so sol.states = []
)e_ops
is not specified: return states at each time pointI think we should also do the same thing.
However, if the keyword argument saveat
is given, we can still store the state even e_ops
is given
Currently, the following code works:
julia> using QuantumToolbox
julia> fock(2, 0; dims = [2, 3])
Quantum Object: type=Ket dims=[2, 3] size=(2,)
2-element Vector{ComplexF64}:
1.0 + 0.0im
0.0 + 0.0im
However, this makes no sense and should raise an error.
I suggest to call
QuantumObject(array; type = Ket, dims = dims)
instead of
QuantumObject(array, Ket, dims)
in fock
function, cause the previous one will do argument check.
Also, we should add additional test_throws
checks for the above invalid cases
The simple code below fails when running the type-inference test
using Test
using QuantumToolbox
A = sprand(ComplexF64, 100, 100, 0.1)
@inferred QuantumObject(A)
return type QuantumObject{SparseMatrixCSC{ComplexF64, Int64}, OperatorQuantumObject} does not match inferred return type Union{QuantumObject{SparseMatrixCSC{ComplexF64, Int64}, BraQuantumObject}, QuantumObject{SparseMatrixCSC{ComplexF64, Int64}, KetQuantumObject}, QuantumObject{SparseMatrixCSC{ComplexF64, Int64}, OperatorQuantumObject}}
I think that we should fix it using multiple dispatch when calling it.
This might be too breaking to be feasible, but I think it is worth considering defining tensor(a, b) = kron(b, a)
(and in general tensor(args...) = kron(reverse(args)...)
), essentially because the definition of kron
is (I believe) only natural for row-oriented languages like C and python, and not column-oriented languages like Julia/fortran/MATLAB, and does not interact nicely with reshape
, permute
, etc in a column-oriented language (as e.g. mentioned here).
I will add that reversing the order like I'm advocating is the convention used in QuantumOptics.jl (see also qojulia/QuantumOptics.jl#177 (comment)).
Let me give some more examples about why it would help. I'll define kron_rev
here for clarity, and use kron
and kron_rev
.
kron_rev(a, b) = kron(b, a)
In section 1.1.3 of Watrous's TQI, he defines the vec
operator:
then:
N = 5
e(i) = (z = zeros(N); z[i] = 1; z) # quick basis vector
E(i,j) = (z = zeros(N, N); z[i, j] = 1; z)
vec(E(1,2)) == kron(e(1), e(2)) # false
vec(E(1,2)) == kron_rev(e(1), e(2)) # true
That is, if ⊗
were kron_rev
, then this identity holds; if not, it does not. This has a bunch of consequences. For example, we have the identity:
and again we need kron_rev
for this to hold:
julia> kron(A₀, A₁) * vec(B) ≈ vec(A₀*B*transpose(A₁))
false
julia> kron_rev(A₀, A₁) * vec(B) ≈ vec(A₀*B*transpose(A₁))
true
This issue also has this example:
julia> A = rand(2,2); B = rand(2,2);
julia> AB = reshape(kron(A, B), (2,2,2,2));
julia> all(AB[a1,b1,a2,b2] ≈ A[a1,a2] * B[b1,b2] for a1=1:2, a2=1:2, b1=1:2, b2=1:2)
false
whereas
julia> AB = reshape(kron_rev(A, B), (2,2,2,2));
julia> all(AB[a1,b1,a2,b2] ≈ A[a1,a2] * B[b1,b2] for a1=1:2, a2=1:2, b1=1:2, b2=1:2)
true
Add plotting tool for Wigner functions
I think that a detailed README
file is needed.
The file should explain the potentialities of this package
And making also some examples such as
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.