Any-Order Magnus Expansion

CHANGELOG.md

@@ -2,7 +2,12 @@ QuantumAnnealing.jl Change Log
 ==============================
 
 ### Staged
-- nothing
+- Add a generic Magnus expansion solver for any order
+- Update hard-coded Magnus expansion solver to support orders 1 through 4
+- Change default solver order from 2 to 4 (breaking)
+- Update d-wave simulation tools to use adaptive solvers (breaking)
+- Reversed coefficient ordering of `get_quadratic_coefficients` (breaking)
+- Remove export of a variety of internal helper functions (breaking)
 
 ### v0.1.0
 - Add variant of `solve_de` with adaptive solve tolerance

src/base.jl

@@ -14,23 +14,22 @@ struct AnnealingSchedule
 end
 
 """
-A short hand AnnealingSchedule constructor that uses the default_initial_state,
+A short hand AnnealingSchedule constructor that uses the initial_state_default,
 which is the most common case for the conventions of this implementation.
 """
-AnnealingSchedule(A,B) = AnnealingSchedule(A, B, default_initial_state)
-
+AnnealingSchedule(A,B) = AnnealingSchedule(A, B, initial_state_default)
 
 #predefining Pauli Matrices
-const IMAT = SparseArrays.sparse([1,2], [1,2], [1.0+0im;1.0+0im])
-const XMAT = SparseArrays.sparse([1,2], [2,1], [1.0+0im;1.0+0im])
-const YMAT = SparseArrays.sparse([1,2], [2,1], [-im;im])
-const ZMAT = SparseArrays.sparse([1,2], [1,2], [1.0+0im;-1.0+0im])
+const _IMAT = SparseArrays.sparse([1,2], [1,2], [1.0+0im;1.0+0im])
+const _XMAT = SparseArrays.sparse([1,2], [2,1], [1.0+0im;1.0+0im])
+const _YMAT = SparseArrays.sparse([1,2], [2,1], [-im;im])
+const _ZMAT = SparseArrays.sparse([1,2], [1,2], [1.0+0im;-1.0+0im])
 
 
 """
 ground state of sum_i A(0) X_i where A(0) > 0 and B(0) = 0
 """
-function default_initial_state(n)
+function initial_state_default(n)
     return complex(foldl(kron,[[1;-1] for i in 1:n]) ./ (2^(n/2)))
 end
 

src/dwave.jl

@@ -4,7 +4,7 @@
 """
 ground state of sum_i A(0)X_i where A(0) < 0 and B(0) = 0
 """
-function default_dwave_initial_state(n)
+function initial_state_default_dwave(n)
     return complex(ones(2^n)./(2^(n/2)))
 end
 
@@ -16,17 +16,15 @@ NOTE: users are strongly encouraged to download annealing schedules for specific
 D-Wave Systems devices and load them using `parse_dwave_annealing_schedule`.
 """
 const AS_DW_QUADRATIC = AnnealingSchedule(
-    function A(s)
+    (s) -> begin
         if s >= 0.69
             return 0
         else
             return (6.366401*((1.449275)^2*s^2 + (-2.898551)*s + 1.0)*(2.0*π))/-2.0
         end
     end,
-    function B(s)
-        return (14.55571*(0.85*s^2 + 0.15*s + 0.0)*(2.0*π))/2.0
-    end,
-    default_dwave_initial_state
+    (s) -> (14.55571*(0.85*s^2 + 0.15*s + 0.0)*(2.0*π))/2.0,
+    initial_state_default_dwave
 )
 
 
@@ -133,7 +131,7 @@ function to take a CSV of DWave annealing schedule values and convert it into
 an annealing schedule usable by the simulator.
 valid values for interpolation are :none, :linear, :quadratic
 """
-function parse_dwave_annealing_schedule(infile; header=1, delim=',', interpolation=:linear, initial_state=default_dwave_initial_state)
+function read_dwave_annealing_schedule(infile; header=1, delim=',', interpolation=:linear, initial_state=initial_state_default_dwave)
     s_values = Float64[]
     a_values = Float64[]
     b_values = Float64[]
@@ -186,17 +184,50 @@ function parse_dwave_annealing_schedule(infile; header=1, delim=',', interpolati
     )
 end
 
