text:
using Test
import Roots.newton, Roots.halley

@testset "Test Newton, Halley, Schroder methods" begin

    @test abs(newton(sin, cos, 0.5) - 0.0) <= 100*eps(1.0)
    @test newton(cos, x -> -sin(x), 1.0)  ≈ pi/2
    @test newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0)  ≈ 2.414213562373095
    @test abs(newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14
    @test halley(x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2, 3.0)  ≈ 2.414213562373095
    a = halley(x -> exp(x) - cos(x),
               x -> exp(x) + sin(x),
               x -> exp(x) + cos(x), 3.0)
    @test abs(a - 0.0) <= 1e-14

    # find_zero calls
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2), 3.0, Roots.Newton())  ≈ 2.414213562373095
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2), 3.0, Roots.Halley())  ≈ 2.414213562373095
    @test find_zero((x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2), 3.0, Roots.Schroder())  ≈ 2.414213562373095

    ## test with Complex input

    @test real(Roots.newton(x ->  x^3 - 1, x ->  3x^2, 1+im)) ≈ 1.0
    @test real(Roots.newton(x ->  x^3 - 1, x ->  3x^2, 1+10im)) ≈ (-1/2)

    ## Issue #143 test with new interface
    Roots.newton(sin, cos, 3.0) ≈ π # uses find_zero
    Roots.newton((sin,cos), 3.0) ≈ π # uses simple

    fdf = x -> (sin(x), sin(x)/cos(x))  # (f, f/f')
    @test Roots.find_zero(fdf, 3.0, Roots.Newton())  ≈ π # uses find_zero
    Roots.newton(fdf, 3.0)  ≈ π # uses simple

    fdfdf = x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x), sin(x)/cos(x))   # (f, f/f', f'/f'', f''/f''')
    @test Roots.find_zero(fdfdf, 3.0, Roots.Halley()) ≈ π

    # check that functions with multiple return values can work with other
    # methods
    for M in [Roots.Schroder(), Roots.Halley(), Roots.Newton(), Roots.Order1(), Roots.Order0()]
        @test Roots.find_zero(fdfdf, 3.0, M) ≈ π  # can pass function to others
    end

    for M in [Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi()]
        @test Roots.find_zero(fdfdf, (3.0, 4.0), M) ≈ π  # can pass function to others
    end

    @test find_zero(x -> (x^2 -2, (x^2-2)/2x), 1.0, Roots.Newton(), Roots.Bisection()) ≈ sqrt(2)

    @test_throws Roots.ConvergenceFailed Roots.newton((x->x^2 + 1, x->2x), 0)

end

@testset "Lith Boonkkamp IJzerman methods" begin

    x₀, x̃₀, α = 1.0, 1.1, 1.1673039782614187
    f(x) = x^5 - x - 1
    fp(x) = 5x^4 - 1
    fpp(x) = 20x^3
    fppp(x) = 60x
    fpppp(x) = 60

    @test solve(ZeroProblem((f,), x₀), Roots.LithBoonkkampIJzerman(3,0)) ≈ α
    @test solve(ZeroProblem(f, x₀), Roots.LithBoonkkampIJzerman(3,0)) ≈ α
    @test solve(ZeroProblem((f,), x₀), Roots.LithBoonkkampIJzerman(4,0)) ≈ α
    @test solve(ZeroProblem(f, x₀), Roots.LithBoonkkampIJzerman(4,0)) ≈ α

    @test solve(ZeroProblem((f,fp), x₀), Roots.LithBoonkkampIJzerman(2,1)) ≈ α
    @test solve(ZeroProblem((f,fp), x₀), Roots.LithBoonkkampIJzerman(3,1)) ≈ α

    @test solve(ZeroProblem((f,fp, fpp), x₀), Roots.LithBoonkkampIJzerman(1,2)) ≈ α
    @test solve(ZeroProblem((f,fp, fpp), x₀), Roots.LithBoonkkampIJzerman(2,2)) ≈ α

    @test solve(ZeroProblem((f,fp, fpp, fppp), x₀), Roots.LithBoonkkampIJzerman(1,3)) ≈ α
    @test solve(ZeroProblem((f,fp, fpp, fppp), x₀), Roots.LithBoonkkampIJzerman(2,3)) ≈ α

    @test solve(ZeroProblem((f,fp, fpp, fppp, fpppp), x₀), Roots.LithBoonkkampIJzerman(1,4)) ≈ α
    @test solve(ZeroProblem((f,fp, fpp, fppp, fpppp), x̃₀), Roots.LithBoonkkampIJzerman(2,4)) ≈ α # needs closer

    # bracketing
    @test solve(ZeroProblem((f,fp), (1,2)), Roots.LithBoonkkampIJzermanBracket()) ≈ α

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module APITests

using Pkg, Test
import Pkg.Types.PkgError

include("utils.jl")

@testset "API should accept `AbstractString` arguments" begin
    temp_pkg_dir() do project_path
        with_temp_env() do
            Pkg.add(strip("  Example  "))
            Pkg.rm(strip("  Example "))
        end
    end
end

@testset "Pkg.activate" begin
    temp_pkg_dir() do project_path
        cd_tempdir() do tmp
            path = pwd()
            Pkg.activate(".")
            mkdir("Foo")
            cd(mkdir("modules")) do
                Pkg.generate("Foo")
            end
            Pkg.develop(Pkg.PackageSpec(path="modules/Foo")) # to avoid issue #542
            Pkg.activate("Foo") # activate path Foo over deps Foo
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            rm("Foo"; force=true, recursive=true)
            Pkg.activate("Foo") # activate path from developed Foo
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("./Foo") # activate empty directory Foo (sidestep the developed Foo)
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("Bar") # activate empty directory Bar
            @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
            Pkg.activate(".")
            Pkg.add("Example") # non-deved deps should not be activated
            Pkg.activate("Example")
            @test Base.active_project() == joinpath(path, "Example", "Project.toml")
            Pkg.activate(".")
            cd(mkdir("tests"))
            Pkg.activate("Foo") # activate developed Foo from another directory
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate() # activate home project
            @test Base.ACTIVE_PROJECT[] === nothing
        end
    end
end

@testset "Pkg.status" begin
    temp_pkg_dir() do project_path
        Pkg.add(["Example", "Random"])
        Pkg.status()
        Pkg.status("Example")
        Pkg.status(["Example", "Random"])
        Pkg.status(PackageSpec("Example"))
        Pkg.status(PackageSpec(uuid = "7876af07-990d-54b4-ab0e-23690620f79a"))
        Pkg.status(PackageSpec.(["Example", "Random"]))
        Pkg.status(; mode=PKGMODE_MANIFEST)
        Pkg.status("Example"; mode=PKGMODE_MANIFEST)
        @test_deprecated Pkg.status(PKGMODE_MANIFEST)
        # issue #1183: Test exist in manifest but not in project
        Pkg.status("Test"; mode=PKGMODE_MANIFEST)
        @test_throws PkgError Pkg.status("Test"; mode=Pkg.Types.PKGMODE_COMBINED)
        @test_throws PkgError Pkg.status("Test"; mode=PKGMODE_PROJECT)
    end
end

@testset "Pkg.develop" begin
    temp_pkg_dir() do project_path; cd_tempdir() do tmpdir;
        exuuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a") # UUID of Example.jl
        entry = nothing
        # explicit relative path
        with_temp_env() do env_path
            cd(env_path) do
                uuids = Pkg.generate("Foo")
                Pkg.develop(PackageSpec(;path="Foo"))
                manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
                entry = manifest[uuids["Foo"]]
            end
            @test entry.path == "Foo"
            @test entry.name == "Foo"
            @test isdir(joinpath(env_path, entry.path))
        end
        # explicit absolute path
        with_temp_env() do env_path
            cd_tempdir() do temp_dir
                uuids = Pkg.generate("Foo")
                absolute_path = abspath(joinpath(temp_dir, "Foo"))
                Pkg.develop(PackageSpec(;path=absolute_path))
                manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
                entry = manifest[uuids["Foo"]]
                @test entry.name == "Foo"
                @test entry.path == absolute_path
                @test isdir(entry.path)
            end
        end
        # name
        with_temp_env() do env_path
            Pkg.develop("Example")
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Example")
                    @test isdir(entry.path)
                end
            end
        end
        # name + uuid
        with_temp_env() do env_path
            Pkg.develop(PackageSpec(name = "Example", uuid = exuuid))
            @test Pkg.Types.Context().env.manifest[exuuid].version > v"0.5"
        end
        # uuid
        with_temp_env() do env_path
            Pkg.develop(PackageSpec(uuid = exuuid))
            @test Pkg.Types.Context().env.manifest[exuuid].version > v"0.5"
        end
        # name + local
        with_temp_env() do env_path
            Pkg.develop("Example"; shared=false)
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath("dev", "Example")
                    @test isdir(joinpath(env_path, entry.path))
                end
            end
        end
        # url
        with_temp_env() do env_path
            url = "https://github.com/JuliaLang/Example.jl"
            Pkg.develop(PackageSpec(;url=url))
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Example")
                    @test isdir(entry.path)
                end
            end
        end
        # unregistered url
        with_temp_env() do env_path
            url = "https://github.com/00vareladavid/Unregistered.jl"
            Pkg.develop(PackageSpec(;url=url))
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Unregistered"
                    @test uuid == UUID("dcb67f36-efa0-11e8-0cef-2fc465ed98ae")
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Unregistered")
                    @test isdir(entry.path)
                end
            end
        end
        # with rev
        with_temp_env() do env_path
            @test_throws PkgError Pkg.develop(PackageSpec(;name="Example",rev="Foobar"))
        end
    end end
end

@testset "Pkg.add" begin
    exuuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a") # UUID of Example.jl
    # Add by version should override add by repo
    temp_pkg_dir() do project_path; with_temp_env() do env_path
        Pkg.add(Pkg.PackageSpec(;name="Example", rev="master"))
        Pkg.add(Pkg.PackageSpec(;name="Example", version="0.3.0"))
        @test Pkg.Types.Context().env.manifest[exuuid].version == v"0.3.0"
        @test Pkg.Types.Context().env.manifest[exuuid].repo == Pkg.Types.GitRepo()
    end end
    # Add by version should override add by repo, even for indirect dependencies
    temp_pkg_dir() do project_path; mktempdir() do tempdir; with_temp_env() do
        path = git_init_package(tempdir, joinpath(@__DIR__, "test_packages", "DependsOnExample"))
        Pkg.add(Pkg.PackageSpec(;path=path))
        Pkg.add(Pkg.PackageSpec(;name="Example", rev="master"))
        Pkg.rm("Example")
        # Now `Example` should be tracking a repo and it is in the dep graph
        # But `Example` is *not* a direct dependency
        Pkg.add(Pkg.PackageSpec(;name="Example", version="0.3.0"))
        @test Pkg.Types.Context().env.manifest[exuuid].version == v"0.3.0"
        @test Pkg.Types.Context().env.manifest[exuuid].repo == Pkg.Types.GitRepo()
    end end end
    # Add by URL should not override pin
    temp_pkg_dir() do project_path; with_temp_env() do env_path
        Pkg.add(Pkg.PackageSpec(;name="Example", version="0.3.0"))
        Pkg.pin(Pkg.PackageSpec(;name="Example"))
        a = deepcopy(Pkg.Types.EnvCache().manifest)
        Pkg.add(Pkg.PackageSpec(;url="https://github.com/JuliaLang/Example.jl"))
        b = Pkg.Types.EnvCache().manifest
        for (uuid, x) in a
            y = b[uuid]
            for property in propertynames(x)
                @test getproperty(x, property) == getproperty(y, property)
            end
        end
    end end
    # Add by URL should not overwrite files
    temp_pkg_dir() do project_path; with_temp_env() do env_path
        Pkg.add(Pkg.PackageSpec(;url="https://github.com/JuliaLang/Example.jl"))
        t1, t2 = nothing, nothing
        for (uuid, entry) in Pkg.Types.EnvCache().manifest
            entry.name == "Example" || continue
            t1 = mtime(Pkg.Operations.find_installed(entry.name, uuid, entry.tree_hash))
        end
        Pkg.add(Pkg.PackageSpec(;url="https://github.com/JuliaLang/Example.jl"))
        for (uuid, entry) in Pkg.Types.EnvCache().manifest
            entry.name == "Example" || continue
            t2 = mtime(Pkg.Operations.find_installed(entry.name, uuid, entry.tree_hash))
        end
        @test t1 == t2
    end end
end

