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)
prob_ode_sin = ODEProblem(ODEFunction(f; analytic=(u0,p,t)->sin(t)), 0.,(0.0,1.0))

f = (du,u,p,t)->du[1]=cos(t)
prob_ode_sin_inplace = ODEProblem(ODEFunction(f;analytic=(u0,p,t)->[sin(t)]), [0.], (0.0,1.0))

f = (u,p,t)->sin(u)
prob_ode_nonlinear = ODEProblem(ODEFunction(f;analytic=(u0,p,t)->2*acot(exp(-t)*cot(0.5))), 1.,(0.,0.5))

f = (du,u,p,t)->du[1]=sin(u[1])
prob_ode_nonlinear_inplace = ODEProblem(ODEFunction(f;analytic=(u0,p,t)->[2*acot(exp(-t)*cot(0.5))]),[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))

@testset "ORK256" begin
  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]
end

@testset "CarpenterKennedy2N54" begin
  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]
end

@testset "SHLDDRK64" begin
  alg = SHLDDRK64()
  alg2 = SHLDDRK64(;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]
end

@testset "DGLDDRK73_C" begin
  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]
end

@testset "DGLDDRK84_C" begin
  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]
end

@testset "DGLDDRK84_F" begin
  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]
end

@testset "NDBLSRK124" begin
  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]
end

@testset "NDBLSRK134" begin
  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]
end

@testset "NDBLSRK144" begin
  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]
end

@testset "CFRLDDRK64" begin
  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]
end

@testset "TSLDDRK74" begin
  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]
end

# 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)
prob_nonlinear_A = ODEProblem(ODEFunction(f;analytic=(u0,p,t)->2*acot(exp(-t)*cot(BigFloat(0.5)))), BigFloat(1.),(BigFloat(0.),BigFloat(0.5)))

f = (du,u,p,t)->du[1]=sin(u[1])
prob_nonlinear_B = ODEProblem(ODEFunction(f;analytic=(u0,p,t)->[2*acot(exp(-t)*cot(BigFloat(0.5)))]),[BigFloat(1.)],(BigFloat(0.),BigFloat(0.5)))
test_problems_nonlinear_BigFloat = [prob_nonlinear_A,prob_nonlinear_B]

@testset "CKLLSRK43_2" begin
  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]
end

@testset "CKLLSRK54_3C" begin
  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]
end

@testset "CKLLSRK95_4S" begin
  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]
end

@testset "CKLLSRK95_4C" begin
  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]
end

@testset "CKLLSRK95_4M" begin
  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]
end

@testset "CKLLSRK54_3C_3R" begin
  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]
end

@testset "CKLLSRK54_3M_3R" begin
  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]
end

@testset "CKLLSRK54_3N_3R" begin
  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]
end

@testset "CKLLSRK85_4C_3R" begin
  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]
end

@testset "CKLLSRK85_4M_3R" begin
  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]
end

@testset "CKLLSRK85_4P_3R" begin
  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]
end

@testset "CKLLSRK54_3N_4R" begin
  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]
end

@testset "CKLLSRK54_3M_4R" begin
  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]
end

@testset "CKLLSRK65_4M_4R" begin
  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]
end

@testset "CKLLSRK85_4FM_4R" begin
  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]
end

@testset "CKLLSRK75_4M_5R" begin
  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]
end

# Methods from Parsani, Ketcheson, Deconinck (2013)

@testset "ParsaniKetchesonDeconinck3S32" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S82" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S53" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S173" begin
  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 .^(6:-1:3)
  for prob in test_problems_nonlinear
    sim = test_convergence(dts, prob, alg)
    @test sim.𝒪est[:final] ≈ OrdinaryDiffEq.alg_order(alg) atol=1
  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]
end

@testset "ParsaniKetchesonDeconinck3S94" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S184" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S105" begin
  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]
end

@testset "ParsaniKetchesonDeconinck3S205" begin
  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]
end

# Methods from Ranocha, Dalcin, Parsani, Ketcheson (2021)

@testset "RDPK3Sp35" begin
  alg = RDPK3Sp35()
  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) <= 6
  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) <= 5
  # 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]
end

@testset "RDPK3Sp49" begin
  alg = RDPK3Sp49()
  dts = 1 ./ 2 .^(5:-1:2)
  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 .^(8:-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) <= 6
  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) <= 5
  # 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]
end

@testset "RDPK3Sp510" begin
  alg = RDPK3Sp510()
  dts = 1 ./ 2 .^(4.5:-1:1.5)
  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) <= 6
  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) <= 5
  # 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]
end

@testset "RDPK3SpFSAL35" begin
  alg = RDPK3SpFSAL35()
  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) <= 6
  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) <= 5
  # 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]
end

