Errow when showing operators containing symbolic variables

[mabuni1998]

When working with symbolic variables, they combine almost seamlessly with the quantum optics operators (nice!). But when printing the operator, an error occurs.

MWE:

using Symbolics
using QuantumOptics

@variables α
ba = FockBasis(2)
a = dense(destroy(ba))
H1 = α * dagger(a) * a

Producing:

ERROR: MethodError: no method matching round(::Num, ::RoundingMode{:Nearest})
This error has been manually thrown, explicitly, so the method may exist but be intentionally marked as unimplemented.

Closest candidates are:
  round(::Real, ::RoundingMode; digits, sigdigits, base)
   @ Base floatfuncs.jl:50
  round(::Type{T}, ::Any, ::RoundingMode) where T>:Missing
   @ Base missing.jl:148
  round(::Type{T}, ::Any, ::RoundingMode) where T
   @ Base rounding.jl:479
  ...

Stacktrace:
  [1] round(x::Num, r::RoundingMode{:Nearest}; digits::Nothing, sigdigits::Nothing, base::Nothing)
    @ Base .\floatfuncs.jl:56
  [2] round(x::Num, r::RoundingMode{:Nearest})
    @ Base .\floatfuncs.jl:50
  [3] _round_invstep(x::Num, invstep::Num, r::RoundingMode{:Nearest})
    @ Base .\floatfuncs.jl:77
  [4] _round_digits(x::Num, r::RoundingMode{:Nearest}, digits::Int32, base::Int64)
    @ Base .\floatfuncs.jl:115
  [5] round(x::Num, r::RoundingMode{:Nearest}; digits::Int32, sigdigits::Nothing, base::Nothing)
    @ Base .\floatfuncs.jl:68
  [6] round(z::Complex{Num}, rr::RoundingMode{:Nearest}, ri::RoundingMode{:Nearest}; kwargs::@Kwargs{digits::Int32})
    @ Base .\complex.jl:1113
  [7] (::Base.Broadcast.var"#29#30"{@Kwargs{digits::Int32}, typeof(round)})(args::Complex{Num})
    @ Base.Broadcast .\broadcast.jl:1306
  [8] _broadcast_getindex_evalf
    @ .\broadcast.jl:673 [inlined]
  [9] _broadcast_getindex
    @ .\broadcast.jl:646 [inlined]
 [10] getindex
    @ .\broadcast.jl:605 [inlined]
 [11] copy(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, Base.Broadcast.var"#29#30"{@Kwargs{digits::Int32}, typeof(round)}, Tuple{Matrix{Complex{Num}}}})
    @ Base.Broadcast .\broadcast.jl:906
 [12] materialize
    @ .\broadcast.jl:867 [inlined]
 [13] show(stream::IOContext{Base.TTY}, x::Operator{FockBasis{Int64}, FockBasis{Int64}, Matrix{Complex{Num}}})
    @ QuantumOpticsBase C:\Users\mabun\.julia\packages\QuantumOpticsBase\bXZS2\src\printing.jl:79
 [14] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, x::Operator{FockBasis{Int64}, FockBasis{Int64}, Matrix{Complex{Num}}})
    @ Base.Multimedia .\multimedia.jl:47
 [15] (::REPL.var"#68#69"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:367
 [16] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:661
 [17] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:353
 [18] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:372
 [19] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [20] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [21] invokelatest
    @ .\essentials.jl:1052 [inlined]
 [22] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:237
 [23] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [24] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [25] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [26] #67
    @ c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:150 [inlined]
 [27] with_logstate(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging .\logging\logging.jl:522
 [28] with_logger
    @ .\logging\logging.jl:632 [inlined]
 [29] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [30] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [31] invokelatest(::Any)
    @ Base .\essentials.jl:1052
 [32] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:34

The error occurs due to the show function for dense operators:

function show(stream::IO, x::DenseOpType)
    summary(stream, x)
    print(stream, "\n")
    if !_std_order
        if !haskey(stream, :compact)
            stream = IOContext(stream, :compact => true)
        end
        Base.print_array(stream, round.(x.data; digits=machineprecorder))
    else
        showarray_stdord(stream, round.(x.data; digits=machineprecorder), x.basis_l.shape, x.basis_r.shape, false, header=false)
    end
end

where round.(x.data; digits=machineprecorder) throws an error when the data is symbolic.

A similar error occurs for a non-dense operator:

a = destroy(ba)
H1 = α * dagger(a) * a

giving:

Closest candidates are:
  round(::Real, ::RoundingMode; digits, sigdigits, base)
   @ Base floatfuncs.jl:50
  round(::Type{T}, ::Any, ::RoundingMode) where T>:Missing
   @ Base missing.jl:148
  round(::Type{T}, ::Any, ::RoundingMode) where T
   @ Base rounding.jl:479
  ...

Stacktrace:
  [1] round(x::Num, r::RoundingMode{:Nearest}; digits::Nothing, sigdigits::Nothing, base::Nothing)
    @ Base .\floatfuncs.jl:56
  [2] round(x::Num, r::RoundingMode{:Nearest})
    @ Base .\floatfuncs.jl:50
  [3] _round_invstep(x::Num, invstep::Num, r::RoundingMode{:Nearest})
    @ Base .\floatfuncs.jl:77
  [4] _round_digits(x::Num, r::RoundingMode{:Nearest}, digits::Int32, base::Int64)
    @ Base .\floatfuncs.jl:115
  [5] round(x::Num, r::RoundingMode{:Nearest}; digits::Int32, sigdigits::Nothing, base::Nothing)
    @ Base .\floatfuncs.jl:68
  [6] round(z::Complex{Num}, rr::RoundingMode{:Nearest}, ri::RoundingMode{:Nearest}; kwargs::@Kwargs{digits::Int32})
    @ Base .\complex.jl:1113
  [7] (::Base.Broadcast.var"#29#30"{@Kwargs{digits::Int32}, typeof(round)})(args::Complex{Num})
    @ Base.Broadcast .\broadcast.jl:1306
  [8] _noshapecheck_map
    @ C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\SparseArrays\src\higherorderfns.jl:180 [inlined]
  [9] copy(bc::Base.Broadcast.Broadcasted{SparseArrays.HigherOrderFns.SparseMatStyle, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, Base.Broadcast.var"#29#30"{@Kwargs{digits::Int32}, typeof(round)}, Tuple{SparseArrays.SparseMatrixCSC{Complex{Num}, Int64}}})
    @ SparseArrays.HigherOrderFns C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\SparseArrays\src\higherorderfns.jl:195
 [10] materialize
    @ .\broadcast.jl:867 [inlined]
 [11] show(stream::IOContext{Base.TTY}, x::Operator{FockBasis{Int64}, FockBasis{Int64}, SparseArrays.SparseMatrixCSC{Complex{Num}, Int64}})
    @ QuantumOpticsBase C:\Users\mabun\.julia\packages\QuantumOpticsBase\bXZS2\src\printing.jl:94
 [12] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, x::Operator{FockBasis{Int64}, FockBasis{Int64}, SparseArrays.SparseMatrixCSC{Complex{Num}, Int64}})
    @ Base.Multimedia .\multimedia.jl:47
 [13] (::REPL.var"#68#69"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:367
 [14] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:661
 [15] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:353
 [16] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\mabun\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:372
 [17] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [18] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [19] invokelatest
    @ .\essentials.jl:1052 [inlined]
 [20] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:237
 [21] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [22] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [23] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [24] #67
    @ c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:150 [inlined]
 [25] with_logstate(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)       
    @ Base.CoreLogging .\logging\logging.jl:522
 [26] with_logger
    @ .\logging\logging.jl:632 [inlined]
 [27] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [28] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [29] invokelatest(::Any)
    @ Base .\essentials.jl:1052
 [30] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\mabun\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\eval.jl:34

[Krastanov]

Hi, Matias! I am surprised that any of this actually works. I expect there will be unexpected breakages in many places.

Making the printing more robust would of course be good, but there probably will be more issues that pop up if this is used. I think there are a few action items here, of increasing cost in time and effort. Listing them here just to keep it logged (not sure anyone will be tackling them near term):

• Could you list in what circumstances you have found Symbolics + QuantumOptics useful together? A couple of usecases would inform how we can think about working on this.
• Have you checked QuantumSymbolics.jl and its express interface used to transform symbolic expressions to numeric expressions (including QuantumOptics)? It might be a better tool for your goals.
• Making printing of these objects more robust.
• Making a bunch of tests (presumably some failing, and marked as broken) of various tasks related to Symbolics + QuantumOptics, just so we do not break more things in the future (if they happen to work by accident right now)

Tagging @apkille as he has done a lot of work on QuantumSymbolics

[apkille]

That's really cool that it works. I remember talking to @jgr-rgb at QNumerics about implementing such a feature for QuantumSymbolics. Perhaps he has some thoughts on this? I don't have any great immediate use cases, but I'll give some thought to it over the next few days. Perhaps being able to calculate symbolic expressions for metrics of states might be useful. This is less on the state-vector side and more for something like Gabs.jl, but analytical calculations involving covariance matrices of two-mode Gaussian bosonic systems tend to be pretty useful (see e.g. this paper and this paper for lots of discussion). Or, maybe using symbolic matrix representations of spin systems for studying like the Ising model or NMR experiments is somehow worthwhile.

Regarding the printing, I'm linking #45 as it's related. In case anyone who's reading this isn't aware, it might be tricky to deal with the printing problem without messing with QOBase.set_printing. Defining a simpler method for Base.show as follows

function Base.show(io::IO, mime::MIME"text/plain", x::DenseOpType)
    summary(io, x)
    print(io, "\n")
    Base.show(io, mime, x.data)
end

gives

julia> α * dagger(a) * a
Operator(dim=3x3)
  basis: Fock(cutoff=2)
3×3 Matrix{Complex{Num}}:
 0.0  0.0   0.0
 0.0    α   0.0
 0.0  0.0  2.0α

but that removes the feature of modifying REPL output with QO types.

[mabuni1998]

Wow, thanks for the input! I was not aware of QuantumSymbolics.jl and will check it out.

As for use-cases of symbolic quantum operators, I have been playing around with symbolically deriving the equations of motion of a master equation and translating the resulting equations into an ODE (either with custom function or using modeling toolkits). This serves two purposes: First of all, auto-diff-related things basically work out of the box. Secondly, the function that returns the differential df(du,u,p,t) is compiled and should be non-allocating even if you have time-dependent operators. Thus, in general, I have observed that df(du,u,p,t) is faster than applying the quantum operator itself. The following demonstrates a very basic example of what I am thinking (note how well symbolic quantum operators work).

# 1. Define symbolic parameters for the system

using QuantumOptics, ModelingToolkit, Symbolics
using ModelingToolkit: t_nounits as t, D_nounits as D
@parameters ω γ  # ω is the frequency, γ is the dissipation rate

# Pauli lowering and raising operators
basis = NLevelBasis(2)
bc = FockBasis(1)
σ₋ = destroy(bc)
σ₊ =  create(bc)

# Hamiltonian for the Jaynes-Cummings model (symbolic ω)
H = ω * (σ₊ + σ₋);

#Define function that returns a symbolic representation of the master equation
function me_symbolic(t,H, J,Jdagger=dagger.(J))
    ρ = symbolic_density(H,t);
    # Define the commutator [H, ρ] and the Lindblad terms (symbolic γ)
    H_comm = H * ρ - ρ * H;
    # Lindblad dissipator term: γ (σ₋ρσ₊ - 0.5 * (σ₊σ₋ρ + ρσ₊σ₋))
    master_eq = -1im * H_comm;
    for (i,j) in enumerate(J)
        master_eq = master_eq + ( j* ρ * Jdagger[i] - 0.5 * (Jdagger[i] * j * ρ + ρ * Jdagger[i] * j));
    end
    
    return ρ ,master_eq
end

ρ_symbolic ,master_eq = me_symbolic(t,H, [γ  * σ₋],[γ * σ₊]);

Now printing

julia> master_eq.data
2×2 Matrix{Complex{Num}}:
                                         -2ρ₁ˏ₂_im(t)*ω + ρ₂ˏ₂_re(t)*(γ^2)      -0.5ρ₁ˏ₂_re(t)*(γ^2) + im*(-ρ₂ˏ₂_re(t)*ω + ρ₁ˏ₁_re(t)*ω - 0.5ρ₁ˏ₂_im(t)*(γ^2))
 -0.5ρ₁ˏ₂_re(t)*(γ^2) + im*(ρ₂ˏ₂_re(t)*ω - ρ₁ˏ₁_re(t)*ω + 0.5ρ₁ˏ₂_im(t)*(γ^2))                                            

Thus, we now have the master equation equations of motion. The density matrix elements are generated using:

Here ρ = symbolic_density(H,t); is a symbolic density matrix that is generated like this:

#The following are functions that help construct the density matrix as a symbolic matrix.
function variable_complex(name,t, idx...)
    if idx[1] == idx[2]
        name_ij_re = Symbol(name, join(Symbolics.map_subscripts.(idx), "ˏ"), "_re")
        vre = Sym{Symbolics.FnType{Tuple{Real}, Real}}(name_ij_re)(t)
        return Complex(Num(setmetadata(vre, Symbolics.VariableSource, (:variables, name_ij_re))))
    elseif idx[1]<idx[2]
        name_ij_re = Symbol(name, join(Symbolics.map_subscripts.(idx), "ˏ"), "_re")
        name_ij_im = Symbol(name, join(Symbolics.map_subscripts.(idx), "ˏ"), "_im")
        vre = Sym{Symbolics.FnType{Tuple{Real}, Real}}(name_ij_re)(t)
        vim = Sym{Symbolics.FnType{Tuple{Real}, Real}}(name_ij_im)(t)
        return Complex(Num(setmetadata(vre, Symbolics.VariableSource, (:variables, name_ij_re))), Num(setmetadata(vim, Symbolics.VariableSource, (:variables, name_ij_im))))
    else
        name_ij_re = Symbol(name, join(Symbolics.map_subscripts.(reverse(idx)), "ˏ"), "_re")
        name_ij_im = Symbol(name, join(Symbolics.map_subscripts.(reverse(idx)), "ˏ"), "_im")
        vre = Sym{Symbolics.FnType{Tuple{Real}, Real}}(name_ij_re)(t)
        vim = Sym{Symbolics.FnType{Tuple{Real}, Real}}(name_ij_im)(t)
        return Complex(Num(setmetadata(vre, Symbolics.VariableSource, (:variables, name_ij_re))), -Num(setmetadata(vim, Symbolics.VariableSource, (:variables, name_ij_im))))
    end
end
function variables_complex(name,t, indices...)
    [variable_complex(name,t, ij...) for ij in Iterators.product(indices...)]
end
function symbolic_density(ρ,t)
    dims = size(ρ)
    Operator(ρ.basis_l,variables_complex(:ρ,t,1:dims[1],1:dims[2])) 
end

We can convert the master equations of motion to a modeling toolkits ODE using:

#Function that returns the master equation as a set of equations compatible with modelling toolkit
function convert_to_ode(ρ,master_eq)
    lhs_diag = zeros(Complex{Num}, size(ρ.data, 1))
    rhs_diag = zeros(Complex{Num}, size(ρ.data, 1))
    lhs_real = zeros(Complex{Num}, size(ρ.data, 1)*(size(ρ.data, 2)-1)÷2)
    rhs_real = zeros(Complex{Num}, size(ρ.data, 1)*(size(ρ.data, 2)-1)÷2)
    lhs_imag = zeros(Complex{Num}, size(ρ.data, 1)*(size(ρ.data, 2)-1)÷2)
    rhs_imag = zeros(Complex{Num}, size(ρ.data, 1)*(size(ρ.data, 2)-1)÷2)
    

    for i in 1:size(ρ.data, 1)  # Loop over rows
        lhs_diag[i] = D(ρ.data[i, i].re)
        rhs_diag[i] = master_eq.data[i, i].re
        for j in i+1:size(ρ.data, 2)  # Loop over columns
            idx = ((j - 2) * (j-1)) ÷ 2 + i 
            lhs_real[idx] = D(ρ.data[i, j].re)
            rhs_real[idx] = master_eq.data[i, j].re
            lhs_imag[idx] = D(ρ.data[i, j].im)
            rhs_imag[idx] = master_eq.data[i, j].im
        end
    end
    ode_eqs = vcat(lhs_diag .~ rhs_diag,lhs_real .~ rhs_real,lhs_imag .~ rhs_imag)
    return ode_eqs 
end


eqs = convert_to_ode(ρ_symbolic,master_eq)

#Build model and solve
@mtkbuild model = ODESystem(eqs, t)
u0=[1.0,0,0,0]
prob = ODEProblem(model, u0, (0.0, 10.0),[1,1])

using DifferentialEquations
solve(prob)

Whether this is actually a good approach, I am not sure. I have tested that autodiff works, and also, initial benchmarks show that this is faster than solving the system using the quantum operators. This is, however, only when everything has been compiled, which is VERY slow for large quantum systems (I have found that using custom code that splits the large odefunction up into smaller bits significantly lowers compile time while not suffering in performance). It seems a general problem with ModelingToolkits.jl when you have a large system of equations.