text:
function createClients(apiUrl, userId, robotId, sessionId)
  client = NavAbilityHttpsClient(apiUrl)
  context = Client(userId,robotId,sessionId)
  return client, context
end

function exampleGraph1D(client, context; doSolve=true)
    variables = [
        Variable("x0", :ContinuousScalar),
        Variable("x1", :ContinuousScalar),
        Variable("x2", :ContinuousScalar),
        Variable("x3", :ContinuousScalar),
    ]
    factors = [
        Factor( "x0f1",
                "Prior",
                ["x0"],
                PriorData( Z=Normal(0, 1) )
        ),
        Factor(
            "x0x1f1",
            "LinearRelative",
            ["x0", "x1"],
            LinearRelativeData(Z=Normal(10, 0.1)),
        ),
        Factor(
            "x1x2f1",
            "Mixture",
            ["x1", "x2"],
                MixtureData(
                    "LinearRelative",
                    (;hypo1=Normal(0, 2),hypo2=Uniform(30, 55)),
                    [0.4, 0.6],
                    2,
                )
        ),
        Factor(
            "x2x3f1",
            "LinearRelative",
            ["x2", "x3"],
            LinearRelativeData(Z=Normal(-50, 1)),
        ),
        Factor(
            "x3x0f1",
            "LinearRelative",
            ["x3", "x0"],
            LinearRelativeData(Z=Normal(40, 1)),
        ),
    ]
    # Variables
    @info "[Fixture] Adding variables, waiting for completion"
    resultIds = Task[]
    for v in variables
      push!(resultIds, addVariable(client, context, v))
    end
    waitForCompletion(client, resultIds; expectedStatuses=["Complete"])

    # Add the factors
    @info "[Fixture] Adding factors, waiting for completion"
    resultIds = Task[]
    for f in factors
      push!(resultIds, addFactor(client, context, f))
    end
    waitForCompletion(client, resultIds; expectedStatuses=["Complete"])

    if doSolve
      @info "[Fixture] solving, waiting for completion"
      resultId = solveSession(client, context)
      waitForCompletion(client, Task[resultIds;]; expectedStatuses=["Complete"])
    end

    # and done
    return (client, context, variables, factors)
end
# vals = getVariables(client, context; detail=SUMMARY) .|> x->x["ppes"][1]["suggested"]

# @pytest.fixture(scope="module")
# async def example_1d_graph_solved(example_1d_graph):
#     """Get the graph after it has been solved.
#     NOTE this changes the graph, so tests need to be defensive.
#     """
#     navability_https_client, client, variables, factors = example_1d_graph
#     logging.info(f"[Fixture] Solving graph, client = {client.dumps()}")
#     requestId = await solveSession(navability_https_client, client)
#     await waitForCompletion(navability_https_client, [requestId], maxSeconds=180)
#     return (navability_https_client, client, variables, factors)

# @pytest.fixture(scope="module")
# async def example_2d_graph(
#     navability_https_client: NavAbilityClient, client_2d: Client
# ):
#     variables = [
#         Variable("x0", VariableType.Pose2.value),
#         Variable("x1", VariableType.Pose2.value),
#         Variable("x2", VariableType.Pose2.value),
#         Variable("l0", VariableType.Point2.value),
#     ]
#     factors = [
#         Factor(
#             "x0f1",
#             "PriorPose2",
#             ["x0"],
#             FactorData(
#                 fnc=PriorPose2(
#                     FullNormal(np.zeros(3), np.diag([0.1, 0.1, 0.1]))
#                 ).dump()  # This is a generator for a PriorPose2
#             ),
#         ),
#         Factor(
#             "x0x1f1",
#             "Pose2Pose2",
#             ["x0", "x1"],
#             FactorData(
#                 fnc=Pose2Pose2(
#                     FullNormal([1, 1, np.pi / 3], np.diag([0.1, 0.1, 0.1]))
#                 ).dump()  # This is a generator for a PriorPose2
#             ),
#         ),
#         Factor(
#             "x1x2f1",
#             "Pose2Pose2",
#             ["x1", "x2"],
#             FactorData(
#                 fnc=Pose2Pose2(
#                     FullNormal([1, 1, np.pi / 3], np.diag([0.1, 0.1, 0.1]))
#                 ).dump()  # This is a generator for a PriorPose2
#             ),
#         ),
#         # TODO: Improve problem setup in future.
#         Factor(
#             "l0f1",
#             "PriorPoint2",
#             ["l0"],
#             FactorData(
#                 fnc=PriorPoint2(FullNormal(np.asarray([5, 0]), np.diag([2, 2]))).dump()
#             ),
#         ),
#         Factor(
#             "x0l0f1",
#             "Point2Point2Range",
#             ["x0", "l0"],
#             FactorData(fnc=Point2Point2Range(Normal(5, 0.1)).dump()),  # Range
#         ),
#         Factor(
#             "x0l0f2",
#             "Pose2Point2BearingRange",
#             ["x0", "l0"],
#             FactorData(
#                 fnc=Pose2Point2BearingRange(
#                     Normal(0, 0.3), Normal(5, 0.1)  # Bearing, range
#                 ).dump()
#             ),
#         ),
#     ]
#     # Variables
#     result_ids = [
#         await addVariable(navability_https_client, client_2d, v) for v in variables
#     ] + [await addFactor(navability_https_client, client_2d, f) for f in factors]

#     logging.info(f"[Fixture] Adding variables and factors, waiting for completion")

#     # Await for only Complete messages, otherwise fail.
#     await waitForCompletion(
#         navability_https_client,
#         result_ids,
#         expectedStatuses=["Complete"],
#         maxSeconds=120,
#     )

#     return (navability_https_client, client_2d, variables, factors)

# @pytest.fixture(scope="module")
# async def example_2d_graph_solved(example_2d_graph):
#     """Get the graph after it has been solved.
#     NOTE this changes the graph, so tests need to be defensive.
#     """
#     navability_https_client, client, variables, factors = example_2d_graph
#     logging.info(f"[Fixture] Solving graph, client = {client.dumps()}")
#     requestId = await solveSession(navability_https_client, client)
#     await waitForCompletion(navability_https_client, [requestId], maxSeconds=180)
#     return (navability_https_client, client, variables, factors)


text:
"""
    abstract DataView{TElem, TData} <: AbstractVector{TElem}

Baseclass for all vector-like views of some data structure.
This allow for example to see some design matrix as a vector of
individual observation-vectors instead of one matrix.
see `MLDataPattern.ObsView` and `MLDataPattern.BatchView` for examples.
"""
abstract type DataView{TElem, TData} <: AbstractVector{TElem} end

"""
    abstract AbstractObsView{TElem, TData} <: DataView{TElem, TData}

Baseclass for all vector-like views of some data structure,
that views it as some form or vector of observations.
see `MLDataPattern.ObsView` for a concrete example.
"""
abstract type AbstractObsView{TElem, TData} <: DataView{TElem, TData} end

"""
    abstract AbstractBatchView{TElem, TData} <: DataView{TElem, TData}

Baseclass for all vector-like views of some data structure,
that views it as some form or vector of equally sized batches.
see `MLDataPattern.BatchView` for a concrete example.
"""
abstract type AbstractBatchView{TElem, TData} <: DataView{TElem, TData} end

# --------------------------------------------------------------------

"""
    abstract DataIterator{TElem,TData}

Baseclass for all types that iterate over a `data` source
in some manner. The total number of observations may or may
not be known or defined and in general there is no contract that
`getobs` or `nobs` has to be supported by the type of `data`.
Furthermore, `length` should be used to query how many elements
the iterator can provide, while `nobs` may return the underlying
true amount of observations available (if known).
see `MLDataPattern.RandomObs`, `MLDataPattern.RandomBatches`
"""
abstract type DataIterator{TElem,TData} end

"""
    abstract ObsIterator{TElem,TData} <: DataIterator{TElem,TData}

Baseclass for all types that iterate over some data source
one observation at a time.
```julia
using MLDataPattern
@assert typeof(RandomObs(X)) <: ObsIterator
for x in RandomObs(X)
    # ...
end
```
see `MLDataPattern.RandomObs`
"""
abstract type ObsIterator{TElem,TData} <: DataIterator{TElem,TData} end

"""
    abstract BatchIterator{TElem,TData} <: DataIterator{TElem,TData}

Baseclass for all types that iterate over of some data source one
batch at a time.
```julia
@assert typeof(RandomBatches(X, size=10)) <: BatchIterator

for x in RandomBatches(X, size=10)
    @assert nobs(x) == 10
    # ...
end
```
see `MLDataPattern.RandomBatches`
"""
abstract type BatchIterator{TElem,TData} <: DataIterator{TElem,TData} end

# just for dispatch for those who care to
const AbstractDataIterator{E,T}  = Union{DataIterator{E,T}, DataView{E,T}}
const AbstractObsIterator{E,T}   = Union{ObsIterator{E,T},  AbstractObsView{E,T}}
const AbstractBatchIterator{E,T} = Union{BatchIterator{E,T},AbstractBatchView{E,T}}

ext: jl


text:
module CirruSepal

  using CirruParser
  export compile, evaluateFile

  function transformProgram(ast)
    exprs = map(transform, ast)
    return Expr(:toplevel, exprs...)
  end

  function transform(x)
    if typeof(x) <: String
      if ismatch(r"^\d+$", x)
        return int(x)
      elseif x == "true"
        return true
      elseif x == "false"
        return false
      elseif ismatch(r"^:", x)
        return string(x[2:end])
      else
        return symbol(x)
      end
    elseif typeof(x) <: Array
      head = x[1]
      exprs = map(transform, x[2:end])
      return Expr(symbol(head), exprs...)
    else
      return :error
    end
  end

  function compile(code)
    tree = CirruParser.pare(code, "")
    transformProgram(tree)
  end

  function evaluateFile(relativeName)
    code = readall(open(relativeName))
    ast = compile(code)
    eval(ast)
  end

end # module


text:
# Test hashing algorithm

using BenchmarkTools