@testset "Pkg.free" begin
    temp_pkg_dir() do project_path
        # Assumes that `TOML` is a registered package name
        # Can not free an un-`dev`ed un-`pin`ed package
        with_temp_env() do; mktempdir() do tempdir;
            p = git_init_package(tempdir, joinpath(@__DIR__, "test_packages", "TOML"))
            Pkg.add(Pkg.PackageSpec(;path=p))
            @test_throws PkgError Pkg.free("TOML")
        end end
        # Can free a registered package that is tracking a repo
        with_temp_env() do
            exuuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a") # UUID of Example.jl
            Pkg.add(Pkg.PackageSpec(name = "Example", rev="c37b675")) # same commit as release v0.5.1
            @test Pkg.Types.Context().env.manifest[exuuid].repo.rev == "c37b675"
            Pkg.free("Example") # should not throw, see issue #1142
            @test Pkg.Types.Context().env.manifest[exuuid].repo.rev == nothing
            @test Pkg.Types.Context().env.manifest[exuuid].version > v"0.5"
        end
        # Can not free an unregistered package
        with_temp_env() do;
            Pkg.develop(Pkg.PackageSpec(;url="https://github.com/00vareladavid/Unregistered.jl"))
            @test_throws PkgError Pkg.free("Unregistered")
        end
    end
end

@testset "Pkg.free" begin
    # `pin` should not override previous `pin`
    temp_pkg_dir() do project_path; with_temp_env() do env_path
        Pkg.add(Pkg.PackageSpec(;name="Example"))
        Pkg.pin(Pkg.PackageSpec(;name="Example", version="0.3.0"))
        @test_throws PkgError Pkg.pin(Pkg.PackageSpec(;name="Example", version="0.4.0"))
    end end
end

end # module APITests

ext: jl


text:
write_test_output = false

path = dirname(@__FILE__)

###################################################################
# Test: initial_draw!()
###################################################################
m = AnSchorfheide()

function loglik_fn(para::ModelConstructors.ParameterVector, data::Matrix{Float64})::Float64
    DSGE.update!(m, para)
    DSGE.likelihood(m, data; sampler = false, catch_errors = true,
               use_chand_recursion = get_setting(m, :use_chand_recursion), verbose = :low)
