diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5166f39..67d5fa5 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/src/base.jl b/src/base.jl
index 546fe08..ec5b533 100644
--- a/src/base.jl
+++ b/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
 
diff --git a/src/dwave.jl b/src/dwave.jl
index 03a9cfd..80a818e 100644
--- a/src/dwave.jl
+++ b/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
diff --git a/src/simulate.jl b/src/simulate.jl
index 1d03368..89e0029 100644
--- a/src/simulate.jl
+++ b/src/simulate.jl
@@ -1,5 +1,4 @@
 
-
 """
 An AnnealingSchedule implementing uniform circular motion with an analytical
 solution on a single qubit
@@ -16,75 +15,71 @@ const AS_QUADRATIC = AnnealingSchedule(s -> (1.0-s)^2, s -> s^2)
 converts an arbitrary function `f(s)` into a quadratic form based on
 interpolation between two extreme points `s0` and `s1`
 """
-function get_function_coefficients(f, s0, s1)
-    smid = (s0+s1)/2.0
+function _get_quadratic_coefficients(f, s0, s1)
+    sm = (s0+s1)/2.0
 
-    b = f.([s0, smid, s1])
-    A = [(s0)^2 (s0) 1
-         (smid)^2 (smid) 1
-         (s1)^2 (s1) 1]
+    b = f.([s0, sm, s1])
+    A = [
+        1 s0 s0^2;
+        1 sm sm^2;
+        1 s1 s1^2
+    ]
 
     x = A\b
 
     return x
 end
 
-function lie_bracket(A, B)
-    return -im * (A*B - B*A)
-end
+"shifts a quadatric function by a value of x"
+function _shift_quadratic_coefficients(x, c0, c1, c2)
+    c0_shift = c0 + c1*x + c2*x^2
+    c1_shift = c1 + 2*c2*x
+    c2_shift = c2
 
-function tensor_sum_single_qubit(mat, n::Int)
-    return sum([foldl(kron,[j == i ? mat : IMAT for i in n:-1:1]) for j in 1:n])
+    return [c0_shift, c1_shift, c2_shift]
 end
 
-function tensor_sum_single_qubit(mat, n::Int, weights::Vector)
-    return sum([foldl(kron,[j == i ? weights[j] * mat : IMAT for i in n:-1:1]) for j in 1:n])
+
+function _lie_bracket(A, B)
+    return -im * (A*B - B*A)
 end
 
-function sum_x(n::Int)
-    return tensor_sum_single_qubit(XMAT, n)
+function _tensor_sum_single_qubit(mat, n::Int)
+    return sum([foldl(kron,[j == i ? mat : _IMAT for i in n:-1:1]) for j in 1:n])
 end
 
-function sum_x(n::Int, w::Vector)
-    return tensor_sum_single_qubit(XMAT, n, w)
+function _tensor_sum_single_qubit(mat, n::Int, weights::Vector)
+    return sum([foldl(kron,[j == i ? weights[j] * mat : _IMAT for i in n:-1:1]) for j in 1:n])
 end
 
-function sum_y(n::Int)
-    return tensor_sum_single_qubit(YMAT, n)
+function _sum_X(n::Int)
+    return _tensor_sum_single_qubit(_XMAT, n)
 end
 
-function sum_y(n::Int, w::Vector)
-    return tensor_sum_single_qubit(YMAT, n, w)
+function _sum_X(n::Int, w::Vector)
+    return _tensor_sum_single_qubit(_XMAT, n, w)
 end
 
-function sum_z(n::Int)
-    return tensor_sum_single_qubit(ZMAT, n)
+function _sum_Y(n::Int)
+    return _tensor_sum_single_qubit(_YMAT, n)
 end
 
-function sum_z(n::Int, w::Vector)
-    return tensor_sum_single_qubit(ZMAT, n, w)
+function _sum_Y(n::Int, w::Vector)
+    return _tensor_sum_single_qubit(_YMAT, n, w)
 end
 
-function zizj_vectorized(n::Int, i::Int, j::Int, J_val)
-    Ivec = [1.0+0im,1.0+0im]
-    Zvec = [1.0+0im,-1.0+0im]
-    matvec = [(k == i || k == j) ? Zvec : Ivec for k in n:-1:1]
-    return J_val * foldl(kron, matvec)
+function _sum_Z(n::Int)
+    return _tensor_sum_single_qubit(_ZMAT, n)
 end
 
-function sum_zizj(n, J::Dict)
-    sum_zs = zeros(2^n)
-    for (i,j) in keys(J)
-        J_val = J[(i,j)]
-        sum_zs +=  zizj_vectorized(n, i, j, J_val)
-    end
-    return SparseArrays.spdiagm(sum_zs)
+function _sum_Z(n::Int, w::Vector)
+    return _tensor_sum_single_qubit(_ZMAT, n, w)
 end
 
-function sum_z_tup(n, tup, w)
+function _kron_Z(n::Int, t, w::Real)
     Ivec = [1;1]
     Zvec = [1;-1]
-    matvec = [k in tup ? Zvec : Ivec for k in n:-1:1]
+    matvec = [k in t ? Zvec : Ivec for k in n:-1:1]
     return SparseArrays.spdiagm(complex(w * foldl(kron, matvec)))
 end
 
