text:
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear

Random.seed!(100)

testTol = 0.25

f = (u,p,t)->cos(t)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = sin(t)
prob_ode_sin = ODEProblem(f, 0.,(0.0,1.0))

f = (du,u,p,t)->du[1]=cos(t)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = [sin(t)]
prob_ode_sin_inplace = ODEProblem(f, [0.], (0.0,1.0))

f = (u,p,t)->sin(u)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = 2*acot(exp(-t)*cot(0.5))
prob_ode_nonlinear = ODEProblem(f, 1.,(0.,0.5))

f = (du,u,p,t)->du[1]=sin(u[1])
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = [2*acot(exp(-t)*cot(0.5))]
prob_ode_nonlinear_inplace = ODEProblem(f,[1.],(0.,0.5))

test_problems_only_time = [prob_ode_sin, prob_ode_sin_inplace]
test_problems_linear = [prob_ode_linear, prob_ode_2Dlinear, prob_ode_bigfloat2Dlinear]
test_problems_nonlinear = [prob_ode_nonlinear, prob_ode_nonlinear_inplace]

# Test the memory usage, cf. #640
# Note: Basically, the size of the integrator should be the size of the cache
# plus the size of the initial condition (stored is integ.sol.prob.u0) if the
# keyword argument `alias_u0` is not set to `true` (default).
# Note: The memory requirements of the 2N methods can be reduced if an assignment
# of the form `tmp = A2end[i]*tmp + dt*f(u, p, t+c2end[i]*dt)` can be carried out
# without saving `f(u, p, t+c2end[i]*dt)` as `k`.
u0_large = rand(10^6)
prob_ode_large = ODEProblem((du,u,p,t)-> du .= u, u0_large, (0.0,1.0))

