Change DifferentialEquations dependence to OrdinaryDiffEq, Fix simula…

…te_de crashes (#50)

Project.toml

@@ -6,15 +6,15 @@ version = "0.2.1"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 CSV = "~0.10, ~0.9, ~0.8, ~0.7"
-DifferentialEquations = "~7, ~6"
+OrdinaryDiffEq = "~6, ~5"
 JSON = "~0.21, ~0.20, ~0.19, ~0.18"
 julia = "^1"
 

src/QuantumAnnealing.jl

@@ -4,7 +4,7 @@ module QuantumAnnealing
     import SparseArrays
     import JSON
     import Printf
-    import DifferentialEquations
+    import OrdinaryDiffEq
 
     include("base.jl")
     include("ising.jl")

src/simulate_de.jl

@@ -1,6 +1,6 @@
 
 """
-simulator that uses the DifferentialEquations.jl differential equation solver
+simulator that uses the OrdinaryDiffEq.jl differential equation solver
 to solve the Schrodinger Equation, rather than the Magnus Expansion.  This
 simulator can run into issues at higher anneal times.  The API is the same as
 for the simulate function.
@@ -40,8 +40,8 @@ function simulate_de(ising_model, annealing_time, annealing_schedule, reltol; ab
     schrod_eq(state, time, s) = -im * time * H(s) * state
 
     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], saveat=[1.0], kwargs...)
+    prob = OrdinaryDiffEq.ODEProblem(schrod_eq, initial_state, s_range, annealing_time)
+    sol = OrdinaryDiffEq.solve(prob, OrdinaryDiffEq.Tsit5(); abstol=abstol, reltol=reltol, saveat=[1.0], kwargs...)
 
 
     state = sol(1)

test/simulate.jl

@@ -459,40 +459,32 @@ 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])
-
         @test isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "1 qubit, csv schedule, analytical solution" begin
         ρ = simulate_de(one_spin_model, 1.0, AS_CIRCULAR_pwl_csv_1000, 1e-6)
-
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-6)
-        @test !isapprox(one_spin_ρ(1.0), ρ, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 1.0, AS_CIRCULAR, 100, 4)
         ρ_de = simulate_de(ising_model, 1.0, AS_CIRCULAR, 1e-6)
-
         @test isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule, long annealing time" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 1000.0, AS_CIRCULAR, 4000, 4)
-        ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-7)
-
+        ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-8)
         @test isapprox(ρ, ρ_de, atol=1e-6)
-        @test !isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule, adaptive tolerance" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 100.0, AS_CIRCULAR, 500, 4)
         ρ_de = simulate_de(ising_model, 100.0, AS_CIRCULAR, silence=true)
-
-        @test isapprox(ρ, ρ_de, atol=1e-7)
-        @test !isapprox(ρ, ρ_de, atol=1e-8)
+        @test isapprox(ρ, ρ_de, atol=1e-6)
     end
 end