From fbef7c8aed214fff950e913765471f7f2f4c65d2 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 19 Jan 2025 20:23:15 -0500 Subject: [PATCH 1/2] Initial implementation --- README.md | 35 ++- examples/README.jl | 34 ++- src/LazyApply/LazyApply.jl | 393 +++++++++++++++++++++++++++++++ src/QuantumOperatorAlgebra.jl | 7 +- src/op.jl | 333 ++++++++++++++++++++++++++ src/ops_itensor.jl | 73 ++++++ src/trotter.jl | 35 +++ test/{basics => }/test_basics.jl | 0 8 files changed, 904 insertions(+), 6 deletions(-) create mode 100644 src/LazyApply/LazyApply.jl create mode 100644 src/op.jl create mode 100644 src/ops_itensor.jl create mode 100644 src/trotter.jl rename test/{basics => }/test_basics.jl (100%) diff --git a/README.md b/README.md index 2aaefc4..87b8e4f 100644 --- a/README.md +++ b/README.md @@ -32,10 +32,39 @@ julia> Pkg.add("QuantumOperatorAlgebra") ## Examples ````julia -using QuantumOperatorAlgebra: QuantumOperatorAlgebra -```` +using QuantumOperatorAlgebra: Op, Prod, Scaled, Sum, coefficient, sites, terms, which_op +using Test: @test + +o1 = Op("X", 1) +o2 = Op("Y", 2) + +@test which_op(o1) == "X" +@test sites(o1) == (1,) + +o = o1 + o2 + +@test o isa Sum{Op} +@test terms(o)[1] == o1 +@test terms(o)[2] == o2 -Examples go here. +o *= 2 + +@test o isa Sum{Scaled{Int,Op}} +@test terms(o)[1] == 2 * o1 +@test terms(o)[2] == 2 * o2 +@test coefficient(terms(o)[1]) == 2 +@test coefficient(terms(o)[2]) == 2 + +o3 = Op("Z", 3) + +o *= o3 + +@test o isa Sum{Scaled{Int,Prod{Op}}} +@test terms(o)[1] == 2 * o1 * o3 +@test terms(o)[2] == 2 * o2 * o3 +@test coefficient(terms(o)[1]) == 2 +@test coefficient(terms(o)[2]) == 2 +```` --- diff --git a/examples/README.jl b/examples/README.jl index ef44cf2..05defec 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -37,5 +37,35 @@ julia> Pkg.add("QuantumOperatorAlgebra") # ## Examples -using QuantumOperatorAlgebra: QuantumOperatorAlgebra -# Examples go here. +using QuantumOperatorAlgebra: Op, Prod, Scaled, Sum, coefficient, sites, terms, which_op +using Test: @test + +o1 = Op("X", 1) +o2 = Op("Y", 2) + +@test which_op(o1) == "X" +@test sites(o1) == (1,) + +o = o1 + o2 + +@test o isa Sum{Op} +@test terms(o)[1] == o1 +@test terms(o)[2] == o2 + +o *= 2 + +@test o isa Sum{Scaled{Int,Op}} +@test terms(o)[1] == 2 * o1 +@test terms(o)[2] == 2 * o2 +@test coefficient(terms(o)[1]) == 2 +@test coefficient(terms(o)[2]) == 2 + +o3 = Op("Z", 3) + +o *= o3 + +@test o isa Sum{Scaled{Int,Prod{Op}}} +@test terms(o)[1] == 2 * o1 * o3 +@test terms(o)[2] == 2 * o2 * o3 +@test coefficient(terms(o)[1]) == 2 +@test coefficient(terms(o)[2]) == 2 diff --git a/src/LazyApply/LazyApply.jl b/src/LazyApply/LazyApply.jl new file mode 100644 index 0000000..58b5274 --- /dev/null +++ b/src/LazyApply/LazyApply.jl @@ -0,0 +1,393 @@ +module LazyApply + +import Base: + ==, + +, + -, + *, + /, + ^, + exp, + adjoint, + copy, + show, + getindex, + length, + isless, + iterate, + firstindex, + lastindex, + keys, + reverse, + size + +export Applied, Scaled, Sum, Prod, Exp, coefficient, argument, expand, materialize, terms + +struct Applied{F,Args<:Tuple,Kwargs<:NamedTuple} + f::F + args::Args + kwargs::Kwargs +end +Applied(f, args::Tuple) = Applied(f, args, (;)) + +materialize(x) = x +function materialize(a::Applied) + return a.f(materialize.(a.args)...; a.kwargs...) +end + +function (a1::Applied == a2::Applied) + return a1.f == a2.f && a1.args == a2.args && a1.kwargs == a2.kwargs +end + +# +# Applied algebra +# + +# Used for dispatch +const Scaled{C<:Number,A} = Applied{typeof(*),Tuple{C,A},NamedTuple{(),Tuple{}}} +const Sum{A} = Applied{typeof(sum),Tuple{Vector{A}},NamedTuple{(),Tuple{}}} +const Prod{A} = Applied{typeof(prod),Tuple{Vector{A}},NamedTuple{(),Tuple{}}} + +# Some convenient empty constructors +Sum{A}() where {A} = Applied(sum, (A[],)) +Prod{A}() where {A} = Applied(prod, (A[],)) + +coefficient(co::Scaled{C}) where {C} = co.args[1] +argument(co::Scaled{C}) where {C} = co.args[2] + +# +# Generic algebra +# + +# 1.3 * Op("X", 1) + 1.3 * Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + 1.3 * Op("X", 3) * Op("X", 4) +function (a1::Scaled{C,A} + a2::Scaled{C,A}) where {C,A} + return Sum{Scaled{C,A}}() + a1 + a2 +end + +function (a1::Prod{A} + a2::Prod{A}) where {A} + return Sum{Prod{A}}() + a1 + a2 +end + +(c::Number * a::Scaled{C}) where {C} = (c * coefficient(a)) * argument(a) +(a::Scaled{C} * c::Number) where {C} = (coefficient(a) * c) * argument(a) + +-(a::Scaled{C}) where {C} = (-one(C) * a) +-(a::Sum) = (-1 * a) +-(a::Prod) = (-1 * a) + +(os::Sum{A} + o::A) where {A} = Applied(sum, (vcat(os.args[1], [o]),)) +(o::A + os::Sum{A}) where {A} = Applied(sum, (vcat([o], os.args[1]),)) + +(a1::Sum{A} - a2::A) where {A} = a1 + (-a2) +(a1::A - a2::Sum{A}) where {A} = a1 + (-a2) + +(a1::Sum{A} - a2::Prod{A}) where {A} = a1 + (-a2) +(a1::Sum{A} - a2::Scaled{C,Prod{A}}) where {C,A} = a1 + (-a2) +(a1::Sum{A} - a2::Sum{Scaled{C,Prod{A}}}) where {C,A} = a1 + (-a2) + +(a1::Prod{A} * a2::A) where {A} = Applied(prod, (vcat(only(a1.args), [a2]),)) +(a1::A * a2::Prod{A}) where {A} = Applied(prod, (vcat([a1], only(a2.args)),)) + +# Fixes ambiguity error with: +# *(a1::Applied, a2::Sum) +# *(os::Prod{A}, o::A) +(a1::Prod{Sum{A}} * a2::Sum{A}) where {A} = Applied(prod, (vcat(only(a1.args), [a2]),)) + +# 1.3 * Op("X", 1) + 1 * Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + 1 * Op("X", 3) +# 1.3 * Op("X", 1) * Op("X", 2) + 1 * Op("X", 3) * Op("X", 4) +function (co1::Scaled{C1,A} + co2::Scaled{C2,A}) where {C1,C2,A} + c1, c2 = promote(coefficient(co1), coefficient(co2)) + return c1 * argument(co1) + c2 * argument(co2) +end + +# (1.3 * Op("X", 1)) * (1.3 * Op("X", 2)) +function (co1::Scaled{C1} * co2::Scaled{C2}) where {C1,C2} + c = coefficient(co1) * coefficient(co2) + o = argument(co1) * argument(co2) + return c * o +end + +function (a1::Prod{A} * a2::Scaled{C,A}) where {C,A} + return coefficient(a2) * (a1 * argument(a2)) +end + +function (a1::Prod{A} + a2::Scaled{C,A}) where {C,A} + return one(C) * a1 + Prod{A}() * a2 +end + +# (Op("X", 1) + Op("X", 2)) + (Op("X", 3) + Op("X", 4)) +# (Op("X", 1) * Op("X", 2) + Op("X", 3) * Op("X", 4)) + (Op("X", 5) * Op("X", 6) + Op("X", 7) * Op("X", 8)) +(a1::Sum{A} + a2::Sum{A}) where {A} = Applied(sum, (vcat(a1.args[1], a2.args[1]),)) +(a1::Sum{A} - a2::Sum{A}) where {A} = a1 + (-a2) + +(a1::Prod{A} * a2::Prod{A}) where {A} = Applied(prod, (vcat(only(a1.args), only(a2.args)),)) + +(os::Sum{Scaled{C,A}} + o::A) where {C,A} = os + one(C) * o +(o::A + os::Sum{Scaled{C,A}}) where {C,A} = one(C) * o + os + +# Op("X", 1) + Op("X", 2) + 1.3 * Op("X", 3) +(os::Sum{A} + co::Scaled{C,A}) where {C,A} = one(C) * os + co + +# 1.3 * Op("X", 1) + (Op("X", 2) + Op("X", 3)) +(co::Scaled{C,A} + os::Sum{A}) where {C,A} = co + one(C) * os + +# 1.3 * (Op("X", 1) + Op("X", 2)) +(c::Number * os::Sum) = Applied(sum, (c * os.args[1],)) + +(a1::Applied * a2::Sum) = Applied(sum, (map(a -> a1 * a, only(a2.args)),)) +(a1::Sum * a2::Applied) = Applied(sum, (map(a -> a * a2, only(a1.args)),)) +(a1::Sum * a2::Sum) = Applied(prod, ([a1, a2],)) + +function _expand(a1::Sum, a2::Sum) + return Applied(sum, (vec([a1[i] * a2[j] for i in 1:length(a1), j in 1:length(a2)]),)) +end + +function expand(a::Prod) + if length(a) == 1 + return a[1] + elseif length(a) ≥ 2 + a12 = _expand(a[1], a[2]) + return expand(Applied(prod, (vcat([a12], a[3:end]),))) + end +end + +# (Op("X", 1) + Op("X", 2)) * 1.3 +(os::Sum * c::Number) = c * os + +# (Op("X", 1) + Op("X", 2)) / 1.3 +(os::Sum / c::Number) = inv(c) * os + +# Promotions +function (co1::Scaled{C,Prod{A}} + co2::Scaled{C,A}) where {C,A} + return co1 + coefficient(co2) * Applied(prod, ([argument(co2)],)) +end + +function (a1::Scaled - a2::Scaled) + return a1 + (-a2) +end + +function (a1::Prod{A} + a2::A) where {A} + return a1 + Applied(prod, ([a2],)) +end + +function (a1::Sum{A} + a2::Prod{A}) where {A} + return Prod{A}() * a1 + a2 +end + +function (a1::Sum{A} + a2::Sum{Scaled{C,Prod{A}}}) where {C,A} + return (one(C) * Prod{A}() * a1) + a2 +end + +function (a1::Prod{A} - a2::A) where {A} + return a1 + (-a2) +end + +function (co1::Sum{Scaled{C,Prod{A}}} + co2::Scaled{C,A}) where {C,A} + return co1 + coefficient(co2) * Applied(prod, ([argument(co2)],)) +end + +function (a1::Sum{Scaled{C1,Prod{A}}} - a2::Scaled{C2,A}) where {C1,C2,A} + return a1 + (-a2) +end + +function (a1::Sum{Scaled{C,Prod{A}}} - a2::Prod{A}) where {C,A} + return a1 + (-a2) +end + +function (a1::Sum{Scaled{C1,Prod{A}}} - a2::Scaled{C2,Prod{A}}) where {C1,C2,A} + return a1 + (-a2) +end + +function (a1::Sum{A} + a2::Scaled{C,Prod{A}}) where {C,A} + return Sum{Scaled{C,Prod{A}}}() + a1 + a2 +end + +function (a1::Sum{Scaled{C1,Prod{A}}} + a2::Scaled{C2,A}) where {C1,C2,A} + C = promote_type(C1, C2) + return one(C) * a1 + one(C) * a2 +end + +# (::Sum{Scaled{Bool,Prod{Op}}} + ::Scaled{Float64,Prod{Op}}) +function (a1::Sum{Scaled{C1,A}} + a2::Scaled{C2,A}) where {C1,C2,A} + C = promote_type(C1, C2) + return one(C) * a1 + one(C) * a2 +end + +# TODO: Is this needed? It seems like: +# +# (a1::Sum{A} + a2::A) +# +# is not being called. +function (a1::Sum{Scaled{C,A}} + a2::Scaled{C,A}) where {C,A} + return Applied(sum, (vcat(only(a1.args), [a2]),)) +end + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::Sum{A}) where {C,A} + a2 = one(C) * a2 + a2 = Prod{A}() * a2 + return a1 + one(C) * Prod{A}() * a2 +end + +function (a1::Sum{Prod{A}} + a2::A) where {A} + return a1 + (Prod{A}() * a2) +end + +function (a1::Sum{Prod{A}} + a2::Scaled{C,A}) where {C,A} + return a1 + (Prod{A}() * a2) +end + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::A) where {C,A} + return a1 + one(C) * a2 +end +(a1::Sum{Scaled{C,Prod{A}}} - a2::A) where {C,A} = a1 + (-a2) + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::Sum{Scaled{C,A}}) where {C,A} + return a1 + (Prod{A}() * a2) +end + +function (o::A + os::Sum{Scaled{C,Prod{A}}}) where {C,A} + return one(C) * o + os +end + +function (a::Sum^n::Int) + r = a + for _ in 2:n + r *= a + end + return r +end + +function (a::Prod^n::Int) + r = a + for _ in 2:n + r *= a + end + return r +end + +exp(a::Applied) = Applied(exp, (a,)) + +const Exp{A} = Applied{typeof(exp),Tuple{A},NamedTuple{(),Tuple{}}} +const Adjoint{A} = Applied{typeof(adjoint),Tuple{A},NamedTuple{(),Tuple{}}} + +argument(a::Exp) = a.args[1] + +(c::Number * e::Exp) = Applied(*, (c, e)) +(e::Exp * c::Number) = c * e +(e1::Exp * e2::Exp) = Applied(prod, ([e1, e2],)) +(e1::Applied * e2::Exp) = Applied(prod, ([e1, e2],)) +(e1::Exp * e2::Applied) = Applied(prod, ([e1, e2],)) + +# Fix ambiguity error with `Applied * Sum`. +# TODO: Share code, define `mul_sum`. +(a1::Exp * a2::Sum) = Applied(sum, (map(a -> a1 * a, only(a2.args)),)) +# Fix ambiguity error with `Sum * Applied`. +# TODO: Share code, define `mul_sum`. +(a1::Sum * a2::Exp) = Applied(sum, (map(a -> a * a2, only(a1.args)),)) + +function reverse(a::Prod) + return Applied(prod, (reverse(only(a.args)),)) +end + +adjoint(a::Prod) = Applied(prod, (map(adjoint, reverse(only(a.args))),)) + +# +# Convenient indexing +# + +getindex(a::Union{Sum,Prod}, I...) = only(a.args)[I...] +iterate(a::Union{Sum,Prod}, args...) = iterate(only(a.args), args...) +size(a::Union{Sum,Prod}) = size(only(a.args)) +length(a::Union{Sum,Prod}) = length(only(a.args)) +firstindex(a::Union{Sum,Prod}) = 1 +lastindex(a::Union{Sum,Prod}) = length(a) +keys(a::Union{Sum,Prod}) = 1:length(a) + +length(a::Scaled{C,<:Sum}) where {C} = length(argument(a)) +length(a::Scaled{C,<:Prod}) where {C} = length(argument(a)) +getindex(a::Scaled{C,<:Sum}, I...) where {C} = getindex(argument(a), I...) +getindex(a::Scaled{C,<:Prod}, I...) where {C} = getindex(argument(a), I...) +lastindex(a::Scaled{C,<:Sum}) where {C} = lastindex(argument(a)) +lastindex(a::Scaled{C,<:Prod}) where {C} = lastindex(argument(a)) + +# +# Functions convenient for OpSum code +# + +terms(a::Union{Sum,Prod}) = only(a.args) +terms(a::Scaled{C,<:Union{Sum,Prod}}) where {C} = terms(argument(a)) +copy(a::Applied) = Applied(deepcopy(a.f), deepcopy(a.args), deepcopy(a.kwargs)) +Sum(a::Vector) = Applied(sum, (a,)) +Prod(a::Vector) = Applied(prod, (a,)) +function isless(a1::Applied{F}, a2::Applied{F}) where {F} + return (isless(a1.args, a2.args) && isless(a1.kwargs, a2.kwargs)) +end + +# +# Printing +# + +function show(io::IO, ::MIME"text/plain", a::Sum) + print(io, "sum(\n") + for n in eachindex(a) + print(io, " ", a[n]) + if n ≠ lastindex(a) + print(io, "\n") + end + end + print(io, "\n)") + return nothing +end +show(io::IO, a::Sum) = show(io, MIME("text/plain"), a) + +function show(io::IO, ::MIME"text/plain", a::Prod) + print(io, "prod(\n") + for n in eachindex(a) + print(io, " ", a[n]) + if n ≠ lastindex(a) + print(io, "\n") + end + end + print(io, "\n)") + return nothing +end +show(io::IO, a::Prod) = show(io, MIME("text/plain"), a) + +function show(io::IO, m::MIME"text/plain", a::Exp) + print(io, a.f, "(") + for n in 1:length(a.args) + print(io, a.args[n]) + if n < length(a.args) + print(io, ", ") + end + end + print(io, ")") + return nothing +end +show(io::IO, a::Exp) = show(io, MIME("text/plain"), a) + +function show(io::IO, m::MIME"text/plain", a::Applied) + print(io, a.f, "(") + for n in eachindex(a.args) + print(io, a.args[n]) + if n < length(a.args) + print(io, ", ") + end + end + if !isempty(a.kwargs) + print(io, "; ") + for n in 1:length(a.kwargs) + print(io, keys(a.kwargs)[n], "=", a.kwargs[n]) + if n < length(a.kwargs) + print(io, ", ") + end + end + end + print(io, ")") + return nothing +end +show(io::IO, a::Applied) = show(io, MIME("text/plain"), a) + +end diff --git a/src/QuantumOperatorAlgebra.jl b/src/QuantumOperatorAlgebra.jl index 3ea2aeb..be630a3 100644 --- a/src/QuantumOperatorAlgebra.jl +++ b/src/QuantumOperatorAlgebra.jl @@ -1,5 +1,10 @@ module QuantumOperatorAlgebra -# Write your package code here. +include("LazyApply/LazyApply.jl") +# Make these available as `QuantumOperatorAlgebra.f`. +using .LazyApply: coefficient, terms + +include("op.jl") +include("trotter.jl") end diff --git a/src/op.jl b/src/op.jl new file mode 100644 index 0000000..17070cf --- /dev/null +++ b/src/op.jl @@ -0,0 +1,333 @@ +using .LazyApply: LazyApply, Exp, Prod, Scaled, Sum, argument, coefficient + +##################################################################################### +# General functionality +# + +# Helper function to split a `Tuple` according to the function `f`. +# For example: +# +# julia> t = (1, "X", 1, 2, "Y", 2, "Z", 4) +# (1, "X", 1, 2, "Y", 2, "Z", 4) +# +# julia> split(x -> x isa AbstractString, t) +# [(1,), ("X", 1, 2), ("Y", 2), ("Z", 4)] +# +function split(f, t::Tuple) + n = findall(f, t) + nsplit = length(n) + 1 + s = Vector{Any}(undef, nsplit) + s[1] = t[1:(first(n) - 1)] + for i in 2:(nsplit - 1) + s[i] = t[n[i - 1]:(n[i] - 1)] + end + s[end] = t[last(n):end] + return s +end + +## XXX: Very long compile times: +## https://github.com/JuliaLang/julia/issues/45545 +## +## julia> using ITensors +## +## julia> @time ITensors.Ops.split(x -> x isa String, ("X", 1)) +## 7.588123 seconds (2.34 M allocations: 100.919 MiB, 1.71% gc time, 100.00% compilation time) +## ((), ("X", 1)) +## +## julia> @time ITensors.Ops.split(x -> x isa String, ("X", 1)) +## 0.042590 seconds (88.59 k allocations: 4.823 MiB, 19.13% gc time, 99.84% compilation time) +## ((), ("X", 1)) +## +## function split(f, t::Tuple) +## n = findall(f, t) +## ti = t[1:(first(n) - 1)] +## ts = ntuple(i -> t[n[i]:(n[i + 1] - 1)], length(n) - 1) +## tf = t[last(n):end] +## return ti, ts..., tf +## end + +struct Op + which_op + sites::Tuple + params::NamedTuple + function Op(which_op, site...; kwargs...) + return new(which_op, site, NamedTuple(kwargs)) + end +end + +which_op(o::Op) = o.which_op +name(o::Op) = which_op(o) +sites(o::Op) = o.sites +site(o::Op) = only(sites(o)) +params(o::Op) = o.params + +function Base.:(==)(o1::Op, o2::Op) + return o1.which_op == o2.which_op && o1.sites == o2.sites && o1.params == o2.params +end + +function Base.hash(o::Op, h::UInt) + return hash(which_op(o), hash(sites(o), hash(params(o), hash(:Op, h)))) +end + +# Version of `isless` defined for matrices +_isless(a, b) = isless(a, b) +_isless(a::AbstractMatrix, b::AbstractMatrix) = isless(hash(a), hash(b)) +_isless(a::AbstractString, b::AbstractMatrix) = true +_isless(a::AbstractMatrix, b::AbstractString) = !_isless(b, a) + +function Base.isless(o1::Op, o2::Op) + if sites(o1) ≠ sites(o2) + return sites(o1) < sites(o2) + end + if which_op(o1) ≠ which_op(o2) + return _isless(which_op(o1), which_op(o2)) + end + return params(o1) < params(o2) +end + +function Base.isless(o1::Prod{Op}, o2::Prod{Op}) + if length(o1) ≠ length(o2) + return length(o1) < length(o2) + end + for n in 1:length(o1) + if o1[n] ≠ o2[n] + return (o1[n] < o2[n]) + end + end + return false +end + +function Base.isless(o1::Scaled{C1,Prod{Op}}, o2::Scaled{C2,Prod{Op}}) where {C1,C2} + if argument(o1) == argument(o2) + if coefficient(o1) ≈ coefficient(o2) + return false + else + c1 = coefficient(o1) + c2 = coefficient(o2) + #"lexicographic" ordering on complex numbers + return real(c1) < real(c2) || (real(c1) ≈ real(c2) && imag(c1) < imag(c2)) + end + end + return argument(o1) < argument(o2) +end + +## function Op(t::Tuple) +## which_op = first(t) +## site_params = Base.tail(t) +## if last(site_params) isa NamedTuple +## site = Base.front(site_params) +## params = last(site_params) +## else +## site = site_params +## params = (;) +## end +## return Op(which_op, site; params...) +## end + +## function Op(t::Tuple{WhichOp,NamedTuple,Vararg}) where {WhichOp} +## params = t[2] +## which_op = t[1] +## sites = t[3:end] +## return Op(which_op, sites...; params...) +## end + +function sites(a::Union{Sum,Prod}) + s = [] + for n in 1:length(a) + s = s ∪ sites(a[n]) + end + return map(identity, s) +end +sites(a::Scaled{C,<:Sum}) where {C} = sites(argument(a)) +sites(a::Scaled{C,<:Prod}) where {C} = sites(argument(a)) + +params(a::Scaled{C,<:Prod}) where {C} = params(only(argument(a))) + +which_op(a::Scaled{C,Op}) where {C} = which_op(argument(a)) +sites(a::Scaled{C,Op}) where {C} = sites(argument(a)) +params(a::Scaled{C,Op}) where {C} = params(argument(a)) + +# +# Op algebra +# + +function Base.convert(::Type{Scaled{C1,Prod{Op}}}, o::Scaled{C2,Prod{Op}}) where {C1,C2} + c = convert(C1, coefficient(o)) + return c * argument(o) +end + +""" +An `OpSum` represents a sum of operator +terms. + +Often it is used to create matrix +product operator (`MPO`) approximation +of the sum of the terms in the `OpSum` oject. +Each term is a product of local operators +specified by names such as `"Sz"` or `"N"`, +times an optional coefficient which +can be real or complex. + +Which local operator names are available +is determined by the function `op` +associated with the `TagType` defined by +special Index tags, such as `"S=1/2"`, `"S=1"`, +`"Fermion"`, and `"Electron"`. +""" +const OpSum{C} = Sum{Scaled{C,Prod{Op}}} + +# This helps with in-place operations +OpSum() = OpSum{ComplexF64}() + +Base.:+(o1::Op, o2::Op) = Applied(sum, ([o1, o2],)) +Base.:*(o1::Op, o2::Op) = Applied(prod, ([o1, o2],)) +Base.:-(o::Op) = -one(Int) * o +Base.:-(o1::Op, o2::Op) = o1 + (-o2) + +Base.:*(c::Number, o::Op) = Applied(*, (c, o)) +Base.:*(o::Op, c::Number) = Applied(*, (c, o)) +Base.:/(o::Op, c::Number) = Applied(*, (inv(c), o)) + +Base.:*(c::Number, o::Prod{Op}) = Applied(*, (c, o)) +Base.:*(o::Prod{Op}, c::Number) = Applied(*, (c, o)) +Base.:/(o::Prod{Op}, c::Number) = Applied(*, (inv(c), o)) + +# 1.3 * Op("X", 1) + Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + Op("X", 3) +Base.:+(co1::Scaled{C}, o2::Op) where {C} = co1 + one(C) * o2 + +# Op("X", 1) + 1.3 * Op("X", 2) +Base.:+(o1::Op, co2::Scaled{C}) where {C} = one(C) * o1 + co2 + +Base.:*(o1::Op, o2::Sum) = Applied(sum, (map(a -> o1 * a, only(o2.args)),)) +Base.:*(o1::Sum, o2::Op) = Applied(sum, (map(a -> a * o2, only(o1.args)),)) + +# 1.3 * Op("X", 1) + Op("X", 2) * Op("X", 3) +# 1.3 * Op("X", 1) * Op("X", 2) + Op("X", 3) * Op("X", 4) +Base.:+(co1::Scaled{C}, o2::Prod{Op}) where {C} = co1 + one(C) * o2 + +# 1.3 * Op("X", 1) * Op("X", 2) +Base.:*(co1::Scaled{C}, o2::Op) where {C} = co1 * (one(C) * o2) + +Base.exp(o::Op) = Applied(exp, (o,)) + +Base.adjoint(o::Op) = Applied(adjoint, (o,)) +Base.adjoint(o::LazyApply.Adjoint{Op}) = only(o.args) + +Base.:*(o1::Exp{Op}, o2::Op) = Applied(prod, ([o1, o2],)) + +# +# Tuple interface +# + +const OpSumLike{C} = Union{ + Sum{Op}, + Sum{Scaled{C,Op}}, + Sum{Prod{Op}}, + Sum{Scaled{C,Prod{Op}}}, + Prod{Op}, + Scaled{C,Prod{Op}}, +} + +const WhichOp = Union{AbstractString,AbstractMatrix{<:Number}} + +# Make a `Scaled{C,Prod{Op}}` from a `Tuple` input, +# for example: +# +# (1.2, "X", 1, "Y", 2) -> 1.2 * Op("X", 1) * Op("Y", 2) +# +function op_term(a::Tuple{Number,Vararg}) + c = first(a) + return c * op_term(Base.tail(a)) +end + +function op_site(which_op, params::NamedTuple, sites...) + return Op(which_op, sites...; params...) +end + +function op_site(which_op, sites_params...) + if last(sites_params) isa NamedTuple + sites = Base.front(sites_params) + params = last(sites_params) + return Op(which_op, sites...; params...) + end + return Op(which_op, sites_params...) +end + +function op_term(a::Tuple{Vararg}) + a_split = split(x -> x isa WhichOp, a) + @assert isempty(first(a_split)) + popfirst!(a_split) + o = op_site(first(a_split)...) + popfirst!(a_split) + for aₙ in a_split + o *= op_site(aₙ...) + end + return o +end + +function Base.:+(o1::OpSumLike, o2::Tuple) + return o1 + op_term(o2) +end + +function Base.:+(o1::Tuple, o2::OpSumLike) + return op_term(o1) + o2 +end + +function Base.:-(o1::OpSumLike, o2::Tuple) + return o1 - op_term(o2) +end + +function Base.:-(o1::Tuple, o2::OpSumLike) + return op_term(o1) - o2 +end + +function Base.:*(o1::OpSumLike, o2::Tuple) + return o1 * op_term(o2) +end + +function Base.:*(o1::Tuple, o2::OpSumLike) + return op_term(o1) * o2 +end + +function Base.show(io::IO, ::MIME"text/plain", o::Op) + print(io, which_op(o)) + print(io, sites(o)) + if !isempty(params(o)) + print(io, params(o)) + end + return nothing +end +Base.show(io::IO, o::Op) = show(io, MIME("text/plain"), o) + +function Base.show(io::IO, ::MIME"text/plain", o::Prod{Op}) + for n in 1:length(o) + print(io, o[n]) + if n < length(o) + print(io, " ") + end + end + return nothing +end +Base.show(io::IO, o::Prod{Op}) = show(io, MIME("text/plain"), o) + +function Base.show( + io::IO, m::MIME"text/plain", o::Scaled{C,O} +) where {C,O<:Union{Op,Prod{Op}}} + c = coefficient(o) + if isreal(c) + c = real(c) + end + print(io, c) + print(io, " ") + show(io, m, argument(o)) + return nothing +end +Base.show(io::IO, o::Scaled{C,Prod{Op}}) where {C} = show(io, MIME("text/plain"), o) + +function Base.show(io::IO, ::MIME"text/plain", o::LazyApply.Adjoint{Op}) + print(io, o') + print(io, "'") + return nothing +end +Base.show(io::IO, o::LazyApply.Adjoint{Op}) = show(io, MIME("text/plain"), o) diff --git a/src/ops_itensor.jl b/src/ops_itensor.jl new file mode 100644 index 0000000..e650376 --- /dev/null +++ b/src/ops_itensor.jl @@ -0,0 +1,73 @@ +using .SiteTypes: SiteTypes, op + +function SiteTypes.op(I::UniformScaling, s::Index...) + return I.λ * op("Id", s...) +end + +function ITensor(o::Op, s::Vector{<:Index}) + return op(o.which_op, map(n -> s[n], o.sites)...; o.params...) +end + +function ITensor(o::Scaled, s::Vector{<:Index}) + c = coefficient(o) + if isreal(c) + c = real(c) + end + return c * ITensor(argument(o), s) +end + +function ITensor(o::Prod, s::Vector{<:Index}) + T = ITensor(true) + for a in o.args[1] + Tₙ = ITensor(a, s) + # TODO: Implement this logic inside `apply` + if hascommoninds(T, Tₙ) + T = T(Tₙ) + else + T *= Tₙ + end + end + return T +end + +function ITensor(o::Sum, s::Vector{<:Index}) + T = ITensor() + for a in o.args[1] + T += ITensor(a, s) + end + return T +end + +function ITensor(o::Exp, s::Vector{<:Index}) + return exp(ITensor(argument(o), s)) +end + +function ITensor(o::LazyApply.Adjoint, s::Vector{<:Index}) + return swapprime(dag(ITensor(o', s)), 0 => 1) +end + +function Sum{ITensor}(o::Sum, s::Vector{<:Index}) + return Applied(sum, (map(oₙ -> ITensor(oₙ, s), only(o.args)),)) +end + +function Prod{ITensor}(o::Prod, s::Vector{<:Index}) + return Applied(prod, (map(oₙ -> ITensor(oₙ, s), only(o.args)),)) +end + +function Prod{ITensor}(o::Scaled{C,Prod{Op}}, s::Vector{<:Index}) where {C} + t = Prod{ITensor}(argument(o), s) + t1 = coefficient(o) * only(t.args)[1] + return Applied(prod, (vcat([t1], only(t.args)[2:end]),)) +end + +function apply(o::Prod{ITensor}, v::ITensor; kwargs...) + ov = v + for oₙ in reverse(only(o.args)) + ov = apply(oₙ, ov; kwargs...) + end + return ov +end + +function (o::Prod{ITensor})(v::ITensor; kwargs...) + return apply(o, v; kwargs...) +end diff --git a/src/trotter.jl b/src/trotter.jl new file mode 100644 index 0000000..9c3b1a5 --- /dev/null +++ b/src/trotter.jl @@ -0,0 +1,35 @@ +using ..LazyApply: Applied, Sum + +abstract type ExpAlgorithm end + +struct Exact <: ExpAlgorithm end + +struct Trotter{Order} <: ExpAlgorithm + nsteps::Int +end +Trotter{Order}() where {Order} = Trotter{Order}(1) +Base.one(::Trotter{Order}) where {Order} = Trotter{Order}(1) + +function Base.exp(o::Sum; alg::ExpAlgorithm=Exact()) + return exp(alg, o) +end + +function Base.exp(::Exact, o::Sum) + return Applied(prod, ([Applied(exp, (o,))],)) +end + +function exp_one_step(trotter::Trotter{1}, o::Sum) + exp_o = Applied(prod, (map(exp, reverse(only(o.args))),)) + return exp_o +end + +function exp_one_step(trotter::Trotter{2}, o::Sum) + exp_o_order_1 = exp_one_step(Trotter{1}(), o / 2) + exp_o = reverse(exp_o_order_1) * exp_o_order_1 + return exp_o +end + +function Base.exp(trotter::Trotter, o::Sum) + expδo = exp_one_step(one(trotter), o / trotter.nsteps) + return expδo^trotter.nsteps +end diff --git a/test/basics/test_basics.jl b/test/test_basics.jl similarity index 100% rename from test/basics/test_basics.jl rename to test/test_basics.jl From 1031ec01db22ece0f8c67d419b1df65169290c93 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 19 Jan 2025 20:42:41 -0500 Subject: [PATCH 2/2] Update README --- README.md | 15 +++++++++++++++ examples/README.jl | 15 +++++++++++++++ 2 files changed, 30 insertions(+) diff --git a/README.md b/README.md index 87b8e4f..3d3c6ed 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,21 @@ [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) +Quantum operator algebra system. This is mostly meant to be used as a backend in [ITensorMPS.jl](https://github.com/ITensor/ITensorMPS.jl) +and [ITensorNetworks.jl](https://github.com/ITensor/ITensorNetworks.jl) for lazily representing operator expressions +that will be turned into quantum circuits and tensor networks. + +See also: +- [ITensorQuantumOperatorDefinitions.jl](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl) for operator definitions +compatible with this system. +- [Yao.jl](https://github.com/QuantumBFS/Yao.jl) +- [Quac.jl](https://github.com/bsc-quantic/Quac.jl) +- [QuantumAlgebra.jl](https://github.com/jfeist/QuantumAlgebra.jl) +- [QuantumCumulants.jl](https://github.com/qojulia/QuantumCumulants.jl) +- [QuantumSymbolics.jl](https://github.com/QuantumSavory/QuantumSymbolics.jl) +- [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl), [QuantumInterface.jl](https://github.com/qojulia/QuantumInterface.jl) +- [QuantumLattices.QuantumOperators](https://github.com/Quantum-Many-Body/QuantumLattices.jl/blob/master/src/QuantumOperators.jl) + ## Installation instructions This package resides in the `ITensor/ITensorRegistry` local registry. diff --git a/examples/README.jl b/examples/README.jl index 05defec..da8caf9 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -6,6 +6,21 @@ # [![Coverage](https://codecov.io/gh/ITensor/QuantumOperatorAlgebra.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/QuantumOperatorAlgebra.jl) # [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) # [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) +# +# Quantum operator algebra system. This is mostly meant to be used as a backend in [ITensorMPS.jl](https://github.com/ITensor/ITensorMPS.jl) +# and [ITensorNetworks.jl](https://github.com/ITensor/ITensorNetworks.jl) for lazily representing operator expressions +# that will be turned into quantum circuits and tensor networks. +# +# See also: +# - [ITensorQuantumOperatorDefinitions.jl](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl) for operator definitions +# compatible with this system. +# - [Yao.jl](https://github.com/QuantumBFS/Yao.jl) +# - [Quac.jl](https://github.com/bsc-quantic/Quac.jl) +# - [QuantumAlgebra.jl](https://github.com/jfeist/QuantumAlgebra.jl) +# - [QuantumCumulants.jl](https://github.com/qojulia/QuantumCumulants.jl) +# - [QuantumSymbolics.jl](https://github.com/QuantumSavory/QuantumSymbolics.jl) +# - [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl), [QuantumInterface.jl](https://github.com/qojulia/QuantumInterface.jl) +# - [QuantumLattices.QuantumOperators](https://github.com/Quantum-Many-Body/QuantumLattices.jl/blob/master/src/QuantumOperators.jl) # ## Installation instructions