From d8aaa6d5aa2db2efa7ccc8cf1b6d91feece8a413 Mon Sep 17 00:00:00 2001 From: Phillip Weinberg Date: Mon, 8 Jan 2024 11:58:42 -0500 Subject: [PATCH] adding interface for DormandPrince.jl as subpackage (#642) * adding interface for DormandPrince.jl as subpackage * adding to CI. * fixing julia compat * adding `SchrodingerProblem` equivilant. * adding schrodinger problem. * adding forward diff test. * fixing helper. * removing uneccesary type. * fixing test. --- .github/workflows/CI.yml | 1 + lib/BloqadeDormandPrince/Project.toml | 30 ++++++++ lib/BloqadeDormandPrince/README.md | 25 +++++++ .../src/BloqadeDormandPrince.jl | 16 +++++ lib/BloqadeDormandPrince/src/impl.jl | 15 ++++ lib/BloqadeDormandPrince/src/types.jl | 53 ++++++++++++++ lib/BloqadeDormandPrince/test/forwarddiff.jl | 29 ++++++++ lib/BloqadeDormandPrince/test/runtests.jl | 69 +++++++++++++++++++ 8 files changed, 238 insertions(+) create mode 100644 lib/BloqadeDormandPrince/Project.toml create mode 100644 lib/BloqadeDormandPrince/README.md create mode 100644 lib/BloqadeDormandPrince/src/BloqadeDormandPrince.jl create mode 100644 lib/BloqadeDormandPrince/src/impl.jl create mode 100644 lib/BloqadeDormandPrince/src/types.jl create mode 100644 lib/BloqadeDormandPrince/test/forwarddiff.jl create mode 100644 lib/BloqadeDormandPrince/test/runtests.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 3675b22de..03255d829 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -45,6 +45,7 @@ jobs: - '.' - 'lib/BloqadeGates' - 'lib/BloqadeExpr' + - 'lib/BloqadeDormandPrince' - 'lib/BloqadeKrylov' - 'lib/BloqadeLattices' - 'lib/BloqadeMIS' diff --git a/lib/BloqadeDormandPrince/Project.toml b/lib/BloqadeDormandPrince/Project.toml new file mode 100644 index 000000000..cfca02b1d --- /dev/null +++ b/lib/BloqadeDormandPrince/Project.toml @@ -0,0 +1,30 @@ +name = "BloqadeDormandPrince" +uuid = "19133132-2f02-44a9-8d8c-5915d2e7e792" +authors = ["Phillip Weinberg"] +version = "0.1.0" + +[deps] +BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2" +DormandPrince = "5e45e72d-22b8-4dd0-9c8b-f96714864bcd" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51" +YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2" + +[compat] +BloqadeExpr = "0.2.1" +DormandPrince = "0.5.1" +Reexport = "1.2" +YaoArrayRegister = "0.9" +YaoSubspaceArrayReg = "0.2" +julia = "1.6" + +[extras] +BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7" +FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "BloqadeWaveforms", "ForwardDiff", "FiniteDifferences", "Random"] diff --git a/lib/BloqadeDormandPrince/README.md b/lib/BloqadeDormandPrince/README.md new file mode 100644 index 000000000..1585ae315 --- /dev/null +++ b/lib/BloqadeDormandPrince/README.md @@ -0,0 +1,25 @@ +# BloqadeDormandPrince + +Bloqade interface for DormandPrince.jl + +## Installation + +

+BloqadeDormandPrince is a   + + + Julia Language + +   package. To install BloqadeDormandPrince, + please open + Julia's interactive session (known as REPL) and press ] + key in the REPL to use the package mode, then type the following command +

+ +```julia +pkg> add BloqadeDormandPrince +``` + +## License + +Apache 2.0 License diff --git a/lib/BloqadeDormandPrince/src/BloqadeDormandPrince.jl b/lib/BloqadeDormandPrince/src/BloqadeDormandPrince.jl new file mode 100644 index 000000000..c43254fee --- /dev/null +++ b/lib/BloqadeDormandPrince/src/BloqadeDormandPrince.jl @@ -0,0 +1,16 @@ +module BloqadeDormandPrince + +using BloqadeExpr: BloqadeExpr, Hamiltonian, nqudits +using YaoArrayRegister: AbstractRegister, ArrayReg, statevec +using YaoSubspaceArrayReg: YaoSubspaceArrayReg +using Reexport: @reexport +@reexport using DormandPrince +using DormandPrince: integrate_core! +using LinearAlgebra: mul! + +export BloqadeDPSolver, register + +include("types.jl") +include("impl.jl") + +end # BloqadeDormandPrince diff --git a/lib/BloqadeDormandPrince/src/impl.jl b/lib/BloqadeDormandPrince/src/impl.jl new file mode 100644 index 000000000..c6592ceb4 --- /dev/null +++ b/lib/BloqadeDormandPrince/src/impl.jl @@ -0,0 +1,15 @@ + +# implement DormandPrince API for BloqadeSolver +DormandPrince.get_current_state(solver::BloqadeDPSolver) = solver.reg +DormandPrince.integrate_core!(solver::BloqadeDPSolver, time::Real) = integrate_core!(solver.dp_solver, time) + +function register(solver::BloqadeDPSolver) + return get_current_state(solver) +end + + +# implement BloqadeExpr API for BloqadeSolver +function BloqadeExpr.emulate!(prob::BloqadeDPSolver) + integrate!(prob.dp_solver, prob.tend) + return prob +end diff --git a/lib/BloqadeDormandPrince/src/types.jl b/lib/BloqadeDormandPrince/src/types.jl new file mode 100644 index 000000000..771bbe2df --- /dev/null +++ b/lib/BloqadeDormandPrince/src/types.jl @@ -0,0 +1,53 @@ + + + +struct BloqadeDPSolver{Reg <: AbstractRegister, T, StateType, F, DPSolverType <: AbstractDPSolver{T, StateType, F}} <: AbstractDPSolver{T, StateType, F} + tend::T + reg::Reg + dp_solver::DPSolverType + function BloqadeDPSolver( + tend::T, + reg, + dp_solver::AbstractDPSolver{T, StateType, F} + ) where {T, StateType, F} + new{typeof(reg), T, StateType, F, typeof(dp_solver)}(tend, reg, dp_solver) + end +end + +promote_tspan(tspan) = (zero(tspan), tspan) +promote_tspan(tspan::Tuple{A, B}) where {A, B} = promote_type(A, B).(tspan) + + +function BloqadeDPSolver(reg::AbstractRegister, tspan, expr; solver_type=DP8Solver, copy_init=true, kw...) + nqudits(reg) == nqudits(expr) || throw(ArgumentError("number of qubits/sites does not match!")) + reg = copy_init ? copy(reg) : reg + tspan = promote_tspan(tspan) + + state = statevec(reg) + space = YaoSubspaceArrayReg.space(reg) + + T = real(eltype(state)) + T = isreal(expr) ? T : Complex{T} + eq = SchrodingerEquation(expr, Hamiltonian(T, expr, space)) + + tspan_type = promote_type(real(eltype(state)), eltype(tspan)) + tstart, tend = tspan_type.(tspan) # promote tspan to T so Dual number works + + solver = solver_type(eq, tstart, state; kw...) + BloqadeDPSolver(tend, reg, solver) +end + + +struct SchrodingerEquation{ExprType,H<:Hamiltonian} + expr::ExprType + hamiltonian::H +end + + +function (eq::SchrodingerEquation)(t::Real, state, dstate) + fill!(dstate, zero(eltype(dstate))) + for (f, term) in zip(eq.hamiltonian.fs, eq.hamiltonian.ts) + mul!(dstate, term, state, -im * f(t), one(t)) + end + return +end diff --git a/lib/BloqadeDormandPrince/test/forwarddiff.jl b/lib/BloqadeDormandPrince/test/forwarddiff.jl new file mode 100644 index 000000000..2e9d73876 --- /dev/null +++ b/lib/BloqadeDormandPrince/test/forwarddiff.jl @@ -0,0 +1,29 @@ +using Test +using Random +using ForwardDiff +using BloqadeDormandPrince +using YaoSubspaceArrayReg +using FiniteDifferences +using YaoArrayRegister +using BloqadeMIS +using BloqadeWaveforms +using BloqadeExpr + +function loss(xs::Vector) + reg = zero_state(Complex{eltype(xs)}, 5) + tspan = (0, 0.1) + atoms = [(i,) for i in 1:5] + Ω = piecewise_constant(clocks = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5], values = xs) + h = rydberg_h(atoms; C = 2π * 109.16, Ω, ϕ = 0.1) + prob = SchrodingerProblem(reg, tspan, h) + emulate!(prob) + return fidelity(reg, product_state(bit"00001")) +end + +Random.seed!(42) +xs = rand(5) + +Δ_ad = ForwardDiff.gradient(loss, xs) +Δ_fd, = FiniteDifferences.grad(central_fdm(5, 1), loss, xs) + +@test Δ_ad ≈ Δ_fd rtol = 1e-2 diff --git a/lib/BloqadeDormandPrince/test/runtests.jl b/lib/BloqadeDormandPrince/test/runtests.jl new file mode 100644 index 000000000..c703b9dc1 --- /dev/null +++ b/lib/BloqadeDormandPrince/test/runtests.jl @@ -0,0 +1,69 @@ +using Test +using BloqadeDormandPrince +using BloqadeExpr +using BloqadeWaveforms +using YaoArrayRegister +using LinearAlgebra + + + + +@testset "Multiple End Time Interface" begin + + times = [0.1, 0.2, 0.3, 0.4] + + # exact solution + function evolution_operator(t::Float64) + ϕ = 2.2 * sin(π * t)^2 + U = zeros(ComplexF64, 2,2) + U[1,1] = 1 / sqrt(2) + U[2,1] = 1 / sqrt(2) + U[2,2] = 1 / sqrt(2) + U[1,2] = -1 / sqrt(2) + + U * diagm(exp.([-im*ϕ, im*ϕ])) * U' + end + + function solution(t) + U = evolution_operator(t) + return U * [1.0, 0.0] + end + + exact_vals = [] + + for time in times + push!(exact_vals, solution(time)) + end + + # equivalent Hamiltonian to exact solution + atoms = [(0,0)] + wf = Waveform(t->2.2*2π*sin(2π*t), duration = 1.3); + h = rydberg_h(atoms; Ω = wf) + + # Bloqade/Dormand Prince multiple end times interface + dp_vals = [] + integrate!(BloqadeDPSolver(zero_state(1), (0, 4), h), times) do t, state + push!(dp_vals, statevec(copy(state)))# register instance returned, get underlying vector out + end + + @test dp_vals ≈ exact_vals + +end + +# test the BloqadeDPSolver no longer holds a copy of the register, it just points +# directly to the inputted register which means it mutates the register +@testset "In-Place Register" begin + atoms = [(i,) for i in 1:5] + h = rydberg_h(atoms; C = 2π * 109.16, Ω = sin, ϕ = cos) + + # DormandPrince.jl + dp_reg = zero_state(5) + bs = BloqadeDPSolver(dp_reg, (0, 10), h; copy_init=false) + integrate!(bs, 1.0) + + # strong equality + @test dp_reg === get_current_state(bs) + @test dp_reg === register(bs) + @test emulate!(bs) === bs +end +