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!

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.

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:

• Generate CHANGELOG.md in root folder
• Write change log generation script in docs/make.jl
• Create change log CI to enforce the update of a changelog file on every pull request (see here as an example)

These are the missing Qobj.attribute functions in qutip.

TODO:

• addressed by #150

  • Move all these functions into a new file qobj/arithmetic_and_attributes.jl
  • For those functions that are not attributes in qutip, move to a new file qobj.functions.jl
  • make tr return real value if the input Qobj is Hermitian

• addressed by #155

  • Make sure all existing functions have docstring
  • generalized dot, same as matrix_element
  • isherm
  • diag
  • expm
  • logm
  • purity
  • proj
  • trans
  • unit

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

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)

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.

The current solver for finding the Floquet steady-state uses the matrix continued fractions (see, e.g., Supplementary Materials of ¹ to obtain an effective Liouvillian, and then compute the steady-state to it.

Another method would solve instead the full linear system at once:

$$ ( \mathcal{L}_0 - i n \omega_d ) \hat{\rho}_n + \mathcal{L}_1 \hat{\rho}_{n-1} + \mathcal{L}_{-1} \hat{\rho}_{n+1} = 0 $$

This is a tridiagonal linear system of the form

$$ \mathbf{A} \cdot \mathbf{b} = 0 $$

where

$$ \mathbf{A} = \begin{pmatrix} \mathcal{L}_0 - i (-n_{\text{max}}) \omega_{\textrm{d}} & \mathcal{L}_{-1} & 0 & \cdots & 0 \\ \mathcal{L}_1 & \mathcal{L}_0 - i (-n_{\text{max}}+1) \omega_{\textrm{d}} & \mathcal{L}_{-1} & \cdots & 0 \\ 0 & \mathcal{L}_1 & \mathcal{L}_0 - i (-n_{\text{max}}+2) \omega_{\textrm{d}} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \mathcal{L}_0 - i n_{\text{max}} \omega_{\textrm{d}} \end{pmatrix} $$

and

$$ \mathbf{b} = \begin{pmatrix} \hat{\rho}_{-n_{\text{max}}} \\ \hat{\rho}_{-n_{\text{max}}+1} \\ \vdots \\ \hat{\rho}_{0} \\ \vdots \\ \hat{\rho}_{n_{\text{max}}-1} \\ \hat{\rho}_{n_{\text{max}}} \end{pmatrix} $$

This will allow to simultaneously obtain all the $\hat{\rho}_n$.

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.

References

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
$$\lVert\frac{\partial |\rho\rangle\rangle}{\partial t}\rVert \leq \textrm{reltol} \times\lVert\frac{\partial |\rho\rangle\rangle}{\partial t}+|\rho\rangle\rangle\rVert$$

or

$$\lVert\frac{\partial |\rho\rangle\rangle}{\partial t}\rVert \leq \textrm{abstol}$$

Implement SciMLOperators.MatrixOperator in all ODE-based solvers

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.

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.

Add @async computations in GPU implementation of the time evolution. Benchmark GPU code with NVIDIA Nsight.

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 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).

Add plotting tool for Wigner functions

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:

• KetQuantumObject
• BraQuantumObject
• OperatorQuantumObject
• SuperOperatorQuantumObject
• OperatorKetQuantumObject
• OperatorBraQuantumObject

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

$$ \ket{\psi} = \ket{\alpha} + \ket{e^{i \frac{2}{3} \pi} \alpha} + \ket{e^{i \frac{4}{3} \pi} \alpha} $$

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.

@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:

• Not sure whether we should change the author in Project.toml as qutip develop team or something else
• Modify the links and authors in docs/make.jl
• Modify the links in docs/src/benchmarks/benchmark_history.md
• Modify LISENCE (should we change to the similar LISENCE as qutip ?)
• Modify all the badges' links in README.md
• Deal with codecov and check whether secrets.CODECOV_TOKEN still works
• Check whether secrets.GITHUB_TOKEN still works
• Check whether secrets.DOCUMENTER_KEY still works
• Change repo link in General
• Check whether Package registration and automatic documentation build for version releases still work
• Check whether Benchmark Track still works
• #40

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

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 only saveat is given, then the arguments save_everystep and dense are false by default. If saveat is given a number, then it will automatically expand to tspan[1]:saveat:tspan[2]. For methods where interpolation is not possible, saveat may be equivalent to tstops. 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.

sol.states

In qutip, the results is stored depends on e_ops:

• if e_ops is specified: don't return states at each time point (even the final state, so sol.states = [])
• if e_ops is not specified: return states at each time point

I 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

There are still some functions (compare to qutip) unsupported yet.

TODO:

• #150

  • Move the existing functions and the following ones into an independent file qobj/states.jl
  • set default value for argument pos in basis and fock as 0

• #164

  • 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

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.

There are still some functions (compare to qutip) unsupported yet.

• #150
  • Move the existing functions and the following ones into an independent file qobj/operators.jl
• #156
  • commutator
  • fdestroy
  • fcreate
• #158
  • jmat
  • spin_Jx
  • spin_Jy
  • spin_Jz
  • spin_Jm
  • spin_Jp
  • spin_J_set
• #169
  • displace
  • num
  • position
  • momentum
  • phase
  • squeeze
  • tunneling
  • qft
• Decide not to implement at this moment
  • 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.

I think that a detailed README file is needed.

The file should explain the potentialities of this package

• Syntax very close to the QuTiP one
• Very fast simulations thanks to Julia
• Easy to extend
• Support for GPUs
• Support for Distributed computing

And making also some examples such as

• The master equation in CPU (showing the similarities to QuTiP) and in the GPU (showing how easy is to perform fast simulations)
• The Distributed computing of quantum trajectories

• Support CUDA functions e.g., cu()
• Multiple dispatch for different type of 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      ⋅    

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

Implement symmetry-decomposed evolution for U(1) strong Liouvillian symmetry. Based on this paper.

https://github.com/JuliaGPU/buildkite

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