@@ -100,64 +95,311 @@ ising_model - ising model represented as a dictionary.  The qubits
 annealing_schedule - The annealing schedule, of the form given by the struct
 s - the imaginary timestep. This should usually be in the range from 0.0-to-1.0
 """
-function transverse_ising_hamiltonian(ising_model::Dict, annealing_schedule::AnnealingSchedule, s::Real)
+function hamiltonian_transverse_ising(ising_model::Dict, annealing_schedule::AnnealingSchedule, s::Real)
     n = _check_ising_model_ids(ising_model)
 
-    x_component = sum_x(n)
+    x_component = _sum_X(n)
     z_component = SparseArrays.spzeros(2^n, 2^n)
     for (tup,w) in ising_model
-        z_component = z_component + sum_z_tup(n, tup, w)
+        z_component += _kron_Z(n, tup, w)
     end
 
     return annealing_schedule.A(s) * x_component + annealing_schedule.B(s) * z_component
 end
 
 
-function integral_1_sched(a_2, a_1, a_0, s0, δ)
-    #TODO: Possibly change to integrate from s0 to s1 to allow
-    #      for a more fine-grained integration with As and Bs
-    #      from spline
-    #δ = s1 - s0
-    left = [(2*(3*s0^2+3*s0*δ+δ^2)) (3*(2*s0+δ)) 6] * δ / 6
-    right = [a_2, a_1, a_0]
-    return (left * right)[1]
+function simulate_fixed_order(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule, steps::Int, order::Int; initial_state=nothing, state_steps=nothing)
+    if steps < 2
+        error("at least two steps are required by simulate_fixed_order, given $(steps)")
+    end
+    if !(1 <= order && order <= 4)
+        error("simulate_fixed_order only supports orders from 1-to-4, given $(order)")
+    end
+
+    n = _check_ising_model_ids(ising_model)
+
+    if initial_state == nothing
+        initial_state = annealing_schedule.init_default(n)
+    end
+
+    track_states = !(state_steps == nothing)
+
+    t0 = 0
+    s0 = 0
+
+    R0 = initial_state * initial_state'
+
+    x_component = _sum_X(n)
+    z_component = SparseArrays.spzeros(2^n, 2^n)
+    for (tup,w) in ising_model
+        z_component = z_component + _kron_Z(n, tup, w)
+    end
+
+    H_parts = _H_parts(x_component, z_component, order)
+
+    s_steps = range(0, 1, length=steps)
+    R_current = R0
+    U = foldl(kron, [_IMAT for i = 1:n])
+
+    if track_states
+        push!(state_steps, R_current)
+    end
+
+    for i in 1:(steps-1)
+        s0 = s_steps[i]
+        s1 = s_steps[i+1]
+        δs = s1 - s0
+
+        Ω_list = _Ω_list(annealing_time, s0, s1, annealing_schedule, H_parts, order)
+
+        #for (i,Ωi) in enumerate(Ω_list)
+        #   println("Ω_$(i)")
+        #   display(Matrix(_hamiltonian_eval2(δs, Ωi)))
+        #end
+
+        U_next = exp(Matrix(sum(Ω_list)))
+        U = U_next * U
+
+        if track_states
+            R_current = U * R0 * U'
+            push!(state_steps, R_current)
+        end
+    end
+
+    return U * R0 * U'
 end
 
-function integral_2_sched(a_2, a_1, a_0, b_2, b_1, b_0, s0, δ)
-    #TODO: see integral_1_sched todo.
-    #returns the magnus expansions double integral for 
-    #(A(s1)B(s2) - A(s2)B(s1))
-    #δ = s0 - s1
-    left = [a_2 a_1 a_0]
-    right = [b_2,b_1,b_0]
-    integ_mat = [(0) (s0^2 + s0*δ + δ^2/5) (2*s0 + δ);
-                 (-(s0^2 + s0*δ + δ^2/5)) (0) (1);
-                 -(2*s0+δ) (-1) (0)] * δ^3 / 6
-    return (left * integ_mat * right)[1]
+
+function _H_parts(x_component, z_component, order::Int)
+    @assert(1 <= order && order <= 4)
+
+    parts = Dict{Any,Any}(
+        (1,) => x_component,
+        (2,) => z_component
+    )
+
+    if order >= 2
+        parts[(2,1)] = _lie_bracket(x_component, z_component)
+    end
+    if order >= 3
+        parts[(3,1)] = _lie_bracket(x_component, parts[(2,1)])
+        parts[(3,2)] = _lie_bracket(parts[(2,1)], z_component)
+    end
+    if order >= 4
+        parts[(4,1)] = _lie_bracket(x_component, parts[(3,1)])
+        parts[(4,2)] = _lie_bracket(x_component, parts[(3,2)])
+        parts[(4,3)] = _lie_bracket(z_component, parts[(3,2)])
+    end
+
+    return parts
 end
 
 
-"""
-Main function for performing quantum annealing simulation via a Magnus Expansion.
-Noise can be simulated by running multiple times with randomized constant fields.
+function _Ω_list(annealing_time::Real, s0::Real, s1::Real, annealing_schedule::AnnealingSchedule, H_parts::Dict, order::Int)
+    @assert(1 <= order && order <= 4)
+    δs = s1 - s0
+    δst = annealing_time*δs
+
+    aqc = _get_quadratic_coefficients(annealing_schedule.A, s0, s1)
+    bqc = _get_quadratic_coefficients(annealing_schedule.B, s0, s1)
+
+    aqc = _shift_quadratic_coefficients(s0, aqc...)
+    bqc = _shift_quadratic_coefficients(s0, bqc...)
+
+    aqc = [aqc[1], δs*aqc[2], δs^2*aqc[3]]
+    bqc = [bqc[1], δs*bqc[2], δs^2*bqc[3]]
+
+    Ω1 = -im*δst*(_integral_1(aqc)*H_parts[(1,)] + _integral_1(bqc)*H_parts[(2,)])
+    Ω_list = [Ω1]
+
+    if order >= 2
+        Ω2 = -im*δst^2/2*(_integral_21(aqc,bqc)*H_parts[(2,1)])
+        push!(Ω_list, Ω2)
+    end
+    if order >= 3
+        Ω3 = -im*δst^3/6*(_integral_31(aqc,bqc)*H_parts[(3,1)] + _integral_31(bqc,aqc)*H_parts[(3,2)])
+        push!(Ω_list, Ω3)
+    end
+    if order >= 4
+        Ω4 = -im*δst^4/6*(_integral_41(aqc,bqc)*H_parts[(4,1)] + _integral_42(aqc,bqc)*H_parts[(4,2)] + _integral_41(bqc,aqc)*H_parts[(4,3)])
+        push!(Ω_list, Ω4)
+    end
+
+    return Ω_list
+end
+
+
+function _integral_1(u::Vector)
+    return u[1] + u[2]/2 + u[3]/3
+end
+
+function _integral_21(u::Vector, v::Vector)
+    return u[2]*v[1]/6 + u[3]*v[1]/6 -
+        u[1]*v[2]/6 + u[3]*v[2]/30 -
+        u[1]*v[3]/6 - u[2]*v[3]/30
+end
+
+function _integral_31(u::Vector, v::Vector)
+    return u[2]^2*v[1]/40 - u[1]*u[3]*v[1]/60 +
+        u[2]*u[3]*v[1]/24 + 5*u[3]^2*v[1]/252 -
+        u[1]*u[2]*v[2]/40 - u[1]*u[3]*v[2]/30 +
+        u[2]*u[3]*v[2]/840 + u[3]^2*v[2]/360 +
+        u[1]^2*v[3]/60 - u[1]*u[2]*v[3]/120 -
+        u[2]^2*v[3]/840 - 5*u[1]*u[3]*v[3]/252 -
+        u[2]*u[3]*v[3]/360
+end
+
+function _integral_41(u::Vector, v::Vector)
+    return -u[1]^2*u[2]*v[1]/120 - u[1]*u[2]^2*v[1]/120 -
+        u[2]^3*v[1]/840 - u[1]^2*u[3]*v[1]/120 -
+        19*u[1]*u[2]*u[3]*v[1]/1260 - u[2]^2*u[3]*v[1]/360 -
+        17*u[1]*u[3]^2*v[1]/2520 - u[2]*u[3]^2*v[1]/504 -
+        u[3]^3*v[1]/2520 + u[1]^3*v[2]/120 +
+        u[1]^2*u[2]*v[2]/120 + u[1]*u[2]^2*v[2]/840 +
+        11*u[1]^2*u[3]*v[2]/2520 - u[1]*u[2]*u[3]*v[2]/1260 -
+        u[2]^2*u[3]*v[2]/2520 - u[1]*u[3]^2*v[2]/720 -
+        u[2]*u[3]^2*v[2]/2016 - u[3]^3*v[2]/7920 +
+        u[1]^3*v[3]/120 + 3*u[1]^2*u[2]*v[3]/280 +
+        u[1]*u[2]^2*v[3]/280 + u[2]^3*v[3]/2520 +
+        17*u[1]^2*u[3]*v[3]/2520 + 17*u[1]*u[2]*u[3]*v[3]/5040 +
+        u[2]^2*u[3]*v[3]/2016 + u[1]*u[3]^2*v[3]/2520 +
+        u[2]*u[3]^2*v[3]/7920
+end
+
+function _integral_42(u::Vector, v::Vector)
+    return u[1]*u[2]*v[1]^2/60 + u[2]^2*v[1]^2/120 +
+        u[1]*u[3]*v[1]^2/60 + 19*u[2]*u[3]*v[1]^2/1260 +
+        17*u[3]^2*v[1]^2/2520 - u[1]^2*v[1]*v[2]/60 +
+        u[2]^2*v[1]*v[2]/420 + 2*u[1]*u[3]*v[1]*v[2]/315 +
+        2*u[2]*u[3]*v[1]*v[2]/315 + 17*u[3]^2*v[1]*v[2]/5040 -
+        u[1]^2*v[2]^2/120 - u[1]*u[2]*v[2]^2/420 +
+        u[1]*u[3]*v[2]^2/1260 + u[2]*u[3]*v[2]^2/1260 +
+        u[3]^2*v[2]^2/2016 - u[1]^2*v[1]*v[3]/60 -
+        2*u[1]*u[2]*v[1]*v[3]/315 - u[2]^2*v[1]*v[3]/1260 +
+        u[2]*u[3]*v[1]*v[3]/1680 + u[3]^2*v[1]*v[3]/1260 -
+        19*u[1]^2*v[2]*v[3]/1260 - 2*u[1]*u[2]*v[2]*v[3]/315 -
+        u[2]^2*v[2]*v[3]/1260 - u[1]*u[3]*v[2]*v[3]/1680 +
+        u[3]^2*v[2]*v[3]/3960 - 17*u[1]^2*v[3]^2/2520 -
+        17*u[1]*u[2]*v[3]^2/5040 - u[2]^2*v[3]^2/2016 -
+        u[1]*u[3]*v[3]^2/1260 - u[2]*u[3]*v[3]^2/3960
+end
+
+
+"computes the n-th Bernoulli number divided by n-th factorial number"
+function _bernoulli_factorial(n)
+    if n == 1
+        return -0.5
+    end
+
+    v = Vector{Rational{BigInt}}(undef, n+1)
+    for m in 0:n
+        v[m + 1] = 1//(m+1)
+        for j in m:-1:1
+            v[j] = j*(v[j] - v[j+1])
+        end
+    end
+
+    return Float64(v[1]/factorial(big(n)))
+end
+
+"given to matricies, applies the commutator operation"
+function _matrix_commutator(a,b)
+    return a*b - b*a
+end
+
+"given two hamiltonians of orders 0-to-(n-1) and 0-to-(m-1), applies the _matrix_commutator operation"
+function _hamiltonian_commutator(h1::Vector, h2::Vector)
+    max_index = (length(h1)-1)+(length(h2)-1)+1
+    zeros = 0*(h1[1]*h2[1])
+    h_prod = [deepcopy(zeros) for i in 1:max_index]
+
+    for (i,m1) in enumerate(h1)
+        for (j,m2) in enumerate(h2)
+            pow = (i-1)+(j-1)
+            h_prod[pow+1] += _matrix_commutator(m1,m2)
+        end
+    end
+
+    return h_prod
+end
+
+"given a hamiltonian of orders 0-to-(n-1), multiplies all orders by a scalar"
+function _hamiltonian_eval(x::Real, h_list::Vector)
+    return sum(h*x^(i-1) for (i,h) in enumerate(h_list))
+end
+
+"given a hamiltonian of orders 0-to-(n-1), multiplies all orders by a scalar"
+function _hamiltonian_scalar(x::Real, h_list::Vector)
+    return [x*h for h in h_list]
+end
+
+"given a hamiltonian of orders 0-to-(n-1), returns an integrated hamiltonian of 0-to-n"
+function _hamiltonian_integrate(h::Vector)
+    ih = [i==1 ? (0*h[1]) : h[i-1] for i in 1:length(h)+1]
+
+    for (i,h) in enumerate(ih)
+        if i != 1
+            ih[i] = h./(i-1)
+        end
+    end
+
+    return ih
+end
+
+
+"given a list of order-based hamiltonians (which are also vectors), returns the sum of these"
+function _hamiltonian_sum(h_list::Vector)
+    max_index = maximum(length(h) for h in h_list)
+    zeros = 0*(h_list[1][1])
+    h_sum = [deepcopy(zeros) for i in 1:max_index]
+
+    for h in h_list
+        for (i,m) in enumerate(h)
+            h_sum[i] += m
+        end
+    end
+
+    return h_sum
+end
+
+
+# https://iopscience.iop.org/article/10.1088/2399-6528/aab291
+function _magnus_generator(H::Vector, order::Int)
+    @assert(order >= 1)
+
+    Ω_list = [_hamiltonian_integrate(H)]
+    S_list = Dict{Tuple{Int64,Int64},Any}()
+
+    for n in 2:order
+        S_list[(1,n)] = _hamiltonian_commutator(Ω_list[n-1], H)
+        for j in 2:n-1
+            S_list[(j,n)] = _hamiltonian_sum([
+                _hamiltonian_commutator(Ω_list[m], S_list[(j-1,n-m)])
+            for m in 1:n-j])
+        end
+
+        Ω_n = _hamiltonian_sum([
+            _hamiltonian_scalar(_bernoulli_factorial(j),_hamiltonian_integrate(S_list[(j,n)]))
+        for j in 1:n-1])
+
+        push!(Ω_list, Ω_n)
+    end
+
+    return Ω_list
+end
 