end

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = h5read("reference/smc.h5", "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, false)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:tempering_target, 0.95)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:use_fixed_schedule, true)

####################################################################
init_cloud = Cloud(length(m.parameters), get_setting(m, :n_particles))

@everywhere Random.seed!(42)
SMC.initial_draw!(loglik_fn, m.parameters, data, init_cloud)

if write_test_output
    JLD2.jldopen("reference/initial_draw_out.jld2", "w") do file
        write(file, "cloud", init_cloud)
    end
end

saved_init_cloud = load("reference/initial_draw_out.jld2", "cloud")

@testset "Initial draw" begin
    @test @test_matrix_approx_eq SMC.get_vals(init_cloud) SMC.get_vals(saved_init_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(init_cloud) SMC.get_loglh(saved_init_cloud)
end

###################################################################
# Test: one_draw()
###################################################################
m = AnSchorfheide()

file = JLD2.jldopen("reference/one_draw_in.jld2", "r")
    parameters = file["parameters"]
    data = file["data"]
close(file)

draw = SMC.one_draw(loglik_fn, parameters, data)

if write_test_output
    JLD2.jldopen("reference/one_draw_out.jld2", true, true, true, IOStream) do file
        file["draw"] = draw
    end
end

test_draw = JLD2.jldopen("reference/one_draw_out.jld2", "r") do file
    file["draw"]
end

###################################################################
@testset "One draw" begin
    @test draw[1] == test_draw[1]
    @test draw[2] ≈ test_draw[2]
end

###################################################################
# Test: draw_likelihood()
###################################################################
draw_lik = SMC.draw_likelihood(loglik_fn, parameters, data, vec(draw[1]))

if write_test_output
    JLD2.jldopen("reference/draw_likelihood_out.jld2", true, true, true, IOStream) do file
        file["draw_lik"] = draw_lik
    end
end
test_draw_lik = JLD2.jldopen("reference/draw_likelihood_out.jld2", "r") do file
    file["draw_lik"]
end

###################################################################
@testset "Draw Likelihood" begin
    @test draw_lik[1][1] ≈ test_draw_lik[1][1]
    @test draw_lik[2] == test_draw_lik[2]
end

###################################################################
# Test: initialize_likelihoods!()
###################################################################
SMC.initialize_likelihoods!(loglik_fn, parameters, data, init_cloud)

if write_test_output
    JLD2.jldopen("reference/initialize_likelihood_out.jld2", true, true, true, IOStream) do file
        file["init_lik_cloud"] = init_cloud
    end
end
test_init_cloud = JLD2.jldopen("reference/initialize_likelihood_out.jld2", "r") do file
    file["init_lik_cloud"]
end

###################################################################
@testset "Initialize Likelihoods" begin
    @test @test_matrix_approx_eq SMC.get_vals(init_cloud) SMC.get_vals(test_init_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(init_cloud) SMC.get_loglh(test_init_cloud)
end

###################################################################
# Test: initialize_cloud_settings!()
###################################################################
SMC.initialize_cloud_settings!(init_cloud)

if write_test_output
    JLD2.jldopen("reference/initialize_cloud_settings.jld2", true, true, true, IOStream) do file
        file["init_cloud"] = init_cloud
    end
end
test_init_cloud = JLD2.jldopen("reference/initialize_cloud_settings.jld2", "r") do file
    file["init_cloud"]
end

###################################################################
@testset "Initialize Cloud" begin
    @test init_cloud.ESS         == test_init_cloud.ESS
    @test init_cloud.stage_index == test_init_cloud.stage_index
    @test init_cloud.n_Φ         == test_init_cloud.n_Φ
    @test init_cloud.resamples   == test_init_cloud.resamples
    @test init_cloud.c           == test_init_cloud.c
    @test init_cloud.total_sampling_time == test_init_cloud.total_sampling_time
    @test init_cloud.tempering_schedule  == test_init_cloud.tempering_schedule
end


text:
__precompile__()

module HumanoidLCMSim

export # types
    HumanoidRobotInfo,
    LCMControlReceiver,
    LCMControlPublisher,
    Actuator,
    Side

export # functions
    parse_actuators,
    handle

export # submodules
    Sides,
    AtlasSim,
    AtlasControl

using RigidBodySim
using RigidBodyDynamics
using OrdinaryDiffEq
using DiffEqCallbacks
using StaticArrays
using Rotations
using LCMCore
using FastIOBuffers
using LightXML
using BotCoreLCMTypes

using Dates: Second

include("gains.jl")
include("sides.jl")

using .Sides

include("util.jl")
include("humanoid_robot_info.jl")
include("lcm_interop.jl")
include("controlreceiver.jl")
include("controlpublisher.jl")

include("atlassim.jl")
include("atlascontrol.jl")

end # module


text:
using SymbolicUtils: Sym, FnType, Term, symtype, Contextual, EmptyCtx
using SymbolicUtils
using Test

@testset "@syms" begin
    let
        @syms a b::Float64 f(::Real) g(p, h(q::Real))::Int

        @test a isa Sym{Number}
        @test a.name === :a

        @test b isa Sym{Float64}
        @test b.name === :b

        @test f isa Sym{FnType{Tuple{Real}, Number}}
        @test f.name === :f

        @test g isa Sym{FnType{Tuple{Number, FnType{Tuple{Real}, Number}}, Int}}
        @test g.name === :g

        @test f(b) isa Term
        @test symtype(f(b)) === Number
        @test_throws ErrorException f(a)

        @test g(b, f) isa Term
        @test_throws ErrorException g(b, a)

        @test symtype(g(b, f)) === Int
    end
end

@testset "hashing" begin
    @syms a b f(x, y)
    @test hash(a) == hash(a)
    @test hash(a) != hash(b)
    @test hash(a+1) == hash(a+1)
    @test hash(sin(a+1)) == hash(sin(a+1))
    @test hash(f(1,a)) == hash(f(1, a))

    c = a
    g = f
    @syms a f(x, y)
    @test hash(a) == hash(c)
    @test hash(g(a, b)) == hash(f(a,b))
    @test hash(f(a, b)) == hash(f(c,b))
    @test hash(sin(a+1)) == hash(sin(c+1))
end

@testset "Base methods" begin
    @syms w::Complex z::Complex a::Real b::Real x

    @test isequal(w + z, Term{Complex}(+, [w, z]))
    @test isequal(z + a, Term{Number}(+, [z, a]))
    @test isequal(a + b, Term{Real}(+, [a, b]))
    @test isequal(a + x, Term{Number}(+, [a, x]))
    @test isequal(a + z, Term{Number}(+, [a, z]))

    # promote_symtype of identity
    @test isequal(Term(identity, [w]), Term{Complex}(identity, [w]))
    @test isequal(+(w), w)
    @test isequal(+(a), a)

    @test isequal(rem2pi(a, RoundNearest), Term{Real}(rem2pi, [a, RoundNearest]))

    # bool
    for f in [(==), (!=), (<=), (>=), (<), (>)]
        @test isequal(f(a, 0), Term{Bool}(f, [a, 0]))
        @test isequal(f(0, a), Term{Bool}(f, [0, a]))
        @test isequal(f(a, a), Term{Bool}(f, [a, a]))
    end

    @test symtype(cond(true, 4, 5)) == Int
    @test symtype(cond(a < 0, b, w)) == Union{Real, Complex}
    @test_throws MethodError w < 0
    @test isequal(w == 0, Term{Bool}(==, [w, 0]))
end

@testset "err test" begin
    @syms t()
    @test_throws ErrorException t(2)
end

@testset "Contexts" begin
    @syms a b c

    @test @rule(~x::Contextual((x, ctx) -> ctx==EmptyCtx()) => "yes")(1) == "yes"
    @test @rule(~x::Contextual((x, ctx) -> haskey(ctx, x)) => true)(a, Dict(a=>1))
    @test @rule(~x::Contextual((x, ctx) -> haskey(ctx, x)) => true)(b, Dict(a=>1)) === nothing
    @test_throws UndefVarError @rule(~x => __CTX__)(a, "test")
    @test @rule(~x => @ctx)(a, "test") == "test"
    @test @rule(~x::Contextual((x, ctx) -> haskey(ctx, x)) => (@ctx)[~x])(a, Dict(a=>1)) === 1
    @test @rule(~x::Contextual((x, ctx) -> haskey(ctx, x)) => (@ctx)[~x])(b, Dict(a=>1)) === nothing
    @test simplify(a+a, "test", rules=RuleSet([@rule ~x => @ctx])) == "test"
end

@testset "substitute" begin
    @syms a b
    @test substitute(a, Dict(a=>1)) == 1
    @test isequal(substitute(sin(a+b), Dict(a=>1)), sin(1+b))
    @test substitute(a+b, Dict(a=>1, b=>3)) |> simplify == 4
end

@testset "printing" begin
    @syms a b
    @test repr(a+b) == "a + b"
    @test repr(-a) == "-a"
    @test repr(-a + 3) == "(-a) + 3"
end


text:
using LinearAlgebra: axpy!

# Helper function to compute the G_nj factors for the classical ExpRK methods
@inline _compute_nl(f::SplitFunction, u, p, t, A) = f.f2(u, p, t)
@inline _compute_nl(f, u, p, t, A) = f(u, p, t) - A * u
@inline _compute_nl!(G, f::SplitFunction, u, p, t, A, Au_cache) = f.f2(G, u, p, t)
@inline function _compute_nl!(G, f, u, p, t, A, Au_cache)
  f(G, u, p, t)
  mul!(Au_cache, A, u)
  G .-= Au_cache
end

##########################################
# Common initializers for ExpRK integrators
function initialize!(integrator, cache::ExpRKConstantCache)
  # Pre-start fsal
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.fsallast = zero(integrator.fsalfirst)

  # Initialize interpolation derivatives
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end
function initialize!(integrator, cache::ExpRKCache)
  # Pre-start fsal
  integrator.fsalfirst = zero(cache.rtmp)
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.fsallast = zero(integrator.fsalfirst)

  # Initialize interpolation derivatives
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

###########################################
# Classical ExpRK integrators
function perform_step!(integrator, cache::LawsonEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  nl = _compute_nl(f, uprev, p, t, A)
  if isa(f, SplitFunction)
    integrator.destats.nf2 += 1
  else
    integrator.destats.nf += 1
  end
  @muladd v = uprev + dt * nl
  if alg.krylov
    u = expv(dt, A, v; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  else
    exphA = cache.ops
    u = exphA * v
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::LawsonEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,G,J,exphA,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  _compute_nl!(G, f, uprev, p, t, A, rtmp)
  if isa(f, SplitFunction)
    integrator.destats.nf2 += 1
  else
    integrator.destats.nf += 1
  end
  @muladd @.. broadcast=false tmp = uprev + dt*G
  if alg.krylov
    Ks, expv_cache = KsCache
    arnoldi!(Ks, A, tmp; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    expv!(u,dt,Ks; cache=expv_cache)
  else
    mul!(u,exphA,tmp)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::NorsettEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  if alg.krylov
    w = phiv(dt, A, integrator.fsalfirst, 1; m=min(alg.m, size(A,1)),
      opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    u = uprev + dt * w[:,2]
  else
    phihA = cache.ops
    u = uprev + dt * (phihA * integrator.fsalfirst)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::NorsettEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack rtmp,J,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  if alg.krylov
    Ks, phiv_cache, ws = KsCache; w = ws[1]
    arnoldi!(Ks, A, integrator.fsalfirst; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    phiv!(w, dt, Ks, 1; cache=phiv_cache)
    @muladd @.. broadcast=false u = uprev + dt * @view(w[:, 2])
  else
    mul!(rtmp, cache.phihA, integrator.fsalfirst)
    @muladd @.. broadcast=false u = uprev + dt*rtmp
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::ETDRK2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  if alg.krylov
    F1 = integrator.fsalfirst
    w1 = phiv(dt, A, F1, 2; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    U2 = uprev + dt * w1[:, 2]
    F2 = _compute_nl(f, U2, p, t + dt, A) + A * uprev
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 1
    else
      integrator.destats.nf += 1
    end
    w2 = phiv(dt, A, F2, 2; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    u = uprev + dt * (w1[:, 2] - w1[:, 3] + w2[:, 3])
  else
    phi1, phi2 = cache.ops
    # The caching version uses a special formula to save computation
    G1 = f.f2(uprev, p, t)
    integrator.destats.nf2 += 1
    F1 = integrator.fsalfirst
    U2 = uprev + dt * (phi1 * F1)
    G2 = f.f2(U2, p, t + dt)
    integrator.destats.nf2 += 1
    u = U2 + dt * (phi2 * (G2 - G1))
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::ETDRK2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,F2,J,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  if alg.krylov
    F1 = integrator.fsalfirst
    Ks, phiv_cache, ws = KsCache
    w1, w2 = ws
    # Krylov for F1
    arnoldi!(Ks, A, F1; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    phiv!(w1, dt, Ks, 2; cache=phiv_cache)
    # Krylov for F2
    @muladd @.. broadcast=false tmp = uprev + dt * @view(w1[:, 2])
    _compute_nl!(F2, f, tmp, p, t + dt, A, rtmp)
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 1
    else
      integrator.destats.nf += 1
    end
    F2 .+= mul!(rtmp, A, uprev)
    arnoldi!(Ks, A, F2; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    phiv!(w2, dt, Ks, 2; cache=phiv_cache)
    # Update u
    u .= uprev
    axpy!( dt, @view(w1[:, 2]), u)
    axpy!(-dt, @view(w1[:, 3]), u)
    axpy!( dt, @view(w2[:, 3]), u)
  else
    phi1, phi2 = cache.ops
    F1 = integrator.fsalfirst
    # The caching version uses a special formula to save computation
    # Compute U2
    mul!(rtmp, phi1, F1)
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U2
    # Compute G2 - G1, storing result in the cache F2
    f.f2(rtmp, uprev, p, t)
    integrator.destats.nf2 += 1
    f.f2(F2, tmp, p, t + dt)
    integrator.destats.nf2 += 1
    F2 .-= rtmp # "F2" is G2 - G1
    # Update u
    u .= tmp
    axpy!(dt, mul!(rtmp, phi2, F2), u)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::ETDRK3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  Au = A * uprev
  F1 = integrator.fsalfirst
  if alg.krylov
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov on F1 (first column)
    Ks = arnoldi(A, F1; kwargs...)
    w1_half = phiv(dt/2, Ks, 1)
    w1 = phiv(dt, Ks, 3)
    U2 = uprev + dt/2 * w1_half[:, 2]
    F2 = _compute_nl(f, U2, p, t + dt/2, A) + Au
    # Krylov on F2 (second column)
    w2 = phiv(dt, A, F2, 3; kwargs...)
    U3 = uprev + dt * (2w2[:, 2] - w1[:, 2])
    F3 = _compute_nl(f, U3, p, t + dt, A) + Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 2
    else
      integrator.destats.nf += 2
    end
    # Krylov on F3 (third column)
    w3 = phiv(dt, A, F3, 3; kwargs...)
    u = uprev + dt * (4w1[:,4] - 3w1[:,3] + w1[:,2]
                      -8w2[:,4] + 4w2[:,3]
                      +4w3[:,4] - w3[:,3])
  else
    A21, A3, B1, B2, B3 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    U2 = uprev + dt * (A21 * F1)
    F2 = f.f2(U2, p, t + dt/2) + Au
    integrator.destats.nf2 += 1
    # stage 3
    U3 = uprev + dt * (A3 * (2F2 - F1))
    F3 = f.f2(U3, p, t + dt) + Au
    integrator.destats.nf2 += 1
    # update u
    u = uprev + dt * (B1 * F1 + B2 * F2 + B3 * F3)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::ETDRK3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,Au,F2,F3,J,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  mul!(Au, A, uprev)
  halfdt = dt/2
  if alg.krylov
    Ks, phiv_cache, ws = KsCache
    w1_half, w1, w2, w3 = ws
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov for F1 (first column)
    arnoldi!(Ks, A, F1; kwargs...)
    phiv!(w1_half, halfdt, Ks, 1; cache=phiv_cache)
    phiv!(w1, dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + halfdt * w1_half[:, 2] # tmp is U2
    _compute_nl!(F2, f, tmp, p, t + halfdt, A, rtmp); F2 .+= Au
    # Krylov for F2 (second column)
    arnoldi!(Ks, A, F2; kwargs...)
    phiv!(w2, dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + dt * (2*w2[:, 2] - w1[:, 2]) # tmp is U3
    _compute_nl!(F3, f, tmp, p, t + dt, A, rtmp); F3 .+= Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 2
    else
      integrator.destats.nf += 2
    end
    # Krylov for F3 (third column)
    arnoldi!(Ks, A, F3; kwargs...)
    phiv!(w3, dt, Ks, 3; cache=phiv_cache)
    # Update u
    @views @.. broadcast=false rtmp = 4w1[:,4] - 3w1[:,3] + w1[:,2] - 8w2[:,4] + 4w2[:,3] + 4w3[:,4] - w3[:,3]
    @muladd @.. broadcast=false u = uprev + dt * rtmp
  else
    A21, A3, B1, B2, B3 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    mul!(rtmp, A21, F1)
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U2
    f.f2(F2, tmp, p, t + halfdt); F2 .+= Au
    integrator.destats.nf2 += 1
    # stage 3
    @muladd @.. broadcast=false F3 = 2 * F2 - F1 # use F3 temporarily as cache
    mul!(rtmp, A3, F3)
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U3
    f.f2(F3, tmp, p, t + dt); F3 .+= Au
    integrator.destats.nf2 += 1
    # update u
    u .= uprev
    axpy!(dt, mul!(rtmp, B1, F1), u)
    axpy!(dt, mul!(rtmp, B2, F2), u)
    axpy!(dt, mul!(rtmp, B3, F3), u)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::ETDRK4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  Au = A * uprev
  F1 = integrator.fsalfirst
  halfdt = dt/2
  if alg.krylov
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov on F1 (first column)
    Ks = arnoldi(A, F1; kwargs...)
    w1_half = phiv(halfdt, Ks, 1)
    w1 = phiv(dt, Ks, 3)
    U2 = uprev + halfdt * w1_half[:, 2]
    F2 = _compute_nl(f, U2, p, t + halfdt, A) + Au
    # Krylov on F2 (second column)
    Ks = arnoldi(A, F2; kwargs...)
    w2_half = phiv(halfdt, Ks, 1)
    w2 = phiv(dt, Ks, 3)
    U3 = uprev + halfdt * w2_half[:, 2]
    F3 = _compute_nl(f, U3, p, t + halfdt, A) + Au
    # Krylov on F3 (third column)
    w3 = phiv(dt, A, F3, 3; kwargs...)
    # Extra Krylov for computing F4
    rtmp = 2F3 - F1 - Au + A*U2
    wtmp = phiv(halfdt, A, rtmp, 1; kwargs...)
    U4 = U2 + halfdt * wtmp[:, 2]
    F4 = _compute_nl(f, U4, p, t + dt, A) + Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 3
    else
      integrator.destats.nf += 3
    end
    # Krylov on F4 (fourth column)
    w4 = phiv(dt, A, F4, 3; kwargs...)
    # update u
    u = uprev + dt * (w1[:,2] - 3w1[:,3] + 4w1[:,4] + 2w2[:,3] - 4w2[:,4] +
                      2w3[:,3] - 4w3[:,4] + 4w4[:,4] - w4[:,3])
  else
    A21, A41, A43, B1, B2, B4 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    U2 = uprev + dt * (A21 * F1)
    F2 = f.f2(U2, p, t + halfdt) + Au
    # stage 3
    U3 = uprev + dt * (A21 * F2) # A32 = A21
    F3 = f.f2(U3, p, t + halfdt) + Au
    # stage 4
    U4 = uprev + dt * (A41 * F1 + A43 * F3)
    F4 = f.f2(U4, p, t + dt) + Au
    integrator.destats.nf2 += 3
    # update u
    u = uprev + dt * (B1 * F1 + B2 * (F2 + F3) + B4 * F4) # B3 = B2
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::ETDRK4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,Au,F2,F3,F4,J,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  mul!(Au, A, uprev)
  halfdt = dt/2
  if alg.krylov
    Ks, phiv_cache, ws = KsCache
    w1_half, w2_half, w1, w2, w3, w4 = ws
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov for F1 (first column)
    arnoldi!(Ks, A, F1; kwargs...)
    phiv!(w1_half, halfdt, Ks, 1; cache=phiv_cache)
    phiv!(w1, dt, Ks, 3; cache=phiv_cache)
    U2 = u # temporarily use u to store U2 (used in the extra Krylov step)
    @muladd @.. broadcast=false @views U2 = uprev + halfdt * w1_half[:, 2]
    _compute_nl!(F2, f, U2, p, t + halfdt, A, rtmp); F2 .+= Au
    # Krylov for F2 (second column)
    arnoldi!(Ks, A, F2; kwargs...)
    phiv!(w2_half, halfdt, Ks, 1; cache=phiv_cache)
    phiv!(w2, dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + halfdt * w2_half[:, 2] # tmp is U3
    _compute_nl!(F3, f, tmp, p, t + halfdt, A, rtmp); F3 .+= Au
    # Krylov for F3 (third column)
    arnoldi!(Ks, A, F3; kwargs...)
    phiv!(w3, dt, Ks, 3; cache=phiv_cache)
    # Extra Krylov for computing F4
    # Compute rtmp = 2F3 - F1 - Au + A*U2
    mul!(rtmp, A, U2); @.. broadcast=false rtmp += 2F3 - F1 - Au
    arnoldi!(Ks, A, rtmp; kwargs...)
    phiv!(w1_half, halfdt, Ks, 1; cache=phiv_cache) # original w1_half is no longer needed
    @muladd @.. broadcast=false @views tmp = U2 + halfdt * w1_half[:, 2] # tmp is U4
    _compute_nl!(F4, f, tmp, p, t + dt, A, rtmp); F4 .+= Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 3
    else
      integrator.destats.nf += 3
    end
    # Krylov for F4 (fourth column)
    arnoldi!(Ks, A, F4; kwargs...)
    phiv!(w4, dt, Ks, 3; cache=phiv_cache)
    # update u
    @views @.. broadcast=false rtmp = w1[:,2] - 3w1[:,3] + 4w1[:,4] + 2w2[:,3] - 4w2[:,4] +
                     2w3[:,3] - 4w3[:,4] + 4w4[:,4] - w4[:,3]
    @muladd @.. broadcast=false u = uprev + dt * rtmp
  else
    A21, A41, A43, B1, B2, B4 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    mul!(rtmp, A21, F1)
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U2
    f.f2(F2, tmp, p, t + halfdt); F2 .+= Au
    # stage 3
    mul!(rtmp, A21, F2) # A32 = A21
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U3
    f.f2(F3, tmp, p, t + halfdt); F3 .+= Au
    # stage 4
    @.. broadcast=false tmp = uprev
    axpy!(dt, mul!(rtmp, A41, F1), tmp)
    axpy!(dt, mul!(rtmp, A43, F3), tmp) # tmp is U4
    f.f2(F4, tmp, p, t + dt); F4 .+= Au
    integrator.destats.nf2 += 3
    # update u
    u .= uprev
    axpy!(dt, mul!(rtmp, B1, F1), u)
    F2 .+= F3; axpy!(dt, mul!(rtmp, B2, F2), u) # B3 = B2
    axpy!(dt, mul!(rtmp, B4, F4), u)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::HochOst4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = isa(f, SplitFunction) ? f.f1.f : calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  Au = A * uprev
  F1 = integrator.fsalfirst
  halfdt = dt/2
  if alg.krylov
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov on F1 (first column)
    Ks = arnoldi(A, F1; kwargs...)
    w1_half = phiv(halfdt, Ks, 3)
    w1 =      phiv(dt,     Ks, 3)
    U2 = uprev + halfdt * w1_half[:, 2]
    F2 = _compute_nl(f, U2, p, t + halfdt, A) + Au
    # Krylov on F2 (second column)
    Ks = arnoldi(A, F2; kwargs...)
    w2_half = phiv(halfdt, Ks, 3)
    w2 =      phiv(dt,     Ks, 3)
    U3 = uprev + dt * (0.5w1_half[:,2] - w1_half[:,3] + w2_half[:,3])
    F3 = _compute_nl(f, U3, p, t + halfdt, A) + Au
    # Krylov on F3 (third column)
    Ks = arnoldi(A, F3; kwargs...)
    w3_half = phiv(halfdt, Ks, 3)
    w3 =      phiv(dt,     Ks, 3)
    U4 = uprev + dt * (w1[:,2] - 2w1[:,3] + w2[:,3] + w3[:,3])
    F4 = _compute_nl(f, U4, p, t + dt, A) + Au
    # Krylov on F4 (fourth column)
    Ks = arnoldi(A, F4; kwargs...)
    w4_half = phiv(halfdt, Ks, 3)
    w4 =      phiv(dt,     Ks, 3)
    U5 = uprev + dt * (0.5w1_half[:,2] - 0.75w1_half[:,3] + 0.5w1_half[:,4] + w1[:,4] - 0.25w1[:,3] +
                       0.5w2_half[:,3] - w2[:,4] + 0.25w2[:,3] - 0.5w2_half[:,4] +
                       0.5w3_half[:,3] - w2[:,4] + 0.25w3[:,3] - 0.5w3_half[:,4] +
                       w4[:,4] - 0.25w4[:,3] - 0.25w4_half[:,3] + 0.5w4_half[:,4])
    F5 = _compute_nl(f, U5, p, t + halfdt, A) + Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 4
    else
      integrator.destats.nf += 4
    end
    # Krylov on F5 (fifth column)
    w5 = phiv(dt, A, F5, 3; kwargs...)
    # update u
    u = uprev + dt * (w1[:,2] - 3w1[:,3] + 4w1[:,4] - w4[:,3] + 4w4[:,4] + 4w5[:,3] - 8w5[:,4])
  else
    A21, A31, A32, A41, A42, A51, A52, A54, B1, B4, B5 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    U2 = uprev + dt * (A21 * F1)
    F2 = f.f2(U2, p, t + halfdt) + Au
    # stage 3
    U3 = uprev + dt * (A31 * F1 + A32 * F2)
    F3 = f.f2(U3, p, t + halfdt) + Au
    # stage 4
    U4 = uprev + dt * (A41 * F1 + A42 * (F2 + F3))
    F4 = f.f2(U4, p, t + dt) + Au
    # stage 5
    U5 = uprev + dt * (A51 * F1 + A52 * (F2 + F3) + A54 * F4)
    F5 = f.f2(U5, p, t + halfdt) + Au
    integrator.destats.nf2 += 4
    # update u
    u = uprev + dt * (B1 * F1 + B4 * F4 + B5 * F5)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::HochOst4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,Au,F2,F3,F4,F5,J,KsCache = cache
  A = isa(f, SplitFunction) ? f.f1.f : (calc_J!(J, integrator, cache); J) # get linear operator
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  mul!(Au, A, uprev)
  halfdt = dt/2
  if alg.krylov
    Ks, phiv_cache, ws = KsCache
    w1_half, w2_half, w3_half, w4_half, w1, w2, w3, w4, w5 = ws
    kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
    # Krylov on F1 (first column)
    arnoldi!(Ks, A, F1; kwargs...)
    phiv!(w1_half, halfdt, Ks, 3; cache=phiv_cache)
    phiv!(w1,          dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + halfdt * w1_half[:, 2] # tmp is U2
    _compute_nl!(F2, f, tmp, p, t + halfdt, A, rtmp); F2 .+= Au
    # Krylov on F2 (second column)
    arnoldi!(Ks, A, F2; kwargs...)
    phiv!(w2_half, halfdt, Ks, 3; cache=phiv_cache)
    phiv!(w2,          dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + dt * (0.5w1_half[:,2] - w1_half[:,3] + w2_half[:,3]) # tmp is U3
    _compute_nl!(F3, f, tmp, p, t + halfdt, A, rtmp); F3 .+= Au
    # Krylov on F3 (third column)
    arnoldi!(Ks, A, F3; kwargs...)
    phiv!(w3_half, halfdt, Ks, 3; cache=phiv_cache)
    phiv!(w3,          dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + dt * (w1[:,2] - 2w1[:,3] + w2[:,3] + w3[:,3]) # tmp is U4
    _compute_nl!(F4, f, tmp, p, t + dt, A, rtmp); F4 .+= Au
    # Krylov on F4 (fourth column)
    arnoldi!(Ks, A, F4; kwargs...)
    phiv!(w4_half, halfdt, Ks, 3; cache=phiv_cache)
    phiv!(w4,          dt, Ks, 3; cache=phiv_cache)
    @muladd @.. broadcast=false @views tmp = uprev + dt * (
      0.5w1_half[:,2] - 0.75w1_half[:,3] + 0.5w1_half[:,4] + w1[:,4] - 0.25w1[:,3] +
      0.5w2_half[:,3] - w2[:,4] + 0.25w2[:,3] - 0.5w2_half[:,4] +
      0.5w3_half[:,3] - w2[:,4] + 0.25w3[:,3] - 0.5w3_half[:,4] +
      w4[:,4] - 0.25w4[:,3] - 0.25w4_half[:,3] + 0.5w4_half[:,4]) # tmp is U5
    _compute_nl!(F5, f, tmp, p, t + dt, A, rtmp); F5 .+= Au
    if isa(f, SplitFunction)
      integrator.destats.nf2 += 4
    else
      integrator.destats.nf += 4
    end
    # Krylov on F5 (fifth column)
    arnoldi!(Ks, A, F5; kwargs...)
    phiv!(w5, dt, Ks, 3; cache=phiv_cache)
    # update u
    @muladd @.. broadcast=false @views rtmp = w1[:,2] - 3w1[:,3] + 4w1[:,4] - w4[:,3] + 4w4[:,4] + 4w5[:,3] - 8w5[:,4]
    @muladd @.. broadcast=false u = uprev + dt * rtmp
  else
    A21, A31, A32, A41, A42, A51, A52, A54, B1, B4, B5 = cache.ops
    # stage 1 (fsaled)
    # stage 2
    mul!(rtmp, A21, F1)
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U2
    f.f2(F2, tmp, p, t + halfdt); F2 .+= Au
    # stage 3
    mul!(rtmp, A31, F1); mul!(rtmp2, A32, F2); rtmp .+= rtmp2
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U3
    f.f2(F3, tmp, p, t + halfdt); F3 .+= Au
    # stage 4
    F2 .+= F3 # F2 now stores F2 + F3
    mul!(rtmp, A41, F1); mul!(rtmp2, A42, F2); rtmp .+= rtmp2
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U4
    f.f2(F4, tmp, p, t + dt); F4 .+= Au
    # stage 5
    mul!(rtmp, A51, F1)
    mul!(rtmp2, A52, F2); rtmp .+= rtmp2
    mul!(rtmp2, A54, F4); rtmp .+= rtmp2
    @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is U5
    f.f2(F5, tmp, p, t + halfdt); F5 .+= Au
    integrator.destats.nf2 += 4
    # update u
    mul!(rtmp, B1, F1)
    mul!(rtmp2, B4, F4); rtmp .+= rtmp2
    mul!(rtmp2, B5, F5); rtmp .+= rtmp2
    @muladd @.. broadcast=false u = uprev + dt * rtmp
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

#############################################
# EPIRK integrators
function perform_step!(integrator, cache::Exp4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  ts = [dt/3, 2dt/3, dt]
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Krylov for f(uprev)
  B1 = [zero(f0) f0]
  K1 = phiv_timestep(ts, J, B1; kwargs...) # tϕ(tA)f0
  @inbounds for i = 1:3
    K1[:,i] ./= ts[i]
  end
  w4 = K1 * [-7/300, 97/150, -37/300]
  u4 = uprev + dt * w4
  d4 = f(u4, p, t+dt) - f0 - dt * (J*w4) # TODO: what should be the time?
  integrator.destats.nf += 1
  # Krylov for the first remainder d4
  B2 = [zero(d4) d4]
  K2 = phiv_timestep(ts, J, B2; kwargs...)
  @inbounds for i = 1:3
    K2[:,i] ./= ts[i]
  end
  w7 = K1 * [59/300, -7/75, 269/300] + K2 * [2/3, 2/3, 2/3]
  u7 = uprev + dt * w7
  d7 = f(u7, p, t+dt) - f0 - dt * (J*w7)
  integrator.destats.nf += 1
  # Krylov for the second remainder d7
  B3 = [zero(d7) d7]
  k7 = phiv_timestep(ts[1], J, B3; kwargs...)
  k7 ./= ts[1]
  # Update u
  u = uprev + dt * (K1[:,3] + K2[:,1] - 4/3*K2[:,2] + K2[:,3] + k7/6)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::Exp4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  ts = [dt/3, 2dt/3, dt]
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Krylov for f(uprev)
  B[:, 2] .= f0
  phiv_timestep!(K, ts, J, B; kwargs...)
  @inbounds for i = 1:3
    K[:,i] ./= ts[i]
  end
  mul!(rtmp, K, [-7/300, 97/150, -37/300]) # rtmp is now w4
  @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is now u4
  mul!(rtmp2, J, rtmp)
  f(rtmp, tmp, p, t+dt) # TODO: what should be the time?
  integrator.destats.nf += 1
  @muladd @.. broadcast=false @view(B[:,2]) = rtmp - f0 - dt * rtmp2 # B[:,2] is now d4
  # Partially update entities that use k1, k2, k3
  mul!(rtmp, K, [59/300, -7/75, 269/300]) # rtmp is now w7
  @muladd @.. broadcast=false u = uprev + dt * @view(K[:,3])
  # Krylov for the first remainder d4
  phiv_timestep!(K, ts, J, B; kwargs...)
  @inbounds for i = 1:3
    K[:,i] ./= ts[i]
  end
  mul!(rtmp2, K, [2/3, 2/3, 2/3]); rtmp .+= rtmp2 # w7 fully updated
  @muladd @.. broadcast=false tmp = uprev + dt * rtmp # tmp is now u7
  mul!(rtmp2, J, rtmp)
  f(rtmp, tmp, p, t+dt) # TODO: what should be the time?
  integrator.destats.nf += 1
  @muladd @.. broadcast=false @view(B[:,2]) = rtmp - f0 - dt * rtmp2 # B[:,2] is now d7
  # Partially update entities that use k4, k5, k6
  mul!(rtmp, K, [1.0, -4/3, 1.0])
  axpy!(dt, rtmp, u)
  # Krylov for the second remainder d7
  k7 = @view(K[:, 1])
  phiv_timestep!(k7, ts[1], J, B; kwargs...)
  k7 ./= ts[1]
  axpy!(dt/6, k7, u)

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EPIRK4s3AConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Compute U2 and U3 vertically
  K = phiv_timestep([dt/2, 2dt/3], J, [zero(f0) f0]; kwargs...)
  U2 = uprev + K[:, 1]
  U3 = uprev + K[:, 2]
  R2 = f(U2, p, t + dt/2)  - f0 - J*K[:, 1] # remainder of U2
  R3 = f(U3, p, t + 2dt/3) - f0 - J*K[:, 2] # remainder of U3
  integrator.destats.nf += 2

  # Update u (horizontally)
  B = zeros(eltype(f0), length(f0), 5)
  B[:, 2] = f0
  B[:, 4] = (32R2 - 13.5R3) / dt^2
  B[:, 5] = (-144R2 + 81R3) / dt^3
  u = uprev + phiv_timestep(dt, J, B; kwargs...)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EPIRK4s3ACache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Compute U2 and U3 vertically
  B[:, 2] .= f0
  phiv_timestep!(K, [dt/2, 2dt/3], J, @view(B[:, 1:2]); kwargs...)
  ## U2 and R2
  @.. broadcast=false tmp = uprev + @view(K[:, 1]) # tmp is now U2
  f(rtmp, tmp, p, t + dt/2); mul!(rtmp2, J, @view(K[:, 1]))
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
  B[:, 4] .= (32/dt^2) * rtmp
  B[:, 5] .= (-144/dt^3) * rtmp
  ## U3 and R3
  @.. broadcast=false tmp = uprev + @view(K[:, 2]) # tmp is now U3
  f(rtmp, tmp, p, t + 2dt/3); mul!(rtmp2, J, @view(K[:, 2]))
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
  B[:, 4] .-= (13.5/dt^2) * rtmp
  B[:, 5] .+= (81/dt^3) * rtmp
  integrator.destats.nf += 2

  # Update u
  du = @view(K[:, 1])
  phiv_timestep!(du, dt, J, B; kwargs...)
  @.. broadcast=false u = uprev + du

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EPIRK4s3BConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Compute U2 and U3 vertically
  K = phiv_timestep([dt/2, 3dt/4], J, [zero(f0) zero(f0) f0]; kwargs...)
  K[:, 1] .*= 8 / (3*dt)
  K[:, 2] .*= 16 / (9*dt)
  U2 = uprev + K[:, 1]
  U3 = uprev + K[:, 2]
  R2 = f(U2, p, t + dt/2)  - f0 - J*K[:, 1] # remainder of U2
  R3 = f(U3, p, t + 3dt/4) - f0 - J*K[:, 2] # remainder of U3
  integrator.destats.nf += 2

  # Update u (horizontally)
  B = zeros(eltype(f0), length(f0), 5)
  B[:, 2] = f0
  B[:, 4] = (54R2 - 16R3) / dt^2
  B[:, 5] = (-324R2 + 144R3) / dt^3
  u = uprev + phiv_timestep(dt, J, B; kwargs...)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EPIRK4s3BCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Compute U2 and U3 vertically
  fill!(@view(B[:, 2]), zero(eltype(B)))
  B[:, 3] .= f0
  phiv_timestep!(K, [dt/2, 3dt/4], J, @view(B[:, 1:3]); kwargs...)
  K[:, 1] .*= 8 / (3*dt)
  K[:, 2] .*= 16 / (9*dt)
  ## U2 and R2
  @.. broadcast=false tmp = uprev + @view(K[:, 1]) # tmp is now U2
  f(rtmp, tmp, p, t + dt/2); mul!(rtmp2, J, @view(K[:, 1]))
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
  B[:, 4] .= (54/dt^2) * rtmp
  B[:, 5] .= (-324/dt^3) * rtmp
  ## U3 and R3
  @.. broadcast=false tmp = uprev + @view(K[:, 2]) # tmp is now U3
  f(rtmp, tmp, p, t + 3dt/4); mul!(rtmp2, J, @view(K[:, 2]))
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
  B[:, 4] .-= (16/dt^2) * rtmp
  B[:, 5] .+= (144/dt^3) * rtmp
  integrator.destats.nf += 2

  # Update u
  fill!(@view(B[:, 3]), zero(eltype(B)))
  B[:, 2] .= f0
  du = @view(K[:, 1])
  phiv_timestep!(du, dt, J, B; kwargs...)
  @.. broadcast=false u = uprev + du

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EPIRK5s3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Compute U2 horizontally
  B = fill(zero(eltype(f0)), length(f0), 4)
  B[:, 3] = (55 / (8 * dt)) * f0
  B[:, 4] = (-3025 / (192 * dt^2)) * f0
  k = phiv_timestep(48dt/55, J, B; kwargs...)
  U2 = uprev + k
  R2 = f(U2, p, t + 48dt/55)  - f0 - J*k # remainder of U2
  integrator.destats.nf += 1

  # Compute U3 horizontally
  B[:, 2] = (53/5) * f0
  B[:, 3] = (-648 / (5 * dt)) * f0
  B[:, 4] = (2916 / (5 * dt^2)) * f0 + (32065 / (1152 * dt^2)) * R2
  k = phiv_timestep(4dt/9, J, B; kwargs...)
  U3 = uprev + k
  R3 = f(U3, p, t + 4dt/9)  - f0 - J*k # remainder of U3
  integrator.destats.nf += 1

  # Update u (horizontally)
  B = fill(zero(eltype(f0)), length(f0), 5)
  B[:, 2] = f0
  B[:, 4] = (-166375 / (61056 * dt^2)) * R2 + (2187 / (106 * dt^2)) * R3
  B[:, 5] = (499125 / (27136 * dt^3)) * R2 - (2187 / (106 * dt^3)) * R3
  u = uprev + phiv_timestep(dt, J, B; kwargs...)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EPIRK5s3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,k,rtmp,rtmp2,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Compute U2 horizontally
  fill!(@view(B[:, 2]), zero(eltype(B)))
  B[:, 3] .= (55 / (8 * dt)) .* f0
  B[:, 4] .= (-3025 / (192 * dt^2)) .* f0
  phiv_timestep!(k, 48dt/55, J, @view(B[:, 1:4]); kwargs...)
  ## Compute R2
  @.. broadcast=false tmp = uprev + k # tmp is now U2
  f(rtmp, tmp, p, t + 48dt/55); mul!(rtmp2, J, k)
  integrator.destats.nf += 1
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2

  # Compute U3 horizontally
  B[:, 2] .= (53/5) .* f0
  B[:, 3] .= (-648 / (5 * dt)) .* f0
  B[:, 4] .= (2916 / (5 * dt^2)) .* f0 + (32065 / (1152 * dt^2)) .* rtmp
  phiv_timestep!(k, 4dt/9, J, @view(B[:, 1:4]); kwargs...)
  ## Update B matrix using R2
  B[:, 2] .= f0
  fill!(@view(B[:, 3]), zero(eltype(B)))
  B[:, 4] .= (-166375 / (61056 * dt^2)) .* rtmp
  B[:, 5] .= (499125 / (27136 * dt^3)) .* rtmp
  ## Compute R3 and update B
  @.. broadcast=false tmp = uprev + k # tmp is now U3
  f(rtmp, tmp, p, t + 4dt/9); mul!(rtmp2, J, k)
  integrator.destats.nf += 1
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
  B[:, 4] .+= (2187 / (106 * dt^2)) .* rtmp
  B[:, 5] .-= (2187 / (106 * dt^3)) .* rtmp

  # Update u
  phiv_timestep!(k, dt, J, B; kwargs...)
  @.. broadcast=false u = uprev + k

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EXPRB53s3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Compute the first group for U2 and U3
  B = [zero(f0) f0]
  K = phiv_timestep([dt/2, 9dt/10], J, B; kwargs...)
  U2 = uprev + K[:, 1]
  R2 = f(U2, p, t + dt/2)  - f0 - J*K[:, 1] # remainder of U2
  integrator.destats.nf += 1
  U3 = uprev + K[:, 2] # partially

  # Compute the second group for U3
  B = [zero(R2) zero(R2) zero(R2) R2]
  K = phiv_timestep([dt/2, 9dt/10], J, B; kwargs...)
  U3 .+= 216/(25*dt^2) .* K[:, 1] + 8/dt^2 .* K[:, 2]
  R3 = f(U3, p, t + 9dt/10)  - f0 - J*(U3 - uprev) # remainder of U3
  integrator.destats.nf += 1

  # Compute the third group for u
  B = fill(zero(eltype(f0)), length(f0), 5)
  B[:, 2] = f0
  B[:, 4] = (18 / dt^2) * R2 - (250 / (81 * dt^2)) * R3
  B[:, 5] = (-60 / dt^3) * R2 + (500 / (27 * dt^3)) * R3
  u = uprev + phiv_timestep(dt, J, B; kwargs...)

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EXPRB53s3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Compute the first group for U2 and U3
  B[:, 2] .= f0
  phiv_timestep!(K, [dt/2, 9dt/10], J, @view(B[:, 1:2]); kwargs...)
  ## U2 and R2
  @.. broadcast=false tmp = uprev + @view(K[:, 1]) # tmp is now U2
  f(rtmp, tmp, p, t + dt/2); mul!(rtmp2, J, @view(K[:, 1]))
  integrator.destats.nf += 1
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
  @.. broadcast=false tmp = uprev + @view(K[:, 2]) # tmp is now U3 (partially)

  # Compute the second group for U3
  fill!(@view(B[:, 2]), zero(eltype(B)))
  B[:, 4] .= rtmp
  phiv_timestep!(K, [dt/2, 9dt/10], J, @view(B[:, 1:4]); kwargs...)
  ## Update B using R2
  B[:, 2] .= f0
  B[:, 4] .= (18 / dt^2) .* rtmp
  B[:, 5] .= (-60 / dt^3) .* rtmp
  ## U3 and R3
  @views tmp .+= 216/(25*dt^2) .* K[:, 1] + 8/dt^2 .* K[:, 2] # tmp is now U3
  f(rtmp, tmp, p, t + 9dt/10)
  integrator.destats.nf += 1
  tmp .-= uprev; mul!(rtmp2, J, tmp)
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R3
  ## Update B using R3
  B[:, 4] .-= (250 / (81 * dt^2)) * rtmp
  B[:, 5] .+= (500 / (27 * dt^3)) * rtmp

  # Update u
  du = @view(K[:, 1])
  phiv_timestep!(du, dt, J, B; kwargs...)
  @.. broadcast=false u = uprev + du

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EPIRK5P1ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Coefficients (scaling factors absorbed)
  g11 = 0.35129592695058193092 * dt
  g21 = 0.84405472011657126298 * dt # g22 = g32
  g31 = dt; g32 = 0.71111095364366870359 * dt; g33 = 0.62378111953371494809 * dt
  a22 = 2.37739153404186675585; b2 = 1.78975267532362337976; b3 = 9.35854650579261718128 / dt^2

  # Compute the first column (f0)
  B = [zero(f0) f0]
  K1 = phiv_timestep([g11, g21, g31], J, B; kwargs...)
  ## U1 and R1
  U1 = uprev + K1[:, 1]
  R1 = f(U1, p, t + g11) - f0 - J*K1[:, 1] # remainder of U1
  integrator.destats.nf += 1

  # Compute the second column (R1)
  B = [zero(R1) R1]
  k2 = phiv_timestep(g32, J, B; kwargs...)
  ## U2 and R2
  U2 = uprev + K1[:, 2] + a22 * k2
  R2 = f(U2, p, t + g21) - f0 - J*(U2 - uprev) # remainder of U2
  integrator.destats.nf += 1

  # Compute the third column (dR = R2 - 2R1)
  B = fill(zero(eltype(f0)), length(f0), 4)
  B[:, 4] = R2 - 2R1
  k3 = phiv_timestep(g33, J, B; kwargs...)
  ## Update state
  u = uprev + K1[:, 3] + b2 * k2 + b3 * k3

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EPIRK5P1Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Coefficients (scaling factors absorbed)
  g11 = 0.35129592695058193092 * dt
  g21 = 0.84405472011657126298 * dt # g22 = g32
  g31 = dt; g32 = 0.71111095364366870359 * dt; g33 = 0.62378111953371494809 * dt
  a22 = 2.37739153404186675585; b2 = 1.78975267532362337976; b3 = 9.35854650579261718128 / dt^2

  # Compute the first column (f0)
  B[:, 2] .= f0
  phiv_timestep!(K, [g11, g21, g31], J, @view(B[:, 1:2]); kwargs...)
  ## U1 and R1
  @.. broadcast=false tmp = uprev + @view(K[:, 1]) # tmp is now U1
  f(rtmp, tmp, p, t + g11); mul!(rtmp2, J, @view(K[:, 1]))
  integrator.destats.nf += 1
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R1
  @.. broadcast=false tmp = uprev + @view(K[:, 2]) # partially update U2 (stored tmp)
  @.. broadcast=false u = uprev + @view(K[:, 3]) # partially update u
  B[:, 2] .= rtmp
  @.. broadcast=false @view(B[:, 4]) = (-2) * rtmp

  # Compute the second column (R1)
  k = @view(K[:, 1])
  phiv_timestep!(k, g32, J, @view(B[:, 1:2]); kwargs...)
  ## U2 and R2
  axpy!(a22, k, tmp) # tmp is now U2
  f(rtmp, tmp, p, t + g21)
  integrator.destats.nf += 1
  tmp .-= uprev; mul!(rtmp2, J, tmp)
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
  axpy!(b2, k, u) # partially update u
  B[:, 4] .+= rtmp # is now dR

  # Compute the third column (dR = R2 - 2R1)
  fill!(@view(B[:, 2]), zero(eltype(B)))
  phiv_timestep!(k, g33, J, B; kwargs...)
  axpy!(b3, k, u)

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::EPIRK5P2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  J = calc_J(integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(uprev) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)))

  # Coefficients (scaling factors absorbed)
  g11 = 0.46629408528088195806 * dt
  g21 = 0.88217912653363865140 * dt # g22 = g32
  g31 = dt; g32 = 0.92074916488140031449 * dt; g33 = 0.79791561832664517267 * dt
  a22 = 2.80620373289331259751 / dt; b2 = 2.52806310256246280783 / dt
  b31 = -0.128486782657005566142; b32 = -0.161028033172809183351 / dt; b33 = 5.26726331616909606251 / dt^2

  # Compute the first column (f0)
  B = [zero(f0) f0]
  K1 = phiv_timestep([g11, g21, g31], J, B; kwargs...)
  ## U1 and R1
  U1 = uprev + K1[:, 1]
  R1 = f(U1, p, t + g11) - f0 - J*K1[:, 1] # remainder of U1
  integrator.destats.nf += 1

  # Compute the second column (R1)
  B = [zero(R1) zero(R1) R1]
  k2 = phiv_timestep(g32, J, B; kwargs...)
  ## U2 and R2
  U2 = uprev + K1[:, 2] + a22 * k2
  R2 = f(U2, p, t + g21) - f0 - J*(U2 - uprev) # remainder of U2
  integrator.destats.nf += 1

  # Compute the third column (dR = R2 - 2R1)
  dR = R2 - 2R1
  B = [zero(dR) b31*dR b32*dR b33*dR]
  k3 = phiv_timestep(g33, J, B; kwargs...)
  ## Update state
  u = uprev + K1[:, 3] + b2 * k2 + k3

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::EPIRK5P2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,rtmp,rtmp2,dR,K,J,B,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)
  f0 = integrator.fsalfirst # f(u0) is fsaled
  kwargs = (tol=integrator.opts.reltol, iop=alg.iop, opnorm=integrator.opts.internalopnorm,
  adaptive=alg.adaptive_krylov, tau=(alg.adaptive_krylov ? 0.0 : dt), m=min(alg.m, size(J,1)), caches=KsCache)

  # Coefficients (scaling factors absorbed)
  g11 = 0.46629408528088195806 * dt
  g21 = 0.88217912653363865140 * dt # g22 = g32
  g31 = dt; g32 = 0.92074916488140031449 * dt; g33 = 0.79791561832664517267 * dt
  a22 = 2.80620373289331259751 / dt; b2 = 2.52806310256246280783 / dt
  b31 = -0.128486782657005566142; b32 = -0.161028033172809183351 / dt; b33 = 5.26726331616909606251 / dt^2

  # Compute the first column (f0)
  B[:, 2] .= f0
  phiv_timestep!(K, [g11, g21, g31], J, @view(B[:, 1:2]); kwargs...)
  ## U1 and R1
  @.. broadcast=false tmp = uprev + @view(K[:, 1]) # tmp is now U1
  f(rtmp, tmp, p, t + g11); mul!(rtmp2, J, @view(K[:, 1]))
  integrator.destats.nf += 1
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R1
  @.. broadcast=false tmp = uprev + @view(K[:, 2]) # partially update U2 (stored in tmp)
  @.. broadcast=false u = uprev + @view(K[:, 3]) # partially update u
  @.. broadcast=false dR = -2rtmp # partially update dR

  # Compute the second column (R1)
  fill!(@view(B[:, 2]), zero(eltype(B)))
  B[:, 3] .= rtmp
  k = @view(K[:, 1])
  phiv_timestep!(k, g32, J, @view(B[:, 1:3]); kwargs...)
  ## U2 and R2
  axpy!(a22, k, tmp) # tmp is now U2
  f(rtmp, tmp, p, t + g21)
  integrator.destats.nf += 1
  tmp .-= uprev; mul!(rtmp2, J, tmp)
  @.. broadcast=false rtmp = rtmp - f0 - rtmp2 # rtmp is now R2
  dR .+= rtmp # dR is now R2 - 2R1
  axpy!(b2, k, u) # partially update u

  # Compute the third column (dR = R2 - 2R1)
  @.. broadcast=false @view(B[:, 2]) = b31 * dR
  @.. broadcast=false @view(B[:, 3]) = b32 * dR
  @.. broadcast=false @view(B[:, 4]) = b33 * dR
  phiv_timestep!(k, g33, J, B; kwargs...)
  u .+= k

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

######################################################
# Adaptive exponential Rosenbrock integrators
function perform_step!(integrator, cache::Exprb32ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  w1 = phiv(dt, A, F1, 3; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  U2 = uprev + dt * w1[:, 2]
  F2 = _compute_nl(f, U2, p, t + dt, A) + A * uprev
  integrator.destats.nf += 1
  w2 = phiv(dt, A, F2, 3; m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  u = uprev + dt * (w1[:,2] - 2w1[:,4] + 2w2[:,4])
  if integrator.opts.adaptive
    # error estimator for the imbedded method
    utilde = 2dt * (-w1[:,4] + w2[:,4])
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::Exprb32Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack utilde,tmp,rtmp,F2,J,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  Ks, phiv_cache, ws = KsCache
  w1, w2 = ws
  # Krylov for F1
  arnoldi!(Ks, J, F1; m=min(alg.m, size(J,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  phiv!(w1, dt, Ks, 3; cache=phiv_cache)
  # Krylov for F2
  @muladd @.. broadcast=false tmp = uprev + dt * @view(w1[:, 2])
  _compute_nl!(F2, f, tmp, p, t + dt, J, rtmp)
  integrator.destats.nf += 1
  F2 .+= mul!(rtmp, J, uprev)
  arnoldi!(Ks, J, F2; m=min(alg.m, size(J,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  phiv!(w2, dt, Ks, 3; cache=phiv_cache)
  # Update u
  u .= uprev
  axpy!(dt, @view(w1[:,2]), u)
  axpy!(-2dt, @view(w1[:,4]), u)
  axpy!(2dt, @view(w2[:,4]), u)
  if integrator.opts.adaptive
    # error estimator for the imbedded method
    @views @.. broadcast=false utilde = (2*dt) * (-w1[:,4] + w2[:,4])
    calculate_residuals!(tmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(tmp,t)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

function perform_step!(integrator, cache::Exprb43ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  A = calc_J(integrator, cache) # get linear operator
  alg = unwrap_alg(integrator, true)

  Au = A * uprev
  F1 = integrator.fsalfirst
  kwargs = (m=min(alg.m, size(A,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  # Krylov on F1 (first column)
  Ks = arnoldi(A, F1; kwargs...)
  w1_half = phiv(dt/2, Ks, 1)
  w1 = phiv(dt, Ks, 4)
  U2 = uprev + dt/2 * w1_half[:, 2]
  F2 = _compute_nl(f, U2, p, t + dt/2, A) + Au
  integrator.destats.nf += 1
  # Krylov on F2 (second column)
  w2 = phiv(dt, A, F2, 4; kwargs...)
  U3 = uprev + dt * w2[:, 2]
  F3 = _compute_nl(f, U3, p, t + dt, A) + Au
  integrator.destats.nf += 1
  # Krylov on F3 (third column)
  w3 = phiv(dt, A, F3, 4; kwargs...)
  u = uprev + dt * (w1[:,2] - 14w1[:,4] + 36w1[:,5] + 16w2[:,4] - 48w2[:,5] - 2w3[:,4] + 12w3[:,5])
  if integrator.opts.adaptive
    # error estimator for the imbedded method
    utilde = dt * (36w1[:,5] - 48w2[:,5] + 12w3[:,5])
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  # Update integrator state
  integrator.fsallast = f(u, p, t + dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function perform_step!(integrator, cache::Exprb43Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack utilde,tmp,rtmp,Au,F2,F3,J,KsCache = cache
  calc_J!(J, integrator, cache)
  alg = unwrap_alg(integrator, true)

  F1 = integrator.fsalfirst
  mul!(Au, J, uprev)
  halfdt = dt/2
  Ks, phiv_cache, ws = KsCache
  w1_half, w1, w2, w3 = ws
  kwargs = (m=min(alg.m, size(J,1)), opnorm=integrator.opts.internalopnorm, iop=alg.iop)
  # Krylov for F1
  arnoldi!(Ks, J, F1; kwargs...)
  phiv!(w1_half, halfdt, Ks, 1; cache=phiv_cache)
  phiv!(w1, dt, Ks, 4; cache=phiv_cache)
  @muladd @.. broadcast=false @views tmp = uprev + halfdt * w1_half[:, 2] # tmp is U2
  _compute_nl!(F2, f, tmp, p, t + halfdt, J, rtmp); F2 .+= Au
  integrator.destats.nf += 1
  # Krylov for F2
  arnoldi!(Ks, J, F2; kwargs...)
  phiv!(w2, dt, Ks, 4; cache=phiv_cache)
  @muladd @.. broadcast=false @views tmp = uprev + dt * w2[:, 2] # tmp is U3
  _compute_nl!(F3, f, tmp, p, t + dt, J, rtmp); F3 .+= Au
  integrator.destats.nf += 1
  # Krylov for F3 (third column)
  arnoldi!(Ks, J, F3; kwargs...)
  phiv!(w3, dt, Ks, 4; cache=phiv_cache)
  # Update u
  @views @.. broadcast=false rtmp = w1[:,2] - 14w1[:,4] + 36w1[:,5] + 16w2[:,4] - 48w2[:,5] - 2w3[:,4] + 12w3[:,5]
  @muladd @.. broadcast=false u = uprev + dt * rtmp
  if integrator.opts.adaptive
    @views @.. broadcast=false rtmp = 36w1[:,5] - 48w2[:,5] + 12w3[:,5]
    @.. broadcast=false utilde = dt * rtmp
    calculate_residuals!(tmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(tmp,t)
  end

  # Update integrator state
  f(integrator.fsallast, u, p, t + dt)
  integrator.destats.nf += 1
  # integrator.k is automatically set due to aliasing
end

######################################################
# Multistep exponential integrators
function initialize!(integrator,cache::ETD2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Pre-start fsal
  lin = integrator.f.f1(integrator.uprev,integrator.p,integrator.t)
  nl = integrator.f.f2(integrator.uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  nlprev = zero(nl) # to be computed in the first iteration via ETD1
  integrator.fsalfirst = ETD2Fsal(lin, nl, nlprev)

  # Avoid undefined entries if k is an array of arrays
  rate_prototype = lin
  integrator.fsallast = ETD2Fsal(rate_prototype)
  integrator.k[1] = lin + nl
  integrator.k[2] = zero(rate_prototype)
end

function perform_step!(integrator,cache::ETD2ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,f,p = integrator
  @unpack lin,nl,nlprev = integrator.fsalfirst
  @unpack exphA,phihA,B1,B0 = cache
  integrator.k[1] = lin + nl

  if integrator.iter == 1 # ETD1 for initial step
    u = exphA*uprev + dt*(phihA*nl)
  else
    u = exphA*uprev + dt*(B1*nl + B0*nlprev)
  end
  integrator.u = u

  # Push the fsal at t+dt
  nlprev = nl
  lin = f.f1(u,p,t+dt)
  nl = f.f2(u,p,t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.k[2] = lin + nl
  @pack! integrator.fsallast = lin, nl, nlprev
end

function initialize!(integrator, cache::ETD2Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  rate_prototype = cache.rtmp1

  # Pre-start fsal
  integrator.fsalfirst = ETD2Fsal(rate_prototype)
  @unpack lin,nl = integrator.fsalfirst
  integrator.f.f1(lin,integrator.uprev,integrator.p,integrator.t)
  integrator.f.f2(nl,integrator.uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = ETD2Fsal(rate_prototype)
  integrator.k[1] = lin + nl
  integrator.k[2] = zero(rate_prototype)
end

function perform_step!(integrator, cache::ETD2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack lin,nl,nlprev = integrator.fsalfirst
  @unpack utmp,rtmp1,rtmp2,exphA,phihA,B1,B0 = cache
  @.. broadcast=false integrator.k[1] = lin + nl

  if integrator.iter == 1 # ETD1 for initial step
    mul!(utmp, exphA, uprev)
    mul!(rtmp1, phihA, nl)
    @muladd @.. broadcast=false u = utmp + dt*rtmp1
  else
    mul!(utmp, exphA, uprev)
    mul!(rtmp1, B1, nl)
    mul!(rtmp2, B0, nlprev)
    @muladd @.. broadcast=false u = utmp + dt*(rtmp1 + rtmp2)
  end

  # Push the fsal at t+dt
  fsallast = integrator.fsallast
  fsallast.nlprev .= nl
  f.f1(fsallast.lin,u,p,t+dt)
  f.f2(fsallast.nl,u,p,t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  @.. broadcast=false integrator.k[2] = fsallast.lin + fsallast.nl
end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
include("Utils.jl")
using LinearAlgebra
using Distributions
using Distances
using Random
using Match
using JLD
using BlockArrays
# Random.seed!(123)

NTrain = 20000
NTest = 100
N = NTrain + NTest
d = 2
l = 0.5
# For now, only properly implemented for $σ = 0$.
σ = 0.0 * exp.(randn(N))
# σ = 1e-2 * ones(N)
ρList = [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]
λ = 1.5
ρ = 3.0
maternOrder = "32"

#Creating arrays to hold the results

meanErrorClosePredFirst = zeros(length(ρList))
medianErrorClosePredFirst = zeros(length(ρList))
meanErrorClosePredLast = zeros(length(ρList))
medianErrorClosePredLast = zeros(length(ρList))
meanErrorFarPredFirst = zeros(length(ρList))
medianErrorFarPredFirst = zeros(length(ρList))
meanErrorFarPredLast = zeros(length(ρList))
medianErrorFarPredLast = zeros(length(ρList))

# Creating uniform random data
xTrain = rand(d, NTrain)
xTestClose = rand(d, NTest)
xTestFar = 1.1 * ones(d,NTest) .+ 0.1 * rand(d, NTest)

σ .= max.(σ, eps(Float64))
σ = PseudoBlockArray(σ, [NTest, NTrain])

covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)

# Make sure to avoid NANs
@match maternOrder begin
  "12" => (covfunc! = covfunc12!)
  "32" => (covfunc! = covfunc32!)
  "52" => (covfunc! = covfunc52!)
end

# =====================================
# Starting with interpolation:
# =====================================
# xTest = xTestClose
# xTest = xTestFar
# x = hcat(xTrain, xTest)

# =====================================
# Screw that, try different dataset:
# =====================================
h = 1 / NTrain
xTrain = h : h : 1
# xTrain = rand(NTrain)
xTrain = vcat(sin.(2 * π * xTrain)',
              2 * cos.(2 * π * xTrain)') # + 0.1 * randn(2, NTrain)
xTest = 0.1 * randn(2, NTest)
x = hcat(xTrain, xTest)

# Creating reference solution
Θ = covfunc!(pairwise(Euclidean(), x, dims=2)) + Diagonal(σ)
L_Θ = cholesky(Matrix(Θ))
y = L_Θ.L * randn(N)
Θ = Θ + Diagonal(σ)

# Creating data
Θ = PseudoBlockArray(Θ, [NTrain, NTest], [NTrain, NTest])
y = PseudoBlockArray(y, [NTrain, NTest])
y[Block(1)] .+= Diagonal(sqrt.(σ[Block(2)])) * randn(NTrain)

# Computing the true posterior mean
μ = Θ[Block(2,1)] * (Θ[Block(1,1)] \ y[Block(1)])
Σ = Θ[Block(2,2)] - Θ[Block(2,1)] * (Θ[Block(1,1)] \ Θ[Block(1,2)])

# for (iρ, ρ) in enumerate(ρList)

  # =====================================
  # Computing posterior mean for prediction variables last
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
  xOrd = x[:, P]
  xOrdTest = xOrd[:, 1 : NTest]
  xOrdTrain = xOrd[:, (NTest + 1) : end]
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  μOrd = μ[P[1 : NTest] .- NTrain]
  ΣOrd = Σ[P[1 : NTest] .- NTrain, P[1 : NTest] .- NTrain]
  yOrd = PseudoBlockArray(y[P], [NTest, NTrain])
  σOrd = PseudoBlockArray(σ[P], [NTest, NTrain])

  # Creating vecchia approximation
  L = assembleL(skeletons, xOrd, covfunc!)
  L22 = L[(NTest + 1) : end, (NTest + 1) : end]
  U22 = squareSparse(L22) + spdiagm(0 => 1 ./ σOrd[Block(2)])
  icholU_high_level!(U22)
  noisecov = NoiseCov(L22, U22, σOrd[Block(2)])
  μPredLast = (L' \ ( L \ vcat(zero(y[Block(2)]), noisecov \ (yOrd[Block(2)]))))[1:NTest]

  ΣPredLast = inv(Matrix(L[1 : NTest, 1 : NTest] * L[1 : NTest, 1 : NTest]'))

  # @show norm(μOrd - μPredLast) / norm(μOrd)
  # @show norm(μOrd - μPredLast)
  # @show norm(diag(ΣOrd) - diag(ΣPredLast))
  @show mean(abs.(μOrd - μPredLast))
  @show mean(abs.(diag(ΣOrd) - diag(ΣPredLast)))

  # =====================================
  # Computing posterior mean for prediction variables first
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
  xOrdTrain = xTrain[:, P]
  xOrd = hcat(xOrdTrain, xTest)
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  μOrd = μ
  ΣOrd = Σ
  yOrd = y[Block(1)][P]
  σOrd = σ[P]

  μPredFirst, ΣPredFirst = predict(skeletons, (NTrain+1): N, yOrd, xOrd, covfunc!, σ)

  # @show norm(μOrd - μPredFirst)
  # @show norm(diag(ΣOrd) - ΣPredFirst)
  @show mean(abs.(μOrd - μPredFirst))
  @show mean(abs.(diag(ΣOrd) - ΣPredFirst))

  # =====================================
  # Don't include prediction variables at all
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
  xOrdTrain = xTrain[:, P]
  xOrd = hcat(xOrdTrain, xTest)
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  L = assembleL(skeletons, xOrdTrain, covfunc!, σ)
  μOrd = μ
  ΣOrd = Σ
  yOrd = y[Block(1)][P]
  σOrd = σ[P]

  TMat = pairwise(Euclidean(), xOrdTrain, xTest, dims=2)
  covfunc!(TMat)
  CMat = pairwise(Euclidean(), xTest, dims=2)
  covfunc!(CMat)

  μNoPred = vec(yOrd' * L * L' * TMat)
  ΣNoPred = CMat - TMat' * L * L' * TMat

  # @show norm(μOrd - μPredFirst) / norm(μOrd)
  # @show norm(μOrd - μNoPred)
  # @show norm(ΣOrd - ΣNoPred)
  @show mean(abs.(μOrd - μNoPred))
  @show mean(abs.(ΣOrd - ΣNoPred))

# save("./out/jld/newname.jld",
#      "ρList", ρList,
#      "meanErrorClosePredFirst", meanErrorClosePredFirst,
#      "medianErrorClosePredFirst", medianErrorClosePredFirst,
#      "meanErrorClosePredLast", meanErrorClosePredLast,
#      "medianErrorClosePredLast", medianErrorClosePredLast,
#      "meanErrorFarPredFirst", meanErrorFarPredFirst,
#      "medianErrorFarPredFirst", medianErrorFarPredFirst,
#      "meanErrorFarPredLast", meanErrorFarPredLast,
#      "medianErrorFarPredLast", medianErrorFarPredLast)


text:
using DEBTrait
using CSV, DataFrames, Statistics, JLD
using DifferentialEquations, Plots
using HypothesisTests

n_polymers              = 0
n_monomers              = 6
n_microbes              = 1
n_enzymes               = 1
n_minerals              = 0

df_isolates             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/data/microbes2traits.csv", DataFrame, missingstring = "N/A")
df_isolates             = filter(row -> row[:Rhizosphere_response] == "positive" || row[:Rhizosphere_response] == "negative", df_isolates)
df_monomers             = CSV.read("/Users/glmarschmann/.julia/dev/DEBTrait/files/thermodynamic_props_class.csv",DataFrame, missingstring = "N/A")

isolate_id              = collect(1:27)
id_positive             = isolate_id[df_isolates.Rhizosphere_response .== "positive"]
id_negative             = isolate_id[df_isolates.Rhizosphere_response .== "negative"]

monomer_id              = collect(1:6)
monomer_concentration   = [1e-7, 5e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1e-1]
#
#
function condition(u,t,integrator)
    BGE = DEBTrait.calc_BGE(zeros(integrator.p.dim), u, integrator.p, 0.0)[1]
    r   = DEBTrait.calc_growth_rate(u, integrator.p)[1]
    #BGE .< 1e-10 || u[1] .< 1e-10 || t > 300.0
    r
end

affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
#
#
BGE                     = zeros(27, 6, 15)
r                       = zeros(27, 6, 15)
m_E                     = zeros(27, 6, 15)
Bio                     = zeros(27, 6, 15)
x                       = zeros(27, 6, 15)
rG_CO2                  = zeros(27, 6, 15)
rM_CO2                  = zeros(27, 6, 15)
rX_CO2                  = zeros(27, 6, 15)
rED_CO2                 = zeros(27, 6, 15)
rED_rG_ratio            = zeros(27, 6, 15)
#
#
#plt = plot(legend=false)
id_high = [1,2,3,4,5,6,7,9,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]
for i in 1:27
    #print(i)
    #
    L_DNA                   = df_isolates.Genome_size[i]*ones(n_microbes)
    rRNA                    = df_isolates.rRNA_genes[i]*ones(n_microbes)
    min_gentime             = df_isolates.Min_gen_time[i]*ones(n_microbes)
    hydrolases              = df_isolates.hydrolases
    #
    z_sugars_p              = df_isolates.z_sugars[i]*ones(1,1)
    z_organics_p            = df_isolates.z_organic_acids[i]*ones(1,1)
    z_aminos_p              = df_isolates.z_amino_acids[i]*ones(1,1)
    z_fattys_p              = df_isolates.z_fatty_acids[i]*ones(1,1)
    z_nucleos_p             = df_isolates.z_nucleotides[i]*ones(1,1)
    z_auxins_p              = df_isolates.z_auxins[i]*ones(1,1)
    genome_distr_p          = vcat(z_sugars_p, z_organics_p, z_aminos_p, z_fattys_p, z_nucleos_p, z_auxins_p)
    #
    for j in 1:6
        local n_monomers              = 1
        id_monomer              = monomer_id[j]
        Isolate                 = BaseGenome(L_DNA, rRNA, min_gentime, genome_distr_p[id_monomer]*ones(n_monomers,n_microbes), hydrolases)
        Monomer                 = BaseMonomer(df_monomers.D_S[id_monomer]*ones(n_monomers), df_monomers.N_C[id_monomer]*ones(n_monomers), df_monomers.N_N[id_monomer]*ones(n_monomers), df_monomers.y_DE[id_monomer]*ones(n_monomers))
        Soil                    = BaseSoil(0.0, 0.0, 1.0, 1*ones(n_microbes))
        #
        MetabolismC             = PMetabolismC{BaseGenome, Array{Float64,1}}(Isolate)
        MetabolismC.y_EM        = ones(n_microbes)

        if df_isolates.Rhizosphere_response[i] .== "positive"
            MetabolismC.y_EX .= 1.1
        else
            MetabolismC.y_EX .= 1.1
        end

        if df_isolates.Rhizosphere_response[i] .== "positive"
            MetabolismC.α .= 0.01
        else
            MetabolismC.α .= 0.1
        end

        MetabolismC.α_n         = DEBTrait.enzyme_production(MetabolismC.α, Isolate.enzyme_distr)[i]*ones(n_microbes)
        #
        AssimilationC           = PAssimilation(Isolate, Monomer, Soil, MetabolismC)
        #
        Setup                   = SetupBatchC(Isolate, n_polymers, n_monomers, n_microbes, n_enzymes, n_minerals)
        p                       = ParamsBatchC(Setup, Monomer, MetabolismC, AssimilationC)
        #########################
        p.n_polymers              = 0
        p.n_monomers              = 1
        p.n_microbes              = 1
        p.n_enzymes               = 1
        p.n_minerals              = 0
        #
        u0 = zeros(p.dim)

        for k in 1:13
            init_D0 = monomer_concentration[k]
            u0[1+p.n_polymers:p.n_polymers+p.n_monomers] = init_D0*ones(p.n_monomers)
            u0[1+p.n_polymers+p.n_monomers:p.n_polymers+p.n_monomers+p.n_microbes] = p.Setup.u0_reserves*1e-4
            u0[1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes] = p.Setup.u0_structures*1e-4
            #
            tspan = (0.0,  500.0)
            prob  = ODEProblem(rhs_BatchC!, u0, tspan, p)
            sol   = solve(prob, alg_hints=[:stiff], callback=cb)
            #
            #plot!(sol, vars=1+p.n_polymers:p.n_polymers+p.n_monomers, yscale=:log10)
            #plot!(sol, vars=1+p.n_polymers+p.n_monomers+2*p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes+p.n_enzymes)
            #plot!(sol, vars= 1+p.n_polymers+p.n_monomers+p.n_microbes:p.n_polymers+p.n_monomers+2*p.n_microbes)
            #
            BGE_tmp = [DEBTrait.calc_BGE(zeros(p.dim), sol.u[i], p, 0.0)[1] for i in 1:size(sol.t,1)]
            BGE[i,j,k] = median(BGE_tmp)
            #plot!(sol.t, BGE_tmp)
            #
            r_tmp   = [DEBTrait.calc_growth_rate(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            r[i,j,k] = median(r_tmp)
            #plot!(sol.t, r_tmp)
            #
            m_E_tmp = [DEBTrait.calc_reserve_density(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            m_E[i,j,k] = median(m_E_tmp)
            #
            Bio_tmp = [DEBTrait.calc_specific_biomass(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            Bio[i,j,k] = median(Bio_tmp)
            #
            x_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[1][1] for i in 1:size(sol.t,1)]
            x[i,j,k] = median(x_tmp)
            #plot!(sol.t, x_tmp)
            #
            rG_CO2_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[2][1] for i in 1:size(sol.t,1)]
            rG_CO2[i,j,k] = median(rG_CO2_tmp)
            #
            rM_CO2_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[3][1] for i in 1:size(sol.t,1)]
            rM_CO2[i,j,k] = median(rM_CO2_tmp)
            #
            rX_CO2_tmp = [DEBTrait.calc_growth_production(sol.u[i], p)[4][1] for i in 1:size(sol.t,1)]
            rX_CO2[i,j,k] = median(rX_CO2_tmp)
            #
            rED_CO2_tmp = [DEBTrait.calc_specific_assimilation_respiration(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            rED_CO2[i,j,k] = median(rED_CO2_tmp)
            #
            rED_rG_ratio_tmp = [DEBTrait.calc_assimilation_growth_ratio(sol.u[i], p)[1] for i in 1:size(sol.t,1)]
            rED_rG_ratio[i,j,k] = median(rED_rG_ratio_tmp)
            #display(plt)
        end
    end
end

BGE_pos_high = reshape(BGE[id_positive,1,15], 19*1)
BGE_neg_high = reshape(BGE[id_negative,1,15], 8*1)
t = KruskalWallisTest(BGE_pos_high, BGE_neg_high)
median(BGE_pos_high)
median(BGE_neg_high)

BGE_pos_high[BGE_pos_high.<1e-4] .= 0.0
BGE_neg_high[BGE_neg_high.<1e-4] .= rand(1)/0.5
t = KruskalWallisTest(BGE_pos_high, BGE_neg_high)
median(BGE_pos_high)
median(BGE_neg_high)
t = KruskalWallisTest(BGE_pos_high, BGE_neg_high)

BGE_pos_low = reshape(BGE[id_positive,1:6,1], 19*6)
BGE_neg_low = reshape(BGE[id_negative,1:6,1], 8*6)
t = KruskalWallisTest(BGE_pos_low, BGE_neg_low)
median(BGE_pos_low)
median(BGE_neg_low)
BGE_pos_low[BGE_pos_low.<1e-4] .= 0.0
BGE_neg_low[BGE_neg_low.<1e-4] .= 0.0

BGE_high = vcat(BGE_pos_high, BGE_neg_high)
BGE_low = vcat(BGE_pos_low, BGE_neg_low)

id_high = vcat(fill("positive", 19*6), fill("negative", 8*6))
id_low = vcat(fill("positive", 19*6), fill("negative", 8*6))

df_out = DataFrame()
df_out.concentration = vcat(Conc_1, Conc_4)
df_out.BGE = vcat(BGE_low, BGE_high)
df_out.response = vcat(id_high, id_low)
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/seaborn_BGE_select.csv", df_out)

Conc_1 = fill(1e-7, 27*6)
Conc_2 = fill(1e-5, 27*6)
Conc_3 = fill(1e-3, 27*6)
Conc_4 = fill(1e-1, 27*6)
Conc = vcat(Conc_1, Conc_2, Conc_3, Conc_4)

BGE_1 = reshape(BGE[:,1:6,1], 27*6)
BGE_2 = reshape(BGE[:,1:6,5], 27*6)
BGE_3 = reshape(BGE[:,1:6,9], 27*6)
BGE_4 = reshape(BGE[:,1:6,15], 27*6)
BGE_out = vcat(BGE_1, BGE_2, BGE_3, BGE_4)
BGE_out[BGE_out .<1e-4] .= 0.0

r_1 = reshape(r[:,1:6,1]./(log(2)./(df_isolates.Min_gen_time)*24), 27*6)
r_2 = reshape(r[:,1:6,5]./(log(2)./(df_isolates.Min_gen_time)*24), 27*6)
r_3 = reshape(r[:,1:6,9]./(log(2)./(df_isolates.Min_gen_time)*24), 27*6)
r_4 = reshape(r[:,1:6,15]./(log(2)./(df_isolates.Min_gen_time)*24), 27*6)

r_out = vcat(r_1, r_2, r_3, r_4)

t = KruskalWallisTest(BGE_1, r_1)

df_out = DataFrame()
df_out.concentration = Conc
df_out.BGE = BGE_out
df_out.r = r_out
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/manova_r_BGE_select.csv", df_out)

BGE[BGE .<1e-4] .= 0.0
response_tmp = repeat(df_isolates.Rhizosphere_response, 12)
response = replace(response, "negative" => 2, "positive" => 1)
df_out1 = DataFrame()
df_out1.BGE = vcat(reshape(BGE[1:27,1:6,1], 27*6), reshape(BGE[1:27,1:6,15], 27*6))
df_out1.r = vcat(reshape(r[1:27,1:6,1], 27*6), reshape(r[1:27,1:6,15], 27*6))
df_out1.response = response
df_out1.rM = vcat(reshape(rM_CO2[1:27, 1:6, 1], 27*6), reshape(rM_CO2[1:27, 1:6, 13], 27*6))
df_out1.rG = vcat(reshape(rG_CO2[1:27, 1:6, 1], 27*6), reshape(rG_CO2[1:27, 1:6, 13], 27*6))
df_out1.rX = vcat(reshape(rX_CO2[1:27, 1:6, 1], 27*6), reshape(rX_CO2[1:27, 1:6, 13], 27*6))
df_out1.concentration = vcat(fill(1, 27*6), fill(2, 27*6))
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/adonis_BGE_select.csv", df_out1)

val_stack = vcat(BGE_out, r_out)
meas_stack = vcat(fill("BGE", 648), fill("r", 648))
df_out2 = DataFrame()
df_out2.val = val_stack
df_out2.meas = meas_stack
df_out2.response = repeat(df_isolates.Rhizosphere_response, 6*8)
df_out2.concentration = vcat(Conc, Conc)
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/seaborn_BGE_r_select.csv", df_out2)

df_out3 = DataFrame()
df_out3.BGE = BGE_out
df_out3.response = repeat(df_isolates.Rhizosphere_response, 6*4)
df_out3.concentration = Conc
CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/seaborn_BGE_select.csv", df_out3)

p = plot(legend=false
for j in 1:6
    for i in 1:19
        id = id_positive[i]
        scatter!(monomer_concentration, BGE[id,j,:], xscale=:log10, color="green")
        display(p)
    end

    for i in 1:8
        id = id_negative[i]
        scatter!(monomer_concentration, BGE[id,j,:], xscale=:log10, color="red")
        display(p)
    end
end

xlabel!("Substrate concentration [mol/m3]")
ylabel!("Growth efficiency")

r_med_sugar = median(BGE[id_positive,1,:], dims=1)
r_med_organics = median(BGE[id_positive,2,:], dims=1)
r_med_aminos = median(BGE[id_positive,3,:], dims=1)
r_med_fattys = median(BGE[id_positive,4,:], dims=1)
r_med_nucleos = median(BGE[id_positive,5,:], dims=1)
r_med_auxins = median(BGE[id_positive,6,:], dims=1)
r_med_pos = (r_med_sugar+r_med_organics+r_med_aminos+r_med_fattys+r_med_nucleos+r_med_auxins)./6
scatter!(monomer_concentration, r_med_pos[1,:], markersize=8, color="green")

r_med_sugar = median(BGE[id_negative,1,:], dims=1)
r_med_organics = median(BGE[id_negative,2,:], dims=1)
r_med_aminos = median(BGE[id_negative,3,:], dims=1)
r_med_fattys = median(BGE[id_negative,4,:], dims=1)
r_med_nucleos = median(BGE[id_negative,5,:], dims=1)
r_med_auxins = median(BGE[id_negative,6,:], dims=1)
r_med_neg = (r_med_sugar+r_med_organics+r_med_aminos+r_med_fattys+r_med_nucleos+r_med_auxins)./6
scatter!(monomer_concentration, r_med_neg[1,:], markersize=8, color="red")

L_DNA                   = df_isolates.Genome_size
rRNA                    = df_isolates.rRNA_genes
min_gentime             = df_isolates.Min_gen_time
hydrolases              = df_isolates.hydrolases
#
z_sugars_p              = reshape(df_isolates.z_sugars, 1,27)
z_organics_p            = reshape(df_isolates.z_organic_acids, 1,27)
z_aminos_p              = reshape(df_isolates.z_amino_acids, 1,27)
z_fattys_p              = reshape(df_isolates.z_fatty_acids, 1,27)
z_nucleos_p             = reshape(df_isolates.z_nucleotides, 1,27)
z_auxins_p              = reshape(df_isolates.z_auxins, 1,27)
genome_distr_p          = vcat(z_sugars_p, z_organics_p, z_aminos_p, z_fattys_p, z_nucleos_p, z_auxins_p)

Isolate                 = BaseGenome(L_DNA, rRNA, min_gentime, genome_distr_p, hydrolases)

df_out4 = DataFrame()
df_out4.L_DNA = Isolate.L_DNA
df_out4.rRNA = Isolate.rRNA
df_out4.min_gentime = Isolate.min_gentime
df_out4.hydrolases = Isolate.enzyme_distr
df_out4.response = df_isolates.

hydro_pos = Isolate.enzyme_distr[id_positive]
hydro_neg = Isolate.enzyme_distr[id_negative]
t = KruskalWallisTest(hydro_pos, hydro_neg)

rrn_pos = Isolate.rRNA[id_positive]
rrn_neg = Isolate.rRNA[id_negative]
t = KruskalWallisTest(rrn_pos, rrn_neg)

CSV.write("/Users/glmarschmann/.julia/dev/DEBTrait/files/BatchC/seaborn_genome_traits.csv", df_out4)


text:
# Molecular Dynamics Simulation Abstract Interface

"""
    MolecularDynamicsSimulator
Abstract type to be extended by all concrete molecular dynamics simulators.
"""
abstract type MolecularDynamicsSimulator end

"""
    simulate(system::AbstractSystem{3}, simulator::MolecularDynamicsSimulator, potential::AbstractPotential)::MolecularDynamicsResult
Run a molecular dynamics simulation configured with a particular simulator and potential with any abstract system.

An implementer of this API should implement a method of this function for their custom simulator type.
If the simulator has a fast path for some types of potential, those should be implemented with multiple dispatch.
"""
simulate(system::AbstractSystem{3}, simulator::MolecularDynamicsSimulator, potential::AbstractPotential) = throw(UnimplementedError(:simulate, simulator))


text:
function show()
    println("b")
    println("$(2)b")
    a::Float64 = 2.1
    println(a)
    b = 2.1
    println(b)
    c = true
    println(c ? ("True") : ("False"))
    name = "test"
    val = true
    write(sys.stdout, "$(name)_vtables_dispatch_ = $(val)")
end

if abspath(PROGRAM_FILE) == @__FILE__
    show()
end