alg = ORK256()
alg2 = ORK256(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CarpenterKennedy2N54()
alg2 = CarpenterKennedy2N54(;williamson_condition=false)
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = HSLDDRK64()
alg2 = HSLDDRK64(;williamson_condition=true)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = DGLDDRK73_C()
alg2 = DGLDDRK73_C(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = DGLDDRK84_C()
alg2 = DGLDDRK84_C(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = DGLDDRK84_F()
alg2 = DGLDDRK84_F(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = NDBLSRK124()
alg2 = NDBLSRK124(;williamson_condition=false)
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = NDBLSRK134()
alg2 = NDBLSRK134(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = NDBLSRK144()
alg2 = NDBLSRK144(;williamson_condition=false)
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
  sim = test_convergence(dts, prob, alg2)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg2, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CFRLDDRK64()
dts = 1 ./ 2 .^(7:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = TSLDDRK74()
dts = 1 ./ 2 .^(8:-1:4)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

println("Methods from Carpenter, Kennedy, Lewis (2000)")

function RemakeNew(p::ODEProblem)
  u1 = @. BigFloat(p.u0)
  tsp1 = @. BigFloat(p.tspan)
  remake(p; u0 = u1, tspan = tsp1)
end

test_problems_only_time_BigFloat = @. RemakeNew(test_problems_only_time)
test_problems_linear_BigFloat = @. RemakeNew(test_problems_linear)
f = (u,p,t)->sin(u)
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = 2*acot(exp(-t)*cot(BigFloat(0.5)))
prob_nonlinear_A = ODEProblem(f, BigFloat(1.),(BigFloat(0.),BigFloat(0.5)))

f = (du,u,p,t)->du[1]=sin(u[1])
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = [2*acot(exp(-t)*cot(BigFloat(0.5)))]
prob_nonlinear_B = ODEProblem(f,[BigFloat(1.)],(BigFloat(0.),BigFloat(0.5)))
test_problems_nonlinear_BigFloat = [prob_nonlinear_A,prob_nonlinear_B]

alg = CKLLSRK43_2()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol    # This scheme has linear order of 4
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3C()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ 1 atol=testTol          # The CI plot is linear but the evaluated order is 1
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK95_4S()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK95_4C()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test_broken sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK95_4M()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 7
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 6
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3C_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3M_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+0.5 atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3N_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK85_4C_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK85_4M_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK85_4P_3R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+2 atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 10
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 9
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3N_4R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK54_3M_4R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+0.5 atol=testTol                              # This scheme has linear orderof 4.5
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK65_4M_4R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK85_4FM_4R()
dts = BigFloat(1) ./ 2 .^(10:-1:6)
for prob in test_problems_only_time_BigFloat
	sim = test_convergence(dts, prob, alg)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) + 1 atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
	sim = test_convergence(dts, prob, alg)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
	sim = test_convergence(dts, prob, alg)
	@test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 12
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 11
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = CKLLSRK75_4M_5R()
dts = BigFloat(1) ./ 2 .^(8:-1:4)
for prob in test_problems_only_time_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear_BigFloat
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, adaptive=false,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 13
integ = init(prob_ode_large, alg, adaptive=true,dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 14
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 13
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

println("Methods from Parsani, Ketcheson, Deconinck (2013)")

alg = ParsaniKetchesonDeconinck3S32()
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # higher order as pure quadrature
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S82()
dts = 1 ./ 2 .^(8:-1:5)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # higher order as pure quadrature
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S53()
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # higher order as pure quadrature
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S173()
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  # higher order as pure quadrature
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg)+1 atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1 ./ 2 .^(10:-1:8)
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S94()
dts = 1 ./ 2 .^(7:-1:3)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S184()
dts = 1 ./ 2 .^(6:-1:2)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1 ./ 2 .^(7:-1:2)
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S105()
dts = 1 ./ 1.95 .^(5:-1:1)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1 ./ 2 .^(5:-1:2)
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1.5 ./ 2 .^(5:-1:2)
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]

alg = ParsaniKetchesonDeconinck3S205()
dts = 1 ./ 1.95 .^(5:-1:1)
for prob in test_problems_only_time
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1 ./ 2 .^(5:-1:2)
for prob in test_problems_linear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
dts = 1.5 ./ 2 .^(5:-1:2)
for prob in test_problems_nonlinear
  sim = test_convergence(dts, prob, alg)
  @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=testTol
end
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5
integ = init(prob_ode_large, alg, dt=1.e-2, save_start=false, save_end=false, save_everystep=false, alias_u0=true)
@test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
# test whether aliasing u0 is bad
new_prob_ode_nonlinear_inplace = ODEProblem(prob_ode_nonlinear_inplace.f,[1.],(0.,0.5))
sol_old = solve(prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false)
sol_new = solve(new_prob_ode_nonlinear_inplace, alg, dt=1.e-4, save_everystep=false, save_start=false, alias_u0=true)
@test sol_old[end] ≈ sol_new[end]


text:
using Random
import Actors: Actor, AddressBook

@testset "Addressing" begin
    null_id = Id(UInt32(0))
    test_actor(state) = Actor(Channel{Any}(), state, null_id)
    rng = MersenneTwister(33)
    book = AddressBook()

    @test_throws BoundsError book[Id(1)]
    @test book.live[1].readers[] == 0
    @test book.live[2].readers[] == 0
    @test book.flips[] == 0
    @test isempty(book)

    a = test_actor(1)
    @test Id(1) == push!(book, a)
    @test book.flips[] == 1
    @test book.live[1].entries[1] == a
    @test book.live[2].entries[1] == a
    @test book[Id(1)] == a
    @test book[a] == [UInt32(1)]
    @test !isempty(book)

    for a in [test_actor(i) for i in 2:1000]
        push!(book, a)
    end

    Threads.@threads for a in [test_actor(i) for i in 1:1000]
        push!(book, a)
    end

    @test book.live[1].entries == book.live[2].entries

    filler = test_actor(:filler)
    Threads.@threads for i in 1:1000
        i % 100 == 0 && push!(book, filler)
        @assert i == book[Id(i)].state
    end

    Threads.@threads for i in 1:lastindex(book)
        a = book[Id(i)]
        if a.state === :filler
            book[Id(i)] = nothing
        end
    end

    @test_throws KeyError book[filler]
    @test book.live[1].entries == book.live[2].entries
    @test book.live[1].readers[] == 0
    @test book.live[2].readers[] == 0

    Threads.@threads for i in rand(rng, 1:2000, 10000)
        if i % 100 == 0
            book[Id(i)] = filler
        end
        @assert book[Id(i)].state != nothing
    end

    @test book.live[1].entries == book.live[2].entries
    @test book.live[1].readers[] == 0
    @test book.live[2].readers[] == 0

    book[Id(1)] = nothing
    @test_throws KeyError book[Id(1)]
end

ext: jl


text:
import .Traceur

function trace(ex, args...)
    :($(showtrace)($(gettrace)(() -> $(esc(ex)); $(map(esc, args)...))))
end

function gettrace(f; kwargs...)
    warns = []
    Traceur.trace(x -> push!(warns, transformwarn(x)), f; kwargs...)
    warns
end

function transformwarn(w)
    meth = Traceur.method(w.call)
    file, path = expandpath(String(meth.file))
    line = max(0, w.line - 1)

    return Dict(:file => file,
                :realpath => path,
                :description => w.message,
                :range => [[line, 0], [line, 999]],
                :provider => "Traceur",
                :severity => "Warning")
end

function showtrace(trace)
    Atom.msg("staticLint", trace)
end


text:
# script to obtain the oscillating atoms using a vanilla verlet Integrator
# The main idea is to build a reference model to study the behavior fo the
# Neural network accelerated Molecular dynamic.

# we use the modelscf to compute the forces, and the velocity verlet algorithm
# to evolve the system.

# in this case we save the evolution of the system in a hd5f file.

push!(LOAD_PATH, "../src/")
using modelscf
using HDF5
using LinearAlgebra

# flag to save the data

# number of MD simulations and
Nit = 3

# getting all the parameters
dx = 0.5;
Nunit = 8;   # number of units
Lat = 10;     # size of the lattice
Ls = Nunit*Lat;
Ns = round(Integer, Ls / dx); # number of discretization points

# using the default values in Lin's code
YukawaK = 0.0100
n_extra = 10;
epsil0 = 10.0;
T_elec = 100.0;

kb = 3.1668e-6;
au2K = 315774.67;
Tbeta = au2K / T_elec;

betamix = 0.5;
mixdim = 10;

Ndist = 1;
Natoms = round(Integer, Nunit / Ndist);

sigma  = ones(Natoms,1)*(1.0);  # insulator
omega  = ones(Natoms,1)*0.03;
Eqdist = ones(Natoms,1)*10.0;
mass   = ones(Natoms,1)*42000.0;
nocc   = ones(Natoms,1)*2;          # number of electrons per atom
Z      = nocc;

function forces(x::Array{Float64,1}, rho0::Array{Float64,2})
    # input
    #       x: vector with the position of the atoms
    # output
    #       f: forces at the center of the atoms

    R = reshape(x, length(x), 1) # we need to make it a two dimensional array
    # creating an atom structure
    atoms = Atoms(Natoms, R, sigma,  omega,  Eqdist, mass, Z, nocc);
    # allocating a Hamiltonian
    ham = Ham(Lat, Nunit, n_extra, dx, atoms,YukawaK, epsil0, Tbeta)

    # total number of occupied orbitals
    Nocc = round(Integer, sum(atoms.nocc) / ham.nspin);

    # setting the options for the scf iteration
    mixOpts = andersonMixOptions(ham.Ns, betamix, mixdim )
    # eigOpts = eigOptions(1.e-10, 1000, "eigs");
    eigOpts = eigOptions(1.e-10, 1000, "lobpcg_sep");
    scfOpts = scfOptions(1.e-8, 300, eigOpts, mixOpts)

    # initialize the potentials within the Hemiltonian, setting H[\rho_0]
    init_pot!(ham, Nocc)

    # for the first iteration
    if norm(rho0) == 0
        rho0[:,:] = ham.rho
    end
    ham.rho = reshape(rho0, Ns, 1)

    (Vnew, err) = update_pot!(ham)
    ham.Vtot = Vnew

    # running the scf iteration
    VtoterrHist = scf!(ham, scfOpts)

    if VtoterrHist[end] > scfOpts.SCFtol
        println("convergence not achieved!! ")
    end

    rho0[:,:] = ham.rho

    println(length(VtoterrHist))

    # we compute the forces
    get_force!(ham)
    # computing the energy
    Vhar = hartree_pot_bc(ham.rho+ham.rhoa,ham);

    # NOTE: ham.Fband is only the band energy here.  The real total energy
    # is calculated using the formula below:
    Etot = ham.Eband + 1/2*sum((ham.rhoa-ham.rho).*Vhar)*dx;

    return ham.atoms.force[:]
end

# Setting the time evolution
dt = 0.01

x0 = zeros(Natoms); # this is defined as an 2D array
for j = 1:Natoms
  x0[j] = (j-0.5)*Lat*Ndist+dx;
end

# initial density
rhoI = zeros(Ns,1)

x0[1] = x0[1] + 1.0

v0 = zeros(Natoms)
x1 = x0 + dt*v0

#### Computing the trajectories
# TODO modify this to take in account the last
(X_traj, v, vdot) = time_evolution(velocity_verlet, x -> forces(x, rhoI), dt, Nit, x0, x1)

Rho    = zeros(Ns, Nit)
Etotal = zeros(1, Nit)
Ftotal = zeros(Natoms, Nit)

# we re-run the DFT computation using the trajectory
for ii = 1:Nit

    # we use the data of the evolved function
    R = reshape(X_traj[ii, :],Natoms,1)

    # creating an atom structure
    atoms = Atoms(Natoms, R, sigma,  omega,  Eqdist, mass, Z, nocc);

    # allocating a Hamiltonian
    ham = Ham(Lat, Nunit, n_extra, dx, atoms,YukawaK, epsil0, Tbeta)

    # total number of occupied orbitals
    Nocc = round(Integer, sum(atoms.nocc) / ham.nspin);

    # initialize the potentials within the Hemiltonian, setting H[\rho_0]
    init_pot!(ham, Nocc)

    # we use the anderson mixing of the potential
    mixOpts = andersonMixOptions(ham.Ns, betamix, mixdim )

    # we use the default options
    eigOpts = eigOptions(1.e-12, 1000, "eigs");
    scfOpts = scfOptions(1.e-10, 3000, eigOpts, mixOpts)

    # running the scf iteration
    @time VtoterrHist = scf!(ham, scfOpts)

    # we compute the forces
    get_force!(ham)

    # computing the energy
    Vhar = hartree_pot_bc(ham.rho+ham.rhoa,ham);

    # NOTE: ham.Fband is only the band energy here.  The real total energy
    # is calculated using the formula below:
    Etot = ham.Eband + 1/2*sum((ham.rhoa-ham.rho).*Vhar)*dx;

    Etotal[ii] = Etot
    Ftotal[:,ii] = ham.atoms.force[:]
    Rho[:,ii] = ham.rho[:];

    println(length(VtoterrHist))
end

print("simulation finished")
Input_str = string("../training/data/Input_KS_MD_scf_", Natoms,"_sigma_", sigma[1],".h5")
Output_str = string("../training/data/Output_KS_MD_scf_", Natoms,"_sigma_", sigma[1],".h5")

# we remoce the files if they are still around
isfile(Output_str) && rm(Output_str)
isfile(Input_str)  && rm(Input_str)

h5write(Input_str, "R", X_traj[1:Nit,:])

h5open(Output_str, "w") do file
    write(file,"Rho", Rho) # alternatively, say "@write file A"
    write(file,"Etotal", Etotal)
    write(file,"Ftotal", Ftotal)
end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
using BinDeps
using Compat
@BinDeps.setup

snappy = library_dependency("libsnappy")

# BinDeps issue: fails on sudo prompt
#provides(AptGet, {
#    "libsnappy-dev" => snappy
#})

provides(Sources, URI("https://github.com/google/snappy/releases/download/1.1.3/snappy-1.1.3.tar.gz"), snappy, unpacked_dir="snappy-1.1.3")
provides(BuildProcess, Autotools(libtarget = dirname(@__FILE__) * "/builds/libsnappy/.libs/libsnappy."*BinDeps.shlib_ext), snappy, os=:Unix)

@osx_only begin
    using Homebrew
    provides(Homebrew.HB, "snappy", snappy, os=:Darwin)
end

@compat @BinDeps.install Dict(:libsnappy => :libsnappy)


text:
# The purpose of this script is to compute the empirical mean squared error of
# various algorithms for line spetral estimation for a varying ambient dimension
# of the data while keeping the number of measurements fixed.
# The performance of classical algorithms for data from N measurements is then
# compared to the performance of matrix completion based methods with N sparse
# measurements randomly sampled from the N2 \geq N ambient dimensions.
# This allows to show how much sparse subsampling can increase the performance
# when the number of samples is kept fixed.
#
# Parameters:
# N: number of measurements
# N2range: vector of number of ambient dimensions
# K: number of sinusoids
# L: number of snapshots
# Niters: number of Monte Carlo iterations
# SNR_db: signal-to-noise-ratio [dB] of the AWGN model
# f_true: location of the frequencies of the sinusoids

include("MUSIC.jl")
include("rootMUSIC.jl")
include("ESPRIT.jl")
include("TLSesprit.jl")
include("MPM.jl")
include("ANM.jl")
include("EMAC.jl")
include("periodogram.jl")
include("findMinWrapDist.jl")
include("findMinWrapDistUnique.jl")
include("WrapDist.jl")

using Plots, Random

N2range = ceil.(Int,range(8,32,length=4))
N = 8
L = 1
SNR_db = 60
Niters = 1
f_true = [0.3 0.31]
K = length(f_true)
MSE = zeros(length(N2range),3)

A(f,N) = exp.(2*pi*1im*kron(f,(0:N-1)))
SNR = 10 .^(SNR_db./10)

for ii = 1:length(N2range)
    N2 = N2range[ii]
    for jj = 1:Niters

        noise = (randn(N,L)+1im*randn(N,L))/ (sqrt(2)*SNR)
        omega = sort(randperm(N2)[1:N])
        S = exp.(2*pi*rand(K,L)*1im)

        Y1 = A(f_true,N)*S + noise
        Y2 = A(f_true,N2)[omega,:]*S + noise

        t0 = time_ns()
        f_MUSIC = MPM(Y1,K)
        t1 = (time_ns()-t0)/1e9
        MinWrapDist = findMinWrapDistUnique(f_true,f_MUSIC)
        MSE[ii,1] += norm(MinWrapDist)^2

        f_EMAC = EMAC(Y2,N2,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_EMAC)
        MSE[ii,2] += norm(MinWrapDist)^2

        f_ANM = ANM(Y2,N2,K,omega,1/SNR)
        MinWrapDist = findMinWrapDistUnique(f_true,f_ANM)
        MSE[ii,3] += norm(MinWrapDist)^2

    end
end

MSE /= Niters

plot(N2range, MSE[:,1],yaxis=:log,tickfontsize = 5,
label="Matrix Pencil Method",
#xlims = (0,1),
#xticks = round.(f_range.-f1; digits=2),
xlabel = "M",
ylabel = "MSE"
)
plot!(N2range,MSE[:,2],label="EMAC")
plot!(N2range,MSE[:,3],label="ANM")
#savefig("/home/thomas/atom-amd64/usr/share/atom/example/mse_N2.png")


text:
"""
    XeLaTeX

The `xelatex` engine.

Command line options `-no-pdf`, `-quiet`, and `-halt-on-error` are passed to the `xelatex` executable, and both the input and output directories are also handled. If you would like to pass additional command line arguments, use the [`latexsvg`](@ref) function.
"""
struct XeLaTeX <: LaTeXEngine end

function runlatex(
    ::XeLaTeX,
    input_file::AbstractString,
    output_path::AbstractString;
    extra_args::Vector{String}=String[],
    runcommand::Bool=true
)
    command = `xelatex $extra_args -no-pdf -output-directory="$output_path" -quiet -halt-on-error "$input_file"`
    if runcommand
        run(command)
    end
    return command
end

dviext(::XeLaTeX) = "xdv"


text:
module CuTextures

export CuTextureArray, CuTexture, CuDeviceTexture, cuda_texture_alias_type

import CUDAdrv
import CUDAnative # TODO: handle CUDA context creation. This is just to create a context already (?)
import CuArrays
import Adapt

include("formats.jl")
include("texturearray.jl")
include("texture.jl")
include("native.jl")

end # module


text:
"""
    c = @sc_native f(args...)
    c = @sc_intel f(args...)
    c = (@sc_llvm f(args...)).native
    c = (@sc_llvm f(args...)).att
    c = (@sc_llvm f(args...)).intel

Native code explore.

```julia
c                  # view code in the REPL
display(c)         # (ditto)
edit(c)            # open the code in editor
print(c)           # print the code
abspath(c)         # file path to the text containing the code
```
"""
(:(@native), :(@intel))

"""
    ShowCode.from_native(code::AbstractString)

Construct a native code explore from a snippet of ASM.

For example,

```julia
a = @sc_native dump_module=true f(...)
b = ShowCode.from_native(string(a))
```

should be roughly equivalent.
"""
function ShowCode.from_native(code::AbstractString; syntax = :att)
    # TODO: guess syntax?
    return CodeNative(String(code), syntax, true, ("f?", "Tuple{?}"), nothing)
end

struct CodeNative <: AbstractCode
    code::String
    syntax::Symbol
    user_dump_module::Bool
    args::Any
    kwargs::Any
    cache::Dict{Symbol,Any}
    abspath::Base.RefValue{Union{Nothing,String}}
end

CodeNative(args::Vararg{Any,5}) =
    CodeNative(args..., Dict{Symbol,Any}(), Ref{Union{Nothing,String}}(nothing))

Base.string(native::CodeNative) = Fields(native).code
Base.print(io::IO, native::CodeNative) = print(io, string(native))

function Base.abspath(native::CodeNative)
    p = Fields(native).abspath[]
    p === nothing || return p
    p = joinpath(mktempdir(prefix = "jl_codeviz_"), "code.s")
    write(p, string(native))
    Fields(native).abspath[] = p
    return p
end

Base.propertynames(native::CodeNative) = (
    # what else?
)

Base.getproperty(native::CodeNative, name::Symbol) =
    begin
        error("Unknown property: ", name)
    end

function Base.summary(io::IO, native::CodeNative)
    f, t = Fields(native).args
    print(io, "CodeNative of ", f, " with ", t)
    return
end

function Base.show(io::IO, ::MIME"text/plain", native::CodeNative)
    @unpack code = Fields(native)
    summary(io, native)
    println(io)
    if get(io, :color, false)
        print_native(io, code)
    else
        print(io, code)
    end
    return
end

macro native(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, sc_native, args)
end

macro intel(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, sc_intel, args)
end

function sc_native(args...; dump_module = false, syntax = :att, kwargs...)
    if dump_module
        return getproperty(
            sc_llvm(args...; dump_module = dump_module, kwargs...),
            syntax,
        )
    end
    @nospecialize
    code = sprint() do io
        @nospecialize
        code_native(io, args...; dump_module = true, syntax = syntax, kwargs...)
    end
    return CodeNative(code, syntax, dump_module, args, kwargs)
end

sc_intel(args...; kwargs...) = sc_native(args...; syntax = :intel, kwargs...)

function CodeNative(llvm::CodeLLVM, syntax::Symbol)
    @unpack user_dump_module, args, kwargs = Fields(llvm)
    CodeNative(
        llvm_to_native(llvm, `--x86-asm-syntax=$syntax`),
        syntax,
        user_dump_module,
        args,
        kwargs,
    )
end

function llvm_to_native(ir, options = ``)
    cmd = getcmd(:llc)
    cmd = `$cmd $options -o=- --filetype=asm -`
    io = IOBuffer()
    write_silently(cmd, string(ir); stdout = io)
    return String(take!(io))
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

using LinearAlgebra

N1 = 100
N2 = div(N1, 2) + 5
A = randn(N1, N2)

x0 = 10 .* randn(N2)
x0[abs.(x0) .< 1f-6] .= 1.0
inds = rand(1:N2, div(N2, 4))
ninds = [i for i=1:N2 if i ∉ inds]
x0[inds] .= 0
b = A*x0

function obj(x)
    fun = .5*norm(A*x - b)^2
    grad = A'*(A*x - b)
    return fun, grad
end

opt = bregman_options(maxIter=200, progTol=0, verbose=2)
sol = bregman(obj, 1 .+ randn(N2), opt)

@show sol.x[inds]
@show x0[inds]
@show sol.x[ninds]
@show x0[ninds]

part_n = i -> norm(sol.x[i] - x0[i])/(norm(x0[i]) + norm(sol.x[i]) + eps(Float64))
part_nz = i -> norm(sol.x[i], 1)/N2
@show part_nz(inds)
@show part_n(ninds)

@test part_nz(inds) < 1f-1
@test part_n(ninds) < 1f-1
@test sol.residual/sol.r_trace[1] < 1f-1


text:
baremodule ConcurrentUtils

export
    # Macros
    @tasklet,
    # Constructors
    Backoff,
    Guard,
    NotSetError,
    OccupiedError,
    Once,
    Promise,
    ReadWriteGuard,
    ReadWriteLock,
    ThreadLocalStorage

export Try, Err, Ok
using Try: Try, Ok, Err

module InternalPrelude
include("prelude.jl")
end  # module InternalPrelude

InternalPrelude.@exported_function race_put!
InternalPrelude.@exported_function race_put_with!
InternalPrelude.@exported_function try_race_fetch
InternalPrelude.@exported_function try_race_put!
InternalPrelude.@exported_function try_race_put_with!

macro tasklet end

struct OccupiedError{T} <: InternalPrelude.Exception
    value::T
end

struct NotSetError <: InternalPrelude.Exception end

#=
InternalPrelude.@exported_function isacquirable
InternalPrelude.@exported_function isacquirable_read
=#

InternalPrelude.@exported_function trylock_read
InternalPrelude.@exported_function lock_read
InternalPrelude.@exported_function unlock_read

abstract type AbstractGuard end
abstract type AbstractReadWriteGuard end
# Should abstract types have `Data` as a parameter? But maybe some implementations may want
# to provide a read-only wrapper type?

struct GenericGuard{Lock,Data} <: AbstractGuard
    lock::Lock
    data::Data
    GenericGuard(lock::Lock, data::Data) where {Lock,Data} = new{Lock,Data}(lock, data)
end

struct GenericReadWriteGuard{Lock,Data} <: AbstractReadWriteGuard
    lock::Lock
    data::Data
    GenericReadWriteGuard(lock::Lock, data::Data) where {Lock,Data} =
        new{Lock,Data}(lock, data)
end

InternalPrelude.@exported_function guardwith
InternalPrelude.@exported_function guarding
InternalPrelude.@exported_function guarding_read
# InternalPrelude.@exported_function read_write_guard

InternalPrelude.@exported_function unsafe_takestorages!

InternalPrelude.@exported_function spinloop
InternalPrelude.@exported_function spinfor

# Maybe:
# * @static_thread_local_storage
# * Copy AsyncFinalizers.SingleReaderDualBag?

"""
    Internal

Internal module that contains main implementations.
"""
module Internal

using Core.Intrinsics: atomic_fence
using Core: OpaqueClosure
using Random: Xoshiro

import UnsafeAtomics: UnsafeAtomics, acq_rel
using ExternalDocstrings: @define_docstrings
using Try: Try, Ok, Err, @?

import ..ConcurrentUtils: @tasklet
using ..ConcurrentUtils:
    AbstractGuard,
    AbstractReadWriteGuard,
    ConcurrentUtils,
    GenericGuard,
    GenericReadWriteGuard,
    NotSetError,
    OccupiedError,
    lock_read,
    race_put!,
    race_put_with!,
    spinfor,
    spinloop,
    try_race_fetch,
    try_race_put!,
    try_race_put_with!,
    trylock_read,
    unlock_read

#=
if isfile(joinpath(@__DIR__, "config.jl"))
    include("config.jl")
else
    include("default-config.jl")
end
=#

include("utils.jl")
include("promise.jl")
include("once.jl")
include("tasklet.jl")
include("thread_local_storage.jl")

# Locks
include("read_write_lock.jl")
include("guards.jl")
include("backoff.jl")

end  # module Internal

const Promise = Internal.Promise
const Once = Internal.Once
const ThreadLocalStorage = Internal.ThreadLocalStorage
const ReadWriteLock = Internal.ReadWriteLock
const Backoff = Internal.Backoff

const Guard = Internal.Guard
const ReadWriteGuard = Internal.ReadWriteGuard

Internal.@define_docstrings

end  # baremodule ConcurrentUtils


text:
module SwapStreams

export SwapStream

const LittleEndian = 0x01020304
const BigEndian = 0x04030201

"""

```jldoctest
julia> s = SwapStream(IOBuffer(), true);

julia> write(s, [1:10...])

julia> seek(s, 0);

julia> read!(s.io, Vector{Int}(undef, 10))
10-element Array{Int64,1}:
  72057594037927936
 144115188075855872
 216172782113783808
 288230376151711744
 360287970189639680
 432345564227567616
 504403158265495552
 576460752303423488
 648518346341351424
 720575940379279360

julia> seek(s, 0);

julia> read!(s, Vector{Int}(undef, 10))
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

```
"""
struct SwapStream{S,IOType<:IO} <: IO
    io::IOType

    SwapStream{S}(io::IOType) where {S,IOType<:IO} = new{S,IOType}(io)
end

SwapStream(io::IO, endianness::UInt32) =
    SwapStream{ENDIAN_BOM != endianness}(io)
SwapStream(io::IO; needswap::Bool=false) = SwapStream{needswap}(io)

Base.seek(s::SwapStream, n::Integer) = seek(s.io, n)
Base.position(s::SwapStream)  = position(s.io)
Base.skip(s::SwapStream, n::Integer) = skip(s.io, n)
Base.eof(s::SwapStream) = eof(s.io)
Base.isreadonly(s::SwapStream) = isreadonly(s.io)
Base.isreadable(s::SwapStream) = isreadable(s.io)
Base.iswritable(s::SwapStream) = iswritable(s.io)
Base.stat(s::SwapStream) = stat(s.io)
Base.close(s::SwapStream) = close(s.io)
Base.isopen(s::SwapStream) = isopen(s.io)
Base.ismarked(s::SwapStream) = ismarked(s.io)
Base.mark(s::SwapStream) = mark(s.io)
Base.unmark(s::SwapStream) = unmark(s.io)
Base.reset(s::SwapStream) = reset(s.io)
Base.seekend(s::SwapStream) = seekend(s.io)

read(s::SwapStream, n::Int) = read(s.io, n)
read(s::SwapStream{true}, x::Type{<:AbstractArray}) = bswap.(read(s.io, x))
read(s::SwapStream{false}, x::Type{<:AbstractArray}) = read(s.io, x)
read(s::SwapStream, x::Type{Int8}) = read(s.io, Int8)

read!(s::SwapStream{true,<:IO}, a::Array{Int8,N}) where N = read!(s.io, a)
read!(s::SwapStream{false,<:IO}, a::Array{Int8,N}) where N = read!(s.io, a)

read!(s::SwapStream{true,<:IO}, a::Array{UInt8,N}) where N = read!(s.io, a)
read!(s::SwapStream{false,<:IO}, a::Array{UInt8,N}) where N = read!(s.io, a)

read!(s::SwapStream{false,IOType}, a::Array{T}) where {IOType<:IO,T} = read!(s.io, a)
read!(s::SwapStream{true,IOType}, a::Array{T}) where {IOType<:IO,T} = bswap.(read!(s.io, a))

function read(s::SwapStream{false},
              T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}}
             )
    return read!(s.io, Ref{T}(0))[]::T
end

function read(s::SwapStream{true},
              T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}}
             )
    return bswap(read!(s.io, Ref{T}(0))[]::T)