-Arguments:
-ising_model - ising model represented as a dictionary.  The qubits
-              and couplings are represented as tuples, and the weights
-              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
-steps - number of iterations for the Magnus Expansion
 
-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(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule, steps::Int; initial_state=nothing, constant_field_x=nothing, constant_field_z=nothing, state_steps=nothing)
+an any-order magnus expansion solver with a fixed number of time steps
+"""
+function simulate_flexible_order(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule, steps::Int, order::Int; initial_state=nothing, constant_field_x=nothing, constant_field_z=nothing, state_steps=nothing)
+    @warn("this any-order magnus expansion solver is not optimized, runtime overheads for high orders are significant", maxlog=1)
     if steps < 2
         error("at least two steps are required by simulate, given $(steps)")
     end
+    if order > 10
+        @warn("magnus expansion orders above 10 can produce numerical stability issues, given $(order)", maxlog=1)
+    end
 
     n = _check_ising_model_ids(ising_model)
 
@@ -180,53 +422,50 @@ function simulate(ising_model::Dict, annealing_time::Real, annealing_schedule::A
 
     R0 = initial_state * initial_state'
 
-    ηs = ones(n)
-    hs = zeros(n)
-
-    x_component = sum_x(n)
+    x_component = _sum_X(n)
     z_component = SparseArrays.spzeros(2^n, 2^n)
     for (tup,w) in ising_model
-        z_component = z_component + sum_z_tup(n, tup, w)
+        z_component = z_component + _kron_Z(n, tup, w)
     end
 
-    xz_bracket = lie_bracket(x_component, z_component)
-
-    constant_component = sum_x(n, constant_field_x) + sum_z(n, constant_field_z)
-    constant_bracket_x = lie_bracket(x_component, constant_component)
-    constant_bracket_z = lie_bracket(z_component, constant_component)
-
     s_steps = range(0, 1, length=steps)
     R_current = R0
-    U = foldl(kron, [IMAT for i = 1:n])
+    U = foldl(kron, [_IMAT for i = 1:n])
 
     if track_states
         push!(state_steps, R_current)
     end
 
+    # explore use of https://github.com/JuliaSymbolics/Symbolics.jl
     for i in 1:(steps-1)
         s0 = s_steps[i]
         s1 = s_steps[i+1]
         δs = s1 - s0
 
-        a_2, a_1, a_0 = get_function_coefficients(annealing_schedule.A, s0, s1)
-        b_2, b_1, b_0 = get_function_coefficients(annealing_schedule.B, s0, s1)
+        aqc = _get_quadratic_coefficients(annealing_schedule.A, s0, s1)
+        bqc = _get_quadratic_coefficients(annealing_schedule.B, s0, s1)
 
-        integA = integral_1_sched(a_2, a_1, a_0, s0, δs)
-        integB = integral_1_sched(b_2, b_1, b_0, s0, δs)
-        integ2 = integral_2_sched(a_2, a_1, a_0, b_2, b_1, b_0, s0, δs)
+        aqc = _shift_quadratic_coefficients(s0, aqc...)
+        bqc = _shift_quadratic_coefficients(s0, bqc...)
 
-        integ2A = a_1*δs^3/6 + a_2*s0*δs^3/3 + a_2*δs^4/6
-        integ2B = b_1*δs^3/6 + b_2*s0*δs^3/3 + b_2*δs^4/6
+        constant_x = _sum_X(n, constant_field_x)
+        constant_z = _sum_Z(n, constant_field_z)
 
-        Ω1Sched = integA * x_component + integB * z_component
-        Ω1Const = δs * constant_component
-        Ω1 = Ω1Sched + Ω1Const
+        H = -im*annealing_time*[
+            aqc[1] * x_component + bqc[1] * z_component + constant_x + constant_z,
+            aqc[2] * x_component + bqc[2] * z_component,
+            aqc[3] * x_component + bqc[3] * z_component,
+        ]
 
-        Ω2Sched = integ2 * xz_bracket 
-        Ω2Const = integ2A * constant_bracket_x + integ2B * constant_bracket_z
-        Ω2 = (Ω2Sched + Ω2Const)/2
+        #Ω_list = _magnus_generator([Matrix(h) for h in H], order)
+        Ω_list = _magnus_generator(H, order)
+        #for (i,Ωi) in enumerate(Ω_list)
+        #   println("Ω_$(i)")
+        #   display(Matrix(_hamiltonian_eval2(δs, Ωi)))
+        #end
+        Ω = sum(_hamiltonian_eval(δs, Ωi) for Ωi in Ω_list)
 
-        U_next = exp(Matrix(-im * (annealing_time*Ω1 + (annealing_time^2)*Ω2)))
+        U_next = exp(Matrix(Ω))
         U = U_next * U
 
         if track_states
@@ -241,22 +480,40 @@ end
 
 
 """
-A simplified interface to the core `simulate` routine that determines a
-suitable number of steps to ensure a high accuracy simulation.
+Main function for performing quantum annealing simulation via a Magnus Expansion (second order).
+Noise can be simulated by running multiple times with randomized constant fields.
+
+Arguments:
+ising_model - ising model represented as a dictionary.  The qubits
+              and couplings are represented as tuples, and the weights
+              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
+
+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)
 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(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule; steps=2, mean_tol=1e-6, max_tol=1e-4, iteration_limit=100, silence=false, state_steps=nothing, kwargs...)
+function simulate(ising_model::Dict, annealing_time::Real, annealing_schedule::AnnealingSchedule; steps=2, order=4, mean_tol=1e-6, max_tol=1e-4, iteration_limit=100, silence=false, state_steps=nothing, kwargs...)
     start_time = time()
     mean_delta = mean_tol + 1.0
     max_delta = max_tol + 1.0
 
+    simulate_method = simulate_fixed_order
+    if order > 4 || haskey(kwargs, :constant_field_x) || haskey(kwargs, :constant_field_z)
+        simulate_method = simulate_flexible_order
+    end
+
     if !silence
         println()
+        println("order: $(order)  method: $(simulate_method)")
         println("iter |  steps  |    max(Δ)    |    mean(Δ)   |")
     end
 
