Inplace exponential

CHANGELOG.md

@@ -3,6 +3,8 @@ QuantumAnnealing.jl Change Log
 
 ### Staged
 - Add Gibbs distribution
+- Fix bug when using kwargs in simulate_de
+- Reduce memory usage in simulate and simulate_de
 
 ### v0.2.0
 - Add a generic Magnus expansion solver for any order

Project.toml

@@ -8,6 +8,7 @@ version = "0.2.0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
@@ -18,9 +19,9 @@ JSON = "~0.21, ~0.20, ~0.19, ~0.18"
 julia = "^1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","LinearAlgebra","Random"]
+test = ["Test", "LinearAlgebra", "Random"]

src/QuantumAnnealing.jl

@@ -1,5 +1,5 @@
 module QuantumAnnealing
-
+    import LinearAlgebra
     import CSV
     import SparseArrays
     import JSON

src/simulate.jl

@@ -144,7 +144,7 @@ function simulate_magnus_optimized(ising_model::Dict, annealing_time::Real, anne
     s_steps = range(0, 1, length=steps)
     R_current = R0
     U = foldl(kron, [_IMAT for i = 1:n])
-
+    U_next = similar(U)
     if track_states
         push!(state_steps, R_current)
     end
@@ -160,7 +160,7 @@ function simulate_magnus_optimized(ising_model::Dict, annealing_time::Real, anne
         #   display(Matrix(Ωi))
         #end
 
-        U_next = exp(Matrix(sum(Ω_list)))
+        U_next = LinearAlgebra.exp!(Matrix(sum(Ω_list)))
         U = U_next * U
 
         if track_states
@@ -471,7 +471,7 @@ function simulate_magnus_generic(ising_model::Dict, annealing_time::Real, anneal
         #   display(Matrix(Ωi))
         #end
 
-        U_next = exp(Matrix(sum(Ω_list)))
+        U_next = LinearAlgebra.exp!(Matrix(sum(Ω_list)))
         U = U_next * U
 
         if track_states

src/simulate_de.jl

@@ -41,7 +41,7 @@ function simulate_de(ising_model, annealing_time, annealing_schedule, reltol; ab
 
     s_range = (0.0, 1.0)
     prob = DifferentialEquations.ODEProblem(schrod_eq, initial_state, s_range, annealing_time)
-    sol = DifferentialEquations.solve(prob, abstol=abstol, reltol=reltol, alg_hints=[:nonstiff], kwargs...)
+    sol = DifferentialEquations.solve(prob; abstol=abstol, reltol=reltol, alg_hints=[:nonstiff], saveat=[1.0], kwargs...)
 
     state = sol(1)
     density = state * state'

test/common.jl

@@ -45,7 +45,7 @@ end
 
 function two_spin_ρ(t; s=1.0)
     s0 = 1/2*cos(π/4*s*sqrt(1+64*t^2/π^2))*sqrt(1-sin(π/2*s)) +
-        (8im*t*sqrt(1-sin(π/2*s)) + sqrt(1+sin(π/2*s)))*sin(π/4*s*sqrt(1+64*t^2/π^2))/(2*sqrt(1+64*t^2/π^2))
+        (8im*t*sqrt(1-sin(π/2*s))/π + sqrt(1+sin(π/2*s)))*sin(π/4*s*sqrt(1+64*t^2/π^2)) / (2*sqrt(1+64*t^2/π^2))
     s1 = -(cos(π/4*s*sqrt(1+64*t^2/π^2))*(1+sin(π/2*s)) + 
         ((-cos(π/2*s) + 8im*t*(1+sin(π/2*s))/π)*(sin(π/4*s*sqrt(1+64*t^2/π^2))))/(sqrt(1+64*t^2/π^2))
         )/(2*sqrt(1+sin(π/2*s)))

test/simulate.jl

@@ -429,6 +429,11 @@ end
         @test isapprox(two_spin_ρ(1.0), ρ)
     end
 
+    @testset "1 qubit, function schedule, analytical solution, with kwargs" begin
+        ρ = simulate_de(one_spin_model, 1.0, AS_CIRCULAR, 1e-6, saveat=[1])
+        @test isapprox(one_spin_ρ(1.0), ρ)
+    end
+
     @testset "1 qubit, function schedule, constant terms" begin
         ρ = simulate_magnus_generic(one_spin_model, 1.0, AS_CIRCULAR, 100, 4, constant_field_x=[1], constant_field_z=[1])
         ρ_de = simulate_de(one_spin_model, 1.0, AS_CIRCULAR, 1e-6, constant_field_x=[1], constant_field_z=[1])