ModelingToolkit? about quantumcumulants.jl HOT 14 CLOSED

Pinging @shashi here.

from quantumcumulants.jl.

I'm not sure how to implement noncommutative algebra in combination with a commutative one in MTK. It seemed to me as if that was exactly what SymbolicUtils was for: defining your own algebra and simplification. I'm all for using MTK for more parts in this package (specifically code generation).

For reference I'll briefly go through everything I need for Qumulants.jl:

1. Operators (or QNumbers) are used to define the system. These form the noncommutative algebra, but you need to be able to combine them with (symbolic) numbers as well.
2. Equations of motion for the operators are derived by resolving fundamental commutation relations. Currently, this is done by converting to SymbolicUtils and simplifying with a custom set of rules.
3. The noncommutative equations of motion are converted to a system of symbolic, commutative ODEs by using the cumulant expansion approach.
4. Numerical code usable with OrdinaryDiffeq is generated.

Now clearly, in step 3. I could directly convert to MTK. This would simply remove step 4. by the much more extensive functionality of MTK. Note that I am planning on doing this in the future. Also, this would replace the custom Parameter implementation in Qumulants by MTK parameters.

However, my problem is with steps 1, and 2. I'd love to leave as much to MTK as possible, but I don't see how I can implement noncommutative algebra in MTK and have it as convenient as it currently is using SymbolicUtils. At the very least, I think I'll need SymbolicUtils to define the commutation relations used for simplification. So from my point of view, I need to use SymbolicUtils.jl in any case, and MTK is just not used yet.

Please let me know if I'm missing something here, or if you see more parts that can be left to MTK.

from quantumcumulants.jl.

I'm not sure how to implement noncommutative algebra in combination with a commutative one in MTK. It seemed to me as if that was exactly what SymbolicUtils was for: defining your own algebra and simplification. I'm all for using MTK for more parts in this package (specifically code generation).

MTK was related to SymbolicUtils, and now with v4 its built on SymbolicUtils, so in theory we should be able to get all of this to work. Anything that we can't get to work we should consider an issue. MTK uses Num{Sym} <: Number so that tracing works as before, but IIUC everything internally is just Sym, so if you use non-commutative Sym it should work as well? That's a question for @shashi . I think we need to get non-commutative algebras working in MTK anyways since we want array/matrix semantics for variables, so let's work together to try and get their. On your end, you'll get a lot of free parallelism and all of that, so that should be a nice feature boost. On our end, I'd think that if we figure out how to make your case work then a lot of these "weirder" cases could similarly be MTK DSLs, and it would be a good example to show something that's particularly difficult to do in other systems.

from quantumcumulants.jl.

tracing

BTW what chris means here is being able to call foo(x::Number) with a symbolic object. I think for the purposes of Qumulants we just don't need that. Also MTK types never store Nums inside them, or restrict inputs to be Num, there's always methods on Terms and Syms, so you should be able to just construct MTK types such as Equations and systems directly with your Q variables!

Once you get to that point let us know if something does not work the way you expect it to! But yes, let's first address JuliaSymbolics/SymbolicUtils.jl#134

from quantumcumulants.jl.

Alright, I guess I was not aware of how extensive the latest changes to MTK were. Just to make sure I'm starting off in the right direction and I understand your suggestions correctly, I drafted up some code which I think has the minimal custom functionality I need on top of SymbolicUtils and MTK. Note that I skipped the definition of the custom rules to shorten the code, so it won't actually run.

using SymbolicUtils
using ModelingToolkit

### QNumber types
abstract type QNumber end
struct Destroy <: QNumber end
struct Create <: QNumber end
Base.one(::Type{<:QNumber}) = 1

### Interface with SymbolicUtils
import SymbolicUtils: islike
islike(::SymbolicUtils.Symbolic{<:QNumber}, ::Type{<:Number}) = true
islike(::SymbolicUtils.Term{<:QNumber}, ::Type{<:Number}) = true
islike(::SymbolicUtils.Symbolic{<:SymbolicUtils.FnType{<:Any,<:QNumber}}, ::Number) = true

# Type promotion
SymbolicUtils.promote_symtype(f, Ts::Type{<:QNumber}...) = promote_type(QNumber,Ts...)
SymbolicUtils.promote_symtype(f, T::Type{<:QNumber}, Ts...) = promote_type(QNumber,T)
for f in [+,-,*,/,^]
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)),
                   T::Type{<:QNumber},
                   S::Type{<:Number}) = QNumber
    @eval SymbolicUtils.promote_symtype(::$(typeof(f)),
                   T::Type{<:Number},
                   S::Type{<:QNumber}) = QNumber
end

### Functions needed for simplification
# Handle noncommutative multiplication
iscommutative(x::SymbolicUtils.Symbolic{<:QNumber}) = false
iscommutative(x::QNumber) = false
iscommutative(x::SymbolicUtils.Symbolic{<:Number}) = true
iscommutative(x::Number) =  true

needs_sorting_nc(x) = (x.f === (*)) && !issorted_nc(x)
function issorted_nc(x)
    args = SymbolicUtils.arguments(x)
    is_c = iscommutative.(args)
    args_c = @views(args[is_c])
    args_nc = @views(args[.!is_c])
    return issorted(is_c, lt=(>)) && SymbolicUtils.issortedₑ(args_c)
end

using SymbolicUtils: <ₑ
function sort_args_nc(x)
    args = SymbolicUtils.arguments(x)
    is_c = iscommutative.(args)
    args_c = sort(args[is_c], lt=(<ₑ))
    args_nc = args[.!is_c]
    return term(*, args_c..., args_nc...)
end

### Rules
[..] # define default_qnumber_simplifier