-    ρ_prev = simulate(ising_model, annealing_time, annealing_schedule, steps; kwargs...)
+    ρ_prev = simulate_method(ising_model, annealing_time, annealing_schedule, steps, order; kwargs...)
 
     iteration = 1
     while mean_delta >= mean_tol || max_delta >= max_tol
@@ -266,7 +523,7 @@ function simulate(ising_model::Dict, annealing_time::Real, annealing_schedule::A
             empty!(state_steps)
         end
 
-        ρ = simulate(ising_model, annealing_time, annealing_schedule, steps; state_steps=state_steps, kwargs...)
+        ρ = simulate_method(ising_model, annealing_time, annealing_schedule, steps, order; state_steps=state_steps, kwargs...)
 
         ρ_delta = abs.(ρ .- ρ_prev)
         mean_delta = sum(ρ_delta)/length(ρ_delta)
@@ -295,4 +552,3 @@ function simulate(ising_model::Dict, annealing_time::Real, annealing_schedule::A
     return ρ_prev
 end
 
-
diff --git a/src/simulate_de.jl b/src/simulate_de.jl
index f027dec..590e8d9 100644
--- a/src/simulate_de.jl
+++ b/src/simulate_de.jl
@@ -33,10 +33,10 @@ function simulate_de(ising_model, annealing_time, annealing_schedule, reltol; ab
         constant_field_z = zeros(n)
     end
 
-    const_x_component = sum_x(n, constant_field_x)
-    const_z_component = sum_z(n, constant_field_z)
+    const_x_component = _sum_X(n, constant_field_x)
+    const_z_component = _sum_Z(n, constant_field_z)
 
-    H(s) = transverse_ising_hamiltonian(ising_model, annealing_schedule, s) + const_x_component + const_z_component
+    H(s) = hamiltonian_transverse_ising(ising_model, annealing_schedule, s) + const_x_component + const_z_component
     schrod_eq(state, time, s) = -im * time * H(s) * state
 
     s_range = (0.0, 1.0)
diff --git a/test/base.jl b/test/base.jl
index da4a1c0..30831b7 100644
--- a/test/base.jl
+++ b/test/base.jl
@@ -130,19 +130,19 @@ end
 @testset "transverse ising hamiltonian" begin
 
     @testset "1 qubit, analytical solution" begin
-        @test all(isapprox(one_spin_H(s), transverse_ising_hamiltonian(one_spin_model, AS_CIRCULAR, s)) for s in s_100)
+        @test all(isapprox(one_spin_H(s), hamiltonian_transverse_ising(one_spin_model, AS_CIRCULAR, s)) for s in s_100)
     end
 
     @testset "2 qubit, analytical solution" begin
-        @test all(isapprox(two_spin_H(s), transverse_ising_hamiltonian(two_spin_model, AS_CIRCULAR, s)) for s in s_100)
+        @test all(isapprox(two_spin_H(s), hamiltonian_transverse_ising(two_spin_model, AS_CIRCULAR, s)) for s in s_100)
     end
 
     @testset "1 qubit, linear schedule" begin
         annealing_schedule = AS_LINEAR
 
-        H_00 = transverse_ising_hamiltonian(one_spin_model, annealing_schedule, 0.0)
-        H_05 = transverse_ising_hamiltonian(one_spin_model, annealing_schedule, 0.5)
-        H_10 = transverse_ising_hamiltonian(one_spin_model, annealing_schedule, 1.0)
+        H_00 = hamiltonian_transverse_ising(one_spin_model, annealing_schedule, 0.0)
+        H_05 = hamiltonian_transverse_ising(one_spin_model, annealing_schedule, 0.5)
+        H_10 = hamiltonian_transverse_ising(one_spin_model, annealing_schedule, 1.0)
 
         @test isapprox(H_00, [0 1; 1 0])
         @test isapprox(H_05, [0.5 0.5; 0.5 -0.5])
@@ -152,9 +152,9 @@ end
     @testset "2 qubit, linear schedule" begin
         annealing_schedule = AS_LINEAR
 
-        H_00 = transverse_ising_hamiltonian(two_spin_model, annealing_schedule, 0.0)
-        H_05 = transverse_ising_hamiltonian(two_spin_model, annealing_schedule, 0.5)
-        H_10 = transverse_ising_hamiltonian(two_spin_model, annealing_schedule, 1.0)
+        H_00 = hamiltonian_transverse_ising(two_spin_model, annealing_schedule, 0.0)
+        H_05 = hamiltonian_transverse_ising(two_spin_model, annealing_schedule, 0.5)
+        H_10 = hamiltonian_transverse_ising(two_spin_model, annealing_schedule, 1.0)
 
         @test isapprox(H_00, [0.0 1.0 1.0 0.0; 1.0  0.0 0.0 1.0; 1.0 0.0  0.0 1.0; 0.0 1.0 1.0 0.0])
         @test isapprox(H_05, [1.0 0.5 0.5 0.0; 0.5 -1.0 0.0 0.5; 0.5 0.0 -1.0 0.5; 0.0 0.5 0.5 1.0])
@@ -198,14 +198,14 @@ end
 
         ss = 0.0:0.001:1.0
         asch = [(0.0,0.0), (1.0,1.0)]
-        mod_asch = dwave_annealing_protocol(annealing_schedule, asch=asch)
+        mod_asch = annealing_protocol_dwave(annealing_schedule, asch=asch)
         @test isapprox(mod_asch.A.(ss), annealing_schedule.A.(ss))
         @test isapprox(mod_asch.B.(ss), annealing_schedule.B.(ss))
 
         ss_reg = [0.0, 0.25, 0.5, 0.5, 0.5, 0.75, 1.0]
         ss_mod = [0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0]
         asch = [(0.0,0.0), (0.4,0.5), (0.6,0.5), (1.0,1.0)]
-        mod_asch = dwave_annealing_protocol(annealing_schedule, asch=asch)
+        mod_asch = annealing_protocol_dwave(annealing_schedule, asch=asch)
         @test isapprox(mod_asch.A.(ss_mod), annealing_schedule.A.(ss_reg))
         @test isapprox(mod_asch.B.(ss_mod), annealing_schedule.B.(ss_reg))
     end
diff --git a/test/common.jl b/test/common.jl
index d6a078d..214332a 100644
--- a/test/common.jl
+++ b/test/common.jl
@@ -62,12 +62,39 @@ end
 s_100 = range(0, 1, length=100)
 s_10000 = range(0, 1, length=10000)
 
-AS_CIRCULAR_pwc_csv_100 = parse_dwave_annealing_schedule("data/trig_sched_100.csv", interpolation=:none, initial_state=default_initial_state)
-AS_CIRCULAR_pwc_csv_1000 = parse_dwave_annealing_schedule("data/trig_sched_1000.csv", interpolation=:none, initial_state=default_initial_state)
+AS_CIRCULAR_pwc_csv_100 = read_dwave_annealing_schedule("data/trig_sched_100.csv", interpolation=:none, initial_state=initial_state_default)
+AS_CIRCULAR_pwc_csv_1000 = read_dwave_annealing_schedule("data/trig_sched_1000.csv", interpolation=:none, initial_state=initial_state_default)
 
-AS_CIRCULAR_pwl_csv_100 = parse_dwave_annealing_schedule("data/trig_sched_100.csv", initial_state=default_initial_state)
-AS_CIRCULAR_pwl_csv_1000 = parse_dwave_annealing_schedule("data/trig_sched_1000.csv", initial_state=default_initial_state)
+AS_CIRCULAR_pwl_csv_100 = read_dwave_annealing_schedule("data/trig_sched_100.csv", initial_state=initial_state_default)
+AS_CIRCULAR_pwl_csv_1000 = read_dwave_annealing_schedule("data/trig_sched_1000.csv", initial_state=initial_state_default)
 
-AS_CIRCULAR_pwq_csv_100 = parse_dwave_annealing_schedule("data/trig_sched_100.csv", interpolation=:quadratic, initial_state=default_initial_state)
-AS_CIRCULAR_pwq_csv_1000 = parse_dwave_annealing_schedule("data/trig_sched_1000.csv", interpolation=:quadratic, initial_state=default_initial_state)
+AS_CIRCULAR_pwq_csv_100 = read_dwave_annealing_schedule("data/trig_sched_100.csv", interpolation=:quadratic, initial_state=initial_state_default)
+AS_CIRCULAR_pwq_csv_1000 = read_dwave_annealing_schedule("data/trig_sched_1000.csv", interpolation=:quadratic, initial_state=initial_state_default)
 
+
+function sum_zizj(n, J::Dict)
+    sum_zs = zeros(2^n)
+    for (i,j) in keys(J)
+        J_val = J[(i,j)]
+        sum_zs += zizj_vectorized(n, i, j, J_val)
+    end
+    return diagm(sum_zs)
+end
+
+function zizj_vectorized(n::Int, i::Int, j::Int, J_val)
+    Ivec = [1.0+0im,1.0+0im]
+    Zvec = [1.0+0im,-1.0+0im]
+    matvec = [(k == i || k == j) ? Zvec : Ivec for k in n:-1:1]
+    return J_val * foldl(kron, matvec)
+end
+
+function ising_hamiltonian(ising_model)
+    n = QuantumAnnealing._check_ising_model_ids(ising_model)
+
+    z_component = zeros(2^n, 2^n)
+    for (tup,w) in ising_model
+        z_component += QuantumAnnealing._kron_Z(n, tup, w)
+    end
+
+    return z_component
+end
diff --git a/test/dwave.jl b/test/dwave.jl
new file mode 100644
index 0000000..638d452
--- /dev/null
+++ b/test/dwave.jl
@@ -0,0 +1,110 @@
+
+@testset "simulate bqpjson" begin
+
+    @testset "1 qubit, function schedule" begin
+        annealing_time = 1000.0
+        annealing_schedule = AS_CIRCULAR
+        #steps = 2500
+
+        ising_intended = Dict((1,) => 1)
+
+        bqpjson_file = "data/bqpjson_1q.json"
+        dwisc_file = "tmp.json"
+
+        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, silence=true)
+
+        dwisc_data = JSON.parsefile(dwisc_file)
+        rm(dwisc_file)
+
+        ρ = simulate(ising_intended, annealing_time, annealing_schedule, silence=true)
+
+        @test dwisc_data["solutions"][1]["prob"] >= 0.99
+        @test isapprox(z_measure_probabilities(ρ)[2], dwisc_data["solutions"][1]["prob"])
+        @test dwisc_data["variable_ids"] == [100]
+    end
+
+    @testset "1 qubit, dwave schedule" begin
+        annealing_time = 1000.0
+        annealing_schedule = AS_DW_QUADRATIC
+
+        ising_intended = Dict((1,) => 1)
+
+        bqpjson_file = "data/bqpjson_1q.json"
+        dwisc_file = "tmp.json"
+
+        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, silence=true)
+
+        dwisc_data = JSON.parsefile(dwisc_file)
+        rm(dwisc_file)
+
+        ρ = simulate(ising_intended, annealing_time, annealing_schedule, silence=true)
+
+        @test dwisc_data["solutions"][1]["prob"] >= 0.99
+        @test isapprox(z_measure_probabilities(ρ)[2], dwisc_data["solutions"][1]["prob"])
+        @test dwisc_data["variable_ids"] == [100]
+    end
+
+    @testset "2 qubit, function schedule" begin
+        annealing_time = 1000.0
+        annealing_schedule = AS_CIRCULAR
+
+        # the ising model that is encoded in bqpjson_2q.json
+        ising_intended = Dict((1,) => -1, (2,) => -1, (1,2) => -1)
+
+        bqpjson_file = "data/bqpjson_2q.json"
+        dwisc_file = "tmp.json"
+
+        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, silence=true)
+
+        dwisc_data = JSON.parsefile(dwisc_file)
+        rm(dwisc_file)
+
+        ρ = simulate(ising_intended, annealing_time, annealing_schedule, silence=true)
+
+        @test dwisc_data["solutions"][1]["prob"] >= 0.99
+        @test isapprox(z_measure_probabilities(ρ)[1], dwisc_data["solutions"][1]["prob"])
+        @test dwisc_data["variable_ids"] == [304, 308]
+    end
+
+    @testset "1 qubit, function schedule, z noise" begin
+        annealing_time = 100.0
+        annealing_schedule = AS_CIRCULAR
+        numshots = 10
+
+        # random numbers between -0.75 and 0.75
+        z_bias = (Random.rand(numshots) .- 0.5)*1.5
+
+        bqpjson_file = "data/bqpjson_1q.json"
+        dwisc_file = "tmp.json"
+
+        simulate_bqpjson_noisy(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, z_bias=z_bias, silence=true)
+
+        dwisc_data = JSON.parsefile(dwisc_file)
+        rm(dwisc_file)
+
+        @test dwisc_data["solutions"][1]["prob"] > 0.70 # true w.h.p.
+        @test dwisc_data["variable_ids"] == [100]
+    end
+
+    @testset "1 qubit, function schedule, x and z noise" begin
+        annealing_time = 100.0
+        annealing_schedule = AS_CIRCULAR
+        numshots = 10
+
+        # random numbers between -0.75 and 0.75
+        x_bias = (Random.rand(numshots) .- 0.5)*1.5
+        z_bias = (Random.rand(numshots) .- 0.5)*1.5
+
+        bqpjson_file = "data/bqpjson_1q.json"
+        dwisc_file = "tmp.json"
+
+        simulate_bqpjson_noisy(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, x_bias=x_bias, z_bias=z_bias, silence=true)
+
+        dwisc_data = JSON.parsefile(dwisc_file)
+        rm(dwisc_file)
+
+        @test dwisc_data["solutions"][1]["prob"] > 0.70 # true w.h.p.
+        @test dwisc_data["variable_ids"] == [100]
+    end
+
+end
diff --git a/test/magnus.jl b/test/magnus.jl
new file mode 100644
index 0000000..c592ae4
--- /dev/null
+++ b/test/magnus.jl
@@ -0,0 +1,67 @@
+@testset "magnus expansion method" begin
+
+    @testset "magnus generator subroutines" begin
+        mc = QuantumAnnealing._matrix_commutator([1 2; 3 4], [5 6; 7 8])
+        @test isapprox([-4 -12; 12 4], mc)
+
+        he = QuantumAnnealing._hamiltonian_eval(2, [[1],[2],[3],[4]])
+        @test isapprox([49], he)
+
+        hs = QuantumAnnealing._hamiltonian_scalar(2, [[1],[2],[3],[4]])
+        @test isapprox([[2],[4],[6],[8]], hs)
+
+        hi = QuantumAnnealing._hamiltonian_integrate([[1],[2],[3],[4]])
+        @test isapprox([[0], [1], [1], [1], [1]], hi)
+
+        hs = QuantumAnnealing._hamiltonian_sum([[1,2,3,4],[5,6]])
+        @test isapprox([6,8,3,4], hs)
+
+        hc = QuantumAnnealing._hamiltonian_commutator([[1 2; 3 4],[1 2; 3 4],[1 2; 3 4],[1 2; 3 4]],[[5 6; 7 8],[9 10; 11 12]])
+        @test isapprox([[-4 -12; 12 4], [-12 -36; 36 12], [-12 -36; 36 12], [-12 -36; 36 12], [-8 -24; 24 8]], hc)
+
+        bernoulli_fact_numbers = [-0.5, 0.08333333333333333, 0.0, -0.001388888888888889, 0.0, 3.306878306878307e-5]
+        for (i,v) in enumerate(bernoulli_fact_numbers)
+            @test isapprox(QuantumAnnealing._bernoulli_factorial(i), v)
+        end
+    end
+
+    @testset "Ω computations, orders 1 to 4" begin
+        ising_model = Dict((1,2) => -1, (1,) => 1/2, (2,) => -5/7)
+        n = 2
+        annealing_time = 2.0
+        annealing_schedule = AS_LINEAR
+        order = 4
+
+        s0 = 0.0
+        s1 = 1.0
+        δs = s1 - s0
+
+        x_component = QuantumAnnealing._sum_X(n)
+        z_component = zeros(2^n, 2^n)
+        for (tup,w) in ising_model
+            z_component += QuantumAnnealing._kron_Z(n, tup, w)
+        end
+
+        H_parts = QuantumAnnealing._H_parts(x_component, z_component, order)
+        Ω_list1 = QuantumAnnealing._Ω_list(annealing_time, s0, s1, annealing_schedule, H_parts, order)
+
+        aqc = QuantumAnnealing._get_quadratic_coefficients(annealing_schedule.A, s0, s1)
+        bqc = QuantumAnnealing._get_quadratic_coefficients(annealing_schedule.B, s0, s1)
+
+        aqc = QuantumAnnealing._shift_quadratic_coefficients(s0, aqc...)
+        bqc = QuantumAnnealing._shift_quadratic_coefficients(s0, bqc...)
+
+        H = -im*annealing_time*[
+            aqc[1] * x_component + bqc[1] * z_component,
+            aqc[2] * x_component + bqc[2] * z_component,
+            aqc[3] * x_component + bqc[3] * z_component,
+        ]
+
+        Ω_list_tmp = QuantumAnnealing._magnus_generator(H, order)
+        Ω_list2 = [QuantumAnnealing._hamiltonian_eval(δs, Ωi) for Ωi in Ω_list_tmp]
+
+        for i in 1:order
+            @test isapprox(Ω_list1[i], Ω_list2[i])
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3327aa8..b2b98dd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,6 +11,10 @@ include("common.jl")
 
     include("base.jl")
 