end

read!(s::SwapStream{true}, ref::Base.RefValue{<:NTuple{N,T}}) where {N,T} = bswap.(read!(s.io, ref)[])
read!(s::SwapStream{false}, ref::Base.RefValue{<:NTuple{N,T}}) where {N,T} = read!(s.io, ref)[]

write(s::SwapStream{true,<:IO}, x::Array) = write(s.io, mappedarray(ntoh, hton, x))
write(s::SwapStream{true,<:IO}, x::T) where T = write(s.io, bswap.(x))
write(s::SwapStream{false,<:IO}, x::T) where T = write(s.io, x)

end


text:
Base.sqrt(f::Dual) = Dual(sqrt(f.val), f.grad ./ sqrt(f.val))
Base.log(f::Dual) = Dual(log(f.val), f.grad ./ f.val)
Base.conj(f::Dual) = Dual(conj(f.val), conj.(f.grad))

max_stars_repo_name: kaosmicadei/AutoError

max_issues_repo_name: kaosmicadei/AutoError

max_forks_repo_name: kaosmicadei/AutoError


text:
using Distributions
using PyPlot
using Roots
using JLD

function tauchen86(μ::AbstractFloat,ρ::AbstractFloat,σ::AbstractFloat,n::Int64,λ::AbstractFloat)
    #output is xgrid, xprob
    # x_t+1 = μ + ρ x_t + σ e_{t+1}, e_t+1 ∼ N(0,1)
    xhi = μ/(1-ρ) + λ*sqrt(σ^2/(1-ρ)^2)
    xlo = μ/(1-ρ) - λ*sqrt(σ^2/(1-ρ)^2)
    xgrid = zeros(n);
    xscale =(xhi-xlo)/(n-1)

    for i=1:n
        xgrid[i] = xlo + xscale*(i-1)
    end
    m=zeros(n-1);
    for i=1:n-1
        m[i] = (xgrid[i]+xgrid[i+1])/2
    end
    xprob = zeros(n,n) # xprob[i,j] = Pr(x_t+1=xgrid[j]|x_t=xgrid[i])
    for i=1:n # this is the state today
        normie = Normal(μ+ρ*xgrid[i],σ)
        normpdf(x) = pdf.(normie,x)
        normcdf(x) = cdf.(normie,x)
        for j=2:n-1
            xprob[i,j] = normcdf(m[j]) - normcdf(m[j-1])
        end
        xprob[i,1] = normcdf(m[1])
        xprob[i,n] = 1 - normcdf(m[n-1])
    end
    xprob = xprob ./ sum(xprob,2) # make sure the rows sum to 1
    return ( xgrid,xprob, xscale )