### Test on example
@parameters t
@derivatives D'~t
@syms a(::Real)::Destroy b(::Real)::Create ω::Real
H = ω*b(t)*a(t)
commutator(a,b) = simplify(a*b - b*a; rewriter=default_qnumber_simplifier())
eq = D(a(t)) ~ commutator(H, a(t)) # Heisenberg(ish) equation

### Conversion to c-number
using SymbolicUtils: Symbolic, arguments, FnType, symtype
average(x::Number) = x
average(x::Sym{T}) where T<:Number = (s = Sym{ModelingToolkit.Parameter{T}}(x.name); Num(s))
average(x::Term{T}) where T<:Number = println(x)
function average(x::Term{<:QNumber})
    f = average(x.f)
    args = []
    for arg in arguments(x)
        push!(args, average(arg))
    end
    return term(f, args...)
end
average(f::Function) = f
average(f::Sym{<:FnType{Tuple{Vararg{T,N}},<:Destroy}}) where {T,N} =
    Sym{FnType{Tuple{Vararg{T,N}},Real}}(f.name)
average(eq::Equation) = Equation(average(eq.lhs), average(eq.rhs))

eq_c = average(eq)
sys = ODESystem([eq_c])

This basically has all the steps I need, though the conversion to c-number equations is actually a bit more complex. So all I need is a proper type system that I can feed into SymbolicUtils Syms. The rest should really work "out of the box" with MTK. I already noticed some issues, however:

• The biggest one is that I actually need support for complex numbers in MTK - the above example worked only because I left out an im that should really be there.
• In the example ω is a parameter. I wanted to use an MTK parameter, but this wraps the result of e.g. ω*a(t) in a Num, which causes assert_like to fail in SymbolicUtils when computing e.g. ω*a(t)*b(t). Not sure what the best way to deal with this would be. Also, I messed up in the conversion so sys.ps is empty.

from quantumcumulants.jl.

Thanks for that summary. We will try to get the complex number support in MTK. That has definitely been on our list.

from quantumcumulants.jl.

@ChrisRackauckas Great, thanks!

I also have another question for @shashi: I need to store some additional information in the symbolics representing my QNumbers, e.g. some integers or symbols. This information needs to be accessible (and mutable) during simplification. What do you think is the best way of achieving this? Thanks!

from quantumcumulants.jl.

would this work for you: JuliaSymbolics/SymbolicUtils.jl#122 ?

from quantumcumulants.jl.

Yes, I think this would be ideal! Another way I can think of doing this is by storing the additional info I need as function arguments when representing qnumbers as symbolic functions. However, this will lead to (some) symbols with ~5-10 arguments. IIUC the simplification will step into each function and try to simplify the arguments, which is unnecessary. So ideally I'd use JuliaSymbolics/SymbolicUtils.jl#122 for better performance.

from quantumcumulants.jl.

I'm finally coming back to this issue. In #12 I implemented a possible (partial) solution, that gives me all the goodies from MTK (jacobians, parallelized functions, etc.). Basically Qumulants just takes the role of defining a set of symbolic equations before an ODESystem is constructed.

It's a bit of a hack so I'm not entirely sure about it. This brings me to a bit of an odd question to @ChrisRackauckas: how bad is it to specifiy an ODESystem where the states don't actually depend on the independent variable? 😅
Because of the way equations are derived in Qumulants, the states are guaranteed to have the same time-parameter and there is always a derivative on the lhs, so I just never bothered putting them there. Additionally, in the end my states are already symbolic functions, e.g. f(x) or f(x*y), so I would have to construct something awkward like f(x)(t) here (not sure what the best way to that is).

They way I did it in #12 everything works, but I had to add some overloads I'm not really comfortable with, see e.g. https://github.com/david-pl/Qumulants.jl/blob/b8df71f39b3b86f25decf04ad09e2822d6fcd7c5/src/diffeq.jl#L2

from quantumcumulants.jl.

Interesting... I guess that works. That's probably not safe though.

from quantumcumulants.jl.

That's probably not safe though.

That's what I figured. I changed it so that a mapping which assigns variables of the form x(t) to my variables is stored along with the equations. When constructing an ODESystem those variables are used. See #14

from quantumcumulants.jl.

Okay everything seems good from my side now. In #15 I revamped the way noncommutative algebra is used: the only "simplifcations" we need there is the application of commutation relations. These are now applied directly whenever q-numbers are multiplied together, which ensures correctness of the resulting ODEs and is also much faster than using term rewriting for that part. All remaining simplification is standard simplification of number terms and is done with Symbolics.simplify. Code generation is done by first converting to an ODESystem so we also have all the nice features from MTK for that.

Just to summarize, here's how one would use the package now:

• Specify your model using q-numbers.
• Derive a bunch of differential equations. These are already c-number equations and are basically just stored as Vector{Equation} along with some extra stuff. So it's easy to work on them with Symbolics.
• You can use functions from Qumulants.jl such as cumulant_expansion to ensure your system of ODEs is complete.
• Finally, to solve numerically you convert to an ODESystem and do all the rest in MTK/OrdinaryDiffEq.

There's really not much of a point to have any additional symbolic functionality work on the noncommutative part because the way to go is always to convert to c-number equations as soon as possible. Still if that need should arise it's very easy to interface with SymbolicUtils (the interface is basically already there but unused) and deal with the rest via Symbolics. So in terms of features and performance there's not much more we can use from Symbolics/MTK. We can still reuse some more code from those libraries though, but I'll open more specific issues for those things.

If you think of anything else, let me know! Thanks for the help so far.

from quantumcumulants.jl.

That's super cool!

from quantumcumulants.jl.