+    include("magnus.jl")
+
     include("simulate.jl")
 
+    include("dwave.jl")
+
 end
diff --git a/test/simulate.jl b/test/simulate.jl
index 19fc055..6adf88a 100644
--- a/test/simulate.jl
+++ b/test/simulate.jl
@@ -1,8 +1,7 @@
 
 @testset "simulate, 1 qubit" begin
-
     @testset "function schedule, default anneal time, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR, 100, 4)
         @test isapprox(one_spin_ρ(1.0), ρ)
     end
 
@@ -12,18 +11,18 @@
     end
 
     @testset "function schedule, fast anneal time, analytical solution" begin
-        ρ = simulate(one_spin_model, 0.5, AS_CIRCULAR, 1000)
+        ρ = simulate_fixed_order(one_spin_model, 0.5, AS_CIRCULAR, 1000, 4)
         @test isapprox(one_spin_ρ(0.5), ρ)
     end
 
     @testset "function schedule, slow anneal time, analytical solution" begin
-        ρ = simulate(one_spin_model, 2.0, AS_CIRCULAR, 1000)
+        ρ = simulate_fixed_order(one_spin_model, 2.0, AS_CIRCULAR, 1000, 4)
         @test isapprox(one_spin_ρ(2.0), ρ)
     end
 
     @testset "fractional field value" begin
         ρ_target = [0.420186+0.0im -0.409634+0.275372im; -0.409634-0.275372im 0.579814+2.77556e-17im]
