Inplace exponential (#37)

* Add inplace matrix exp in sim, fix kwargs bug and reduce mem use in de

* Update Changelog

* Fix analytic two spin rho

* Add tests for analytic density matrices

* Remove LinearAlgebra from Project.toml, restore Test location in Project.toml

CHANGELOG.md

@@ -2,6 +2,7 @@ QuantumAnnealing.jl Change Log
 ==============================
 
 ### Staged
+- Reduce memory usage in simulate and simulate_de
 - Add Gibbs distribution (#34)
 - Fix kwargs bug in simulate_de (#39)
 - Fix analytic two spin model for non-integer s (#40)

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"
 
@@ -19,8 +20,8 @@ julia = "^1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+
 [targets]
-test = ["Test","LinearAlgebra","Random"]
+test = ["Test", "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,8 @@ 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/base.jl

@@ -167,6 +167,31 @@ end
     end
 end
 
+@testset "Analytic Density Matrices" begin
+    @testset "One Spin" begin
+        one_spin_s_0 = one_spin_ρ(2, s=0)
+        one_spin_s_0123 = one_spin_ρ(2, s=0.123)
+        one_spin_s_05 = one_spin_ρ(2, s=0.5)
+        one_spin_s_1 = one_spin_ρ(2, s=1)
+
+        @test isapprox(sum([one_spin_s_0[i,i] for i in 1:2]), 1)
+        @test isapprox(sum([one_spin_s_0123[i,i] for i in 1:2]), 1)
+        @test isapprox(sum([one_spin_s_05[i,i] for i in 1:2]), 1)
+        @test isapprox(sum([one_spin_s_1[i,i] for i in 1:2]), 1)
+    end
+
+    @testset "Two Spin" begin
+        two_spin_s_0 = two_spin_ρ(2, s=0)
+        two_spin_s_0123 = two_spin_ρ(2, s=0.123)
+        two_spin_s_05 = two_spin_ρ(2, s=0.5)
+        two_spin_s_1 = two_spin_ρ(2, s=1)
+
+        @test isapprox(sum([two_spin_s_0[i,i] for i in 1:4]), 1)
+        @test isapprox(sum([two_spin_s_0123[i,i] for i in 1:4]), 1)
+        @test isapprox(sum([two_spin_s_05[i,i] for i in 1:4]), 1)
+        @test isapprox(sum([two_spin_s_1[i,i] for i in 1:4]), 1)
+    end
+end
 
 @testset "csv annealing schedules" begin
 

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

@@ -451,6 +451,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])