@testset "RDPK3SpFSAL49" begin
  alg = RDPK3SpFSAL49()
  dts = 1 ./ 2 .^(5:-1:2)
  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 .^(8:-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) <= 6
  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) <= 5
  # 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]
end

@testset "RDPK3SpFSAL510" begin
  alg = RDPK3SpFSAL510()
  dts = 1 ./ 2 .^(4.5:-1:1.5)
  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) <= 6
  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) <= 5
  # 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]
end


text:
# QP example
#
#    minimize 2 x^2 + y^2 + xy + x + y
#
#       s.t.  x, y >= 0
#             x + y = 1
#
#    solution: (0.25, 0.75), objv = 1.875
#

using Gurobi

env = Gurobi.Env()

model = Gurobi.Model(env, "qp_02")

add_cvars!(model, [1., 1.], 0., Inf)
update_model!(model)

add_qpterms!(model, [1, 1, 2], [1, 2, 2], [2., 1., 1.])
add_constr!(model, [1., 1.], '=', 1.)
update_model!(model)

println(model)

optimize(model)

println("sol = $(get_solution(model))")
println("obj = $(get_objval(model))")

ext: jl


text:
using OneDimensionalCorrelation
using Test

@testset "Elements" begin
    e1 = Element1D(0, 1.2)
    e2 = Element1D(1.2, 2.4)

    @test (e1 < e2) == false
    @test (e2 > e1) == false
    @test e1 <= e2
    @test e2 >= e1

    p = 64
    eg1 = ElementGrid(0, 1.2, p)
    eg2 = ElementGrid(1.2, 2.4, p)

    @test (eg1 < eg2) == false
    @test (eg2 > eg1) == false
    @test eg1 <= eg2
    @test eg2 >= eg1

    eg = ElementGrid(0, 2π, p)

    # Integration accuracy
    @test isapprox( sum(sin.(eg) .* get_weight(eg)), 0.0; atol=1E-12 )

    # Derivative accuracy
    @test all( derivative_matrix(eg) * sin.(eg) ≈ cos.(eg) )

end