# Dictionary for approach 1
dict1 = Dict((1:10^6) .=> zip(rand(1:10, 1_000_000),
                           rand(1:1_000_000, 1_000_000),
                           rand(1:1_000_000, 1_000_000),
                           rand(1:1_000_000, 1_000_000),
                           rand(1:1_000_000, 1_000_000)))

@benchmark $dict1[1000]
@benchmark $dict1[100]
@benchmark $dict1[14203]
# This is consistently about 25 ns per lookup
25e-9*1_000_000 # -> about 25 ms to lookup entire dictionary

# Customish solution
# https://discourse.julialang.org/t/compiling-to-branch-table/16599
const CODES = 1:10^6

function findcode(x)
    index = findfirst(isequal(x), CODES)
    index = index == nothing ? 0 : index
end
# This does _not_ scale well.

# Here's a lookup table version (from same url)
lookuptable(v) = (a = zeros(Int, maximum(v)); for (i,x) in enumerate(v); a[x]=i; end; a)
const LUT = lookuptable(CODES)

function findcode_lookuptable(x)
    x < 1 || x > length(LUT) ? 0 : LUT[x]
end

@benchmark findcode_lookuptable(1323)

singledict = Dict(1:10^6 .=> 1:10^6)
@benchmark $singledict[1323]
# The lookuptable approache averages about 10 ns and is faster than a regular dictionary

# With 5 indices...
function lookuptable1(v)
    a = zeros(NTuple{5,Int}, maximum(v))
    for x in eachindex(v)
        a[x] = (rand(1:10), rand(1:1_000_000), rand(1:1_000_000), rand(1:1_000_000), rand(1:1_000_000))
    end
    return a
end
const LUT1 = lookuptable(CODES)

function findcode_lookuptable1(x)
    x < 1 || x > length(LUT1) ? 0 : LUT1[x]
end

@benchmark findcode_lookuptable1(10232)
# This is still faster than a dictionary

# What if the dict is constant?
const dict1c = Dict((1:10^6) .=> zip(rand(1:10, 1_000_000),
                                    rand(1:1_000_000, 1_000_000),
                                    rand(1:1_000_000, 1_000_000),
                                    rand(1:1_000_000, 1_000_000),
                                    rand(1:1_000_000, 1_000_000)))

@benchmark $dict1c[1000]  # why is this faster without $?
# without, this is about 15ns

# TODO: How does size of tuple scale wrt time?


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    function tα_qubit(β, ψ1, ψ2, fα, f)
        2 + 2*β - cos(ψ1) - cos(ψ2) - 2*β*cos(π*fα)*cos(2*π*f + π*fα - ψ1- ψ2)
    end
    ψ1 = ψ2 = LinRange(0, 4*π, 100)
    z = [tα_qubit(.61, x, y, 0.2, .1) for x in ψ1, y in ψ2]

    fig = Figure(resolution = (470, 400))
    ax = Axis(fig, aspect = 1, xlabel = "ψ1", ylabel = "ψ2")
    clines = contour!(ψ1, ψ2, z,colormap = :plasma,levels = 20,linewidth = 1.5)
    limits!(ax, 0, 4π, 0, 4π)
    cbar  = Colorbar(fig, clines, label ="α-q", height = Relative(3.55/4))
    fig[1, 1] = ax
    fig[1, 2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "contourQubit.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
"""
    struct Johnson <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use the [Johnson algorithm](https://en.wikipedia.org/wiki/Johnson%27s_algorithm).
No additional configuration parameters are specified or required.

### Implementation Notes
`Johnson` supports the following shortest-path functionality:
- non-negative distance matrices / weights
- all-pairs shortest paths

### Performance
Complexity: O(|V|*|E|)
"""
struct Johnson <: ShortestPathAlgorithm end

struct JohnsonResult{T, U<:Integer} <: ShortestPathResult
    parents::Matrix{U}
    dists::Matrix{T}
end

function shortest_paths(g::AbstractGraph{U}, distmx::AbstractMatrix{T}, ::Johnson) where {T, U<:Integer}
    nvg = nv(g)
    type_distmx = typeof(distmx)
    #Change when parallel implementation of Bellman Ford available
    wt_transform = LightGraphs.Experimental.ShortestPaths.dists(shortest_paths(g, vertices(g), distmx, BellmanFord()))

    if !type_distmx.mutable && type_distmx !=  LightGraphs.DefaultDistance
        distmx = sparse(distmx) #Change reference, not value
    end

    #Weight transform not needed if all weights are positive.
    if type_distmx !=  LightGraphs.DefaultDistance
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[src(e)] - wt_transform[dst(e)]
        end
    end

    dists = Matrix{T}(undef, nvg, nvg)
    parents = Matrix{U}(undef, nvg, nvg)
    for v in vertices(g)
        dijk_state = dijkstra_shortest_paths(g, v, distmx)
        dists[v, :] = dijk_state.dists
        parents[v, :] = dijk_state.parents
    end

    broadcast!(-, dists, dists, wt_transform)
    for v in vertices(g)
        dists[:, v] .+= wt_transform[v] #Vertical traversal prefered
    end

    if type_distmx.mutable
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[dst(e)] - wt_transform[src(e)]
        end
    end

    return JohnsonResult(parents, dists)
end

shortest_paths(g::AbstractGraph, alg::Johnson) = shortest_paths(g, weights(g), alg)

function paths(s::JohnsonResult{T, U}, v::Integer) where {T, U <: Integer}
    pathinfo = s.parents[v, :]
    paths = Vector{Vector{U}}()
    for i in 1:length(pathinfo)
        if (i == v) || (s.dists[v, i] == typemax(T))
            push!(paths, Vector{U}())
        else
            path = Vector{U}()
            currpathindex = U(i)
            while currpathindex != 0
                push!(path, currpathindex)
                currpathindex = pathinfo[currpathindex]
            end
            push!(paths, reverse(path))
        end
    end
    return paths
end

paths(s::JohnsonResult) = [paths(s, v) for v in 1:size(s.parents, 1)]
paths(st::JohnsonResult, s::Integer, d::Integer) = paths(st, s)[d]


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2),),
    ((2, 2), :SOCExpPSD),
    ]
insts["fast"] = [
    ((:motzkin2,),),
    ((:motzkin2,), :SOCExpPSD),
    ((:motzkin3,),),
    ((:CS16ex8_13,),),
    ((:CS16ex8_14,),),
    ((:CS16ex18,),),
    ((:CS16ex12,),),
    ((:CS16ex13,),),
    ((:MCW19ex1_mod,),),
    ((:MCW19ex8,),),
    ((:MCW19ex8,), :SOCExpPSD),
    ((3, 2),),
    ((3, 2), :SOCExpPSD),
    ((6, 6),),
    ((20, 3),),
    ((20, 3), :SOCExpPSD),
    ]
insts["various"] = [
    ((:motzkin2,),),
    ((:motzkin3,),),
    ((:CS16ex8_13,),),
    ((:CS16ex8_14,),),
    ((:CS16ex18,),),
    ((:CS16ex12,),),
    ((:CS16ex13,),),
    ((:MCW19ex1_mod,),),
    ((:MCW19ex8,), nothing, relaxed_tols),
    ((6, 3),),
    ((6, 3), :SOCExpPSD),
    ((20, 3), nothing, relaxed_tols),
    ((20, 3), :SOCExpPSD, relaxed_tols),
    ]
return (SignomialMinJuMP, insts)


text:
function initialize!(integrator, cache::Union{Kvaerno3ConstantCache,
                                              KenCarp3ConstantCache,
                                              Kvaerno4ConstantCache,
                                              KenCarp4ConstantCache,
                                              Kvaerno5ConstantCache,
                                              KenCarp5ConstantCache})
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function initialize!(integrator, cache::Union{Kvaerno3Cache,
                                              KenCarp3Cache,
                                              Kvaerno4Cache,
                                              KenCarp4Cache,
                                              Kvaerno5Cache,
                                              KenCarp5Cache})
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::Kvaerno3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32 = cache.tab
  alg = unwrap_alg(integrator, true)

  # calculate W
  update_W!(integrator, cache, γ*dt, repeat_step)

  # FSAL Step 1
  nlsolver.z = z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Add extrapolation for guess
  nlsolver.z = z₂ = z₁

  nlsolver.tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  # Guess is from Hermite derivative on z₁ and z₂
  nlsolver.z = z₃ = α31*z₁ + α32*z₂

  nlsolver.tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  nlsolver.z = z₄ = a31*z₁ + a32*z₂ + γ*z₃ # use yhat as prediction

  nlsolver.tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = 1
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₄

  ################################### Finalize

  if integrator.opts.adaptive
    tmp = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \_vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  integrator.fsallast = z₄./dt
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::Kvaerno3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,k,b,tmp,atmp,nlsolver = cache
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  # FSAL Step 1
  @.. z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Add extrapolation for guess
  @.. z₂ = z₁
  nlsolver.z = z₂

  @.. tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  # Guess is from Hermite derivative on z₁ and z₂
  @.. z₃ = α31*z₁ + α32*z₂
  nlsolver.z = z₃

  @.. tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(cache) <: Kvaerno3Cache
    @.. z₄ = a31*z₁ + a32*z₂ + γ*z₃ # use yhat as prediction
  elseif typeof(cache) <: KenCarp3Cache
    @unpack α41,α42 = cache.tab
    @.. z₄ = α41*z₁ + α42*z₂
  end
  nlsolver.z = z₄

  @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = 1
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₄

  ################################### Finalize

  if integrator.opts.adaptive
    @.. dz = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @.. integrator.fsallast = z₄/dt
end

