Adaptive Tolerance Variant of simulate_de

CHANGELOG.md

@@ -2,6 +2,7 @@ QuantumAnnealing.jl Change Log
 ==============================
 
 ### Staged
+- Add variant of `solve_de` with adaptive solve tolerance
 - Add tools for working with classical Ising models (#8)
 - Add data processing tools (#6)
 

src/simulate_de.jl

@@ -11,13 +11,14 @@ ising_model - ising model represented as a dictionary.  The qubits
               are numbers.
               For Example: im = Dict((1,) => 1, (2,) => 0.5, (1,2) => 2)
 annealing_schedule - The annealing schedule, of the form given by the struct
+reltol - the relative tolerance that will be used in the ODE simulation
 
 Parameters:
 initial_state - Initial state vector. Defaults to uniform superposition state on n qubits
 constant_field_x - vector of constant biases in the X basis on each qubit. Default is zeros(n)
 constant_field_z - vector of constant biases in the Z basis on each qubit. Default is zeros(n)
 """
-function simulate_de(ising_model, annealing_time, annealing_schedule; initial_state=nothing, constant_field_x=nothing, constant_field_z=nothing)
+function simulate_de(ising_model, annealing_time, annealing_schedule, reltol; abstol=reltol/1e3, initial_state=nothing, constant_field_x=nothing, constant_field_z=nothing, kwargs...)
     n = _check_ising_model_ids(ising_model)
 
     if initial_state == nothing
@@ -40,11 +41,65 @@ function simulate_de(ising_model, annealing_time, annealing_schedule; initial_st
 
     s_range = (0.0, 1.0)
     prob = DifferentialEquations.ODEProblem(schrod_eq, initial_state, s_range, annealing_time)
-    sol = DifferentialEquations.solve(prob)
+    sol = DifferentialEquations.solve(prob, abstol=abstol, reltol=reltol, alg_hints=[:nonstiff], kwargs...)
 
     state = sol(1)
     density = state * state'
 
     return density
 end
 
+
+
+"""
+A simplified interface to the `simulate_de` routine that determines a
+suitable tolerance values to ensure a high accuracy simulation.
+The parameters `mean_tol` and `max_tol` specify the desired simulation accuracy.
+The `silence` parameter can be used to suppress the progress log.
+"""
+function simulate_de(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule; reltol=1e-4, mean_tol=1e-6, max_tol=1e-4, iteration_limit=100, silence=false, kwargs...)
+    start_time = time()
+    mean_delta = mean_tol + 1.0
+    max_delta = max_tol + 1.0
+
+    if !silence
+        println()
+        println("iter |    reltol    |    max(Δ)    |    mean(Δ)   |")
+    end
+
+    ρ_prev = simulate_de(ising_model, annealing_time, annealing_schedule, reltol; kwargs...)
+
+    iteration = 1
+    while mean_delta >= mean_tol || max_delta >= max_tol
+        reltol /= 10.0
+
+        ρ = simulate_de(ising_model, annealing_time, annealing_schedule, reltol;  kwargs...)
+
+        ρ_delta = abs.(ρ .- ρ_prev)
+        mean_delta = sum(ρ_delta)/length(ρ_delta)
+        max_delta = maximum(ρ_delta)
+
+        !silence && Printf.@printf("%4d | %e | %e | %e |\n", iteration, reltol, max_delta, mean_delta)
+
+        ρ_prev = ρ
+        iteration += 1
+        if iteration > iteration_limit
+            error("iteration limit reached in simulate function without reaching convergence criteria")
+        end
+    end
+
+    if !silence
+        println("")
+        println("\033[1mconverged\033[0m")
+        Printf.@printf("   iterations........: %d\n", iteration-1)
+        Printf.@printf("   simulation reltol.: %e\n", reltol)
+        Printf.@printf("   maximum difference: %e <= %e\n", max_delta, max_tol)
+        Printf.@printf("   mean difference...: %e <= %e\n", mean_delta, mean_tol)
+        Printf.@printf("   runtime (seconds).: %f\n", time()-start_time)
+        println("")
+    end
+
+    return ρ_prev
+end
+
+

test/simulate.jl

@@ -399,30 +399,57 @@ end
 end
 
 
+
 @testset "simulate_de" begin
 
     @testset "1 qubit, function schedule, analytical solution" begin
-        ρ = simulate_de(single_spin_model, 1.0, AS_CIRCULAR)
+        ρ = simulate_de(single_spin_model, 1.0, AS_CIRCULAR, 1e-6)
 
         @test isapprox(single_spin_analytic_ρ, ρ)
     end
 
     @testset "1 qubit, function schedule, constant terms" begin
-        ρ = simulate(single_spin_model, 1.0, AS_CIRCULAR, constant_field_x = [1], constant_field_z = [1])
-        ρ_de = simulate_de(single_spin_model, 1.0, AS_CIRCULAR, constant_field_x = [1], constant_field_z = [1])
+        ρ = simulate(single_spin_model, 1.0, AS_CIRCULAR, 100, constant_field_x=[1], constant_field_z=[1])
+        ρ_de = simulate_de(single_spin_model, 1.0, AS_CIRCULAR, 1e-6, constant_field_x=[1], constant_field_z=[1])
 
-        @test isapprox(ρ, ρ_de, atol=1e-7)
-        @test !isapprox(ρ, ρ_de, atol=1e-8)
+        @test isapprox(ρ, ρ_de, atol=1e-8)
+        @test !isapprox(ρ, ρ_de, atol=1e-9)
     end
 
     @testset "1 qubit, csv schedule, analytical solution" begin
-        ρ = simulate_de(single_spin_model, 1.0, AS_CIRCULAR_pwl_csv_1000)
-        simulation_prob = z_measure_probabilities(ρ)[1]
+        ρ = simulate_de(single_spin_model, 1.0, AS_CIRCULAR_pwl_csv_1000, 1e-6)
 
         @test isapprox(single_spin_analytic_ρ, ρ, atol=1e-6)
         @test !isapprox(single_spin_analytic_ρ, ρ, atol=1e-7)
     end
 
+    @testset "2 qubit, function schedule" begin
+        ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
+        ρ = simulate(ising_model, 1.0, AS_CIRCULAR, 100)
+        ρ_de = simulate_de(ising_model, 1.0, AS_CIRCULAR, 1e-6)
+
+        @test isapprox(ρ, ρ_de, atol=1e-8)
+        @test !isapprox(ρ, ρ_de, atol=1e-9)
+    end
+
+    @testset "2 qubit, function schedule, long annealing time" begin
+        ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
+        ρ = simulate(ising_model, 1000.0, AS_CIRCULAR, 4000)
+        ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-7)
+
+        @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(ising_model, 100.0, AS_CIRCULAR, 1000)
+        ρ_de = simulate_de(ising_model, 100.0, AS_CIRCULAR, silence=true)
+
+        @test isapprox(ρ, ρ_de, atol=1e-6)
+        @test !isapprox(ρ, ρ_de, atol=1e-7)
+    end
+
 end