end

function trunc_lognpdf(x, UP, LO, mu, sig)
    # truncated log normal pdf
    logn = LogNormal(mu, sig)
    return (x-LO.>=-eps()) .* (x-UP.<=eps()).*pdf.(logn,x)/(cdf.(logn,UP)-cdf.(logn,LO))
end

# set parameters
# PARAMETERS THAT AFFECT STEADY STATE:
# note, we have log utility, so no CRRA parameter
# and we have a quarterly model
r = 0.01                        # 4 percent steady state annual real interest rate
α = 0.3                         # capital share
H = 1.0                         # aggregate hours worked
δ  = 0.03                       # depreciation
sH_over_sL = 1.2 / 0.8
pLH = 0.1                       # prob of going from low to high persistent skill
pHL = 0.1                      # prob of going fom
γ = 0.#0.004                       # TFP growth
GoverY = 0.01#18                        # steady state govt spending/gdp
g = 1/(1-GoverY)                # steady state value of g_t, where G_t/Y_t = (1-(1/g_t))
η = 0.                         # borrowing constraint (normalized by TFP)
In    = 0.443993816237631       # normalizing constant for the mollifier
na = 300                         # cash on hand ditn grid points - set to smaller numbers for debugging
ns = 2# 5                         # skill distribution grid points
zlo =  1. / 3.                      # second income shock to mollify actual income
zhi = 2. - zlo                       # upper bound on this process
BoverY = 0.26