@muladd function perform_step!(integrator, cache::KenCarp3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32,ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,ebtilde1,ebtilde2,ebtilde3,ebtilde4 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  # calculate W
  update_W!(integrator, cache, γdt, repeat_step)

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    z₁ = dt.*f(uprev, p, t)
  else
    # FSAL Step 1
    z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Add extrapolation for guess
  nlsolver.z = z₂ = z₁

  nlsolver.tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    k1 = dt*integrator.fsalfirst - z₁
    nlsolver.tmp += ea21*k1
  end

  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  if typeof(integrator.f) <: SplitFunction
    z₃ = z₂
    u = nlsolver.tmp + γ*z₂
    k2 = dt*f2(u,p,t + 2γdt)
    integrator.destats.nf2 += 1
    tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
  else
    # Guess is from Hermite derivative on z₁ and z₂
    z₃ = α31*z₁ + α32*z₂
    tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃
  nlsolver.tmp = tmp
  nlsolver.c = c3

  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ = z₂
    u = nlsolver.tmp + γ*z₃
    k3 = dt*f2( u,p,t+c3*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3
  else
    @unpack α41,α42 = cache.tab
    z₄ = α41*z₁ + α42*z₂
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  end
  nlsolver.z = z₄
  nlsolver.c = 1
  nlsolver.tmp = tmp

  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₄
  if typeof(integrator.f) <: SplitFunction
    k4 = dt*f2(u,p,t+dt)
    integrator.destats.nf2 += 1
    u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + eb1*k1 + eb2*k2 + eb3*k3 + eb4*k4
  end

  ################################### Finalize

  if integrator.opts.adaptive
    if typeof(integrator.f) <: SplitFunction
      tmp = btilde1*z₁  + btilde2*z₂  + btilde3*z₃ + btilde4*z₄ + ebtilde1*k1 + ebtilde2*k2 + ebtilde3*k3 + ebtilde4*k4
    else
      tmp = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄
    end
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.k[1] = integrator.fsalfirst
    integrator.fsallast = integrator.f(u, p, t+dt)
    integrator.k[2] = integrator.fsallast
  else
    integrator.fsallast = z₄./dt
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
  end
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::KenCarp3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,k1,k2,k3,k4,k,b,tmp,atmp,nlsolver = cache
  @unpack γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4 = cache.tab
  @unpack ebtilde1,ebtilde2,ebtilde3,ebtilde4 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    if !repeat_step && !integrator.last_stepfail
      f(z₁, integrator.uprev, p, integrator.t)
      z₁ .*= dt
    end
  else
    # FSAL Step 1
    @.. z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Add extrapolation for guess
  @.. z₂ = z₁
  nlsolver.z = z₂

  @.. tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    @.. k1 = dt*integrator.fsalfirst - z₁
    @.. tmp += ea21*k1
  end

  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  if typeof(integrator.f) <: SplitFunction
    z₃ .= z₂
    @.. u = tmp + γ*z₂
    f2(k2, u, p, t + 2γdt); k2 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a31*z₁[i] + a32*z₂[i] + ea31*k1[i] + ea32*k2[i]
    end
  else
    # Guess is from Hermite derivative on z₁ and z₂
    @.. z₃ = α31*z₁ + α32*z₂
    @.. tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃

  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ .= z₂
    @.. u = tmp + γ*z₃
    f2( k3, u,p,t+c3*dt); k3 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a41*z₁[i] + a42*z₂[i] + a43*z₃[i] + ea41*k1[i] + ea42*k2[i] + ea43*k3[i]
    end
  else
    @unpack α41,α42 = cache.tab
    @.. z₄ = α41*z₁ + α42*z₂
    @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  end
  nlsolver.z = z₄

  nlsolver.c = 1
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₄
  if typeof(integrator.f) <: SplitFunction
    f2( k4, u,p,t+dt); k4 .*= dt
    integrator.destats.nf2 += 1
    #@.. u = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄ + eb1*k1 + eb2*k2 + eb3*k3 + eb4*k4
    for i in eachindex(u)
      u[i] = uprev[i] + a41*z₁[i] + a42*z₂[i] + a43*z₃[i] + γ*z₄[i] + eb1*k1[i] + eb2*k2[i] + eb3*k3[i] + eb4*k4[i]
    end
  end

  ################################### Finalize

  if integrator.opts.adaptive
    if typeof(integrator.f) <: SplitFunction
      #@.. dz = btilde1*z₁  + btilde2*z₂  + btilde3*z₃ + btilde4*z₄ + ebtilde1*k1 + ebtilde2*k2 + ebtilde3*k3 + ebtilde4*k4
      for i in eachindex(dz)
        @inbounds dz[i] = btilde1*z₁[i]  + btilde2*z₂[i]  + btilde3*z₃[i] + btilde4*z₄[i] + ebtilde1*k1[i] + ebtilde2*k2[i] + ebtilde3*k3[i] + ebtilde4*k4[i]
      end
    else
      @.. dz = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄
    end
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.f(integrator.fsallast,u,p,t+dt)
  else
    @.. integrator.fsallast = z₄/dt
  end
end

@muladd function perform_step!(integrator, cache::Kvaerno4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,c3,c4 = cache.tab
  @unpack α21,α31,α32,α41,α42 = cache.tab
  @unpack btilde1,btilde2,btilde3,btilde4,btilde5 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt

  # calculate W
  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Add extrapolation choice
  nlsolver.z = z₂ = zero(u)

  nlsolver.tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  nlsolver.z = z₃ = α31*z₁ + α32*z₂

  nlsolver.tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  nlsolver.z = z₄ = α41*z₁ + α42*z₂

  nlsolver.tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  # Use yhat2 for prediction
  nlsolver.z = z₅ = a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄

  nlsolver.tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = 1
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₅

  ################################### Finalize

  if integrator.opts.adaptive
    tmp = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  integrator.fsallast = z₅./dt
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::Kvaerno4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,z₅,k,b,tmp,atmp,nlsolver = cache
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,c3,c4 = cache.tab
  @unpack α21,α31,α32,α41,α42 = cache.tab
  @unpack btilde1,btilde2,btilde3,btilde4,btilde5 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  @.. z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Allow other choices here
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂

  @.. tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  @.. z₃ = α31*z₁ + α32*z₂
  nlsolver.z = z₃

  @.. tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  # Use constant z prediction
  @.. z₄ = α41*z₁ + α42*z₂
  nlsolver.z = z₄

  @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  # Use yhat prediction
  @.. z₅ = a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄
  nlsolver.z = z₅

  @.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = 1
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₅

  ################################### Finalize

  if integrator.opts.adaptive
    @.. dz = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @.. integrator.fsallast = z₅/dt
end

@muladd function perform_step!(integrator, cache::KenCarp4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,c3,c4,c5 = cache.tab
  @unpack α31,α32,α41,α42,α51,α52,α53,α54,α61,α62,α63,α64,α65 = cache.tab
  @unpack btilde1,btilde3,btilde4,btilde5,btilde6 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea42,ea43,ea51,ea52,ea53,ea54,ea61,ea62,ea63,ea64,ea65 = cache.tab
  @unpack eb1,eb3,eb4,eb5,eb6 = cache.tab
  @unpack ebtilde1,ebtilde3,ebtilde4,ebtilde5,ebtilde6 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  # calculate W
  update_W!(integrator, cache, γdt, repeat_step)

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    z₁ = dt.*f(uprev, p, t)
  else
    # FSAL Step 1
    z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Add extrapolation choice
  nlsolver.z = z₂ = zero(u)

  tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    k1 = dt*integrator.fsalfirst - z₁
    tmp += ea21*k1
  end
  nlsolver.tmp = tmp
  nlsolver.c = γ

  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  if typeof(integrator.f) <: SplitFunction
    z₃ = z₂
    u = nlsolver.tmp + γ*z₂
    k2 = dt*f2(u,p,t+2γdt)
    integrator.destats.nf2 += 1
    tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
  else
    # Guess is from Hermite derivative on z₁ and z₂
    z₃ = α31*z₁ + α32*z₂
    tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃
  nlsolver.tmp = tmp
  nlsolver.c = c3

  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ = z₂
    u = nlsolver.tmp + γ*z₃
    k3 = dt*f2( u,p,t+c3*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3
  else
    z₄ = α41*z₁ + α42*z₂
    tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  end
  nlsolver.z = z₄
  nlsolver.tmp = tmp
  nlsolver.c = c4

  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  if typeof(integrator.f) <: SplitFunction
    z₅ = z₄
    u = nlsolver.tmp + γ*z₄
    k4 = dt*f2( u,p,t+c4*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄ + ea51*k1 + ea52*k2 + ea53*k3 + ea54*k4
  else
    z₅ = α51*z₁ + α52*z₂ + α53*z₃ + α54*z₄
    tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  end
  nlsolver.z = z₅
  nlsolver.tmp = tmp
  nlsolver.c = c5

  u = nlsolver.tmp + γ*z₅
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  if typeof(integrator.f) <: SplitFunction
    z₆ = z₅
    u = nlsolver.tmp + γ*z₅
    k5 = dt*f2( u,p,t+c5*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + ea61*k1 + ea62*k2 + ea63*k3 + ea64*k4 + ea65*k5
  else
    z₆ = α61*z₁ + α62*z₂ + α63*z₃ + α64*z₄ + α65*z₅
    tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
  end
  nlsolver.z = z₆
  nlsolver.tmp = tmp
  nlsolver.c = 1

  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₆
  if typeof(integrator.f) <: SplitFunction
    k6 = dt*f2(u,p,t+dt)
    integrator.destats.nf2 += 1
    u = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + γ*z₆ + eb1*k1 + eb3*k3 + eb4*k4 + eb5*k5 + eb6*k6
  end

  ################################### Finalize

  if integrator.opts.adaptive
    if typeof(integrator.f) <: SplitFunction
      tmp = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + ebtilde1*k1 + ebtilde3*k3 + ebtilde4*k4 + ebtilde5*k5 + ebtilde6*k6
    else
      tmp = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆
    end
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.k[1] = integrator.fsalfirst
    integrator.fsallast = integrator.f(u, p, t+dt)
    integrator.k[2] = integrator.fsallast
  else
    integrator.fsallast = z₆./dt
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
  end
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::KenCarp4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,z₅,z₆,k,b,tmp,atmp,nlsolver = cache
  @unpack k1,k2,k3,k4,k5,k6 = cache
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,c3,c4,c5 = cache.tab
  @unpack α31,α32,α41,α42,α51,α52,α53,α54,α61,α62,α63,α64,α65 = cache.tab
  @unpack btilde1,btilde3,btilde4,btilde5,btilde6 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea42,ea43,ea51,ea52,ea53,ea54,ea61,ea62,ea63,ea64,ea65 = cache.tab
  @unpack eb1,eb3,eb4,eb5,eb6 = cache.tab
  @unpack ebtilde1,ebtilde3,ebtilde4,ebtilde5,ebtilde6 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    if !repeat_step && !integrator.last_stepfail
      f(z₁, integrator.uprev, p, integrator.t)
      z₁ .*= dt
    end
  else
    # FSAL Step 1
    @.. z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Allow other choices here
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂

  @.. tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    @.. k1 = dt*integrator.fsalfirst - z₁
    @.. tmp += ea21*k1
  end

  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  tstep = t + c3*dt

  if typeof(integrator.f) <: SplitFunction
    z₃ .= z₂
    @.. u = tmp + γ*z₂
    f2(k2, u, p, t + 2γdt); k2 .*= dt
    integrator.destats.nf2 += 1
    # @.. tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a31*z₁[i] + a32*z₂[i] + ea31*k1[i] + ea32*k2[i]
    end
  else
    # Guess is from Hermite derivative on z₁ and z₂
    @.. z₃ = α31*z₁ + α32*z₂
    @.. tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃

  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ .= z₂
    @.. u = tmp + γ*z₃
    f2( k3, u,p,t+c3*dt); k3 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃ + ea41*k1 + ea42*k2 + ea43*k3
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a41*z₁[i] + a42*z₂[i] + a43*z₃[i] + ea41*k1[i] + ea42*k2[i] + ea43*k3[i]
    end
  else
    @.. z₄ = α41*z₁ + α42*z₂
    @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  end
  nlsolver.z = z₄

  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  tstep = t + c5*dt

  if typeof(integrator.f) <: SplitFunction
    z₅ .= z₄
    @.. u = tmp + γ*z₄
    f2( k4, u,p,t+c4*dt); k4 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄ + ea51*k1 + ea52*k2 + ea53*k3 + ea54*k4
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a51*z₁[i] + a52*z₂[i] + a53*z₃[i] + a54*z₄[i] + ea51*k1[i] + ea52*k2[i] + ea53*k3[i] + ea54*k4[i]
    end
  else
    @.. z₅ = α51*z₁ + α52*z₂ + α53*z₃ + α54*z₄
    #@.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a51*z₁[i] + a52*z₂[i] + a53*z₃[i] + a54*z₄[i]
    end
  end
  nlsolver.z = z₅

  nlsolver.c = c5
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  if typeof(integrator.f) <: SplitFunction
    z₆ .= z₅
    @.. u = tmp + γ*z₅
    f2( k5, u,p,t+c5*dt); k5 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + ea61*k1 + ea62*k2 + ea63*k3 + ea64*k4 + ea65*k5
    for i in eachindex(tmp)
      @inbounds tmp[i] = uprev[i] + a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i] + ea61*k1[i] + ea62*k2[i] + ea63*k3[i] + ea64*k4[i] + ea65*k5[i]
    end
  else
    #@.. z₆ = α61*z₁ + α62*z₂ + α63*z₃ + α64*z₄ + α65*z₅
    #@.. tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
    @tight_loop_macros for i in eachindex(u)
      @inbounds z₆[i] = α61*z₁[i] + α62*z₂[i] + α63*z₃[i] + α64*z₄[i] + α65*z₅[i]
      @inbounds tmp[i] = uprev[i] + a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i]
    end
  end
  nlsolver.z = z₆

  nlsolver.c = 1
  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₆
  if typeof(integrator.f) <: SplitFunction
    f2( k6, u,p,t+dt); k6 .*= dt
    integrator.destats.nf2 += 1
    @.. u = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + γ*z₆ + eb1*k1 + eb3*k3 + eb4*k4 + eb5*k5 + eb6*k6
    for i in eachindex(u)
      u[i] = uprev[i] + a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i] + γ*z₆[i] + eb1*k1[i] + eb3*k3[i] + eb4*k4[i] + eb5*k5[i] + eb6*k6[i]
    end
  end

  ################################### Finalize

  if integrator.opts.adaptive
    if typeof(integrator.f) <: SplitFunction
      #@.. dz = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + ebtilde1*k1 + ebtilde3*k3 + ebtilde4*k4 + ebtilde5*k5 + ebtilde6*k6
      for i in eachindex(u)
        @inbounds dz[i] = btilde1*z₁[i] + btilde3*z₃[i] + btilde4*z₄[i] + btilde5*z₅[i] + btilde6*z₆[i] + ebtilde1*k1[i] + ebtilde3*k3[i] + ebtilde4*k4[i] + ebtilde5*k5[i] + ebtilde6*k6[i]
      end
    else
      # @.. dz = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆
      @tight_loop_macros for i in eachindex(u)
        @inbounds dz[i] = btilde1*z₁[i] + btilde3*z₃[i] + btilde4*z₄[i] + btilde5*z₅[i] + btilde6*z₆[i]
      end
    end

    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.f(integrator.fsallast,u,p,t+dt)
  else
    @.. integrator.fsallast = z₆/dt
  end
end

@muladd function perform_step!(integrator, cache::Kvaerno5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,c3,c4,c5,c6 = cache.tab
  @unpack btilde1,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  @unpack α31,α32,α41,α42,α43,α51,α52,α53,α61,α62,α63 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt

  # calculate W
  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Add extrapolation choice
  nlsolver.z = z₂ = zero(u)

  nlsolver.tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  nlsolver.z = z₃ = α31*z₁ + α32*z₂

  nlsolver.tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  nlsolver.z = z₄ = α41*z₁ + α42*z₂ + α43*z₃

  nlsolver.tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  nlsolver.z = z₅ = α51*z₁ + α52*z₂ + α53*z₃

  nlsolver.tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = c5
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  nlsolver.z = z₆ = α61*z₁ + α62*z₂ + α63*z₃

  nlsolver.tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
  nlsolver.c = c6
  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 7

  # Prediction from embedding
  nlsolver.z = z₇ = a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + γ*z₆

  nlsolver.tmp = uprev + a71*z₁ + a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆
  nlsolver.c = 1
  z₇ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₇

  ################################### Finalize

  if integrator.opts.adaptive
    tmp = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  integrator.fsallast = z₇./dt
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::Kvaerno5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,z₅,z₆,z₇,k,b,tmp,atmp,nlsolver = cache
  @unpack γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,c3,c4,c5,c6 = cache.tab
  @unpack btilde1,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  @unpack α31,α32,α41,α42,α43,α51,α52,α53,α61,α62,α63 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  @.. z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Allow other choices here
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂

  @.. tmp = uprev + γ*z₁
  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  @.. z₃ = α31*z₁ + α32*z₂
  nlsolver.z = z₃

  @.. tmp = uprev + a31*z₁ + a32*z₂
  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  # Use constant z prediction
  @.. z₄ = α41*z₁ + α42*z₂ + α43*z₃
  nlsolver.z = z₄

  @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  @.. z₅ = α51*z₁ + α52*z₂ + α53*z₃
  nlsolver.z = z₅

  @.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = c5
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  @.. z₆ = α61*z₁ + α62*z₂ + α63*z₃
  nlsolver.z = z₆

  @.. tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
  nlsolver.c = c6
  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 7

  # Prediction is embedded method
  # @.. z₇ = a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + γ*z₆
  @tight_loop_macros for i in eachindex(u)
    @inbounds z₇[i] = a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i] + γ*z₆[i]
  end
  nlsolver.z = z₇

  # @.. tmp = uprev + a71*z₁ + a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆
  @tight_loop_macros for i in eachindex(u)
    @inbounds tmp[i] = uprev[i] + a71*z₁[i] + a73*z₃[i] + a74*z₄[i] + a75*z₅[i] + a76*z₆[i]
  end
  nlsolver.c = 1
  z₇ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₇

  ################################### Finalize

  if integrator.opts.adaptive
    # @.. dz = btilde1*z₁ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇
    @tight_loop_macros for i in eachindex(u)
      @inbounds dz[i] = btilde1*z₁[i] + btilde3*z₃[i] + btilde4*z₄[i] + btilde5*z₅[i] + btilde6*z₆[i] + btilde7*z₇[i]
    end
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @.. integrator.fsallast = z₇/dt
end

@muladd function perform_step!(integrator, cache::KenCarp5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  nlsolver = cache.nlsolver
  @unpack γ,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a84,a85,a86,a87,c3,c4,c5,c6,c7 = cache.tab
  @unpack α31,α32,α41,α42,α51,α52,α61,α62,α71,α72,α73,α74,α75,α81,α82,α83,α84,α85 = cache.tab
  @unpack btilde1,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea43,ea51,ea53,ea54,ea61,ea63,ea64,ea65 = cache.tab
  @unpack ea71,ea73,ea74,ea75,ea76,ea81,ea83,ea84,ea85,ea86,ea87 = cache.tab
  @unpack eb1,eb4,eb5,eb6,eb7,eb8 = cache.tab
  @unpack ebtilde1,ebtilde4,ebtilde5,ebtilde6,ebtilde7,ebtilde8 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  # calculate W
  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    z₁ = dt.*f(uprev, p, t)
  else
    # FSAL Step 1
    z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Add extrapolation choice
  nlsolver.z = z₂ = zero(u)

  tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    k1 = dt*integrator.fsalfirst - z₁
    tmp += ea21*k1
  end
  nlsolver.tmp = tmp
  nlsolver.c = γ

  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 3

  tstep = t + c3*dt

  if typeof(integrator.f) <: SplitFunction
    z₃ = z₂
    u = nlsolver.tmp + γ*z₂
    k2 = dt*f2(u,p,t+2γdt)
    integrator.destats.nf2 += 1
    tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
  else
    # Guess is from Hermite derivative on z₁ and z₂
    z₃ = α31*z₁ + α32*z₂
    tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃
  nlsolver.c = c3
  nlsolver.tmp = tmp

  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ = z₂
    u = nlsolver.tmp + γ*z₃
    k3 = dt*f2( u,p,t+c3*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a41*z₁ + a43*z₃ + ea41*k1 + ea43*k3
  else
    z₄ = α41*z₁ + α42*z₂
    tmp = uprev + a41*z₁ + a43*z₃
  end
  nlsolver.z = z₄
  nlsolver.c = c4
  nlsolver.tmp = tmp

  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  if typeof(integrator.f) <: SplitFunction
    z₅ = z₂
    u = nlsolver.tmp + γ*z₄
    k4 = dt*f2( u,p,t+c4*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a51*z₁ + a53*z₃ + a54*z₄ + ea51*k1 + ea53*k3 + ea54*k4
  else
    z₅ = α51*z₁ + α52*z₂
    tmp = uprev + a51*z₁ + a53*z₃ + a54*z₄
  end
  nlsolver.z = z₅
  nlsolver.c = c5
  nlsolver.tmp = tmp

  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  if typeof(integrator.f) <: SplitFunction
    z₆ = z₃
    u = nlsolver.tmp + γ*z₅
    k5 = dt*f2( u,p,t+c5*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + ea61*k1 + ea63*k3 + ea64*k4 + ea65*k5
  else
    z₆ = α61*z₁ + α62*z₂
    tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
  end
  nlsolver.z = z₆
  nlsolver.c = c6
  nlsolver.tmp = tmp

  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 7

  if typeof(integrator.f) <: SplitFunction
    z₇ = z₂
    u = nlsolver.tmp + γ*z₆
    k6 = dt*f2( u,p,t+c6*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a71*z₁ +  a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆ + ea71*k1 + ea73*k3 + ea74*k4 + ea75*k5 + ea76*k6
  else
    z₇ = α71*z₁ + α72*z₂ + α73*z₃ + α74*z₄ + α75*z₅
    tmp = uprev + a71*z₁ +  a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆
  end
  nlsolver.z = z₇
  nlsolver.c = c7
  nlsolver.tmp = tmp

  z₇ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 8

  if typeof(integrator.f) <: SplitFunction
    z₈ = z₅
    u = nlsolver.tmp + γ*z₇
    k7 = dt*f2( u,p,t+c7*dt)
    integrator.destats.nf2 += 1
    tmp = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇ + ea81*k1 + ea83*k3 + ea84*k4 + ea85*k5 + ea86*k6 + ea87*k7
  else
    z₈ = α81*z₁ + α82*z₂ + α83*z₃ + α84*z₄ + α85*z₅
    tmp = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
  end
  nlsolver.z = z₈
  nlsolver.c = 1
  nlsolver.tmp = tmp

  z₈ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₈
  if typeof(integrator.f) <: SplitFunction
    k8 = dt*f2( u,p,t+dt)
    integrator.destats.nf2 += 1
    u = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇ + γ*z₈ + eb1*k1 + eb4*k4 + eb5*k5 + eb6*k6 + eb7*k7 + eb8*k8
  end

  ################################### Finalize

  if integrator.opts.adaptive
    if typeof(integrator.f) <: SplitFunction
      tmp = btilde1*z₁ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈ + ebtilde1*k1 + ebtilde4*k4 + ebtilde5*k5 + ebtilde6*k6 + ebtilde7*k7 + ebtilde8*k8
    else
      tmp = btilde1*z₁ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈
    end
    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      est = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    else
      est = tmp
    end
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.k[1] = integrator.fsalfirst
    integrator.fsallast = integrator.f(u, p, t+dt)
    integrator.k[2] = integrator.fsallast
  else
    integrator.fsallast = z₈./dt
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
  end
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::KenCarp5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,p = integrator
  @unpack dz,z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,k,b,tmp,atmp,nlsolver = cache
  @unpack k1,k2,k3,k4,k5,k6,k7,k8 = cache
  @unpack γ,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a84,a85,a86,a87,c3,c4,c5,c6,c7 = cache.tab
  @unpack α31,α32,α41,α42,α51,α52,α61,α62,α71,α72,α73,α74,α75,α81,α82,α83,α84,α85 = cache.tab
  @unpack btilde1,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
  @unpack ea21,ea31,ea32,ea41,ea43,ea51,ea53,ea54,ea61,ea63,ea64,ea65 = cache.tab
  @unpack ea71,ea73,ea74,ea75,ea76,ea81,ea83,ea84,ea85,ea86,ea87 = cache.tab
  @unpack eb1,eb4,eb5,eb6,eb7,eb8 = cache.tab
  @unpack ebtilde1,ebtilde4,ebtilde5,ebtilde6,ebtilde7,ebtilde8 = cache.tab
  alg = unwrap_alg(integrator, true)

  if typeof(integrator.f) <: SplitFunction
    f = integrator.f.f1
    f2 = integrator.f.f2
  else
    f = integrator.f
  end

  # precalculations
  γdt = γ*dt

  update_W!(integrator, cache, γdt, repeat_step)

  ##### Step 1

  if typeof(integrator.f) <: SplitFunction
    # Explicit tableau is not FSAL
    # Make this not compute on repeat
    if !repeat_step && !integrator.last_stepfail
      f(z₁, integrator.uprev, p, integrator.t)
      z₁ .*= dt
    end
  else
    # FSAL Step 1
    @.. z₁ = dt*integrator.fsalfirst
  end

  ##### Step 2

  # TODO: Allow other choices here
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂

  tstep = t + 2*γdt
  @.. tmp = uprev + γ*z₁

  if typeof(integrator.f) <: SplitFunction
    # This assumes the implicit part is cheaper than the explicit part
    @.. k1 = dt*integrator.fsalfirst - z₁
    @.. tmp += ea21*k1
  end

  nlsolver.c = γ
  z₂ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  set_new_W!(nlsolver, false)

  ################################## Solve Step 3

  tstep = t + c3*dt

  if typeof(integrator.f) <: SplitFunction
    z₃ .= z₂
    @.. u = tmp + γ*z₂
    f2(k2, u, p, t+2γdt); k2 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a31*z₁ + a32*z₂ + ea31*k1 + ea32*k2
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a31*z₁[i] + a32*z₂[i] + ea31*k1[i] + ea32*k2[i]
    end
  else
    # Guess is from Hermite derivative on z₁ and z₂
    @.. z₃ = a31*z₁ + α32*z₂
    @.. tmp = uprev + a31*z₁ + a32*z₂
  end
  nlsolver.z = z₃

  nlsolver.c = c3
  z₃ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 4

  if typeof(integrator.f) <: SplitFunction
    z₄ .= z₃
    @.. u = tmp + γ*z₃
    f2( k3, u,p,t+c3*dt); k3 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a41*z₁ + a43*z₃ + ea41*k1 + ea43*k3
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a41*z₁[i] + a43*z₃[i] + ea41*k1[i] + ea43*k3[i]
    end
  else
    @.. z₄ = α41*z₁ + α42*z₂
    @.. tmp = uprev + a41*z₁ + a43*z₃
  end
  nlsolver.z = z₄

  nlsolver.c = c4
  z₄ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 5

  if typeof(integrator.f) <: SplitFunction
    z₅ .= z₂
    @.. u = tmp + γ*z₄
    f2( k4, u,p,t+c4*dt); k4 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a51*z₁ + a53*z₃ + a54*z₄ + ea51*k1 + ea53*k3 + ea54*k4
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a51*z₁[i] + a53*z₃[i] + a54*z₄[i] + ea51*k1[i] + ea53*k3[i] + ea54*k4[i]
    end
  else
    @.. z₅ = α51*z₁ + α52*z₂
    @.. tmp = uprev + a51*z₁ + a53*z₃ + a54*z₄
  end
  nlsolver.z = z₅

  nlsolver.c = c5
  z₅ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 6

  if typeof(integrator.f) <: SplitFunction
    z₆ .= z₃
    @.. u = tmp + γ*z₅
    f2( k5, u,p,t+c5*dt); k5 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅ + ea61*k1 + ea63*k3 + ea64*k4 + ea65*k5
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i] + ea61*k1[i] + ea63*k3[i] + ea64*k4[i] + ea65*k5[i]
    end
  else
    @.. z₆ = α61*z₁ + α62*z₂
    #@.. tmp = uprev + a61*z₁ + a63*z₃ + a64*z₄ + a65*z₅
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a61*z₁[i] + a63*z₃[i] + a64*z₄[i] + a65*z₅[i]
    end
  end
  nlsolver.z = z₆

  nlsolver.c = c6
  z₆ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 7

  if typeof(integrator.f) <: SplitFunction
    z₇ .= z₂
    @.. u = tmp + γ*z₆
    f2( k6, u,p,t+c6*dt); k6 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a71*z₁ +  a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆ + ea71*k1 + ea73*k3 + ea74*k4 + ea75*k5 + ea76*k6
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a71*z₁[i] +  a73*z₃[i] + a74*z₄[i] + a75*z₅[i] + a76*z₆[i] + ea71*k1[i] + ea73*k3[i] + ea74*k4[i] + ea75*k5[i] + ea76*k6[i]
    end
  else
    #@.. z₇ = α71*z₁ + α72*z₂ + α73*z₃ + α74*z₄ + α75*z₅
    @tight_loop_macros for i in eachindex(u)
      @inbounds z₇[i] = α71*z₁[i] + α72*z₂[i] + α73*z₃[i] + α74*z₄[i] + α75*z₅[i]
    end
    #@.. tmp = uprev + a71*z₁ + a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆
    @tight_loop_macros for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a71*z₁[i] + a73*z₃[i] + a74*z₄[i] + a75*z₅[i] + a76*z₆[i]
    end
  end
  nlsolver.z = z₇

  nlsolver.c = c7
  z₇ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 8

  if typeof(integrator.f) <: SplitFunction
    z₈ .= z₅
    @.. u = tmp + γ*z₇
    f2( k7, u,p,t+c7*dt); k7 .*= dt
    integrator.destats.nf2 += 1
    #@.. tmp = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇ + ea81*k1 + ea83*k3 + ea84*k4 + ea85*k5 + ea86*k6 + ea87*k7
    for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a81*z₁[i] + a84*z₄[i] + a85*z₅[i] + a86*z₆[i] + a87*z₇[i] + ea81*k1[i] + ea83*k3[i] + ea84*k4[i] + ea85*k5[i] + ea86*k6[i] + ea87*k7[i]
    end
  else
    #@.. z₈ = α81*z₁ + α82*z₂ + α83*z₃ + α84*z₄ + α85*z₅
    @tight_loop_macros for i in eachindex(u)
      @inbounds z₈[i] = α81*z₁[i] + α82*z₂[i] + α83*z₃[i] + α84*z₄[i] + α85*z₅[i]
    end
    #@.. tmp = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
    @tight_loop_macros for i in eachindex(u)
      @inbounds tmp[i] = uprev[i] + a81*z₁[i] + a84*z₄[i] + a85*z₅[i] + a86*z₆[i] + a87*z₇[i]
    end
  end
  nlsolver.z = z₈

  nlsolver.c = 1
  z₈ = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₈
  if typeof(integrator.f) <: SplitFunction
    f2( k8, u,p,t+dt); k8 .*= dt
    integrator.destats.nf += 1
    # @.. u = uprev + a81*z₁ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇ + γ*z₈ + eb1*k1 + eb4*k4 + eb5*k5 + eb6*k6 + eb7*k7 + eb8*k8
    for i in eachindex(u)
      @inbounds u[i] = uprev[i] + a81*z₁[i] + a84*z₄[i] + a85*z₅[i] + a86*z₆[i] + a87*z₇[i] + γ*z₈[i] + eb1*k1[i] + eb4*k4[i] + eb5*k5[i] + eb6*k6[i] + eb7*k7[i] + eb8*k8[i]
    end
  end

  ################################### Finalize

  if integrator.opts.adaptive

    if typeof(integrator.f) <: SplitFunction
      #@.. dz =  btilde1*z₁ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈ + ebtilde1*k1 + ebtilde4*k4 + ebtilde5*k5 + ebtilde6*k6 + ebtilde7*k7 + ebtilde8*k8
      for i in eachindex(u)
        dz[i] =  btilde1*z₁[i] + btilde4*z₄[i] + btilde5*z₅[i] + btilde6*z₆[i] + btilde7*z₇[i] + btilde8*z₈[i] + ebtilde1*k1[i] + ebtilde4*k4[i] + ebtilde5*k5[i] + ebtilde6*k6[i] + ebtilde7*k7[i] + ebtilde8*k8[i]
      end
    else
      # @.. dz = btilde1*z₁ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈
      @tight_loop_macros for i in eachindex(u)
        @inbounds dz[i] = btilde1*z₁[i] + btilde4*z₄[i] + btilde5*z₅[i] + btilde6*z₆[i] + btilde7*z₇[i] + btilde8*z₈[i]
      end

    end

    if isnewton(nlsolver) && alg.smooth_est # From Shampine
      integrator.destats.nsolve += 1
      cache.linsolve(vec(tmp),get_W(nlsolver),vec(dz),false)
    else
      tmp .= dz
    end
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  if typeof(integrator.f) <: SplitFunction
    integrator.f(integrator.fsallast,u,p,t+dt)
  else
    @.. integrator.fsallast = z₈/dt
  end
end


text:
using OptimBase
using Base.Test
import Base.summary

immutable FakeOptimizer <: Optimizer
end

@testset "summary" begin
    a_fake_optimizer = FakeOptimizer()
    @test_throws ErrorException summary(a_fake_optimizer)

    Base.summary(::FakeOptimizer) = "Fake Optimizer"
# add back when moving to v0.6
#    @test_nowarn summary(a_fake_optimizer)
end

@testset "optimization results" begin
    Base.summary(::FakeOptimizer) = "Fake Optimizer"
    uor = UnivariateOptimizationResults(FakeOptimizer(),
                                  0.1,
                                  0.2,
                                  0.3,
                                  0.0,
                                  10,
                                  false,
                                  true,
                                  0.1,
                                  0.2,
                                  OptimizationTrace{FakeOptimizer}(),
                                  32)
    @test OptimBase.summary(uor) == "Fake Optimizer"
    @test OptimBase.rel_tol(uor) == 0.1
    @test OptimBase.abs_tol(uor) == 0.2
    @test OptimBase.f_calls(uor) == 32
    #=
    type MultivariateOptimizationResults{O<:Optimizer,T,N,M} <: OptimizationResults
        method::O
        initial_x::Array{T,N}
        minimizer::Array{T,N}
        minimum::T
        iterations::Int
        iteration_converged::Bool
        x_converged::Bool
        x_tol::Float64
        x_residual::Float64
        f_converged::Bool
        f_tol::Float64
        f_residual::Float64
        g_converged::Bool
        g_tol::Float64
        g_residual::Float64
        f_increased::Bool
        trace::OptimizationTrace{M}
        f_calls::Int
        g_calls::Int
        h_calls::Int
    end
    =#
end


text:
# coding: utf-8

# In[1]:

#counting sort is kinda different from other sorts in this repo
#it works the best when there is large quantity of duplicates
#when the maximum value is larger than the length of list
#it is not economical to use counting sort
#counting sort is counting how many times each element has repeatedly showed up
#we create a list of the size of maximum value
#each index of the list represents the value of elements
#each value of the list represents the frequency of elements
#in the end we concatenate each element and its frequency
#simplest sort among all

# In[2]:

function counting_sort(arr)

    maxval=maximum(arr)

    #first we get the maximum value
    #if maximum value is larger than the length of the list
    #we would create a list with a huge space complexity
    #its not worth it
    #we can use a better sort
    if maxval>length(arr)
        printstyled("choose other sorting techniques",color=:red)
        return
    end

    #we create a list to store the frequency of each element
    #note that we use maximum value +1
    #as zero is also taken into consideration
    freq=zeros(maxval+1)

    #this is the calculation of the frequency
    #use i+1 becuz julia index starts from 1
    for i in arr
        freq[i+1]+=1
    end

    #lets clear the original list
    #we add each element and its frequency back to the list
    output=[]

    for i in 1:length(freq)

        #using i-1 prior to i+1 in the lines above
        #use Int32 to ensure freq[i] is int32
        temp=fill(i-1,Int32(freq[i]))
        append!(output,temp)

    end

    return output

end

# In[3]:

#the time complexity for counting sort is o(n+k)
#k is determined by the maximum value
#its easy to see that counting sort easily beats o(n^2)
#however its space complexity is really large

# In[4]:

for _ in 1:100

    #usually we form a non-duplicate list
    #for counting sort, we need massive amount of duplicates
    test_arr=rand(0:100,1000)

    if !(counting_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: counting sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: counting sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: counting sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
struct StringOptions
        nmax::Int
        tol::Float64
        verbose::Bool
        save_trajectory::Bool
        print_iters::Int
end

"""
    StringOptions(;nmax = 10^3, tol = 1e-6, verbose = false, save_evolution = true)
"""
function StringOptions(;nmax = 10^3, tol = 1e-6, verbose = false, save_trajectory = true, print_iters = 10^3)
    return StringOptions(nmax, tol, verbose, save_trajectory, print_iters)
end

struct SaddleOptions
        nmax::Int
        tol::Float64
        verbose::Bool
        save_trajectory::Bool
        print_iters::Int
end
"""
    SaddleOptions(;nmax = 10^3, tol = 1e-6, verbose = false, save_trajectory = true, print_iters = 10^3)
"""
function SaddleOptions(;nmax = 10^3, tol = 1e-6, verbose = false, save_trajectory = true, print_iters = 10^3)
    return SaddleOptions(nmax, tol, verbose, save_trajectory, print_iters)
end

"""
`periodic_dist`: Compute the periodic distance between values

### Fields
* `u` - First point on the periodic interval
* `v` - Second point on the periodic interval
### Optional Fields
* `Π = 1.0` - Period of the interval, [0, Π)
"""
function periodic_dist(u,v; Π = 1.0)
    return min(abs(u-v), Π- abs(u-v))
end

"""
`linear_string` - Construct the linear interpolant string between x₀ and x₁ with
a given number of images (inclusive).

### Fields
* `x₀` - Initial image on the string
* `x₁` - Final image on the string
* `n_images` - Total number of images on the string, including x₀ and x₁
"""
function linear_string(x₀, x₁, n_images)

    X = typeof(x₀)[];
    x = similar(x₀);

    for j in 0:n_images-1
        t = Float64(j/(n_images-1));
        @. x = (1. - t) * x₀ + t * x₁;
        push!(X, copy(x));
    end
    return X

end

"""
`upwind_tangent` - Estimate the tangent vector from a string's interior using
an upwinding method

### Fields
* `U` - Array of three images, in sequence, along the string
* `V` - Potential energy defining the landscape
"""

function upwind_tangent(U, V)
    V_vals = V.(U);
    if (V_vals[1]<V_vals[2] && V_vals[2]< V_vals[3])
        τ = (U[3] .- U[2])/norm(U[3] .- U[2]);
    elseif (V_vals[1]>V_vals[2] && V_vals[2]> V_vals[3])
        τ = (U[2] .- U[1])/norm(U[2] .- U[1]);
    else
        τ = (U[3] .- U[1])/norm(U[3] .- U[1]);
    end

    return τ
end


text:
##### Beginning of file

module Utils # Begin submodule Snapshots.Utils

__precompile__(true)

include("git.jl")

include("misc.jl")

end # End submodule Snapshots.Utils

##### End of file


text:
module ClTools

include("visualization.jl")
include("config.jl")

end # module


text:
"""
Converts Fahrenheit or Centigrade to absolute temperature
temperature = temperature to Convert
temperatureUnit = (c)entigrade, (f)ahrenheit, or (a)bsolute
"""
function absolute_temperature(temperature,temperatureUnit)
    if lowercase(temperatureUnit) == "c"
        return round(temperature + 273)
    elseif lowercase(temperatureUnit) == "f"
        return round((temperature + 459) / 1.8)
    elseif lowercase(temperatureUnit) == "a"
        return temperature
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+b*coth(c+d*x)^2)^5, x) == :(x*(a+b)^5+-1//9*b^5*d^-1*coth(c+d*x)^9+-1*b*d^-1*(b^4+5*a^4+5*a*b^3+10*b*a^3+10*a^2*b^2)*coth(c+d*x)+-1//3*b^2*d^-1*coth(c+d*x)^3*(b^3+10*a^3+5*a*b^2+10*b*a^2)+-1//5*b^3*d^-1*coth(c+d*x)^5*(b^2+10*a^2+5*a*b)+-1//7*b^4*d^-1*coth(c+d*x)^7*(b+5a))
@test integrate((a+b*coth(c+d*x)^2)^4, x) == :(x*(a+b)^4+-1//7*b^4*d^-1*coth(c+d*x)^7+-1//3*b^2*d^-1*coth(c+d*x)^3*(b^2+6*a^2+4*a*b)+-1//5*b^3*d^-1*coth(c+d*x)^5*(b+4a)+-1*b*d^-1*(b+2a)*(b^2+2*a^2+2*a*b)*coth(c+d*x))
@test integrate((a+b*coth(c+d*x)^2)^3, x) == :(x*(a+b)^3+-1//5*b^3*d^-1*coth(c+d*x)^5+-1*b*d^-1*(b^2+3*a^2+3*a*b)*coth(c+d*x)+-1//3*b^2*d^-1*coth(c+d*x)^3*(b+3a))
@test integrate((a+b*coth(c+d*x)^2)^2, x) == :(x*(a+b)^2+-1//3*b^2*d^-1*coth(c+d*x)^3+-1*b*d^-1*(b+2a)*coth(c+d*x))
@test integrate((a+b*coth(c+d*x)^2)^-1, x) == :(x*(a+b)^-1+-1*a^-1//2*b^(1/2)*d^-1*(a+b)^-1*arctan(a^(1/2)*b^-1//2*tanh(c+d*x)))
@test integrate((a+b*coth(c+d*x)^2)^-2, x) == :(x*(a+b)^-2+(1/2)*b*a^-1*d^-1*(a+b)^-1*(a+b*coth(c+d*x)^2)^-1*coth(c+d*x)+-1//2*a^-3//2*b^(1/2)*d^-1*(a+b)^-2*(b+3a)*arctan(a^(1/2)*b^-1//2*tanh(c+d*x)))
@test integrate((a+b*coth(c+d*x)^2)^-3, x) == :(x*(a+b)^-3+-1//8*a^-5//2*b^(1/2)*d^-1*(a+b)^-3*(3*b^2+15*a^2+10*a*b)*arctan(a^(1/2)*b^-1//2*tanh(c+d*x))+1//4*b*a^-1*d^-1*(a+b)^-1*(a+b*coth(c+d*x)^2)^-2*coth(c+d*x)+1//8*b*a^-2*d^-1*(a+b)^-2*(a+b*coth(c+d*x)^2)^-1*(3b+7a)*coth(c+d*x))
@test integrate((a+b*coth(c+d*x)^2)^-4, x) == :(x*(a+b)^-4+-1//16*a^-7//2*b^(1/2)*d^-1*(a+b)^-4*(5*b^3+35*a^3+21*a*b^2+35*b*a^2)*arctan(a^(1/2)*b^-1//2*tanh(c+d*x))+1//6*b*a^-1*d^-1*(a+b)^-1*(a+b*coth(c+d*x)^2)^-3*coth(c+d*x)+1//16*b*a^-3*d^-1*(a+b)^-3*(a+b*coth(c+d*x)^2)^-1*(5*b^2+19*a^2+16*a*b)*coth(c+d*x)+1//24*b*a^-2*d^-1*(a+b)^-2*(a+b*coth(c+d*x)^2)^-2*(5b+11a)*coth(c+d*x))
@test integrate((1+-2*coth(x)^2)^-1, x) == :(-1x+2^(1/2)*arctanh((1/2)*2^(1/2)*tanh(x)))
@test integrate((1+-1*coth(x)^2)^(1/2), x) == :(arcsin(coth(x)))
@test integrate((-1+coth(x)^2)^(1/2), x) == :(-1*arctanh((csch(x)^2)^-1//2*coth(x)))
@test integrate((1+-1*coth(x)^2)^3//2, x) == :((1/2)*arcsin(coth(x))+(1/2)*(-1*csch(x)^2)^(1/2)*coth(x))
@test integrate((-1+coth(x)^2)^3//2, x) == :((1/2)*arctanh((csch(x)^2)^-1//2*coth(x))+-1//2*(csch(x)^2)^(1/2)*coth(x))
@test integrate((1+-1*coth(x)^2)^-1//2, x) == :((-1*csch(x)^2)^-1//2*coth(x))
@test integrate((-1+coth(x)^2)^-1//2, x) == :((csch(x)^2)^-1//2*coth(x))
@test integrate((a+b*coth(x)^2)^(1/2)*coth(x)^3, x) == :(-1*(a+b*coth(x)^2)^(1/2)+(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1//3*b^-1*(a+b*coth(x)^2)^3//2)
@test integrate((a+b*coth(x)^2)^(1/2)*coth(x)^2, x) == :((a+b)^(1/2)*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1//2*(a+b*coth(x)^2)^(1/2)*coth(x)+-1//2*b^-1//2*(a+2b)*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^(1/2)*coth(x), x) == :(-1*(a+b*coth(x)^2)^(1/2)+(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2)))
@test integrate((a+b*coth(x)^2)^(1/2), x) == :((a+b)^(1/2)*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*b^(1/2)*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^(1/2)*tanh(x), x) == :((a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*a^(1/2)*arctanh(a^-1//2*(a+b*coth(x)^2)^(1/2)))
@test integrate((a+b*coth(x)^2)^(1/2)*tanh(x)^2, x) == :((a+b)^(1/2)*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*(a+b*coth(x)^2)^(1/2)*tanh(x))
@test integrate((a+b*coth(x)^2)^3//2*coth(x)^3, x) == :(-1//3*(a+b*coth(x)^2)^3//2+(a+b)^3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*(a+b*coth(x)^2)^(1/2)*(a+b)+-1//5*b^-1*(a+b*coth(x)^2)^5//2)
@test integrate((a+b*coth(x)^2)^3//2*coth(x)^2, x) == :((a+b)^3//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*(a+b*coth(x)^2)^(1/2)*((1/2)*b+5//8*a)*coth(x)+-1//4*b*(a+b*coth(x)^2)^(1/2)*coth(x)^3+-1//8*b^-1//2*(3*a^2+8*b^2+12*a*b)*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^3//2*coth(x), x) == :(-1//3*(a+b*coth(x)^2)^3//2+(a+b)^3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*(a+b*coth(x)^2)^(1/2)*(a+b))
@test integrate((a+b*coth(x)^2)^3//2, x) == :((a+b)^3//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1//2*b*(a+b*coth(x)^2)^(1/2)*coth(x)+-1//2*b^(1/2)*(2b+3a)*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^3//2*tanh(x), x) == :((a+b)^3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*b*(a+b*coth(x)^2)^(1/2)+-1*a^3//2*arctanh(a^-1//2*(a+b*coth(x)^2)^(1/2)))
@test integrate((a+b*coth(x)^2)^3//2*tanh(x)^2, x) == :((a+b)^3//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*b^3//2*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*a*(a+b*coth(x)^2)^(1/2)*tanh(x))
@test integrate((1+coth(x)^2)^(1/2), x) == :(-1*arcsinh(coth(x))+2^(1/2)*arctanh(2^(1/2)*(1+coth(x)^2)^-1//2*coth(x)))
@test integrate((-1+-1*coth(x)^2)^(1/2), x) == :(-1*2^(1/2)*arctan(2^(1/2)*(-1+-1*coth(x)^2)^-1//2*coth(x))+arctan((-1+-1*coth(x)^2)^-1//2*coth(x)))
@test integrate((1+coth(x)^2)^3//2, x) == :(-5//2*arcsinh(coth(x))+2*2^(1/2)*arctanh(2^(1/2)*(1+coth(x)^2)^-1//2*coth(x))+-1//2*(1+coth(x)^2)^(1/2)*coth(x))
@test integrate((-1+-1*coth(x)^2)^3//2, x) == :(-5//2*arctan((-1+-1*coth(x)^2)^-1//2*coth(x))+(1/2)*(-1+-1*coth(x)^2)^(1/2)*coth(x)+2*2^(1/2)*arctan(2^(1/2)*(-1+-1*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^-1//2*coth(x)^3, x) == :((a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*b^-1*(a+b*coth(x)^2)^(1/2))
@test integrate((a+b*coth(x)^2)^-1//2*coth(x)^2, x) == :((a+b)^-1//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*b^-1//2*arctanh(b^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^-1//2*coth(x), x) == :((a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2)))
@test integrate((a+b*coth(x)^2)^-1//2, x) == :((a+b)^-1//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x)))
@test integrate((a+b*coth(x)^2)^-1//2*tanh(x), x) == :((a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*a^-1//2*arctanh(a^-1//2*(a+b*coth(x)^2)^(1/2)))
@test integrate((a+b*coth(x)^2)^-1//2*tanh(x)^2, x) == :((a+b)^-1//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*a^-1*(a+b*coth(x)^2)^(1/2)*tanh(x))
@test integrate((a+b*coth(x)^2)^-3//2*coth(x)^3, x) == :((a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+a*b^-1*(a+b)^-1*(a+b*coth(x)^2)^-1//2)
@test integrate((a+b*coth(x)^2)^-3//2*coth(x)^2, x) == :((a+b)^-3//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*(a+b)^-1*(a+b*coth(x)^2)^-1//2*coth(x))
@test integrate((a+b*coth(x)^2)^-3//2*coth(x), x) == :((a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*(a+b)^-1*(a+b*coth(x)^2)^-1//2)
@test integrate((a+b*coth(x)^2)^-3//2*tanh(x), x) == :((a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*a^-3//2*arctanh(a^-1//2*(a+b*coth(x)^2)^(1/2))+b*a^-1*(a+b)^-1*(a+b*coth(x)^2)^-1//2)
@test integrate((a+b*coth(x)^2)^-3//2*tanh(x)^2, x) == :((a+b)^-3//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+b*a^-1*(a+b)^-1*(a+b*coth(x)^2)^-1//2*tanh(x)+-1*a^-2*(a+b)^-1*(a+b*coth(x)^2)^(1/2)*(a+2b)*tanh(x))
@test integrate((a+b*coth(x)^2)^-5//2*coth(x)^3, x) == :((a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*(a+b)^-2*(a+b*coth(x)^2)^-1//2+1//3*a*b^-1*(a+b)^-1*(a+b*coth(x)^2)^-3//2)
@test integrate((a+b*coth(x)^2)^-5//2*coth(x)^2, x) == :((a+b)^-5//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+-1*(a+b*coth(x)^2)^-3//2*(3a+3b)^-1*coth(x)+-1//3*a^-1*(a+b)^-2*(a+b*coth(x)^2)^-1//2*(-1b+2a)*coth(x))
@test integrate((a+b*coth(x)^2)^-5//2*coth(x), x) == :((a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*(a+b)^-2*(a+b*coth(x)^2)^-1//2+-1*(a+b*coth(x)^2)^-3//2*(3a+3b)^-1)
@test integrate((a+b*coth(x)^2)^-5//2*tanh(x), x) == :((a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*coth(x)^2)^(1/2))+-1*a^-5//2*arctanh(a^-1//2*(a+b*coth(x)^2)^(1/2))+1//3*b*a^-1*(a+b)^-1*(a+b*coth(x)^2)^-3//2+b*a^-2*(a+b)^-2*(a+b*coth(x)^2)^-1//2*(b+2a))
@test integrate((a+b*coth(x)^2)^-5//2*tanh(x)^2, x) == :((a+b)^-5//2*arctanh((a+b)^(1/2)*(a+b*coth(x)^2)^-1//2*coth(x))+1//3*b*a^-1*(a+b)^-1*(a+b*coth(x)^2)^-3//2*tanh(x)+-1//3*a^-3*(a+b)^-2*(a+b*coth(x)^2)^(1/2)*(a+4b)*(2b+3a)*tanh(x)+1//3*b*a^-2*(a+b)^-2*(a+b*coth(x)^2)^-1//2*(4b+7a)*tanh(x))
@test integrate((1+coth(x)^2)^-1//2, x) == :((1/2)*2^(1/2)*arctanh(2^(1/2)*(1+coth(x)^2)^-1//2*coth(x)))
@test integrate((-1+-1*coth(x)^2)^-1//2, x) == :((1/2)*2^(1/2)*arctan(2^(1/2)*(-1+-1*coth(x)^2)^-1//2*coth(x)))
@test integrate((1+coth(x)^3)^-1, x) == :((1/2)*x+-1*(6+6*coth(x))^-1+-2//9*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2*coth(x))))
@test integrate((a+b*coth(x)^4)^(1/2)*coth(x), x) == :(-1//2*(a+b*coth(x)^4)^(1/2)+(1/2)*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*coth(x)^4)^-1//2*(a+b*coth(x)^2))+-1//2*b^(1/2)*arctanh(b^(1/2)*(a+b*coth(x)^4)^-1//2*coth(x)^2))
@test integrate((a+b*coth(x)^4)^-1//2*coth(x), x) == :((1/2)*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*coth(x)^4)^-1//2*(a+b*coth(x)^2)))
@test integrate((a+b*coth(x)^4)^-3//2*coth(x), x) == :((1/2)*(a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*coth(x)^4)^-1//2*(a+b*coth(x)^2))+-1//2*a^-1*(a+b)^-1*(a+b*coth(x)^4)^-1//2*(a+-1*b*coth(x)^2))

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/6.4.7 (d hyper)^m (a+b (c coth)^n)^p.jl


text:
using GLMakie, LinearAlgebra, ForwardDiff
x = -2:0.25:2
y = -1:0.25:1

zz(x,y) = x*exp(-x^2 - y^2)
g(x,y, z) = [-x*exp(-x^2 - y^2) + z]
J(xx, yy, zz) = ForwardDiff.jacobian(x -> g(x[1], x[2], x[3]), [xx,yy,zz])

ps = [Point3f(x, y, zz(x,y)) for x in x for y in y]
normals = [Point3f(J(x, y, zz(x,y))...) for x in x for y in y]
lengths = norm.(normals)

set_theme!(theme_dark())
fig, ax, = surface(x, y, zz.(x, y'); transparency = true,
    axis = (;type = Axis3, aspect = :data, perspectiveness = 0.5,
    elevation = 0.34, azimuth = 2.11),
    figure = (;resolution = (1200,600)))
wireframe!(ax, x, y,  zz.(x, y'); color = :white, transparency = true,
    linewidth = 0.25)
arrows!(ax, ps, normals./7; color=lengths, colormap = Reverse(:thermal),
    linewidth = 0.02, arrowsize = Vec3f(0.08, 0.08, 0.08))
hidedecorations!(ax; grid = false)
save("arrows_normals.png", fig)
fig


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

# Current version
version = "0.6.0"

# Get current operating system
osqp_platform =
if Sys.islinux()
    "linux"
elseif Sys.isapple()
    "mac"
elseif Sys.iswindows()
    "windows"
else
    error("Platform not supported!")
end

subfolder="osqp-$version-$osqp_platform$(Sys.WORD_SIZE)"
products = [
        LibraryProduct(joinpath(prefix, subfolder, "lib"), ["libosqp"], :osqp),
    ]

# Provide binaries for each operating system

function get_hash_code(platform)
    archive_name = "osqp-$version-$platform.tar.gz"
    readlines(download("https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.sha256"))[1]
end

function get_url(platform)
    archive_name = "osqp-$version-$platform"
    return "https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.tar.gz"
end

download_info = Dict(
    Linux(:x86_64) => (get_url("linux64"), get_hash_code("linux64")),
    MacOS(:x86_64) => (get_url("mac64"), get_hash_code("mac64")),
    Windows(:x86_64) => (get_url("windows64"), get_hash_code("windows64")),
    Windows(:i686) => (get_url("windows32"), get_hash_code("windows32")),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]

    # Check if this build.jl is providing new versions of the binaries, and
    # if so, ovewrite the current binaries even if they were installed by the user
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key_abi())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using GeodesicRendering
using AccretionGeometry
using CarterBoyerLindquist
using AccretionFormulae
using StaticArrays
using Printf

using Plots
gr()

function temperature(m, gp, max_time; kwargs...)
    g = AccretionFormulae.redshift(m, gp, max_time; kwargs...)

    M = m.M
    a_star = m.a
    temperature = AccretionFormulae.observed_temperature(gp.u[2], a_star, M, g)

end

function radius(m, gp, max_time; kwargs...)
    gp.u[2]
end

function temperature_render(;
    mass = 1,
    spin = 0.998,
    obs_angle = 85.0,
    disc_angle = 90.0,
    tolerance = 1e-8,
    dtmax = 1000.0,
    size_multiplier::Int64 = 1,
    resolution = 400,
    fov = 3.0,
    η = 0.1,
    η_phys = 0.1,
    edd_ratio = 0.1,
    edd_ratio_phys = 0.1,
)

    m = CarterMethodBL(M = 1.0, a = spin)
    M = m.M

    # observer position
    u = [0.0, 1000.0, deg2rad(obs_angle), 0.0]
    R_isco = AccretionFormulae.r_isco(m.a, m.M)

    # disc
    d = GeometricThinDisc(R_isco + 1, 50.0, deg2rad(disc_angle))

    # cache the render
    cache = @time prerendergeodesics(
        m,
        u,
        2000.0,
        d,
        fov_factor = fov * size_multiplier,
        abstol = tolerance,
        reltol = tolerance,
        image_width = 350 * size_multiplier,
        image_height = 350 * size_multiplier,
        dtmax = dtmax,
    )

    # create and compose the PointFunctions
    temperature_vf = (
        PointFunction(temperature) ∘
        FilterPointFunction((m, gp, max_time; kwargs...) -> gp.u[2] > R_isco, NaN) ∘
        ConstPointFunctions.filter_early_term
    )

    radius_vf = (
        PointFunction(radius) ∘
        FilterPointFunction((m, gp, max_time; kwargs...) -> gp.u[2] > R_isco, NaN) ∘
        ConstPointFunctions.filter_early_term
    )

    radius_img = GeodesicRendering.apply(radius_vf, cache)
    temperature_img = GeodesicRendering.apply(temperature_vf, cache)

    # correcting for physical mass
    M_phys = mass * 1.99e30
    r_isco_phys = AccretionFormulae.r_isco(M_phys, 0.998)
    r_g_phys = M_phys

    numerators =
        AccretionFormulae.mass_scale_fraction.(
            M_phys,
            η_phys,
            edd_ratio_phys,
            r_isco_phys,
            r_g_phys,
            radius_img * M_phys,
        )
    denominators =
        AccretionFormulae.mass_scale_fraction.(M, η, edd_ratio, R_isco, M, radius_img)

    fractions = numerators ./ denominators
    temperature_img .*= fractions

    # # autoscale
    # scale = maximum(filter(!isnan,temperature_img))
    # scale = floor(log(10, scale))
    # scale = 10^scale

    # scaling image
    scale = 1e7
    scalestr = @sprintf "%.E" scale
    exponent =  Int32(floor(log(10, scale)))
    new_img = reverse(temperature_img, dims = 1)
    new_img ./= scale

    hmap = heatmap(
        new_img,
        aspect_ratio = 1.0,
        size = (resolution * 3 / 2, resolution),
        clim = (0, 3),
    )
    # contour(new_img, aspect_ratio=1.0, size=(resolution*3/2, resolution), clim=(0,3))
    # title!("Temperature Scale = $scalestr, Mass = $mass M_☼, Obs Angle = $obs_angle")
    # title = "Temperature Scale = $scalestr K, Mass = $mass M_☼, Obs Angle = $obs_angle"
    title = "Temperature Scale = \$10^{$exponent}\$ K, Mass = $mass \$\\mathrm{M}_{\\odot}\$, Obs Angle = $obs_angle\$^{\\circ}\$"
    # "")
    return hmap, cache, title, new_img
end

# hmap, cache, title = temperature_render(obs_angle = 62.5, mass = 10, resolution=1080)
# # hmap, cache, title = temperature_render(
# #                                         mass=10,
# #                                         spin=0.998,
# #                                         obs_angle=0.01,
# #                                         tolerance=1e-12,
# #                                         size_multiplier=6,
# #                                         dtmax=1,
# #                                         resolution=2000
# #                                         )
# title!(title)
# display(hmap)