+
+"""
+Function to modify an existing annealing schedule to use a customized
+annealing schedule (asch).  These parameters are the same as those
+used in a dwisc call or a dwave schedule.
+Inputs:
+annealing_schedule - annealing_schedule
+
+Parameters:
+asch - This is the annealing-schedule parameter.  This is a list of tuples of the form
+       [(s₀,s_effective₀), (s₀,s_effective₁), ..., (sₙ,s_effectiveₙ)].
+"""
+function annealing_protocol_dwave(annealing_schedule::AnnealingSchedule; asch=[(0,0) (1,1)])
+    asch_slopes = zeros(length(asch)-1)
+    for i in 1:(length(asch)-1)
+        s0,s_eff_0 = asch[i]
+        s1,s_eff_1 = asch[i+1]
+        asch_slopes[i] = (s_eff_1 - s_eff_0)/(s1 - s0)
+    end
+
+    #branchless piecewise function using linear interpolation from y = m*(x-x0) + y0
+    function asch_func(s)
+        return sum([(asch_slopes[i]*(s-asch[i][1]) + asch[i][2]) * (asch[i][1] <= s < asch[i+1][1]) for i = 1:(length(asch)-1)]) + ((s == asch[end][1])*asch[end][2])
+    end
+
+    new_annealing_schedule = AnnealingSchedule(
+        s -> annealing_schedule.A(asch_func(s)),
+        s -> annealing_schedule.B(asch_func(s))
+    )
+    return new_annealing_schedule
+end
+
+
 """
 function that allows for simulation from a bqpjson data file
 """
-function simulate_bqpjson(infile, outfile, annealing_time, annealing_schedule, steps; simulated_num_reads=1e17, scale=1.0)
+function simulate_bqpjson(infile, outfile, annealing_time, annealing_schedule; simulated_num_reads=1e17, scale=1.0, kwargs...)
     ising_model, qubit_ids = read_bqpjson(infile)
     n = length(qubit_ids)
     for (k,v) in ising_model
         ising_model[k] = v*scale
     end
 
-    ρ = simulate(ising_model, annealing_time, annealing_schedule, steps)
+    ρ = simulate(ising_model, annealing_time, annealing_schedule; kwargs...)
 
     write_dwisc(outfile, ρ, ising_model, qubit_ids, simulated_num_reads=simulated_num_reads, annealing_time=annealing_time)
 end
@@ -206,7 +237,7 @@ end
 function that allows for simulation with x and z noise from a bqpjson data file.
 The `x_bias` and `z_bias` parameters provide vectors of noise realizations.
 """
-function simulate_noisy_bqpjson(infile, outfile, annealing_time, annealing_schedule, steps; simulated_num_reads=1e17, scale=1.0, x_bias::Vector{<:Any}=[], z_bias::Vector{<:Any}=[])
+function simulate_bqpjson_noisy(infile, outfile, annealing_time, annealing_schedule; simulated_num_reads=1e17, scale=1.0, x_bias::Vector{<:Any}=[], z_bias::Vector{<:Any}=[], kwargs...)
     if length(x_bias) > 0 && length(z_bias) > 0 && length(x_bias) != length(z_bias)
         error("x_bias and z_bias require the same number of parameters given, $(length(x_bias)) and $(length(z_bias)) respectively")
     end
@@ -233,7 +264,7 @@ function simulate_noisy_bqpjson(infile, outfile, annealing_time, annealing_sched
         x_field = x_bias[shot]
         z_field = z_bias[shot]
 
-        ρ = simulate(ising_model, annealing_time, annealing_schedule, steps, constant_field_x=[x_field], constant_field_z=[z_field])
+        ρ = simulate(ising_model, annealing_time, annealing_schedule, constant_field_x=[x_field], constant_field_z=[z_field]; kwargs...)
 
         accumulator = accumulator + ρ
     end
@@ -313,34 +344,3 @@ function write_dwisc(outfile::String, ρ, ising_model, qubit_ids; simulated_num_
         write(io, json_string)
     end
 end
-
-"""
-Function to modify an existing annealing schedule to use a customized
-annealing schedule (asch).  These parameters are the same as those
-used in a dwisc call or a dwave schedule.
-Inputs:
-annealing_schedule - annealing_schedule
-
-Parameters:
-asch - This is the annealing-schedule parameter.  This is a list of tuples of the form
-       [(s₀,s_effective₀), (s₀,s_effective₁), ..., (sₙ,s_effectiveₙ)].
-"""
-function dwave_annealing_protocol(annealing_schedule::AnnealingSchedule; asch=[(0,0) (1,1)])
-    asch_slopes = zeros(length(asch)-1)
-    for i in 1:(length(asch)-1)
-        s0,s_eff_0 = asch[i]
-        s1,s_eff_1 = asch[i+1]
-        asch_slopes[i] = (s_eff_1 - s_eff_0)/(s1 - s0)
-    end
-
-    #branchless piecewise function using linear interpolation from y = m*(x-x0) + y0
-    function asch_func(s)
-        return sum([(asch_slopes[i]*(s-asch[i][1]) + asch[i][2]) * (asch[i][1] <= s < asch[i+1][1]) for i = 1:(length(asch)-1)]) + ((s == asch[end][1])*asch[end][2])
-    end
-
-    new_annealing_schedule = AnnealingSchedule(
-        s -> annealing_schedule.A(asch_func(s)),
-        s -> annealing_schedule.B(asch_func(s))
-    )
-    return new_annealing_schedule
-end