# AGGREGATE SHOCKS
ρB = 0.5 # persistence of discount factor shock
ρG = 0.5 # persistence of govt spending shock
ρZ = 0.5 # persistence of tfp growth shock
ρμ = 0.5 # persistence of investment shock
ρlamw = 0.5 # persistence of wage mkup shock
ρlamf = 0.5 # persistence of price mkup shock
ρmon = 0.5 # persistence of mon policy shock

# OTHER PARAMETERS THAT AFFECT DYNAMICS
spp  = 4.                        # second derivative of investment adjustment cost
lamw = 1.5                       # wage markup
ϕh   = 2.                        # inverse frisch elasticity
Φw   = 10.                       # rotemberg cost for wages
lamf = 1.5                       # price markup
Φp   = 1.                       # rotemberg cost for prices
ρR   = 0.75                      # persistence in taylor rule
ψπ   = 1.5                      # weight on inflation in taylor rule
ψy   = 0.5                       # weight on output growth in taylor rule

# OPTIONS
TRUNCATE      = true # default: true
mindens       = 1e-8
RESCALE_XWTS  = true # default: true
SIMPLE_TAYLOR = false # default: false

# solve for steady state of aggregate scalar variables
# here I am assuming firms are subsidized in ss so that real mc = 1
Rk = r + δ                               # rental rate on capital
ω = (α^(α/(1-α)))*(1-α)*Rk^(-α/(1-α))    # real wage
kl = (α/(1-α))*(ω/Rk)*ℯ^γ                # capital/labor ratio
k = kl*H                                 # capital
x = (1-(1-δ)*ℯ^(-γ))*k                   # investment
y = (ℯ^(-α*γ))*(k^α)*H^(1-α)             # gdp
bg = BoverY*y                            # govt debt
T = Rk*k*ℯ^(-γ) - x - (1-(1/g))*y - ((1+r)*ℯ^(-γ)-1)*bg        # net transfer to hhs