-        ρ = simulate(Dict((1,) => 0.5), 1.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(Dict((1,) => 0.5), 1.0, AS_CIRCULAR, 100, 4)
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -31,7 +30,7 @@
 
     @testset "field value above 1.0" begin
         ρ_target = [0.291065-2.77556e-17im 0.114524+0.43958im; 0.114524-0.43958im 0.708935+5.55112e-17im]
-        ρ = simulate(Dict((1,) => 1.5), 1.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(Dict((1,) => 1.5), 1.0, AS_CIRCULAR, 100, 4)
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -39,7 +38,7 @@
 
     @testset "function schedule, constant terms" begin
         ρ_target = [0.0578906+1.38778e-17im -0.165069-0.165202im; -0.165069+0.165202im 0.942109+0.0im]
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, 100, constant_field_x = [1], constant_field_z = [1])
+        ρ = simulate_flexible_order(one_spin_model, 1.0, AS_CIRCULAR, 100, 4, constant_field_x = [1], constant_field_z = [1])
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -47,7 +46,7 @@
 
     @testset "function schedule, different initial state" begin
         ρ_target = [0.326006+0.0im -0.413095-0.221536im; -0.413095+0.221536im 0.673994+0.0im]
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, 100, initial_state = [0,1])
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR, 100, 4, initial_state = [0,1])
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -55,7 +54,7 @@
 
     @testset "function schedule (AS_LINEAR), default anneal time" begin
         ρ_target = [0.422382+2.77556e-17im -0.278818+0.40772im; -0.278818-0.40772im 0.577618+2.77556e-17im]
-        ρ = simulate(one_spin_model, 1.0, AS_LINEAR, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_LINEAR, 100, 4)
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -63,7 +62,7 @@
 
     @testset "function schedule (AS_QUADRATIC), default anneal time" begin
         ρ_target = [0.489037+0.0im -0.393381+0.308433im; -0.393381-0.308433im 0.510963+5.55112e-17im]
-        ρ = simulate(one_spin_model, 1.0, AS_QUADRATIC, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_QUADRATIC, 100, 4)
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -71,7 +70,7 @@
 
     @testset "function schedule (AS_DW_QUADRATIC), default anneal time" begin
         ρ_target = [0.0162536+0.0im 0.121897-0.0336245im; 0.121897+0.0336245im  0.983746+2.77556e-17im]
-        ρ = simulate(one_spin_model, 1.0, AS_DW_QUADRATIC, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_DW_QUADRATIC, 100, 4)
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -79,7 +78,7 @@
 
 
     @testset "function schedule, too few steps, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, 10)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR, 10, 4)
 
         # NOTE, atol required due to too few iterations
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-5)
@@ -87,31 +86,31 @@
     end
 
     @testset "csv schedule pwq, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR_pwq_csv_1000, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR_pwq_csv_1000, 100, 4)
 
         @test isapprox(one_spin_ρ(1.0), ρ)
     end
 
     @testset "csv schedule pwq, low resolution, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR_pwq_csv_100, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR_pwq_csv_100, 100, 4)
 
-        # NOTE, atol required due to pwq approximation in the schedule file
+        # NOTE, atol required due to pwq approx_IMATion in the schedule file
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-6)
         @test !isapprox(one_spin_ρ(1.0), ρ, atol=1e-7)
     end
 
     @testset "csv schedule pwl, low resolution, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR_pwl_csv_100, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR_pwl_csv_100, 100, 4)
 
-        # NOTE, atol required due to pwl approximation in the schedule file
+        # NOTE, atol required due to pwl approx_IMATion in the schedule file
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-4)
         @test !isapprox(one_spin_ρ(1.0), ρ, atol=1e-5)
     end
 
     @testset "csv schedule pwc, low resolution, analytical solution" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR_pwc_csv_100, 100)
+        ρ = simulate_fixed_order(one_spin_model, 1.0, AS_CIRCULAR_pwc_csv_100, 100, 4)
 
-        # NOTE, atol required due to pwc approximation in the schedule file
+        # NOTE, atol required due to pwc approx_IMATion in the schedule file
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-2)
         @test !isapprox(one_spin_ρ(1.0), ρ, atol=1e-3)
     end
@@ -120,9 +119,9 @@
         ρ_target = [0.88389+0.0im -0.197484-0.252247im; -0.197484+0.252247im 0.11611+0.0im]
 
         asch = [(0.0, 1.0), (0.5, 0.5), (1.0, 1.0)]
-        mod_asch = dwave_annealing_protocol(AS_LINEAR, asch=asch)
+        mod_asch = annealing_protocol_dwave(AS_LINEAR, asch=asch)
 
-        ρ = simulate(Dict((1,) => 0.1), 5.0, mod_asch, 100, initial_state = [0,1])
+        ρ = simulate_fixed_order(Dict((1,) => 0.1), 5.0, mod_asch, 100, 4, initial_state = [0,1])
 
         # NOTE, atol required due to too few digits in target
         @test isapprox(ρ_target, ρ, atol=1e-6)
@@ -131,11 +130,11 @@
     @testset "collect probability trajectory, nonadaptive" begin
         ρ_list = []
         steps=100
-        ρ = simulate(one_spin_model, 1, AS_CIRCULAR, steps, state_steps=ρ_list)
+        ρ = simulate_fixed_order(one_spin_model, 1, AS_CIRCULAR, steps, 4, state_steps=ρ_list)
 
         @test isapprox(one_spin_ρ(1.0), ρ)
         @test length(ρ_list) == steps
