jfeist / QuantumAlgebra.jl

Quantum operator algebra in Julia
MIT License
60 stars 10 forks source link

^(-1) and / operation on param() #15

Closed kfkq closed 1 year ago

kfkq commented 1 year ago

is there any reason why arithmetic operation ^(-1) and / is not implemented for param? i want to try to implement 1 / U expression, but i cannot do this since 1/ param(:U,'r'), param(:U,'r')^(-1) is not implemented yet.

jfeist commented 1 year ago

The main reason is simply that no one has needed it before. It would be nice to have for sure. However, implementing it would be quite a bit of work, as the code does not support arbitrary exponents (the exponents shown in the output are just "sugar" for expressions where a parameter or operator is repeated several times). For now, one possible workaround would be to define a parameter Pr"Uinv" that represents the inverse of U, and to clean up the final expressions "manually". Depending on what you use the package for and how long the expressions are, you might just be able to do it by hand, or it would be possible to write a short function that removes any occurrences of U * Uinv in the resulting expressions (to do so, you'd have to use the internals of the representation, there is no public API for this, unfortunately). I have done similar things, I could look for a short snippet that should do this.

kfkq commented 1 year ago

yes short code snippet that automatically removes "U" and "Uinv" will be very helpful to me. I am trying to do canonical transfomation, identifying each coefficient with "U" "Uinv" manually will be very tedious for hundreds terms.

jfeist commented 1 year ago

Hi, sorry for the long delay, didn't have time to get around to this before. Here's a short function that should do what you want:

function filter_xxinv(A::QuExpr,x::QuExpr,xinv::QuExpr)
    # extract the Param type from the QuExpr
    xparam = only(first(only(x.terms)).params)
    xinvparam = only(first(only(xinv.terms)).params)
    x == param(xparam.name,'r') || error("input expr x must be a single real parameter, got x = $x")
    xinv == param(xinvparam.name,'r') || error("input expr xinv must be a single real parameter, got xinv = $xinv")

    # make a new QuantumAlgebra expression
    An = QuExpr()
    # loop over all terms in the previous expression
    for (t,s) in A.terms
        npos = count(==(xparam),t.params)
        ninv = count(==(xinvparam),t.params)
        if npos==ninv==0
            QuantumAlgebra._add_sum_term!(An,t,s)
        else
            newparams = filter(p->(p != xparam)&&(p != xinvparam), t.params)
            for i = 1:(npos-ninv) # does nothing if npos <= nneg
                push!(newparams,xparam)
            end
            for i = 1:(ninv-npos) # does nothing if nneg <= npos
                push!(newparams,xinvparam)
            end
            tn = QuantumAlgebra.QuTerm(t.nsuminds,t.δs,sort!(newparams),t.expvals,t.corrs,t.bares)
            QuantumAlgebra._add_sum_term!(An,tn,s)
        end
    end
    An
end

I tested it with

U = Pr"U"
Uinv = Pr"Uinv"
ex = U * Pr"m" * a()' + Uinv * a() + a()'a()
ex3 = ex * ex * ex
display(ex3)
display(filter_xxinv(ex3,U,Uinv))

and it seems to work as expected.

kfkq commented 1 year ago

thank you so much, it works. i close this issue.