function persistent_skill_process(sH_over_sL::AbstractFloat, pLH::AbstractFloat, pHL::AbstractFloat, ns::Int)
    f1 = [[1-pLH pLH];[pHL 1-pHL]] # f1[i,j] is prob of going from i to j
    ss_skill_distr = [pHL/(pLH+pHL); pLH/(pLH+pHL)]
    slo = 1. / (ss_skill_distr'*[1;sH_over_sL])
    sgrid = slo*[1;sH_over_sL]
    sscale = sgrid[2] - sgrid[1]
    swts     = (sscale/ns)*ones(ns) #quadrature weights
    f = f1 ./ repmat(swts',ns,1)
    return (f, sgrid, swts)
end
(f, sgrid, swts) = persistent_skill_process(sH_over_sL, pLH, pHL, ns)

function cash_grid(sgrid::AbstractArray, ω::AbstractFloat, H::AbstractFloat, r::AbstractFloat, η::AbstractFloat, γ::AbstractFloat,
                   T::AbstractFloat, zlo::AbstractFloat, na::Int)
    smin = minimum(sgrid)*zlo                           # lowest possible skill
    alo_ss = ω*smin*H - (1+r)*η*ℯ^(-γ) + T + sgrid[1]*ω*H*0.05       # lowest possible cash on hand in ss

    alo = alo_ss                    # lower bound on cash on hand - could be < alo_ss
    ahi = max(alo*2,alo+5.)                      # upper bound on cash on hand
    ascale = (ahi-alo)              # size of w grids

    # make grids
    agrid    = collect(linspace(alo,ahi,na)) #Evenly spaced grid
    awts = (ascale/na)*ones(na)          #quadrature weights
    return (agrid, awts)
end
smin = minimum(sgrid)*zlo                           # lowest possible skill
alo_ss = ω*smin*H - (1+r)*η*ℯ^(-γ) + T + sgrid[1]*ω*H*0.05       # lowest possible cash on hand in ss

alo = alo_ss                    # lower bound on cash on hand - could be < alo_ss
ahi = max(alo*2,alo+5.)                      # upper bound on cash on hand
ascale = (ahi-alo)              # size of w grids

(agrid, awts) = cash_grid(sgrid, ω, H, r, η, γ, T, zlo, na)

function mollifier(z::AbstractFloat,ehi::AbstractFloat,elo::AbstractFloat)
    # mollifier function
    In = 0.443993816237631
    if z<ehi && z>elo
        temp = -1.0 + 2.0*(z-elo)/(ehi-elo)
        out = (2.0/(ehi-elo))*exp(-1.0/(1.0-temp*temp))/In
    else
        out = 0.0
    end

    return out
end

zgrid = linspace(zlo,zhi,na)
zwts = (zhi-zlo)/na
sumz=0.
for i=1:na
    sumz += mollifier(zgrid[i],zhi,zlo)*zwts
end

# experiment with a different q
ne = 100
σe = 0.01
(legrid, fe, sscale) = tauchen86(0.0,0.0,σe,ne,2.0)
egrid = exp.(legrid)
eprob = fe[1,:]

function convoluted_q(x::AbstractFloat, zhi::AbstractFloat, zlo::AbstractFloat, ne::Int, egrid::Vector{Float64}, eprob::Vector{Float64})
    sumne = 0.
    for i=1:ne
        sumne += mollifier(x - egrid[i], zhi, zlo)*eprob[i]
    end
    return sumne
end

qfunction(x) = mollifier(x, zhi, zlo)#/sumz # this ensures the mollifer sums to 1, at least on this particular grid
#qfunction(x) = convoluted_q(x, zhi, zlo, ne, egrid, eprob)

function dmollifier(x::AbstractFloat, ehi::AbstractFloat, elo::AbstractFloat)
    In = 0.443993816237631
    if x<ehi && x>elo
        temp = (-1.0 + 2.0*(x-elo)/(ehi-elo))
        out  = -(2*temp ./ ((1 - temp.^2).^2)) .* (2/(ehi-elo)).*mollifier(x, ehi, elo)
    else
        out = 0.0
    end
    return out
end

qp(z) = dmollifier(z,zhi,zlo)

Win = 1 ./ (5. + 0.02*(repeat(agrid,ns)-5))
Win = 2*ones(na*ns)/(ahi+alo)
aswts = kron(swts,awts)

function sspolicy(na::Int, ns::Int, β::AbstractFloat, R::AbstractFloat,
                     ω::AbstractFloat,
                     H::AbstractFloat,
                     η::AbstractFloat,
                     T::AbstractFloat,
                     γ::AbstractFloat,
                     qfunction::Function,
                     agrid::Vector{Float64},
                     sgrid::Vector{Float64},
                     aswts::Vector{Float64},
                     Win::Vector{Float64},
                     f::Array{Float64,2},damp::Float64 = 0.5, dist::Float64 = 1.,
                         tol::Float64 = 1e-4, maxit::Int64 = 500)
    n = na*ns
    c  = zeros(n)      # consumption
    bp = zeros(n)      # savings
    counter = 1
    Wout = copy(Win)
    while dist>tol && counter<maxit # for debugging
        # compute c(w) given guess for Win = β*R*E[u'(c_{t+1})]
        for iss in 1:ns
            for ia in 1:na
                c[na*(iss-1)+ia] =  min(1/Win[na*(iss-1)+ia],agrid[ia]+η)
            end
        end
        bp = repeat(agrid,ns) - c  # compute bp(w) given guess for Win
        Wout = parameterized_expectations(na,ns,β,R,ω,H,η,T,γ,qfunction,agrid,sgrid,aswts,c,bp,f)
        dist = maximum(abs.(Wout-Win))
        #println(dist)
        if mod(counter,25)==1
            println(dist)
            #plot(Wout)
        end
        Win = damp*Wout + (1.0-damp)*Win
        counter += 1
    end
    if counter==maxit
        warn("Euler iteration did not converge")
    end
    tr = kolmogorov_fwd(na,ns,ω,H,η,T,γ,qfunction,agrid,sgrid,bp,f)
    return (c,bp,Wout,tr)
end

function parameterized_expectations(na::Int,ns::Int, β::AbstractFloat, R::AbstractFloat, ω::AbstractFloat,
                                    H::AbstractFloat, η::AbstractFloat, T::AbstractFloat, γ::AbstractFloat,
                                    qfunction::Function,
                                    agrid::Vector{Float64}, sgrid::Vector{Float64},
                                    aswts::Vector{Float64}, c::Vector{Float64},
                                    bp::Vector{Float64}, f::Array{Float64,2})
    l_out = zeros(na*ns)
    for iss=1:ns
        for ia=1:na
            sumn = 0.0
            for isp=1:ns
                for iap=1:na
                    sumn += (aswts[na*(isp-1)+iap]/c[na*(isp-1)+iap])*qfunction((agrid[iap] - R*(ℯ^(-γ))*bp[na*(iss-1)+ia] - T)/(ω*H*sgrid[isp]))*f[iss,isp] ./ sgrid[isp]
                end
            end
            l_out[na*(iss-1)+ia] = (β*R*(ℯ^(-γ))/ω*H)*sumn
        end
    end
    return l_out
end

function kolmogorov_fwd(na::Int, ns::Int, ω::AbstractFloat,
                        H::AbstractFloat, η::AbstractFloat, T::AbstractFloat, γ::AbstractFloat,
                        qfunction::Function,agrid::Vector{Float64},
                        sgrid::Vector{Float64}, bp::Vector{Float64},f::Array{Float64,2})
    tr = zeros(na*ns,na*ns)
    for iss=1:ns
        for ia=1:na
            for isp=1:ns
                for iap=1:na
                    tr[na*(isp-1)+iap,na*(iss-1)+ia] = qfunction((agrid[iap] - R*(ℯ^(-γ))*bp[na*(iss-1)+ia] - T)/(ω*H*sgrid[isp]))*f[iss,isp] ./ (ω*H*sgrid[isp])
                end
            end
        end
    end
    return tr
end

# just to test if code returns
R = 1+r
βlo = 0.5*(ℯ^γ)/R # excess should be -ve
βhi = (ℯ^γ)/R # excess should be +ve

β = (βlo+βhi)/2.0

function findss(na::Int, ns::Int, βlo::AbstractFloat,
                βhi::AbstractFloat, R::AbstractFloat,
                     ω::AbstractFloat,
                     H::AbstractFloat,
                     η::AbstractFloat,
                     T::AbstractFloat,
                     γ::AbstractFloat,
                     bg::AbstractFloat,
                     qfunction::Function,
                     agrid::Vector{Float64},
                     sgrid::Vector{Float64},
                     aswts::Vector{Float64},
                     Win::Vector{Float64},
                     f::Array{Float64,2},
                     excess::AbstractFloat = 5000., tol::AbstractFloat = 1e-4,
                     maxit::Int64 = 20)
    counter=1
    n       = ns*na
    c = zeros(n)
    bp = zeros(n)
    KF = zeros(n,n)
    m = zeros(n)
    report = zeros(maxit,2)
    #Win0=copy(Win)
    while abs(excess)>tol && counter<maxit # clearing markets
        β = (βlo+βhi)/2.0
        Win_guess = ones(n) # or Win
        (c, bp, Win, KF) = sspolicy(na, ns, β, R, ω, H, η, T, γ, qfunction, agrid, sgrid,
                                       aswts, Win_guess, f)

        LPMKF=aswts[1]*KF
        # find eigenvalue closest to 1
        (D,V) = eig(LPMKF)
        order_D = sortperm(abs.(D), rev = true)
        V = V[:,order_D]
        D = D[order_D]
        if abs(D[1]-1)>2e-1 # that's the tolerance we are allowing
            warn("your eigenvalue is too far from 1, something is wrong")
        end
        m = real(V[:,1]) #Pick the eigen vecor associated with the largest eigenvalue and moving it back to values
        m = m/(aswts'*m) #Scale of eigenvectors not determinate: rescale to integrate to exactly 1
        excess = (aswts'*(m .* bp))[1] - bg  #compute excess supply of savings, which is a fn of w
                # bisection
        println([βlo β βhi])
        println(excess)
        report[counter,1] = β # save the history of guesses and excess supply so we can check it's upward-sloping and monotonic
        report[counter,2] = excess
        if excess>0
            βhi=β
        elseif excess<0
            βlo = β
        end
        counter += 1
    end
    return (Win, c, m, β, report)
end

tic()
(ell, c, m, β, report) = findss(na, ns, βlo, βhi, R, ω, H, η, T, γ, bg, qfunction, agrid, sgrid,aswts, Win, f)
toc()


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v03.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")
include("EnergiesT.jl")
include("PotentialsT.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

include("andersonmix.jl")
include("KS_solve_SCF_andersonmix.jl")

include("pulaymix.jl")
include("KS_solve_SCF_pulaymix.jl")

include("get_ub_lb_lanczos.jl")
include("KS_solve_ChebySCF.jl")
include("chebyfilt.jl")
include("norm_matrix_induced.jl")

function test_main( ecutwfc_Ry::Float64; method="SCF" )

    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( 0.5*ecutwfc_Ry, LatVecs )

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Ns = pw.Ns
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("ecutwfc (Ha) = %f\n", ecutwfc_Ry*0.5)
    @printf("Ns = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # C atom
    Zatm = 6.0;
    Nstates = 3
    Focc = [2.0, 2.0, 2.0]

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = calc_strfact( Xpos, 1, [1], pw.gvec.G )
    E_nn = calc_ewald( pw, Sf, Xpos, 1, [1], [Zatm] )

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi*Zatm/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    # Need to sum up over Nspecies for more than one species type
    # We simply need reshape because we only have one species type here.
    V_ionic = reshape( V_ionic, (Npoints) )

    if method == "CG"
        psi, Energies, Potentials = KS_solve_Emin_cg(
            pw, V_ionic, Focc, Nstates, NiterMax=1000, E_NN=E_nn )
        #
        Y = ortho_gram_schmidt(psi)
        mu = Y' * op_H( pw, Potentials, Y )
        evals, evecs = eig(mu)
        psi = Y*evecs
    elseif method == "SCF"
        Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates, E_NN=E_nn )
    elseif method == "ChebySCF"
        Energies, Potentials, psi, evals = KS_solve_ChebySCF( pw, V_ionic, Focc, Nstates, β=0.8, E_NN=E_nn )
    elseif method == "SCF_andersonmix"
        Energies, Potentials, psi, evals = KS_solve_SCF_andersonmix( pw, V_ionic, Focc, Nstates, β=0.1, E_NN=E_nn )
    elseif method == "SCF_pulaymix"
        Energies, Potentials, psi, evals = KS_solve_SCF_pulaymix( pw, V_ionic, Focc, Nstates, β=0.1, E_NN=E_nn)
    else
        println("ERROR: unknown method: ", method)
        exit()
    end

    for st = 1:Nstates
        @printf("State # %d, Energy = %f\n", st, real(evals[st]))
    end

    print_Energies(Energies)

end

#@time test_main( 35.0, method="SCF" )
#@time test_main( 35.0, method="CG" )
#@time test_main( 35.0, method="ChebySCF" )
@time test_main( 35.0, method="SCF_pulaymix" )
@time test_main( 35.0, method="SCF_andersonmix" )