-        @test isapprox(ρ_list[1], (default_initial_state(1) * default_initial_state(1)'))
+        @test isapprox(ρ_list[1], (initial_state_default(1) * initial_state_default(1)'))
         @test isapprox(ρ_list[steps], one_spin_ρ(1.0))
     end
 
@@ -145,42 +144,80 @@
 
         @test isapprox(one_spin_ρ(1.0), ρ)
         @test length(ρ_list) == 64
-        @test isapprox(ρ_list[1], (default_initial_state(1) * default_initial_state(1)'))
+        @test isapprox(ρ_list[1], (initial_state_default(1) * initial_state_default(1)'))
         @test isapprox(ρ_list[64], one_spin_ρ(1.0))
     end
-
 end
 
 
 @testset "simulate, 2 qubit" begin
-
     @testset "function schedule, default anneal time, analytical solution" begin
-        ρ = simulate(two_spin_model, 1.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(two_spin_model, 1.0, AS_CIRCULAR, 100, 4)
         @test isapprox(two_spin_ρ(1.0), ρ)
     end
 
     @testset "function schedule, fast anneal time, analytical solution" begin
-        ρ = simulate(two_spin_model, 0.5, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(two_spin_model, 0.5, AS_CIRCULAR, 100, 4)
         @test isapprox(two_spin_ρ(0.5), ρ)
     end
 
     @testset "function schedule, slow anneal time, analytical solution" begin
-        ρ = simulate(two_spin_model, 2.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(two_spin_model, 2.0, AS_CIRCULAR, 100, 4)
         @test isapprox(two_spin_ρ(2.0), ρ)
     end
+end
+
+
+@testset "simulate, any-order magnus expansion" begin
+    @testset "1 qubit, adaptive, orders 1 to 8" begin
+        for i in 1:8
+            ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, order=i, max_tol=1e-9, silence=true)
+            @test isapprox(one_spin_ρ(1.0), ρ)
+        end
+    end
+
+    @testset "2 qubit, adaptive, orders 1 to 8" begin
+        for i in 1:8
+            ρ = simulate(two_spin_model, 1.0, AS_CIRCULAR, order=i, max_tol=1e-9, silence=true)
+            @test isapprox(two_spin_ρ(1.0), ρ)
+        end
+    end
 
+    @testset "5 qubit, fixed_order solver, from 1 to 4" begin
+        ising_model = Dict((1,) => -1, (1,2) => -1, (1,3) => -1, (1,4) => -1, (1,5) => -1, (2,3) => 1, (4,5) => 1)
+
+        for i in 1:4
+            ρ_target = simulate_fixed_order(ising_model, 2.0, AS_CIRCULAR, 2, i)
+            ρ = simulate_flexible_order(ising_model, 2.0, AS_CIRCULAR, 2, i)
+            @test isapprox(ρ_target, ρ)
+        end
+
+    end
+
+    @testset "5 qubit, hardcoded forth order solver, function schedules" begin
+        ising_model = Dict((1,) => -1, (1,2) => -1, (1,3) => -1, (1,4) => -1, (1,5) => -1, (2,3) => 1, (4,5) => 1)
+
+        ρ_target = simulate_fixed_order(ising_model, 2.0, AS_LINEAR, 2, 4)
+        ρ = simulate_flexible_order(ising_model, 2.0, AS_LINEAR, 2, 4)
+        @test isapprox(ρ_target, ρ)
+
+        ρ_target = simulate_fixed_order(ising_model, 2.0, AS_QUADRATIC, 2, 4)
+        ρ = simulate_flexible_order(ising_model, 2.0, AS_QUADRATIC, 2, 4)
+        @test isapprox(ρ_target, ρ)
+
+        ρ_target = simulate_fixed_order(ising_model, 2.0, AS_DW_QUADRATIC, 2, 4)
+        ρ = simulate_flexible_order(ising_model, 2.0, AS_DW_QUADRATIC, 2, 4)
+        @test isapprox(ρ_target, ρ)
+    end
 end
 
 
 @testset "simulate, multi-qubit" begin
 
     @testset "2 qubit, function schedules (AS_CIRCULAR, AS_LINEAR, AS_QUADRATIC), near adiabatic limit" begin
-        n = 2
-        h = ones(n)
-        J = Dict((1,2) => -1)
         ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
@@ -188,19 +225,19 @@ end
         annealing_time = 100.0
         steps = 100
 
-        ρ = simulate(ising_model, annealing_time, AS_CIRCULAR, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_CIRCULAR, steps, 4)
         @test real(ρ[4,4]) >= 0.9999
 
         @test isapprox(min_ρ, ρ, atol=1e-2)
         @test !isapprox(min_ρ, ρ, atol=1e-3)
 
-        ρ = simulate(ising_model, annealing_time, AS_LINEAR, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_LINEAR, steps, 4)
         @test real(ρ[4,4]) >= 0.9999
 
         @test isapprox(min_ρ, ρ, atol=1e-2)
         @test !isapprox(min_ρ, ρ, atol=1e-3)
 
-        ρ = simulate(ising_model, annealing_time, AS_QUADRATIC, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_QUADRATIC, steps, 4)
         @test real(ρ[4,4]) >= 0.9999
 
         @test isapprox(min_ρ, ρ, atol=1e-3)
@@ -208,12 +245,9 @@ end
     end
 
     @testset "2 qubit, function schedules (AS_DW_QUADRATIC), near adiabatic limit" begin
-        n = 2
-        h = ones(n)
-        J = Dict((1,2) => -1)
         ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
@@ -221,7 +255,7 @@ end
         annealing_time = 10.0
         steps = 1000
 
-        ρ = simulate(ising_model, annealing_time, AS_DW_QUADRATIC, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_DW_QUADRATIC, steps, 4)
         @test real(ρ[4,4]) >= 0.9999
 
         @test isapprox(min_ρ, ρ, atol=1e-3)
@@ -229,17 +263,14 @@ end
     end
 
     @testset "2 qubit, csv schedule, near adiabatic limit" begin
-        n = 2
-        h = ones(2)
-        J = Dict((1,2) => -1)
         ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
 
-        ρ = simulate(ising_model, 100.0, AS_CIRCULAR_pwl_csv_1000, 100)
+        ρ = simulate_fixed_order(ising_model, 100.0, AS_CIRCULAR_pwl_csv_1000, 100, 4)
 
         @test real(ρ[4,4]) >= 0.9999
 
@@ -248,17 +279,14 @@ end
     end
 
     @testset "2 qubit, function schedule, fractional values, near adiabatic limit" begin
-        n = 2
-        h = ones(2)
-        J = Dict((1,2) => -0.76)
         ising_model = Dict((1,) => 0.5, (1,2) => -0.75)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
 
-        ρ = simulate(ising_model, 50.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(ising_model, 50.0, AS_CIRCULAR, 100, 4)
 
         # NOTE, Non-1.0 due to annealing time not being in the fill adiabatic limit
         @test real(tr(ρ*min_ρ)) >= 0.999
@@ -268,73 +296,68 @@ end
     end
 
     @testset "2 qubit probability trajectory, nonadaptive" begin
-        n = 2
-        h = ones(n)
-        J = Dict((1,2) => -1)
         ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
 
+        n = QuantumAnnealing._check_ising_model_ids(ising_model)
+
         annealing_time = 10.0
         steps = 1000
 
         ρ_list = []
-        ρ = simulate(ising_model, annealing_time, AS_DW_QUADRATIC, steps, state_steps=ρ_list)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_DW_QUADRATIC, steps, 4, state_steps=ρ_list)
         @test real(ρ[4,4]) >= 0.9999
         @test length(ρ_list) == steps
-        @test isapprox(ρ_list[1], (default_dwave_initial_state(n) * default_dwave_initial_state(n)'))
+        @test isapprox(ρ_list[1], (initial_state_default_dwave(n) * initial_state_default_dwave(n)'))
         @test real(ρ_list[steps][4,4]) >= 0.9999
         @test isapprox(min_ρ, ρ, atol=1e-3)
         @test !isapprox(min_ρ, ρ, atol=1e-4)
     end
 
     @testset "2 qubit probability trajectory, adaptive" begin
-        n = 2
-        h = ones(n)
-        J = Dict((1,2) => -1)
         ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
         min_vec=evecs[:,1]
         min_ρ = min_vec * min_vec'
 
+        n = QuantumAnnealing._check_ising_model_ids(ising_model)
+
         annealing_time = 10.0
 
         ρ_list = []
         ρ = simulate(ising_model, annealing_time, AS_DW_QUADRATIC, silence=true, state_steps=ρ_list)
         @test real(ρ[4,4]) >= 0.9999
-        @test length(ρ_list) == 1024
-        @test isapprox(ρ_list[1], (default_dwave_initial_state(n) * default_dwave_initial_state(n)'))
-        @test real(ρ_list[1024][4,4]) >= 0.9999
+        @test length(ρ_list) == 512
+        @test isapprox(ρ_list[1], (initial_state_default_dwave(n) * initial_state_default_dwave(n)'))
+        @test real(ρ_list[end][4,4]) >= 0.9999
         @test isapprox(min_ρ, ρ, atol=1e-3)
         @test !isapprox(min_ρ, ρ, atol=1e-4)
     end
 
     @testset "3 qubit, degenerate, function schedules (AS_CIRCULAR, AS_QUADRATIC), near adiabatic limit" begin
         # ring of disagrees => 6 ground states
-        n = 3
-        h = zeros(n)
-        J = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
-        ising_model = J
+        ising_model = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
 
         annealing_time = 100.0
-        steps = 100
+        steps = 500
 
-        ρ = simulate(ising_model, annealing_time, AS_CIRCULAR, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_CIRCULAR, steps, 4)
         for i = 1:6
             min_vec = evecs[:,i]
             min_ρ = min_vec * min_vec'
-            @test isapprox(real(tr(ρ*min_ρ)), 1.0/6.0, atol=1e-8)
+            @test isapprox(real(tr(ρ*min_ρ)), 1.0/6.0, atol=1e-7)
         end
 
-        ρ = simulate(ising_model, annealing_time, AS_QUADRATIC, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_QUADRATIC, steps, 4)
         for i = 1:6
             min_vec = evecs[:,i]
             min_ρ = min_vec * min_vec'
@@ -344,43 +367,34 @@ end
 
     @testset "3 qubit, degenerate, function schedules (AS_DW_QUADRATIC), near adiabatic limit" begin
         # ring of disagrees => 6 ground states
-        n = 3
-        h = zeros(n)
-        J = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
-        ising_model = J
+        ising_model = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
 
-        H = sum_z(n,h) + sum_zizj(n,J)
+        H = ising_hamiltonian(ising_model)
         evals,evecs = eigen(Matrix(H))
 
         annealing_time = 100.0
-        steps = 100
+        steps = 4000
 
-        ρ = simulate(ising_model, annealing_time, AS_DW_QUADRATIC, steps)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_DW_QUADRATIC, steps, 4)
         for i = 1:6
             min_vec = evecs[:,i]
             min_ρ = min_vec * min_vec'
-            @test isapprox(real(tr(ρ*min_ρ)), 1.0/6.0, atol=1e-3)
-            @test !isapprox(real(tr(ρ*min_ρ)), 1.0/6.0, atol=1e-4)
+            @test isapprox(real(tr(ρ*min_ρ)), 1.0/6.0)
         end
     end
 
     @testset "3 qubit, csv schedule, near adiabatic limit" begin
         #ring of disagrees => 6 ground states
-        n = 3
-        h = zeros(n)
-        J = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
-        ising_model = J
-
+        ising_model = Dict((1,2) => 1, (1,3) => 1, (2,3) => 1)
         annealing_time = 100.0
-        steps = 100
+        steps = 500
 
-        ρ_target = simulate(ising_model, annealing_time, AS_CIRCULAR, steps)
-        ρ = simulate(ising_model, annealing_time, AS_CIRCULAR_pwq_csv_1000, steps)
+        ρ_target = simulate_fixed_order(ising_model, annealing_time, AS_CIRCULAR, steps, 4)
+        ρ = simulate_fixed_order(ising_model, annealing_time, AS_CIRCULAR_pwq_csv_1000, steps, 4)
 
         @test isapprox(ρ_target, ρ, atol=1e-7)
         @test !isapprox(ρ_target, ρ, atol=1e-8)
     end
-
 end
 
 
@@ -391,7 +405,7 @@ end
         redirect_stdout(io)
         redirect_stderr(io)
 
-        ρ = simulate(two_spin_model, 1.0, AS_CIRCULAR, 100)
+        ρ = simulate_fixed_order(two_spin_model, 1.0, AS_CIRCULAR, 100, 4)
         print_z_state_probabilities(ρ)
         print_z_state_probabilities(ρ, sort=true)
         print_z_state_probabilities(ρ, sort=true, limit=2)
@@ -416,11 +430,11 @@ end
     end
 
     @testset "1 qubit, function schedule, constant terms" begin
-        ρ = simulate(one_spin_model, 1.0, AS_CIRCULAR, 100, constant_field_x=[1], constant_field_z=[1])
+        ρ = simulate_flexible_order(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-8)
-        @test !isapprox(ρ, ρ_de, atol=1e-9)
+        @test isapprox(ρ, ρ_de, atol=1e-9)
+        @test !isapprox(ρ, ρ_de, atol=1e-10)
     end
 
     @testset "1 qubit, csv schedule, analytical solution" begin
@@ -432,16 +446,16 @@ 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)
+        ρ = simulate_fixed_order(ising_model, 1.0, AS_CIRCULAR, 100, 4)
         ρ_de = simulate_de(ising_model, 1.0, AS_CIRCULAR, 1e-6)
 
-        @test isapprox(ρ, ρ_de, atol=1e-8)
-        @test !isapprox(ρ, ρ_de, atol=1e-9)
+        @test isapprox(ρ, ρ_de, atol=1e-9)
+        @test !isapprox(ρ, ρ_de, atol=1e-10)
     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)
+        ρ = simulate_fixed_order(ising_model, 1000.0, AS_CIRCULAR, 4000, 4)
         ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-7)
 
         @test isapprox(ρ, ρ_de, atol=1e-6)
@@ -450,127 +464,10 @@ 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)
+        ρ = simulate_fixed_order(ising_model, 100.0, AS_CIRCULAR, 500, 4)
         ρ_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
-
-
-@testset "simulate bqpjson" begin
-
-    @testset "1 qubit, function schedule" begin
-        annealing_time = 10000.0
-        annealing_schedule = AS_CIRCULAR
-        steps = 1000
-
-        ising_intended = Dict((1,) => 1)
-
-        bqpjson_file = "data/bqpjson_1q.json"
-        dwisc_file = "tmp.json"
-
-        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, steps)
-
-        dwisc_data = JSON.parsefile(dwisc_file)
-        rm(dwisc_file)
-
-        ρ = simulate(ising_intended, annealing_time, annealing_schedule, steps)
-
-        @test dwisc_data["solutions"][1]["prob"] >= 0.99
-        @test isapprox(z_measure_probabilities(ρ)[2], dwisc_data["solutions"][1]["prob"])
-        @test dwisc_data["variable_ids"] == [100]
-    end
-
-    @testset "1 qubit, dwave schedule" begin
-        annealing_time = 10000.0
-        annealing_schedule = AS_DW_QUADRATIC
-        steps = 1000
-
-        ising_intended = Dict((1,) => 1)
-
-        bqpjson_file = "data/bqpjson_1q.json"
-        dwisc_file = "tmp.json"
-
-        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, steps)
-
-        dwisc_data = JSON.parsefile(dwisc_file)
-        rm(dwisc_file)
-
-        ρ = simulate(ising_intended, annealing_time, annealing_schedule, steps)
-
-        @test dwisc_data["solutions"][1]["prob"] >= 0.99
-        @test isapprox(z_measure_probabilities(ρ)[2], dwisc_data["solutions"][1]["prob"])
-        @test dwisc_data["variable_ids"] == [100]
-    end
-
-    @testset "2 qubit, function schedule" begin
-        annealing_time = 10000.0
-        annealing_schedule = AS_CIRCULAR
-        steps = 1000
-
-        # the ising model that is encoded in bqpjson_2q.json
-        ising_intended = Dict((1,) => -1, (2,) => -1, (1,2) => -1)
-
-        bqpjson_file = "data/bqpjson_2q.json"
-        dwisc_file = "tmp.json"
-
-        simulate_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, steps)
-
-        dwisc_data = JSON.parsefile(dwisc_file)
-        rm(dwisc_file)
-
-        ρ = simulate(ising_intended, annealing_time, annealing_schedule, steps)
-
-        @test dwisc_data["solutions"][1]["prob"] >= 0.99
-        @test isapprox(z_measure_probabilities(ρ)[1], dwisc_data["solutions"][1]["prob"])
-        @test dwisc_data["variable_ids"] == [304, 308]
+        @test isapprox(ρ, ρ_de, atol=1e-7)
+        @test !isapprox(ρ, ρ_de, atol=1e-8)
     end
-
-
-    @testset "1 qubit, function schedule, z noise" begin
-        annealing_time = 100.0
-        annealing_schedule = AS_CIRCULAR
-        steps = 100
-        numshots = 10
-
-        # random numbers between -0.75 and 0.75
-        z_bias = (Random.rand(numshots) .- 0.5)*1.5
-
-        bqpjson_file = "data/bqpjson_1q.json"
-        dwisc_file = "tmp.json"
-
-        simulate_noisy_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, steps, z_bias=z_bias)
-
-        dwisc_data = JSON.parsefile(dwisc_file)
-        rm(dwisc_file)
-
-        @test dwisc_data["solutions"][1]["prob"] > 0.70 # true w.h.p.
-        @test dwisc_data["variable_ids"] == [100]
-    end
-
-    @testset "1 qubit, function schedule, x and z noise" begin
-        annealing_time = 100.0
-        annealing_schedule = AS_CIRCULAR
-        steps = 100
-        numshots = 10
-
-        # random numbers between -0.75 and 0.75
-        x_bias = (Random.rand(numshots) .- 0.5)*1.5
-        z_bias = (Random.rand(numshots) .- 0.5)*1.5
-
-        bqpjson_file = "data/bqpjson_1q.json"
-        dwisc_file = "tmp.json"
-
-        simulate_noisy_bqpjson(bqpjson_file, dwisc_file, annealing_time, annealing_schedule, steps, x_bias=x_bias, z_bias=z_bias)
-
-        dwisc_data = JSON.parsefile(dwisc_file)
-        rm(dwisc_file)
-
-        @test dwisc_data["solutions"][1]["prob"] > 0.70 # true w.h.p.
-        @test dwisc_data["variable_ids"] == [100]
-    end
-
 end