@testset "Basis Gauss-Lagrange" begin
    b = Basis(-π, π, 2, 64)

    # Symmetric range
    @test b[begin] == -b[end]

    w = get_weight(b)
    ∇ = derivative_matrix(b)

    # Integral
    @test isapprox( w' *  sin.(b), 0.0; atol=1E-12 )

    # Derivative
    @test ∇ * sin.(b) ≈ cos.(b)

end

@testset "Basis Gauss-Lobatto" begin
    b = BasisLobatto(-2, 2, 3, 32)

    # Symmetric range
    @test b[begin] == -b[end]

    w = get_weight(b)
    ∇ = derivative_matrix(b)

    # Integral
    @test isapprox( w' *  sin.(b), 0.0; atol=1E-12 )

    # Derivative
    @test  ∇ * sin.(b) ≈ cos.(b)

end

@testset "Integrals" begin
    b = BasisLobatto(-5, 5, 2, 64)

    ψ = exp.(-b.^2)
    I = get_identity(b)

    # Norm of Gaussian
    @test bracket(b, ψ, I, ψ) ≈ sqrt(π/2)
    @test bracket(b, ψ, ψ) ≈ sqrt(π/2)

    # Particle in box states
    ϕ1 = particle_in_box(b, 1)
    ϕ4 = particle_in_box(b, 4)
    @test bracket(b, ϕ1, ϕ1) ≈ 1
    @test bracket(b, ϕ4, ϕ4) ≈ 1
    @test bracket(b, ϕ1, ϕ4) + 1 ≈ 1

    # Particle in box Energy
    E(b, n) = 0.5 * n^2 * π^2 / (get_length(b))^2

    d = derivative_matrix(b)
    e_kin = -0.5 * d^2

    # Kinetic energy accuracy
    abs( bracket(b, ϕ1, e_kin, ϕ1) - E(b, 1) ) < 1E-9
    abs( bracket(b, ϕ4, e_kin, ϕ4) - E(b, 4) ) < 1E-9
end

@testset "Hartree Fock" begin
    b = BasisLobatto(-5, 5, 2, 64)
    orbitals = initial_orbitals(b)
    C = coulomb_matrix(b, orbitals)
    K = exchange_matrix(b, orbitals)
    F = fock_matrix(b)

    # Matrices are Hermitian
    @test C[2,1] ≈ C[2,1]
    @test K[2,1] ≈ K[2,1]
    @test F[2,1] ≈ F[2,1]

    V(x) = -2exp(-0.5x^2)
    E, eo, orbitals = solve_hartree_fock(b, V)
    @test E < -2.2123101
end

@testset "CI" begin
    b = BasisLobatto(-10, 10, 6, 12)
    V(x) = -2exp(-0.5x^2)
    E, eo, orbitals = solve_hartree_fock(b, V)

    eb = block_ci(b, orbitals, V)
    @test E > eb["energy"][1]
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
using Compat
using Test
using MPI

@test !MPI.Initialized()
MPI.Init()
@test MPI.Initialized()

@test MPI.Comm(MPI.CComm(MPI.COMM_WORLD)).val == MPI.COMM_WORLD.val

if !Sys.iswindows()
    MPI.finalize_atexit()
    @test MPI.FINALIZE_ATEXIT[]
end

@test !MPI.Finalized()
MPI.Finalize()
@test MPI.Finalized()


text:
num_unlabel_samples = [100,1000,10000,100000]

for i in num_unlabel_samples
    Mat_Label, labels, Mat_Unlabel = loadCircleData(i)
    print("data generated")
    @time label_propagation(Mat_Label, Mat_Unlabel, labels, kernel_type = "knn", knn_num_neighbors = 10, max_iter = 100)
end


text:
"""
    AttentionalDiffusion(;ν1=5.0, ν2=4.0, α=1.0, z=α*.5, θ=.6, σ=.05, Δ=.005)

Constructor for attentional diffusion model object.

# Keywords

- `ν1=5.0`: relative decision value for alternative 1
- `ν2=4.0`: relative decision value for alternative 2
- `α=1.0`: evidence threshold
- `z=0.0`: initial evidence
- `θ=.3`: bias towards attended alternative (lower indicates more bias)
- `σ=.02`: standard deviation of noise in evidence accumulation
- `Δ=.0004`: constant of evidence accumulation speed (evidence per ms)

# References

Krajbich, I., Armel, C., & Rangel, A. (2010). Visual fixations and the computation and comparison of
value in simple choice. Nature neuroscience, 13(10), 1292-1298.

"""
@concrete struct AttentionalDiffusion <: SequentialSamplingModel
    ν1
    ν2
    α
    z
    θ
    σ
    Δ
end

function AttentionalDiffusion(;ν1=5.0, ν2=4.0, α=1.0, z=0.0, θ=.3, σ=.02, Δ=.0004)
    return AttentionalDiffusion(ν1, ν2, α, z, θ, σ, Δ)
end

Broadcast.broadcastable(x::AttentionalDiffusion) = Ref(x)

"""
    rand(dist::AttentionalDiffusion, n_sim::Int, fixation, args...; kwargs...)

Generate `n_sim` simulated trials from the attention diffusion model.

# Arguments

- `dist`: an attentional diffusion model object
- `n_sim::Int`: the number of simulated trials
- `fixation`: a function of the visual fixation process which returns 1 for alternative
    and 2 for alternative 2
- `args...`: optional positional arguments for the `fixation` function

# Keywords

- `kwargs...`: optional keyword arguments for the `fixation` function
"""
function rand(dist::AttentionalDiffusion, n_sim::Int, fixation, args...; kwargs...)
    rts = [Vector{Int}() for _ in 1:2]
    for sim in 1:n_sim
       choice,rt = _rand(dist, ()->fixation(args...; kwargs...))
       push!(rts[choice], rt)
    end
    return rts
end

"""
    rand(dist::AttentionalDiffusion, fixation, args...; kwargs...)

Generate a single simulated trial from the attention diffusion model.

# Arguments

- `dist`: an attentional diffusion model object
- `fixation`: a function of the visual fixation process which returns 1 for alternative
    and 2 for alternative 2
- `args...`: optional positional arguments for the `fixation` function

# Keywords

- `kwargs...`: optional keyword arguments for the `fixation` function
"""
function rand(dist::AttentionalDiffusion, fixation, args...; kwargs...)
    return _rand(dist, ()->fixation(args...; kwargs...))
end

function _rand(dist::AttentionalDiffusion, fixation)
    @unpack α,z,σ,ν1,ν2,θ,Δ = dist
    t = 0
    v = z
    while abs(v) < α
        t += 1
        option = fixation()
        if option == 1
            v += Δ * (ν1 - θ * ν2) + noise(σ)
        elseif option == 2
            v -= Δ * (ν2 - θ * ν1) + noise(σ)
        else
            v += noise(σ)
        end
    end
    choice = (v < α) + 1
    return choice,t
end

noise(σ) = rand(Normal(0, σ))


text:
function twodLinForm( Ff, test, wg )
#= """
  twodLinForm - routine to compute \int{ f*test }

  Author: Jeff Borggaard, Virginia Tech
  Version: 1.0

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    F = twodLinForm( Ff, test, wg )
  ```

  Arguments:
  - `Ff`: Function values at the Gauss points
  - `test`: Matrix of test functions evaluated at the Gauss points
            (dim: n_gauss, n_dof)
  - `wg`: Row vector of Gauss weights

  Output argument:
  - `F`: A numerical approximation to the integrals of phi * f over an element
""" =#

# n_dof = size(test,2)

# F = zeros(n_dof,1)
# for j=1:n_dof
#   F(j) = test(:,j)' * ( wg .* Ff )
# end
# F = sum(test'*(wg.*Ff),2)

  F = test'*(wg.*Ff)

  return F

end


text:
# lets test long approxConv chains

using Test
using IncrementalInference
using Statistics

@testset "zdim size of factor is correct" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousEuclid{3})
addVariable!(fg, :x1, ContinuousEuclid{3})

lr3 =  LinearRelative(MvNormal([1;0;0.0],diagm(ones(3))))
addFactor!(fg, [:x0; :x1], lr3, graphinit=false)

ccw = IIF._getCCW(fg, :x0x1f1)

@test IIF._getZDim(ccw) == 3

res = testFactorResidualBinary( lr3, ContinuousEuclid{3}, ContinuousEuclid{3},
                                zeros(3),[0;0;1.0], ([0;0;0.5],))
#

@test (sum(abs.(res)) - 0.5) < 1e-10

##

end

@testset "approxConv basic and chains" begin

##

fg = generateCanonicalFG_Kaess()

# from a prior to neighbor
pts = approxConv(fg, :x1f1, :x1)

@test Statistics.mean(pts) |> abs  < 0.3
@test 0.5 < Statistics.std(pts) < 1.5

# set a value in graph to start things off
initManual!(fg, :x1, pts)

# legacy case where relative to neighbor
pts = approxConv(fg, :x1x2f1, :x2)

@test Statistics.mean(pts) |> abs  < 0.4
@test 0.7 < Statistics.std(pts) < 2

# along a chain of variables
pts = approxConv(fg, :x1, :x3)

@test Statistics.mean(pts) |> abs  < 0.5
@test 1.3 < Statistics.std(pts) < 2.5

# from a prior down the chain of variables
pts = approxConv(fg, :x1f1, :l2)

@test Statistics.mean(pts) |> abs  < 0.65
@test 1.6 < Statistics.std(pts) < 3.0

##

end

#


text:
include("euler/euler.jl")
using .Primes: prime_sieve
import StatsBase: countmap

function compute(n::Int)::Int
    sieve = prime_sieve(10 ^ (n - 1))
    for (number, is_prime) ∈ enumerate(sieve)
        if !is_prime || number < 10 continue end
        most_common = findmax(countmap(collect(string(number)[1:end - 1])))
        if most_common[1] == 1 && parse(Int, most_common[2]) > 9 - n continue end
        if sum(sieve[map(x->parse(Int, replace(string(number), most_common[2] => x)), parse(Int, most_common[2]):9)]) ≥ n return number end
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# using Cairo
# using Gadfly
# using DataFrames
# using HypergeometricFunctions
# using SpecialFunctions
# using CategoricalArrays
# using CSV

# function weight(j, s, xi)
#     exp(
#         logabsbinomial(j + s - 2, j - 1)[1] -
#         log(j) +
#         s * log(xi) +
#         (j - 1) * log(1 - xi) +
#         log(_₂F₁(j + s - 1, 1, j + 1, 1 - xi))
#     )
# end

# df_j = DataFrame(j = 1:20)
# df_xi = DataFrame(xi = [0.2, 0.4, 0.6])
# df_s = DataFrame(s = [1, 2, 3, 4])
# df = crossjoin(df_j, df_xi, df_s)
# df[!, :wj] = weight.(df[!, :j], df[!, :s], df[!, :xi])
# df[!, :s] = CategoricalArray(df[!, :s])
# CSV.write("weights.csv", df)

# import Pkg;
# Pkg.activate()
# Pkg.add([
#     "AlgebraOfGraphics",
#     "CategoricalArrays",
#     "CairoMakie",
#     "DataFrames",
#     "HypergeometricFunctions",
#     "LaTeXStrings",
#     "SpecialFunctions"
# ])
# Pkg.activate("")
using AlgebraOfGraphics, CategoricalArrays, CairoMakie, DataFrames
using LaTeXStrings, HypergeometricFunctions, SpecialFunctions
set_aog_theme!()
update_theme!(font = "Arial", fontsize = 28)
theme_minimal
function weight(j, s, xi)
    exp(
        logabsbinomial(j + s - 2, j - 1)[1] -
        log(j) +
        s * log(xi) +
        (j - 1) * log(1 - xi) +
        log(_₂F₁(j + s - 1, 1, j + 1, 1 - xi))
    )
end

df_j = DataFrame(j = 1:20);
df_ϕ = DataFrame(ϕ = [0.1, 0.3, 0.5]);
df_s = DataFrame(s = [1, 2, 3, 4]);
df = crossjoin(df_j, df_ϕ, df_s);
df[!, :w] = weight.(df[!, :j], df[!, :s], df[!, :ϕ]);
df[!, :s] = string.(df[!, :s]);
df[!, :ϕ] = string.(df[!, :ϕ]);
df[!, :ϕ] = LaTeXString.(L"\varphi(x_i \cdot β) = " .* df[!, :ϕ]);
plt = data(df) *
    visual(Lines) *
    mapping(:j, :w, col = :ϕ, color = :s)
axis = (xlabel = L"j", ylabel = L"w_j(x)", color = L"s");
fg = draw(plt; axis, figure = (resolution = (900, 450), px_per_unit = 8))
fg
save("figures/final/lineplot-dgsb-weights.svg", fg, px_per_unit = 3)


text:
include("./LatticeQCD.jl")
using .LatticeQCD
using Random
using Dates
#using JLD

#using GFlops

if length(ARGS) == 0
    error("""
    Use input file:
    Like,
    julia run.jl parameters.jl
    """)
end

function runtest()
    run_LQCD(ARGS[1])
    #parameters = parameterloading()
    #run_LQCD(parameters)
end

@time runtest()


text:
module LinearTests

using Interpolations, Base.Test

for D in (Constant, Linear), G in (OnCell, OnGrid)
    ## 1D
    a = rand(5)
    knots = (collect(linspace(1,length(a),length(a))),)
    itp = @inferred(interpolate(knots, a, Gridded(D())))
    @inferred(getindex(itp, 2))
    @inferred(getindex(itp, CartesianIndex((2,))))
    for i = 2:length(a)-1
        @test_approx_eq itp[i] a[i]
        @test_approx_eq itp[CartesianIndex((i,))] a[i]
    end
    @inferred(getindex(itp, knots...))
    @test_approx_eq itp[knots...] a
    # compare scalar indexing and vector indexing
    x = knots[1]+0.1
    v = itp[x]
    for i = 1:length(x)
        @test_approx_eq v[i] itp[x[i]]
    end
    # check the fallback vector indexing
    x = [2.3,2.2]   # non-increasing order
    v = itp[x]
    for i = 1:length(x)
        @test_approx_eq v[i] itp[x[i]]
    end
    # compare against BSpline
    itpb = @inferred(interpolate(a, BSpline(D()), G()))
    for x in linspace(1.1,4.9,101)
        @test_approx_eq itp[x] itpb[x]
    end

    ## 2D
    A = rand(6,5)
    knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2))))
    itp = @inferred(interpolate(knots, A, Gridded(D())))
    @inferred(getindex(itp, 2, 2))
    @inferred(getindex(itp, CartesianIndex((2,2))))
    for j = 2:size(A,2)-1, i = 2:size(A,1)-1
        @test_approx_eq itp[i,j] A[i,j]
        @test_approx_eq itp[CartesianIndex((i,j))] A[i,j]
    end
    @test_approx_eq itp[knots...] A
    @inferred(getindex(itp, knots...))
    # compare scalar indexing and vector indexing
    x, y = knots[1]+0.1, knots[2]+0.6
    v = itp[x,y]
    for j = 1:length(y), i = 1:length(x)
        @test_approx_eq v[i,j] itp[x[i],y[j]]
    end
    # check the fallback vector indexing
    x = [2.3,2.2]   # non-increasing order
    y = [3.5,2.8]
    v = itp[x,y]
    for j = 1:length(y), i = 1:length(x)
        @test_approx_eq v[i,j] itp[x[i],y[j]]
    end
    # compare against BSpline
    itpb = @inferred(interpolate(A, BSpline(D()), G()))
    for y in linspace(1.1,5.9,101), x in linspace(1.1,4.9,101)
        @test_approx_eq itp[x,y] itpb[x,y]
    end

    A = rand(8,20)
    knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
    itp = interpolate(knots, A, Gridded(D()))
    @test_approx_eq itp[4,1.2] A[2,6]
end

end


text:
using Test
using DataFrames
using BrainFlowML
using DSP
using Statistics
using IterTools

function get_gesture(file_name)
    test_filepath = joinpath(BrainFlowML.testdata_path, file_name)
    return BrainFlowML.load_gesture(test_filepath)
end

@testset "BrainFlowML" begin

    @testset "BioData iteration" begin
        bio_data = get_gesture("left_gesture.csv")
        @test bio_data isa BrainFlowML.BioData
        itr = each_channel(bio_data)
        @test itr isa BrainFlowML.ChannelIterator
        @test first(itr) === view(bio_data.raw, :, 1)
        collected_itr = collect(itr)
        @test length(collected_itr) == 8
        @test collected_itr isa Array{Array{Float64,1},1}
    end

    @testset "BioData append" begin
        nchannels = 4
        nsamples = 16

        A = rand(Float64, nsamples, nchannels)
        labels = ones(Int, nsamples)
        bio_data1 = BrainFlowML.BioData(A, 1:nchannels, 600, labels)

        A = rand(Float64, nsamples, nchannels)
        labels = 2*ones(Int, nsamples)
        bio_data2 = BrainFlowML.BioData(A, 1:nchannels, 600, labels)

        append!(bio_data1, bio_data2)
        @test size(bio_data1.raw) == (nsamples*2, nchannels)
        @test length(bio_data1.labels) == nsamples*2
        @test bio_data1.labels[end] == 2
    end

    @testset "widen range" begin
        v = [false, false, false, true, false, false, false]
        BrainFlowML.widen_range!(v, 2)
        @test v == [false, true, true, true, true, true, false]
        BrainFlowML.widen_range!(v, 3) # push it over the edges
        @test v == [true, true, true, true, true, true, true]

        v = [true, false, false]
        BrainFlowML.widen_range!(v, 1)
        @test v == [true, true, false]
    end

    @testset "partition into training data" begin
        nsamples = 10
        sample_size = 4
        step_size = 2
        nchannels = 4
        all_partitions = collect(partition(1:nsamples, sample_size, step_size))
        partition3 = [x for x in all_partitions[3]]
        expected_partitions = length(all_partitions)

        A = rand(Int, nsamples, nchannels)
        X = BrainFlowML.partition_samples(A, sample_size, step_size)
        @test size(X) == (sample_size, nchannels, expected_partitions)
        @test X[:,:,3] == A[partition3, :]

        v = rand(Int, nsamples)
        y = BrainFlowML.partition_labels(v, sample_size, step_size)
        @test length(y) == expected_partitions
        @test y[3] == v[partition3[end]]

        bio_data = BrainFlowML.BioData(A, 1:nchannels, step_size)
        X2, y = BrainFlowML.partition_samples(bio_data, sample_size , step_size)
        @test X2 == X
        @test length(y) == expected_partitions
    end

    @testset "detrending" begin
        A = rand(Float64, 10, 5) .+ 1.0
        BrainFlowML.detrend!(A, 1:4)
        m = [mean(chan) for chan in each_channel(A)]
        for idx = 1:4
            @test isapprox(m[idx], 0, atol=1e-14)
        end
        @test m[5] > 0.5
    end

    @testset "Labeling gestures with DSP" begin
        bio_data = get_gesture("left_gesture.csv")

        BrainFlowML.detrend!(bio_data)
        a_channel = first(each_channel(bio_data))
        @test isapprox(mean(a_channel), 0.0, atol=1e-10)

        h = BrainFlowML.smooth_envelope(a_channel)
        @test h isa AbstractVector

        smooth_data = BrainFlowML.smooth_envelope(bio_data)
        @test smooth_data isa BrainFlowML.BioData

        gesture = 2
        BrainFlowML.label_gesture!(bio_data; gesture=gesture)
        @test gesture ∈ bio_data.labels

        # for manual verification
        # using Plots
        # summed_data = BrainFlowML.sum_channels(smooth_data)
        # plot(summed_data)
        # plot!(bio_data.labels.*maximum(summed_data)/gesture)
    end

    @testset "load labeled gestures" begin
        file_names = ("left_gesture.csv", "right_gesture.csv", "fist_gesture.csv", "spread_gesture.csv")
        file_paths = [joinpath(BrainFlowML.testdata_path, name) for name in file_names]
        bio_data = BrainFlowML.load_labeled_gestures(file_paths)
        @test unique(bio_data.labels) == [0, 1, 2, 3, 4]
    end

    @testset "preprocessing & prediction" begin
        using JLD2
        using DecisionTree

        test_file = joinpath(BrainFlowML.testdata_path, "fist_gesture.csv")
        bio_data = BrainFlowML.load_labeled_gesture(test_file, 1)

        gesture_start = findfirst(x -> x > 0, diff(bio_data.labels))
        gesture_end = findfirst(x -> x < 0, diff(bio_data.labels))
        halfway_gesture = Int( ceil(gesture_start + (gesture_end-gesture_start)/2) )

        modeled_sample_size = 128

        sample_slice = halfway_gesture:halfway_gesture+modeled_sample_size-1
        A = bio_data.raw[sample_slice, :]

        nchannels = length(bio_data.channels)
        v = BrainFlowML.preprocess_brainflow_data(A, bio_data.channels) # slow function, takes a millisecond due to smoothening
        @test length(v) == nchannels*modeled_sample_size

        test_model = joinpath(BrainFlowML.testdata_path, "model_file.jld2")
        @load test_model model

        prediction = apply_forest(model, v')
        @test prediction == [3]
    end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [
    ((2,), :SOCExpPSD),
    ]
insts["fast"] = [
    ((3,), :SOCExpPSD),
    ]
insts["various"] = [
    ((3,), :SOCExpPSD),
    ((4,), :SOCExpPSD, relaxed_tols),
    ((5,), :SOCExpPSD, relaxed_tols),
    ]
return (LotkaVolterraJuMP, insts)


text:
using Agents.Pathfinding

moore = Pathfinding.moore_neighborhood(2)
vonneumann = Pathfinding.vonneumann_neighborhood(2)
@testset "constructors" begin
    cost = GridSpace((5, 5); pathfinder = Pathfinding.Pathfinder()).pathfinder.cost_metric
    @test typeof(cost) <: DirectDistance{2}
    @test cost.direction_costs == [10, 14]
    cost =
        GridSpace(
            (5, 5);
            pathfinder = Pathfinding.Pathfinder(cost_metric = DirectDistance{2}()),
        ).pathfinder.cost_metric
    @test_throws AssertionError GridSpace(
        (5, 5);
        pathfinder = Pathfinding.Pathfinder(cost_metric = DirectDistance{2}([1])),
    )
    @test_throws AssertionError GridSpace(
        (5, 5);
        pathfinder = Pathfinding.Pathfinder(
            diagonal_movement = false,
            cost_metric = DirectDistance{2}([]),
        ),
    )
    @test typeof(cost) <: DirectDistance{2}
    @test cost.direction_costs == [10, 14]
    cost =
        GridSpace(
            (5, 5);
            pathfinder = Pathfinding.Pathfinder(cost_metric = MaxDistance{2}()),
        ).pathfinder.cost_metric
    @test typeof(cost) <: MaxDistance{2}
    @test_throws MethodError GridSpace(
        (5, 5);
        pathfinder = Pathfinding.Pathfinder(cost_metric = HeightMap),
    )
    @test_throws AssertionError GridSpace(
        (5, 5);
        pathfinder = Pathfinding.Pathfinder(cost_metric = HeightMap([1 1])),
    )
    cost =
        GridSpace(
            (5, 5);
            pathfinder = Pathfinding.Pathfinder(cost_metric = HeightMap(fill(1, 5, 5))),
        ).pathfinder.cost_metric
    @test typeof(cost) <: HeightMap{2}
    @test typeof(cost.base_metric) <: DirectDistance{2}
    @test cost.hmap == fill(1, 5, 5)
    cost =
        GridSpace(
            (5, 5);
            pathfinder = Pathfinding.Pathfinder(
                cost_metric = HeightMap(fill(1, 5, 5), MaxDistance{2}()),
            ),
        ).pathfinder.cost_metric
    @test typeof(cost) <: HeightMap{2}
    @test typeof(cost.base_metric) <: MaxDistance{2}
    @test cost.hmap == fill(1, 5, 5)
    hmap = zeros(Int, 1, 1, 1)
    @test_throws MethodError GridSpace(
        (5, 5);
        pathfinder = Pathfinding.Pathfinder(cost_metric = HeightMap(hmap)),
    )

    space = GridSpace((5, 5); pathfinder = Pathfinding.Pathfinder())
    model = ABM(Agent3, space)
    a = add_agent!((5, 2), model, 654.5)
    @test is_stationary(a, model)
    set_target!(a, (1, 3), model)
    @test !is_stationary(a, model)
    @test length(model.space.pathfinder.agent_paths) == 1
    kill_agent!(a, model)
    @test length(model.space.pathfinder.agent_paths) == 0
    @test heightmap(model) === nothing

    hmap = fill(1, 5, 5)
    space = GridSpace((5, 5); pathfinder = Pathfinding.Pathfinder(cost_metric = HeightMap(hmap)))
    model = ABM(Agent3, space)
    @test heightmap(model) == hmap
end

@testset "metrics" begin
    pfinder_2d_np_m = Pathfinding.AStar{2,false,true}(
        Dict(),
        (10, 10),
        copy(moore),
        0.0,
        trues(10, 10),
        DirectDistance{2}(),
    )
    pfinder_2d_np_nm = Pathfinding.AStar{2,false,false}(
        Dict(),
        (10, 10),
        copy(vonneumann),
        0.0,
        trues(10, 10),
        DirectDistance{2}(),
    )
    pfinder_2d_p_m = Pathfinding.AStar{2,true,true}(
        Dict(),
        (10, 10),
        copy(moore),
        0.0,
        trues(10, 10),
        DirectDistance{2}(),
    )
    pfinder_2d_p_nm = Pathfinding.AStar{2,true,false}(
        Dict(),
        (10, 10),
        copy(vonneumann),
        0.0,
        trues(10, 10),
        DirectDistance{2}(),
    )
    hmap = fill(0, 10, 10)
    hmap[:, 6] .= 100
    hmap[1, 6] = 0

    @test delta_cost(pfinder_2d_np_m, DirectDistance{2}(), (1, 1), (4, 6)) == 62
    @test delta_cost(pfinder_2d_p_m, DirectDistance{2}(), (1, 1), (8, 6)) == 62
    @test delta_cost(pfinder_2d_np_nm, DirectDistance{2}(), (1, 1), (4, 6)) == 80
    @test delta_cost(pfinder_2d_p_nm, DirectDistance{2}(), (1, 1), (8, 6)) == 80

    @test delta_cost(pfinder_2d_np_m, MaxDistance{2}(), (1, 1), (4, 6)) == 5
    @test delta_cost(pfinder_2d_p_m, MaxDistance{2}(), (1, 1), (8, 6)) == 5
    @test delta_cost(pfinder_2d_np_nm, MaxDistance{2}(), (1, 1), (4, 6)) == 5
    @test delta_cost(pfinder_2d_p_nm, MaxDistance{2}(), (1, 1), (8, 6)) == 5

    @test delta_cost(pfinder_2d_np_m, HeightMap(hmap), (1, 1), (4, 6)) == 162
    @test delta_cost(pfinder_2d_p_m, HeightMap(hmap), (1, 1), (8, 6)) == 162
    @test delta_cost(pfinder_2d_np_nm, HeightMap(hmap), (1, 1), (4, 6)) == 180
    @test delta_cost(pfinder_2d_p_nm, HeightMap(hmap), (1, 1), (8, 6)) == 180
end

@testset "pathing" begin
    wlk = trues(7, 6)
    wlk[2:7, 1] .= false
    wlk[7, 3:6] .= false
    wlk[[2:4; 6], 4] .= false
    wlk[2:5, 5] .= false
    wlk[2, 2] = false
    wlk[4, 3] = false
    wlk[5, 3] = false

    pfinder_2d_np_m = Pathfinding.AStar{2,false,true}(
        Dict(),
        (7, 6),
        copy(moore),
        0.0,
        wlk,
        DirectDistance{2}(),
    )
    pfinder_2d_np_nm = Pathfinding.AStar{2,false,false}(
        Dict(),
        (7, 6),
        copy(vonneumann),
        0.0,
        wlk,
        DirectDistance{2}(),
    )
    pfinder_2d_p_m = Pathfinding.AStar{2,true,true}(
        Dict(),
        (7, 6),
        copy(moore),
        0.0,
        wlk,
        DirectDistance{2}(),
    )
    pfinder_2d_p_nm = Pathfinding.AStar{2,true,false}(
        Dict(),
        (7, 6),
        copy(vonneumann),
        0.0,
        wlk,
        DirectDistance{2}(),
    )

    p = collect(Pathfinding.find_path(pfinder_2d_np_m, (1, 1), (6, 6)))
    @test p == [(1, 2), (1, 3), (1, 4), (1, 5), (2, 6), (3, 6), (4, 6), (5, 6), (6, 6)]
    p = collect(Pathfinding.find_path(pfinder_2d_np_nm, (1, 1), (6, 6)))
    @test p ==
          [(1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (2, 6), (3, 6), (4, 6), (5, 6), (6, 6)]
    p = collect(Pathfinding.find_path(pfinder_2d_p_m, (1, 1), (6, 6)))
    @test p == [(2, 6), (3, 6), (4, 6), (5, 6), (6, 6)]
    p = collect(Pathfinding.find_path(pfinder_2d_p_nm, (1, 1), (6, 6)))
    @test p == [(1, 6), (2, 6), (3, 6), (4, 6), (5, 6), (6, 6)]
end


text:
module DevTools

include("get_initial_compat.jl")
include("explicit_import.jl")
include("get_single_src.jl")

end # module


text:
module NoteSequences

using MIDI
using DataStructures

include("constants.jl")
include("instruments.jl")
include("notesequence.jl")
include("utils.jl")
include("PerformanceRepr/PerformanceRepr.jl")
include("MelodyRepr/MelodyRepr.jl")

using .PerformanceRepr
export Performance, PerformanceEvent

using .MelodyRepr
export Melody

export encode_event, decode_event

end


text:
using .CUDA

function set_device!()
    lcomm = MPI.Comm_split_type(COMM_WORLD, MPI.MPI_COMM_TYPE_SHARED, MPI.Comm_rank(COMM_WORLD))
    CUDA.device!(MPI.Comm_rank(lcomm) % length(devices()))
end

function sync()
    synchronize()
end