text:
module aoc_20

export read_input, edges, Tile
export make_edges_dict, get_out_edge, get_outer_tiles

struct Tile
    number::Int
    image::BitMatrix
    edges::Vector{Int}
end

const tile_r = r"Tile (\d+):"

function read_input(filename)
    lines = readlines(filename)
    tile_number = -1
    buffer = Vector{BitVector}()
    tiles = Dict{Int,Tile}()
    for line in lines
        if isempty(line)
            image = permutedims( hcat(buffer...) )
            tiles[ tile_number ] = Tile( tile_number, image, edges(image))
            tile_number = -1
            empty!(buffer)
            continue
        end
        if tile_number == -1
            m = match(tile_r, line)
            if isnothing(m) || length(m.captures) < 1
                error("Could not read tile from \"$line\"")
            else
                tile_number = parse(Int, m.captures[1])
                continue
            end
        else
            push!(buffer, Vector{Char}(line) .== '#')
        end
    end
    if tile_number != -1
        image = permutedims( hcat(buffer...) )
        tiles[ tile_number ] = Tile( tile_number, image, edges(image))
    end
    tiles
end

export edge_order
const edge_order = [:top, :bottom, :left, :right, :top_flipped, :bottom_flipped, :left_flipped, :right_flipped]

function edges(image::BitMatrix)
    a = image[1,:][:] # top
    b = image[end,:][:] # bottom
    c = image[:,1][:] # left
    d = image[:,end][:] # right
    e = [a, b, c, d]
    e = append!(e, reverse.(e)) # flipped
    edge2int.(e)
end

function make_edges_dict( tiles::Dict{Int,Tile} )
    edge2tile = Dict{Int,Vector{Int}}()
    for (k,v) in tiles, e in v.edges
        v = get( edge2tile, e, Vector{Int}() )
        edge2tile[ e ] = push!(v, k)
    end
    edge2tile
end

const bits = 0:9
edge2int(b::BitVector) = sum(b .<< bits)

function place_tiles( tiles::Dict{Int,Tile} )
    edge2tile = make_edges_dict( tiles )
    outer_edge = get_outer_edge( edge2tile )
    outer_tiles = get_outer_tiles( outer_edge )
    placed_tiles = Int[]
    n_tiles = length(tiles)
    tiles_per_edge = round(Int, sqrt(n_tiles) )
    image = zeros(Int, tiles_per_edge, tiles_per_edge)
    # Wait, do I need to do this? I just need the corners
    # Yes for part 2

    corners = get_corners( tiles )

    # Start with a corner
    image[1] = first( corners )
    get_neighbors(t, tiles=tiles, edge2tile=edge2tile) = unique(first.(filter.(x->x!=t,filter(x->length(x)==2,[ edge2tile[ e ] for e in tiles[ t ].edges ]))))

    # Build out neighbors
    neighbors = get_neighbors( image[1] )
    image[1,2] = neighbors[1]
    image[2,1] = neighbors[2]
    push!(placed_tiles, image[1])
    append!(placed_tiles, neighbors)

    indices = CartesianIndices(image)
    #last_ind = indices[2,1]
    outer_edge_itr = vcat( indices[3:end,1], indices[end,2:end], indices[end-1:-1:1,end], indices[1, end-1:-1:3] )
    # outer_tiles = filter!(t->!in(t[1],placed_tiles), filter(t->in(t[1],outer_tiles), tiles) )

    is_unplaced_outer_tile(t::Int) = in(t, outer_tiles) & !in(t, placed_tiles)

    last_tile = image[2,1]
    for ind in outer_edge_itr
        #last_tile = image[ last_ind ]
        neighbor = first( filter!(is_unplaced_outer_tile, get_neighbors(last_tile)) )

        image[ ind ] = neighbor
        push!(placed_tiles, neighbor)

        last_tile = neighbor

    end

    is_unplaced(t) = !in(t, placed_tiles)
    itr = indices[2:end-1,2:end-1][:]
    for ind in itr
        above = image[ ind - CartesianIndex(1,0) ]
        left = image[ ind - CartesianIndex(0,1) ]
        neighbor = first( filter(is_unplaced, intersect( get_neighbors(above), get_neighbors(left) ) ) )

        image[ ind ] = neighbor
        push!(placed_tiles, neighbor)
    end

    image

end

function assemble_image( tiles::Dict{Int,Tile}, image = place_tiles( tiles ))
    tile_image = [tiles[ i ] for i in image]
    n_tiles = length(tiles)
    tiles_per_edge = round(Int, sqrt(n_tiles) )

    actual_image = falses( size( first(tiles).second.image ) .*  tiles_per_edge )

    tn, tm = size(tile_image[1].image)

    actual_image[1:tn,1:tm] = orient_tile_down_right( tile_image[1,1], tile_image[2,1], tile_image[1,2] )

    indices = CartesianIndices(image)

    prev_edge = edge2int(actual_image[tn, 1:tm])

    for ind in indices[2:end,1]
        current = tile_image[ ind ]

        eo = edge_order[ findfirst(current.edges .== prev_edge ) ]
        if eo == :left_flipped
            I = rotr90(current.image)
        elseif eo == :bottom_flipped
            I = rot180(current.image)
        elseif eo == :top
            I = current.image
        elseif eo == :right
            I = rotl90(current.image)
        elseif eo == :left
            I = reverse( rotr90(current.image), dims = 2)
        elseif eo == :bottom
            I = reverse( rot180(current.image), dims = 2)
        elseif eo == :top_flipped
            I = reverse( current.image, dims = 2)
        elseif eo == :right_flipped
            I = reverse( rotl90(current.image), dims = 2 )
        else
            error("Orientation")
        end
        #print(eo)
        actual_image[tn*(ind[1]-1) .+ (1:tn), tm*(ind[2]-1) .+ (1:tm)] = I
        prev_edge = edge2int(I[end,:])
    end

    for row = 1:tiles_per_edge
        prev_edge = edge2int( actual_image[(row-1)*tn .+ (1:tn), tm])
    for ind in indices[row,2:end]
        current = tile_image[ ind ]

        eo = edge_order[ findfirst(current.edges .== prev_edge ) ]
        if eo == :bottom
            I = rotr90(current.image)
        elseif eo == :right_flipped
            I = rot180(current.image)
        elseif eo == :left
            I = current.image
        elseif eo == :top_flipped
            I = rotl90(current.image)
        elseif eo == :bottom_flipped
            I = reverse( rotr90(current.image), dims = 1)
        elseif eo == :right
            I = reverse( rot180(current.image), dims = 1)
        elseif eo == :left_flipped
            I = reverse( current.image, dims = 1)
        elseif eo == :top
            I = reverse( rotl90(current.image), dims = 1)
        else
            error("Orientation")
        end
        #print(eo)
        actual_image[tn*(ind[1]-1) .+ (1:tn), tm*(ind[2]-1) .+ (1:tm)] = I
        prev_edge = edge2int(I[:,end])

    end
    end

    ind = filter(i->mod1(i,10)!=10 && mod1(i,10)!=1,1:size(actual_image)[1] )

    actual_image[ind, ind]
end

function orient_tile_down_right( current, down, right)
    down_edges = intersect( down.edges, current.edges )
    right_edges = intersect( right.edges, current.edges )

    #println(current)
    #println(down)
    #println(right)

    #println(down_edges)
    #println(right_edges)

    down_ori = edge_order[ mod1( findfirst(in.(current.edges,first(down_edges))), 4) ]
    right_ori = edge_order[ mod1( findfirst(in.(current.edges,first(right_edges))), 4) ]

    #println(down_ori)
    #println(right_ori)

    if down_ori == :bottom && right_ori == :right
        return current.image
    elseif down_ori == :right && right_ori == :bottom
        return permutedims(current.image)
    elseif down_ori == :right && right_ori == :top
        return rotr90(current.image)
    elseif down_ori == :left && right_ori == :top
        return permutedims( rot180( current.image ) )
    end

    # more needed
end

using IterTools

function get_outer_edge( tiles::Dict{Int,Tile})
    edge2tile = make_edges_dict( tiles )
    get_outer_edge( edge2tile )
end
get_outer_edge( edge2tile::Dict{Int,Vector{Int}} ) = filter(p->length(p.second) == 1, edge2tile)

get_outer_tiles( outer_edge::Dict{Int,Vector{Int}}) = unique(first.(values(outer_edge)))

function get_corners( tiles::Dict{Int,Tile})
    v = get_outer_edge( tiles )
    g = groupby(p->p.second,sort( collect( v ), by=p->first(p.second)))
    g = filter!(g->length(g)==4, collect(g))
    corners = first.(last.(first.( g )))
    corners
end

function part1(tiles = read_input("input.txt"))
    prod( aoc_20.get_corners( tiles ) )
    # 23497974998093
end

# Part 2

function load_pattern(filename = "pattern.txt")
   permutedims(hcat(Vector{Char}.(readlines(filename))...)) .== '#'
end

using ImageFiltering
function part2(tiles = aoc_20.read_input("input.txt"), pattern = aoc_20.load_pattern())
    actual_image = aoc_20.assemble_image(tiles)
    n_pixels_pattern = sum(pattern)
    # Find how many patterns there are across flips and rotations
    total = 0
    for i=1:4
        total += sum(imfilter(Float64,actual_image,rotr90(pattern,i)) .≈ n_pixels_pattern)
    end
    pattern = permutedims(pattern)
    for i=1:4
        total += sum(imfilter(Float64,actual_image,rotr90(pattern,i)) .≈ n_pixels_pattern)
    end
    total
    sum(actual_image) - total*n_pixels_pattern
    # 2256
end

function __init__()
    println("Initializing Advent of Code Day 20 Solution by Mark Kittisopikul: aoc_20.jl")
    if !isinteractive()
        run()
    end
end

function run()
    if isempty(ARGS)
        files = ["input.txt"]
    else
        files = ARGS
    end
    for file in files
        @show file
        input = read_input(file)
        @show part1( input )
        @show part2( input )
    end
end

end

#tiles = aoc_20.read_input("input.txt")
#aoc_20.assemble_image( tiles )


text:
module Synthesize

using Random, Distributions

export musicFromNotes
export Music
export default, flute, flute2, organ

# instrument profiles (just a bunch of random numbers, needs revising)
default = Dict{Rational, Float64}(1 => 1)
flute = Dict{Rational, Float64}(1//2 => 1, 1 => 10^0.7, 2 => 10^0.6, 3 => 10^0.6, 4 => 10^-0.1, 5 => 1, 6 => 10^-0.3, 7 => 10^-0.7, 8 => 10^-0.7)
flute2 = Dict{Rational, Float64}(1//2 => 2^-8, 1 => 2^-7, 2 => 2^-1.4, 3 => 2^-6, 4 => 2^-1.8, 5 => 2^-6.7, 6 => 2^-3.3, 7 => 2^-6.6, 8 => 2^-3.4, 10 => 2^-4.6, 12 => 2^-5.1, 14 => 2^-5.9, 16 => 2^-5.8, 18 => 2^-6.4, 26 => 2^-7)
organ = Dict{Rational, Float64}(1//16 => 0.5, 1//4 => 2, 1//2 => 0.2, 1 => 6, 2 => 0.2, 3 => 0.1, 4 => 4, 5 => 0.1, 6 => 0.1, 8 => 2, 12 => 1, 16 => 0.5, 20 => 0.2)

# generate sound with given frequency (Hz), duration (s), and instrument profile
function genSound(freq :: Real, duration :: Real, harmonics :: Dict{<:Rational, <:Real}) :: Array{Float64,1}
    data = Array{Float64, 1}(undef, Int64(floor(duration * 44100)))
    Random.seed!(123)
    d = Normal()
    multiplier = repeat([1.0], length(data))
    # first and last n signals are decaying to avoid sudden bumps in signal
    decayLength = 500
    multiplier[1:decayLength] = (1:decayLength) / decayLength
    multiplier[end+1-decayLength:end] = (decayLength:-1:1) / decayLength
    t = (1:length(data)) / 44100
    data = 1000 * sum(h -> h[2] * sin.(h[1] * freq * 2 * π * t), harmonics)
    noise = (freq != 0 ? 20 : 0) * rand(d, length(data))
    return (data + noise) .* multiplier
end

# generate given duration of silence
silence(duration :: Real) = repeat([0], Int64(floor(duration * 44100)))

# list of common note frequencies
freqlist = Dict("0" => 0, "C" => 440.0 * 2^(-9/12), "C#" => 440.0 * 2^(-8/12), "D" => 440.0 * 2^(-7/12),
            "D#" => 440.0 * 2^(-6/12), "E" => 440.0 * 2^(-5/12), "F" => 440.0 * 2^(-4/12),
            "F#" => 440.0 * 2^(-3/12), "G" => 440.0 * 2^(-2/12), "G#" => 440.0 * 2^(-1/12),
            "A" => 440.0, "A#" => 440.0 * 2^(1/12), "H" => 440.0 * 2^(2/12))

struct Music
    notes :: String
    speed :: Number
end

# generate a music from list of notes in a string, with given instrument profile
function musicFromNotes(music :: Music, harmonics :: Dict{<:Rational, <:Real}) :: Array{Float64,1}
    notelist = split(music.notes)
    sound = []
    for note in notelist
        duration = 1
        notedata = split(note, "-")
        if length(notedata) > 1
            duration = parse(Float64, notedata[2])
        end
        octave = count(isequal('\''), notedata[1])
        notedata[1] = split(notedata[1], "'")[1]
        # append sound to current signal
        sound = vcat(sound, genSound(freqlist[notedata[1]] * 2^octave, duration*60/music.speed, harmonics))
        # append a short silence after note
        sound = vcat(sound, silence(0.05*duration*60/music.speed))
    end
    return sound
end

end

ext: jl

max_stars_repo_name: agostonsipos/synthesizer

max_issues_repo_name: agostonsipos/synthesizer

max_forks_repo_name: agostonsipos/synthesizer


text:
using Libdl

const libnames = ["libplasma", "plasma"]

found = false

global const libname = find_library(libnames, [])
if libname != nothing
    dlopen(libname)
else
    error("Plasma library not loaded: You currently need to build plasma locally")
end


text:
using NiLang, NiLang.AD
using Test, Random

@testset "softmax_crossentropy" begin
    Random.seed!(2)
    x = randn(10)
    x0 = copy(x)
    p = randn(10); p=p./maximum(p)
    res = _sce(x, p)
    imax = 0
    Z = 0.0
    out = 0.0
    xmax = 0.0
    x_ = x
    p_ = p
    @instr softmax_cross_entropy(x_, p_, imax, xmax, Z, out)
    @show Z
    @test isapprox(imax, argmax(x0), atol=1e-8)
    @test isapprox(out, res[], atol=1e-8)
    @instr (~softmax_cross_entropy)(x_, p_, imax, xmax, Z, out)
    args = x_, p_, imax, xmax, Z, out
    @test check_inv(softmax_cross_entropy, args)
    args = x_, p_, imax, xmax, Z, Loss(out)
    @test check_grad(softmax_cross_entropy, args; verbose=true)
end


text:
if Sys.isapple()
    function _font_paths()
        [
            "/Library/Fonts", # Additional fonts that can be used by all users. This is generally where fonts go if they are to be used by other applications.
            joinpath(homedir(), "Library/Fonts"), # Fonts specific to each user.
            "/Network/Library/Fonts", # Fonts shared for users on a network
        ]
    end
elseif Sys.iswindows()
    _font_paths() = [joinpath(get(ENV, "SYSTEMROOT", "C:\\Windows"), "Fonts")]
else
    function add_recursive(result, path)
        for p in readdir(path)
            pabs = joinpath(path, p)
            if isdir(pabs)
                push!(result, pabs)
                add_recursive(result, pabs)
            end
        end
    end
    function _font_paths()
        result = String[]
        for p in ("/usr/share/fonts", joinpath(homedir(), "/.fonts"), "/usr/local/share/fonts",)
            if isdir(p)
                add_recursive(result, p)
            end
        end
        result
    end
end

family_name(x::String) = replace(lowercase(x), ' ' => "") # normalize

function family_name(x::FTFont)
    family_name(x.family_name)
end

function style_name(x::FTFont)
    lowercase(x.style_name)
end

function match_font(face::FTFont, name, italic, bold)
    fname = family_name(face)
    sname = style_name(face)
    italic = italic == (sname == "italic")
    bold = bold == (sname == "bold")
    perfect_match = (fname == name) && italic && bold
    fuzzy_match = occursin(name, fname)
    score = fuzzy_match + bold + italic
    return perfect_match, fuzzy_match, score
end

function try_load(fpath)
    try
        return FTFont(fpath)
    catch e
        return nothing
    end
end

function findfont(
        name::String;
        italic::Bool=false, bold::Bool=false, additional_fonts::String=""
    )
    font_folders = copy(fontpaths())
    normalized_name = family_name(name)
    isempty(additional_fonts) || pushfirst!(font_folders, additional_fonts)
    candidates = Pair{FTFont, Int}[]
    for folder in font_folders
        for font in readdir(folder)
            fpath = joinpath(folder, font)
            face = try_load(fpath)
            face === nothing && continue
            perfect_match, fuzzy_match, score = match_font(
                face, normalized_name, italic, bold
            )
            perfect_match && return face
            if fuzzy_match
                push!(candidates, face => score)
            else
                finalize(face) # help gc a bit!
            end
        end
    end
    if !isempty(candidates)
        sort!(candidates; by=last)
        final_candidate = pop!(candidates)
        foreach(x-> finalize(x[1]), candidates)
        return final_candidate[1]
    end
    return nothing
end


text:
using Enrichr
using Test

@testset "Enrichr.jl" begin
    # Write your tests here.
end


text:
using Optim
# using FlightMechanics
# using FlightMechanics.Models

export steady_state_trim

mutable struct Trimmer
    ac::Aircraft
    aerostate::AeroState
    state::State
    env::Environment
    fcs::FCS
    turn_rate::Number
    gamma::Number
end

"""
    steady_state_trim(ac::Aircraft, fcs::FCS, env::Environment, tas::Number,
        pos::Position, psi::Number, gamma::Number, turn_rate::Number,
        α0::Number, β0::Number; show_trace=false)

Find a steady state flight condition. Steady flight is defined as flight where
the aircraft's linear and angular velocity vectors are constant in a body-fixed
reference frame (ie. body frame or wind frame).

Three different conditions can be sought:
- Steady level longitudinal flight: bank angle μ=0 and flight-path angle γ=0.
- Steady level turns: turn_rate≠0 and flight-path angle γ=0.
- Steady longitudinal climbs or descents: γ≠0 and bank angle μ=0.
- Steady turning climbs or descents: γ≠0 and bank angle μ≠0.

# Inputs
ac: aircraft to be trimmed.
fcs: flight control system. Controls given by `get_controls_trimmer` will be
  used to trim the aircraft while the rest will remain constant.
env: environment variables (atmospheric values, wind and gravity).
tas: true airspeed [m/s].
pos: position of the aircraft.
psi: heading of the aircraft [rad].
gamma: aerodynamic rate of climb of the aircraft [rad].
turn_rate: heading rate of change [rad/s].
α0, β0: Inital AOA and AOS for trimmer algorithm.
show_trace: show trimming information and algorithm trace. Optional, default=true.

# Returns
ac: trimmed aircraft.
aerostate: trimmed aerostate.
state: trimmed state.
fcs: trimmed FCS.

# References

See section 2.5 in [1] for the definition of steady-state flight condition.
See section 3.4 in [1] for the algorithm description.

- [1] Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation:
 dynamics, controls design, and autonomous systems. John Wiley & Sons.
 (page 41, formula 1.4-23)
"""
function steady_state_trim(ac::Aircraft, fcs::FCS, env::Environment,
    tas::Number, pos::Position, psi::Number, gamma::Number, turn_rate::Number,
    α0::Number, β0::Number; show_trace=false)

    alpha0 = α0
    beta0 = β0

    phi = coordinated_turn_bank(turn_rate, alpha0, beta0, tas, gamma)
    theta = climb_theta(gamma, alpha0, beta0, phi)
    att  = Attitude(psi, theta, phi)

    p, q, r = turn_rate_angular_velocity(turn_rate, theta, phi)
    ang_vel = [p, q, r]

    accel = [0., 0., 0.]
    ang_accel = [0., 0., 0.]

    state, aerostate = generate_state_aerostate(pos, att, tas, alpha0, beta0, env, ang_vel,
                                       accel, ang_accel)

    # Ensure that environment is calculated at the position given to the trimmer
    env = calculate_environment(env, get_position(state))

    # Store every necessary variable in the trimmer
    trimmer = Trimmer(ac, aerostate, state, env, fcs, turn_rate, gamma)

    # # Varibles in the trimming loop are alpha, beta, and controls.
    trim_vars0 = [alpha0, beta0] # append not fixed controls
    # lower_bounds = [-15 * DEG2RAD, -15 * DEG2RAD]
    # upper_bounds = [ 15 * DEG2RAD,  15 * DEG2RAD]
    #
    for (value, range)=zip(get_controls_trimmer(fcs), get_controls_ranges_trimmer(fcs))
    #     min, max = range
        val = value
        append!(trim_vars0, val)
    #     append!(lower_bounds, min)
    #     append!(upper_bounds, max)
    end

    # Wrapper for trim_cost_function with trimmer
    trimming_function(x) = trim_cost_function(x, trimmer)

    # Trim
    result = optimize(trimming_function, trim_vars0,
                      Optim.Options(
                        g_tol=1e-25,
                        iterations=5000,
                        show_trace=show_trace, show_every=100
                        );
                      )
    if show_trace
        println(result)
        println(Optim.minimizer(result))
    end

    # Return trimmed variables
    trimmer.ac, trimmer.aerostate, trimmer.state, trimmer.fcs
end

function trim_cost_function(trimming_variables, trimmer::Trimmer)

    # alpha, beta, tas, height are known so aero can be created
    # alpha and beta are given by trimming_variables
    # tas is fixed for the trim
    alpha, beta = trimming_variables[1:2]

    tr = trimmer.turn_rate
    tas = get_tas(trimmer.aerostate)
    gamma = trimmer.gamma
    psi = get_euler_angles(trimmer.state)[1]

    # Impose constrains
    phi = coordinated_turn_bank(tr, alpha, beta, tas, gamma)
    theta = climb_theta(gamma, alpha, beta, phi)
    att = Attitude(psi, theta, phi)
    p, q, r = turn_rate_angular_velocity(tr, theta, phi)
    ang_vel = [p, q, r]
    pos = get_position(trimmer.state)
    aerostate = trimmer.aerostate
    accel = [0., 0., 0.]
    ang_accel = [0., 0., 0.]

    env = trimmer.env

    state, aerostate = generate_state_aerostate(pos, att, tas, alpha, beta, env, ang_vel,
                                       accel, ang_accel)

    env = calculate_environment(trimmer.env, get_position(trimmer.state))

    # Set trimmer attributes
    trimmer.state = state
    trimmer.aerostate = aerostate
    trimmer.env = env

    fcs = trimmer.fcs
    ac = trimmer.ac
    grav = env.grav
    # Some controls may be fixed and the rest of them are given in trimming_variables
    set_controls_trimmer(fcs, trimming_variables[3:end]...)

    trimmer.ac = calculate_aircraft(ac, fcs, aerostate, state, grav; consume_fuel=false)
    pfm = trimmer.ac.pfm
    mass_props = get_mass_props(trimmer.ac)
    mass = mass_props.mass
    inertia = mass_props.inertia

    state = get_sixdof_euler_fixed_mass_state(trimmer.state)

    r = six_dof_euler_fixed_mass(state, mass, inertia, pfm.forces, pfm.moments)[1:6]

    return sum(r.^2)
end


text:
function step(rng, model, sampler::Nested; kwargs...)
    # Initialize particles
    # us are in unit space, vs are in prior space
    us, vs, logl = init_particles(rng, model, sampler)

    # Find least likely point
    logl_dead, idx_dead = findmin(logl)
    u_dead = us[:, idx_dead]
    v_dead = vs[:, idx_dead]

    # update weight using trapezoidal rule
    logvol = log1mexp(-1 / sampler.nactive)
    logwt = logl_dead + logvol

    # sample a new live point without bounds
    point = rand(rng, eltype(us), sampler.ndims)
    bound = Bounds.fit(Bounds.NoBounds, us)
    proposal = Proposals.Rejection()
    u, v, ll, nc = proposal(rng, v_dead, logl_dead, bound, model)

    us[:, idx_dead] .= u
    vs[:, idx_dead] .= v
    logl[idx_dead] = ll

    ncall = since_update = nc

    # update evidence and information
    logz = logwt
    h = logl_dead - logz
    logzerr = sqrt(h / sampler.nactive)
    logvol -= 1 / sampler.nactive

    sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead)
    state = (it = 1, ncall = ncall, us = us, vs = vs, logl = logl, logl_dead = logl_dead,
             logz = logz, logzerr = logzerr, h = h, logvol = logvol,
             since_update = since_update, has_bounds = false, active_bound = nothing)

    return sample, state
end

function step(rng, model, sampler, state; kwargs...)
    ## Update bounds
    pointvol = exp(state.logvol) / sampler.nactive
    # check if ready for first update
    if !state.has_bounds && state.ncall > sampler.min_ncall && state.it / state.ncall < sampler.min_eff
        @debug "First update: it=$(state.it), ncall=$(state.ncall), eff=$(state.it / state.ncall)"
        active_bound = Bounds.scale!(Bounds.fit(sampler.bounds, state.us, pointvol=pointvol), sampler.enlarge)
        since_update = 0
        has_bounds = true
    # if accepted first update, is it time to update again?
    elseif iszero(state.since_update % sampler.update_interval)
        @debug "Updating bounds: it=$(state.it), ncall=$(state.ncall), eff=$(state.it / state.ncall)"
        active_bound = Bounds.scale!(Bounds.fit(sampler.bounds, state.us, pointvol=pointvol), sampler.enlarge)
        since_update = 0
        has_bounds = true
    else
        active_bound = state.active_bound
        since_update = state.since_update + 1
        has_bounds = state.has_bounds
    end

    ## Replace least-likely active point
    # Find least likely point
    logl_dead, idx_dead = findmin(state.logl)
    u_dead = state.us[:, idx_dead]
    v_dead = state.vs[:, idx_dead]

    # sample a new live point using bounds and proposal
    if has_bounds
        point, bound = rand_live(rng, active_bound, state.us)
        if isnothing(bound)
            # live point not inside active bounds: refit them
            active_bound = Bounds.scale!(Bounds.fit(sampler.bounds, state.us, pointvol=pointvol), sampler.enlarge)
            since_update = 0
            point, bound = rand_live(rng, active_bound, point[:, :])
        end
        u, v, logl, nc = sampler.proposal(rng, point, logl_dead, bound, model)
    else
        point = rand(rng, eltype(state.us), sampler.ndims)
        bound = Bounds.fit(Bounds.NoBounds, state.us)
        proposal = Proposals.Rejection()
        u, v, logl, nc = proposal(rng, point, logl_dead, bound, model)
    end

    state.us[:, idx_dead] .= u
    state.vs[:, idx_dead] .= v
    state.logl[idx_dead] = logl

    it = state.it + 1
    ncall = state.ncall + nc
    since_update += nc

    # update weight
    logwt = state.logvol + logl_dead

    # update evidence and information
    logz = logaddexp(state.logz, logwt)
    h = (exp(logwt - logz) * logl_dead +
         exp(state.logz - logz) * (state.h + state.logz) - logz)
    logzerr = h ≥ 0 ? sqrt(h / sampler.nactive) : NaN
    logvol = state.logvol - 1 / sampler.nactive

    ## prepare returns
    sample = (u = u_dead, v = v_dead, logwt = logwt, logl = logl_dead)
    state = (it = it, ncall = ncall, us = state.us, vs = state.vs, logl = state.logl, logl_dead = logl_dead,
             logz = logz, logzerr = logzerr, h = h, logvol = logvol,
             since_update = since_update, has_bounds = has_bounds, active_bound = active_bound)

    return sample, state
end

function bundle_samples(samples,
        model::AbstractModel,
        sampler::Nested,
        state,
        ::Type{Chains};
        add_live=true,
        param_names=missing,
        check_wsum=true,
        kwargs...)

    if add_live
        samples, state = add_live_points(samples, model, sampler, state)
    end
    vals = mapreduce(t -> hcat(t.v..., exp(t.logwt - state.logz)), vcat, samples)

    if check_wsum
        wsum = sum(vals[:, end, 1])
        err = !iszero(state.logzerr) ? 3 * state.logzerr : 1e-3
        isapprox(wsum, 1, atol=err) || @warn "Weights sum to $wsum instead of 1; possible bug"
    end

    # Parameter names
    if param_names === missing
        param_names = ["Parameter $i" for i in 1:length(vals[1, :]) - 1]
    end
    push!(param_names, "weights")

    return Chains(vals, param_names, Dict(:internals => ["weights"]), evidence=state.logz), state
end

function bundle_samples(samples,
        model::AbstractModel,
        sampler::Nested,
        state,
        ::Type{Array};
        add_live=true,
        check_wsum=true,
        kwargs...)

    if add_live
        samples, state = add_live_points(samples, model, sampler, state)
    end

    vals = mapreduce(t -> hcat(t.v..., exp(t.logwt - state.logz)), vcat, samples)

    if check_wsum
        wsum = sum(vals[:, end])
        err = !iszero(state.logzerr) ? 3 * state.logzerr : 1e-3
        isapprox(wsum, 1, atol=err) || @warn "Weights sum to $wsum instead of 1; possible bug"
    end

    return vals, state
end

## Helpers

init_particles(rng, ndims, nactive, model) =
    init_particles(rng, Float64, ndims, nactive, model)

init_particles(rng, model, sampler) =
    init_particles(rng, sampler.ndims, sampler.nactive, model)

# loop and fill arrays, checking validity of points
# will retry 100 times before erroring
function init_particles(rng, T, ndims, nactive, model)
    us = rand(rng, T, ndims, nactive)
    vs_and_logl = mapslices(
        Base.Fix1(prior_transform_and_loglikelihood, model), us;
        dims=1
    )
    vs = mapreduce(first, hcat, vs_and_logl)
    logl = dropdims(map(last, vs_and_logl), dims=1)

    ntries = 1
    while true
        any(isfinite, logl) && break
        rand!(rng, us)
        vs_and_logl .= mapslices(
            Base.Fix1(prior_transform_and_loglikelihood, model), us;
            dims=1
        )
        vs .= mapreduce(first, hcat, vs_and_logl)
        map!(last, logl, vs_and_logl)
        ntries += 1
        ntries > 100 && error("After 100 attempts, could not initialize any live points with finite loglikelihood. Please check your prior transform and loglikelihood methods.")
    end

    # force -Inf to be a finite but small number to keep estimators from breaking
    @. logl[logl == -Inf] = -1e300

    return us, vs, logl
end

# add remaining live points to `samples`
function add_live_points(samples, model, sampler, state)
    logvol = -state.it / sampler.nactive - log(sampler.nactive)

    prev_logz = state.logz
    prev_h = state.h

    local logl, logz, h, logzerr
    N = length(samples)

    @inbounds for (i, idx) in enumerate(eachindex(state.logl))
        # get new point
        u = state.us[:, idx]
        v = state.vs[:, idx]
        logl = state.logl[idx]

        # update sampler
        logwt = logvol + logl
        logz = logaddexp(prev_logz, logwt)
        h = (exp(logwt - logz) * logl +
             exp(prev_logz - logz) * (prev_h + prev_logz) - logz)
        logzerr = sqrt(h / sampler.nactive)

        prev_logz = logz
        prev_h = h

        sample = (u = u, v = v, logwt = logwt, logl = logl)
        save!!(samples, sample, N + i, model, sampler)
    end

    state = (it = state.it + sampler.nactive, us = state.us, vs = state.vs, logl = logl,
            logz = logz, logzerr = logzerr, logvol = logvol,
            since_update = state.since_update, has_bounds = state.has_bounds, active_bound = state.active_bound)
    return samples, state
end


text:
##################Utility Functions ###############################

#helper function  to get the model matrix from a dataframe given a formula
#IN: A dataframe and formula object
#OUT: the model matrix
function generateX(::Type{M}, df::T, f::FormulaTerm)::M where
  {M<:AbstractMatrix, T <: AbstractDataFrame}

  #*************WARNING WARNING WARNING WARNING WARNING
  #replace the code below w/ the two lines above
  #once the missing-value coding of numerical variables is fixed see
  #https://github.com/JuliaStats/StatsModels.jl/issues/145

  #f = apply_schema(f, schema(f,df), StatisticalModel)
  #m = modelcols(f.rhs, df)

  #below there was some CUDA code- not sure if it is needed, but not working with CUDA now in any case
  #(M<:CuArray) && CUDA.allowscalar(true)
  m::M = ModelMatrix(ModelFrame(f, df)).m
  #(M<:CuArray) && CUDA.allowscalar(false)

  return m
end

#Same as above but allows for an expression
function generateX(::Type{M}, df::T, rhs::V)::M where
  {M<:AbstractMatrix, T <: AbstractDataFrame, V <: FMExpr}

  #special case of no covariates
  if isnothing(rhs) #(rhs == Symbol(""))
    return ones(Float64,size(df,1),1)
  end

  #WARNING: This is a HACK! Replace with `nothing` or something else (See StatsModels github)
  #or just redo the whole thing in the GLM framework
  #Will be able to replace the LHS w/0 in subsequent versions
  #lhs_hack::Symbol = propertynames(df)[1]
  #f::FormulaTerm = @eval(@formula($lhs_hack ~ $rhs))
  f::FormulaTerm = @eval(@formula(0 ~ $rhs))

  return generateX(M, df, f)
end

#helper function to create a one-sided formula given an expression
#IN: an expression and dataframe
#OUT: A one-sided formula object
#WARNING: Broken! if everything works drop this
#=function get1SidedFormula(RHS::T)::FormulaTerm where
  {T <: FMExpr}
  return @eval(@formula(nothing ~ $RHS)) #WARNING: THIS IS A HACK!!!!
end=#

#helper function  to get the list of symbols in an expression
#IN: A dataframe and a string object, typically representing a formula
#OUT: a vector of symbols in the string object
function symbolsFromExpr(inexpr::V)::Vector{Symbol} where {V <: FMExpr}
  isnothing(inexpr) && return Vector{Symbol}() #special case with no focal var
  symstrings::Vector{String} = split(string(inexpr),
    ['|','=','~',' ','+','*','&',')','(','-'], keepempty=false)

  filter!(s::String->typeof(Meta.parse(s))≠Int,symstrings)

  #convert the strings to symbols
  return Symbol.(symstrings)
end

#drops missings form a dataframe
function cloakmissings(::Type{M}, df::T, syms::Vector{Symbol})::SubDataFrame where {
    T<:AbstractDataFrame, M<:AbstractMatrix}
  local sdf::SubDataFrame

  #CUDA code disabled- might need to reneable it at some point
  #(M <: CuArray) && CUDA.allowscalar(true) #WARNING- this is expensive!
  sdf = view(df,completecases(df[!,syms]),:)
  #(M <: CuArray) && CUDA.allowscalar(false)

  return sdf
end

#helper function for the above, de-nulls entire dataframe
#IN: a subdataframe #OUT: subdataframe less the null entries
function cloakmissings(::Type{M}, df::T)::SubDataFrame where {T<:AbstractDataFrame, M<:AbstractMatrix}
  return cloakmissings(M, df, propertynames(df))
end

#####################FMQR Type###########################
#=This is a storage object specific to the QR optimizations
meant to hold the QR factorization and a handle for X
Values: X matrix that is decomposed, Q matrix, and R matrix
The X matrix has dimensions of NxK. R is the inverse, stored
for computational efficiency.=#

struct FMQR{M<:AbstractMatrix}
  Q::M
  R::M

  N::Int
  K::Int
  Rinv::M
end

#main constructor for the FMQR object
#IN: X Matrix
#OUT: FMQR object
function FMQR(::Type{T}, X::M)::FMQR where {T<:AbstractMatrix, M<:AbstractMatrix}
  local Q::M
  local R::M
  local Rinv::M
  local QRcompact

  if M == T
    QRcompact = qr(X)
  else
    QRcompact = qr(T(X))
  end

  R = QRcompact.R

  #hack to get around a bug involving scalar operations when pulling the Q matrix
  #hack disabled
  #Q = T<:CuArray ? CuMatrix(QRcompact.Q) : QRcompact.Q
  Q = QRcompact.Q

  try
    #hack below disabled
    #Rinv = T<:CuArray ? pinv(Matrix(R)) : pinv(R)
    Rinv = pinv(R)

  catch
    println("X matrix:")
    display(X)

    println("R matrix:")
    display(R)
    error("Inversion of R matrix failed!")
  end
  return FMQR{M}(Q,
      R,
      size(X,1),
      size(X,2),
      Rinv)
end

#default constructor, enables legacy compatability
function FMQR(X::M)::FMQR{M} where {M<:AbstractMatrix}
  return FMQR(M, X)
end

const FMClusters = Vector{Vector{C} where C<:FMData}

#####################FMLM Model Type#########################
#This is designed to hold the minimal info for a linear model
#Plus the associated FMQR object

struct FMLM{M<:AbstractMatrix, V<:AbstractVector} <: FMModel
  xqr::FMQR #handle to the xqr object
  X::M #the data and x matrix (handle)
  Y::V #Y data

  N::Int #Number of data points (rows)
  K::Int #Number of data columns
  dof::Int
  β::V #beta coefficient
  ε::V #residuals

  Xnames::Vector{String} #names of the X variables
  Yname::String #names of the Y variables

  clusters::FMClusters  #holds the clusters

end

function FMLM(X::M, Y::V; clusters::FMClusters = FMClusters([]),
        Xnames::Vector{String} = string.(:X, 1:size(X,2)),
        Yname::String = "Y", dof = size(X,1)-size(X,2),
        qrtype::Type = M)::FMLM where {
          M<:AbstractMatrix, V<:AbstractVector, C<:FMData}
    local xqr::FMQR

    xqr = FMQR(qrtype, X)

    β::V = Vector{Float64}(undef, xqr.K)
    coeff!(xqr,Y,β)
    ε::V = Vector{Float64}(undef, xqr.N)
    resid!(X,Y,β,ε)

    return FMLM{M,V}(xqr, X, Y, xqr.N, xqr.K, dof, β, ε, Xnames, Yname, clusters)
end

#=Constructor and helper method which gets required info from DataFrame
IN: The source dataframe, a formula expression for the RHS,
    the dependent variable as a symbol, optionally names for X columns,
    Y and a switch to elliminate null values
OUT: A FMLM object
NOTE: RHS expression must be ordered with factors last, otherwise
    the factor names will be incorrect =#
function FMLM(df::AbstractDataFrame,  Xexpr::FMExpr, Ysym::Symbol,
    ::Type{M} = Matrix{Float64}, ::Type{V} = Vector{Float64};
    Xnames::Vector{String}=Vector{String}(), Yname::String="Y",
    containsmissings::Bool=true, withinsym::NSymbol = nothing,
    clustersyms::C= nothing,
    qrtype::Type = M, checkwithin::Bool = false)::FMLM  where {
     M<:AbstractMatrix, V<:AbstractVector, C<:Union{NSymbol, Vector{Symbol}}}

    local Xmodelmatrix::M
    local Y::V
    local clusters::FMClusters = []

    (C <: Vector) && (length(clustersyms)>2) && error("Max of 2 clusters supported!!!")

    #get the list of symbols from the expression for X
    Xsyms::Vector{Symbol} = symbolsFromExpr(Xexpr)

    #println("flag2.2")
    #form an array of all symbols for which we want to drop nulls
    sarray::Vector{Symbol} = [Xsyms; Ysym]
    if withinsym ≠ nothing && withinsym ∉ sarray
      sarray = [sarray; withinsym]
    end
    if !isnothing(clustersyms)
      for clustersym ∈ [clustersyms;]
        (clustersym ∉ sarray) && (sarray = [sarray; clustersym])
      end
    end

    sdf::SubDataFrame = view(df[!, sarray], 1:size(df,1), :)
    #println("flag1")
    if containsmissings    # if we are going to drop nulls
      sdf = cloakmissings(M, sdf)
    end

    #get the within and clustering symbols
    if withinsym ≠ nothing
      within::Vector{W where W <: FMData}  = Vector{W where W <: FMData}(sdf[!, withinsym])
    end

    #println("flag4.2")
    if !isnothing(clustersyms)
      for clustersym ∈ [clustersyms;]
        push!(clusters, Vector(sdf[!, clustersym]))
      end
    end

    #println("flag2")
    #drop the columns which arn't going into the model matrix
    sdf = view(sdf, 1:size(sdf,1), [Xsyms; Ysym])
    #println("flag3")

    #get the factor expanded model matrix (Adds the dummies)
    Xmodelmatrix = generateX(M, sdf, Xexpr)
    Y = sdf[!, Ysym]

    #println("flag41")

    #this is a check to make sure the factor expansion follows other columns
    if length(Xsyms) > 1
        for i ∈ 2:length(Xsyms) #all types
            if typeof(sdf[!, Xsyms[i-1]])<:CategoricalArray &&
                !(typeof(sdf[!, Xsyms[i]])<:CategoricalArray)
                @warn "Type of $(typeof(sdf[i])) is preceeded by type of
                    factor. Column names are likely MISALIGNED. Put factors
                    after numerical variables in RHS expressions."
            end
        end
    end

    #Now assign the names as appropriate
    Xnamesfull::Vector{String} = Vector{String}(undef, size(Xmodelmatrix,2))
    for i ∈ 1:length(Xnamesfull)
        Xnamesfull[i] = i≤length(Xnames) ? Xnames[i] : string(:X,i)
    end

    #get the appropriate matrices and construct the FMLM object

    if (withinsym != nothing) && (!checkwithin)
        return FMLM(Xmodelmatrix, Y,
          within, clusters = clusters,
            Xnames=Xnamesfull, Yname=Yname, qrtype=qrtype)
    elseif (withinsym != nothing) && (checkwithin)
      #  println("flag3b")
        return FMLMaltwithin(Xmodelmatrix,
            Y, within, clusters = clusters,
            Xnames=Xnamesfull, Yname=Yname, qrtype=qrtype)
    else
      #println("flag3b")
      return FMLM(Xmodelmatrix,
          Y, clusters = clusters,
          Xnames=Xnamesfull, Yname=Yname, qrtype=qrtype)
    end

end

#transforms the data using the within estimator, including an adjustment to dof
function FMLM(X::M, Y::V, within::AbstractVector{W};
    clusters::FMClusters = FMClusters([]),
    Xnames::Vector{String} = (string).(:X, 1:size(X,2)),
    Yname::String = "Y",
    qrtype::Type = M)::FMLM where {
      W<:FMData, M<:AbstractMatrix, V<:AbstractVector}

    Xmat::Matrix{Float64} = X
    Yvec::Vector{Float64} = Y
    Xstart::Int = 1

    if minimum(Xmat[:,1]) ≈ maximum(Xmat[:,1]) ≈ 1.0
      @warn("intercept with within effects detected- problem is not well posed
        answer may be wrong!!! (it might be ok, but why take the risk?)
        remediate by using +0 in the rhs regression expression")
      Xstart = 2
    end

    #we get the unique values and make a two way table
    #tVars::Vector{T} = unique(tFI)
    ivars::Vector{W} = unique(within)
    iN::Int = length(ivars)
    itable::Dict{W,Int} = Dict{W,Int}(ivars[i] => i for i::Int ∈ 1:iN)

    K::Int = size(X,2)
    N::Int = size(X,1)

    #now create the codified versions
    #Xi::Matrix{Float64} = Matrix{Float64}(length(ivars),size(X,2))
    withincode::Vector{Int} = ((x::W)->itable[x]).(within)

    meanXY::Matrix{Float64} = zeros(iN,K+1)
    NXY::Vector{Int} = zeros(iN)

    #println("$N,$K,$iN,$(size(NXY))")
    @fastmath for r::Int ∈ 1:N
      meanXY[withincode[r],K+1] += Yvec[r]
      NXY[withincode[r]] += 1.0
    end

    @fastmath for c::Int ∈Xstart:K, r ∈1:N
      meanXY[withincode[r],c] += Xmat[r,c]
    end

    meanXY ./= NXY

    @fastmath for c::Int ∈Xstart:K, r::Int ∈ 1:N
      Xmat[r,c] -= meanXY[withincode[r],c]
    end

    @fastmath for r::Int ∈ 1:N
      Yvec[r] -= meanXY[withincode[r],K+1]
    end

    X = Xmat
    Y = Yvec

    return FMLM(X, Y,
      clusters=clusters, Xnames=Xnames, Yname=Yname,
      dof = N - iN - K, qrtype=qrtype)
end

#this alternate method is easier to understand, but as of 1.2 it is 20% slower
#I expect the speed of this version to increase over time
function FMLMaltwithin(X::M, Y::V, within::Vector{W};
    clusters::FMClusters = FMClusters([]),
    Xnames::Vector{String} = (string).(:X, 1:size(X,2)),
    Yname::String = "Y",
    qrtype::Type = M)::FMLM where {
      W<:FMData, M<:AbstractMatrix, V<:AbstractVector}

    #println("flag4")
    #we get the unique values and make a two way table
    #tVars::Vector{T} = unique(tFI)
    ivars::Vector{W} = unique(within)
    iN::Int = length(ivars)

    #perposefully do not use a copy command, since we are looking to replace original elements
    Xmat::Matrix{Float64} = X
    Yvec::Vector{Float64} = Y

    K::Int = size(X,2)
    N::Int = size(X,1)

    firstcolintercept::Bool = false
    if minimum(Xmat[:,1]) ≈ maximum(Xmat[:,1]) ≈ 1.0
      @warn("intercept with within effects detected- problem is not well posed
        answer may be wrong!!! (it might be ok, but why take the risk?)
        remediate by using +0 in the rhs regression expression")
      firstcolintercept = true
    end

    #println("$N,$K,$iN,$(size(NXY))")
    #warning change these if multi-threading
    meanX::Matrix{Float64} = Matrix{Float64}(undef, 1,K) #WARNING:
    local sX::SubArray
    local sY::SubArray
    iindex::Vector{Bool} = Vector{Bool}(undef, length(within))

    @fastmath for (i::Int, iv::W) ∈ enumerate(ivars)

      iindex .= within.==ivars[i]
      sX=  view(Xmat, iindex, :)
      sY =  view(Yvec, iindex)

      meanY::Float64 = mean(sY)
      meanX .= mean(sX, dims=1) #This is actually a row vector

      if firstcolintercept
        meanX[1] = 0.0 #WARNING: Assumes first column is the intercept
      end

      sY .-= meanY
      sX .-= meanX
    end

    X = Xmat
    Y = Yvec

    return FMLM(X, Y,
      clusters=clusters, Xnames=Xnames, Yname=Yname,
      dof = N - iN - K, qrtype=qrtype)
end

#this function calculates the Pearson correlation coefficient of the predicted values
#relative to the realized values
getR(args...; keyargs...) = error("use R instead")
R(lin::FMLM)::Float64 = cor(lin.Y, lin.Y .- lin.ε)

getR²(args...; keyargs...) = error("use R² instead")

function R²(lin::FMLM; adjusted::Bool = false)::Float64
  R²::Float64 = R(lin)^2
  R² = adjusted ? R² - (1-R²)*(lin.N-lin.dof+1)/lin.dof : R²

  return R²
end

#fucntion for getting regression coefficients
getTerm(args...; keyargs...) = error("use βterm instead")
βterm(lm::FMLM, s::String) = (lm.β[findfirst(lm.Xnames, s)])::Float64
βterm(lm::FMLM, s::Symbol) = βterm(lm, string(s))

#########################Project!##############################
#=gets the projection matrix

IN: THe FMQR decomposition and the memory for  the projection matrix
OUT: Writes and returns the projection matrix
NOTE: This method will allocate the projection matrix if no second
argument is supplied. =#
function project!(xqr::FMQR{M}, P::M=M(undef, xqr.N, xqr.N)
  )::Nothing where M<:AbstractMatrix

  T::Type = eltype(M)
  BLAS.gemm!('N','T',T(1.0),xqr.Q,xqr.Q,T(0.0),P)
  #use the BLAS library for fast matrix algebra
  return nothing

end

#=This version gets only the diagonal of the projection matrix
IN: THe FMQR decomposition and the memory for  the projection diagonal
OUT: Writes and returns the projection matrix=#
function project!(xqr::FMQR{M}, P::V)::Nothing where {M<:AbstractMatrix, V<:AbstractVector}
  #get the diagonal quickly
  #equivlenet to diag(Q*Q')
  P .= vec(sum(xqr.Q .^ 2, dims=2))

  return nothing

end

#=This version runs either of the above methods after generating
the xqr object given an independent variable matrix
IN: Independent variable X matrix and the memory for  the projection matrix
or vector for the diagonal
OUT: Writes and returns the projection matrix=#
function project!(X::M,
  P::AbstractArray=M(undef, size(X,1),size(X,1)))::Nothing where M<:AbstractMatrix

  return project!(FMQR(M, X),P)
end

#=simple test function to verify the projection optimizations
output should be identical to the matrix function above
IN: Independent variable X matrix, memory for the projection matrix
OUT: Writes and returns teh projection matrix=#
projectSlow!(args...; keyargs...) = error("Use projectslow! instead")
function projectslow!(X::M, P::AbstractArray) where M<:AbstractMatrix
  XXinv = M(Matrix(X' * X)\I)
  P .= X * XXinv * X'
  return P
end

#########################resid!#########################

#=Main method to get residuals
IN: X matrix the RHS variables, LHS dependent variable,
a coefficient vector, optionally the memory for the residual vector
OUT: Writes and returns the residual vector=#
getResid!(args...; keyargs...) = error("use resid! instead")
function resid!(X::M, Y::V,
  β::V, ε::V=similar(Y)) where {M<:AbstractMatrix, V<:AbstractVector}

  T::Type = eltype(V)
  return BLAS.axpy!(T(1.0), Y, BLAS.gemv!('N',T(-1.0), X, β,T(0.0),ε)) #ε=Y-Xβ
end

#########################getCoef!###########################
# Gets the coefficients in a regression

#=Main method to get the regression coefficients
IN: QR decomposition of X matrix the RHS variables, LHS dependent variable,
optionally the memory for the coefficient vector
OUT: Writes and returns the coefficient vector=#

getCoeff!(args...; keyargs...) = error("Use coeff! instead")

function coeff!(xqr::FMQR{M}, Y::V,
  β::V = V(undef, xqr.K))::V where  {M<:AbstractMatrix, V<:AbstractVector}

  #println("got here 5")
  #println("xqr.Rinv size: $(size(xqr.Rinv)) xqr.Q size:$(size(xqr.Q))")
  RinvQt::M = BLAS.gemm('N', 'T', xqr.Rinv, xqr.Q)
  #println("got here 6")

  if M<:Matrix{Float64}
    return  BLAS.gemv!('N',1.0,RinvQt,Y,0.0,β)
  else
    T::Type = eltype(M)
    return  BLAS.gemv!('N',T(1.0),RinvQt,Y,T(0.0),β)
  end
end

max_stars_repo_name: clintonTE/Finometrics

max_issues_repo_name: clintonTE/Finometrics

max_forks_repo_name: clintonTE/Finometrics


text:
"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern6Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 9 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,tmp = cache
    @.. tmp = uprev+dt*(a21*k1)
    f(k2,tmp,p,t+c1*dt)
    @.. tmp = uprev+dt*(a31*k1+a32*k2)
    f(k3,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a41*k1+a43*k3)
    f(k4,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a51*k1+a53*k3+a54*k4)
    f(k5,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a61*k1+a63*k3+a64*k4+a65*k5)
    f(k6,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a71*k1+a73*k3+a74*k4+a75*k5+a76*k6)
    f(k7,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
    f(k8,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a91*k1+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
    f(k9,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
  end
  if (allow_calc_end && length(k)< 12) || force_calc_end # Have not added the extra stages yet
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @.. tmp = uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9])
    f(rtmp,tmp,p,t+c10*dt); copyat_or_push!(k,10,rtmp)
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10])
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11])
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
  end
  nothing
end

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern7Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,tmp = cache
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a021*k1)
    f(k2,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a031*k1+a032*k2)
    f(k3,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a041*k1+a043*k3)
    f(k4,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a051*k1+a053*k3+a054*k4)
    f(k5,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a061*k1+a063*k3+a064*k4+a065*k5)
    f(k6,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a071*k1+a073*k3+a074*k4+a075*k5+a076*k6)
    f(k7,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a081*k1+a083*k3+a084*k4+a085*k5+a086*k6+a087*k7)
    f(k8,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a091*k1+a093*k3+a094*k4+a095*k5+a096*k6+a097*k7+a098*k8)
    f(k9,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a101*k1+a103*k3+a104*k4+a105*k5+a106*k6+a107*k7)
    f(k10,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
  end
  if (allow_calc_end && length(k)< 16) || force_calc_end # Have not added the extra stages yet
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9])
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11])
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
    @.. tmp = uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12])
    f(rtmp,tmp,p,t+c13*dt); copyat_or_push!(k,13,rtmp)
    @.. tmp = uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13])
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp=  uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13])
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13])
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
  end
  nothing
end

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern8Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) <13 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,tmp = cache
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a0201*k1)
    f(k2,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
    f(k3,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
    f(k4,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
    f(k5,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
    f(k6,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
    f(k7,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a0801*k1+a0804*k4+a0805*k5+a0806*k6+a0807*k7)
    f(k8,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a0901*k1+a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)
    f(k9,tmp,p,t+c9*dt)
    @.. tmp = uprev+dt*(a1001*k1+a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
    f(k10,tmp,p,t+c10*dt)
    @.. tmp = uprev+dt*(a1101*k1+a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
    f(k11,tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1201*k1+a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
    f(k12,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a1301*k1+a1304*k4+a1305*k5+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10)
    f(k13,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
    copyat_or_push!(k,11,k11)
    copyat_or_push!(k,12,k12)
    copyat_or_push!(k,13,k13)
  end
  if (allow_calc_end && length(k)< 21) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab
    @unpack tmp = cache
    @.. tmp = uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12])
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp = uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14])
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15])
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
    @.. tmp = uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16])
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,17,rtmp)
    @.. tmp = uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17])
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,18,rtmp)
    @.. tmp = uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17])
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,19,rtmp)
    @.. tmp = uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17])
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,20,rtmp)
    @.. tmp = uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17])
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,21,rtmp)
  end
  nothing
end

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern9Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,tmp = cache
    uidx = eachindex(uprev)
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a0201*k1)
    f(k2,tmp,p,t+c1*dt)
    @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
    f(k3,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
    f(k4,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
    f(k5,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
    f(k6,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
    f(k7,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a0801*k1+a0806*k6+a0807*k7)
    f(k8,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a0901*k1+a0906*k6+a0907*k7+a0908*k8)
    f(k9,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a1001*k1+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
    f(k10,tmp,p,t+c9*dt)
    @.. tmp = uprev+dt*(a1101*k1+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
    f(k11,tmp,p,t+c10*dt)
    @.. tmp = uprev+dt*(a1201*k1+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
    f(k12,tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1301*k1+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10+a1311*k11+a1312*k12)
    f(k13,tmp,p,t+c12*dt)
    @.. tmp = uprev+dt*(a1401*k1+a1406*k6+a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13)
    f(k14,tmp,p,t+c13*dt)
    @.. tmp = uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8+a1509*k9+a1510*k10+a1511*k11+a1512*k12+a1513*k13+a1514*k14)
    f(k15,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9+a1610*k10+a1611*k11+a1612*k12+a1613*k13)
    f(k16,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k8)
    copyat_or_push!(k,3,k9)
    copyat_or_push!(k,4,k10)
    copyat_or_push!(k,5,k11)
    copyat_or_push!(k,6,k12)
    copyat_or_push!(k,7,k13)
    copyat_or_push!(k,8,k14)
    copyat_or_push!(k,9,k15)
    copyat_or_push!(k,10,k16)
  end
  if (allow_calc_end && length(k)< 20) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @unpack tmp = cache
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab
    @.. tmp = uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9])
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11])
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,12,rtmp)
    @.. tmp = uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12])
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,13,rtmp)
    @.. tmp = uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13])
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp = uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14])
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15])
    f(rtmp,tmp,p,t+c22*dt); copyat_or_push!(k,16,rtmp)
    @.. tmp = uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15])
    f(rtmp,tmp,p,t+c23*dt); copyat_or_push!(k,17,rtmp)
    @.. tmp  = uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15])
    f(rtmp,tmp,p,t+c24*dt); copyat_or_push!(k,18,rtmp)
    @.. tmp = uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15])
    f(rtmp,tmp,p,t+c25*dt); copyat_or_push!(k,19,rtmp)
    @.. tmp = uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15])
    f(rtmp,tmp,p,t+c26*dt); copyat_or_push!(k,20,rtmp)
 end
  nothing
end

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern6ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 9 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98 = cache
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a21*k[1]),p,t+c1*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a31*k[1]+a32*k[2]),p,t+c2*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a41*k[1]       +a43*k[3]),p,t+c3*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a51*k[1]       +a53*k[3]+a54*k[4]),p,t+c4*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a61*k[1]       +a63*k[3]+a64*k[4]+a65*k[5]),p,t+c5*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a71*k[1]       +a73*k[3]+a74*k[4]+a75*k[5]+a76*k[6]),p,t+c6*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a81*k[1]       +a83*k[3]+a84*k[4]+a85*k[5]+a86*k[6]+a87*k[7]),p,t+dt))
    copyat_or_push!(k,9,f(uprev+dt*(a91*k[1]+a94*k[4]+a95*k[5]+a96*k[6]+a97*k[7]+a98*k[8]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 12) || force_calc_end # Have not added the extra stages yet
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache
    copyat_or_push!(k,10,f(uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9]),p,t+c10*dt))
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11]),p,t+c12*dt))
  end
  nothing
end

#=
"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern6Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 9 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,tmp = cache
    uidx = eachindex(uprev)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a21*k1[i])
    end
    f(k2,tmp,p,t+c1*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    f(k3,tmp,p,t+c2*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a41*k1[i]+a43*k3[i])
    end
    f(k4,tmp,p,t+c3*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a51*k1[i]+a53*k3[i]+a54*k4[i])
    end
    f(k5,tmp,p,t+c4*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a61*k1[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    f(k6,tmp,p,t+c5*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
    end
    f(k7,tmp,p,t+c6*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a81*k1[i]+a83*k3[i]+a84*k4[i]+a85*k5[i]+a86*k6[i]+a87*k7[i])
    end
    f(k8,tmp,p,t+dt)
    @tight_loop_macros for i in uidx
      tmp[i]=uprev[i]+dt*(a91*k1[i]+a94*k4[i]+a95*k5[i]+a96*k6[i]+a97*k7[i]+a98*k8[i])
    end
    f(k9,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
  end
  if (allow_calc_end && length(k)< 12) || force_calc_end # Have not added the extra stages yet
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1001*k[1][i]+a1004*k[4][i]+a1005*k[5][i]+a1006*k[6][i]+a1007*k[7][i]+a1008*k[8][i]+a1009*k[9][i])
    end
    f(rtmp,tmp,p,t+c10*dt); copyat_or_push!(k,10,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1101*k[1][i]+a1104*k[4][i]+a1105*k[5][i]+a1106*k[6][i]+a1107*k[7][i]+a1108*k[8][i]+a1109*k[9][i]+a1110*k[10][i])
    end
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1201*k[1][i]+a1204*k[4][i]+a1205*k[5][i]+a1206*k[6][i]+a1207*k[7][i]+a1208*k[8][i]+a1209*k[9][i]+a1210*k[10][i]+a1211*k[11][i])
    end
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
  end
  nothing
end
=#

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern7ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107 = cache
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a021*k[1]),p,t+c2*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a031*k[1]+a032*k[2]),p,t+c3*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a041*k[1]       +a043*k[3]),p,t+c4*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a051*k[1]       +a053*k[3]+a054*k[4]),p,t+c5*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a061*k[1]       +a063*k[3]+a064*k[4]+a065*k[5]),p,t+c6*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a071*k[1]       +a073*k[3]+a074*k[4]+a075*k[5]+a076*k[6]),p,t+c7*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a081*k[1]       +a083*k[3]+a084*k[4]+a085*k[5]+a086*k[6]+a087*k[7]),p,t+c8*dt))
    copyat_or_push!(k,9,f(uprev+dt*(a091*k[1]          +a093*k[3]+a094*k[4]+a095*k[5]+a096*k[6]+a097*k[7]+a098*k[8]),p,t+dt))
    copyat_or_push!(k,10,f(uprev+dt*(a101*k[1]          +a103*k[3]+a104*k[4]+a105*k[5]+a106*k[6]+a107*k[7]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 16) || force_calc_end # Have not added the extra stages yet
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11]),p,t+c12*dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12]),p,t+c13*dt))
    copyat_or_push!(k,14,f(uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13]),p,t+c14*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13]),p,t+c15*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13]),p,t+c16*dt))
  end
  nothing
end

#=
"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern7Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,tmp = cache
    uidx = eachindex(uprev)
    f(k1,uprev,p,t)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a021*k1[i])
    end
    f(k2,tmp,p,t+c2*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a031*k1[i]+a032*k2[i])
    end
    f(k3,tmp,p,t+c3*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a041*k1[i]+a043*k3[i])
    end
    f(k4,tmp,p,t+c4*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a051*k1[i]+a053*k3[i]+a054*k4[i])
    end
    f(k5,tmp,p,t+c5*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a061*k1[i]+a063*k3[i]+a064*k4[i]+a065*k5[i])
    end
    f(k6,tmp,p,t+c6*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a071*k1[i]+a073*k3[i]+a074*k4[i]+a075*k5[i]+a076*k6[i])
    end
    f(k7,tmp,p,t+c7*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a081*k1[i]+a083*k3[i]+a084*k4[i]+a085*k5[i]+a086*k6[i]+a087*k7[i])
    end
    f(k8,tmp,p,t+c8*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a091*k1[i]+a093*k3[i]+a094*k4[i]+a095*k5[i]+a096*k6[i]+a097*k7[i]+a098*k8[i])
    end
    f(k9,tmp,p,t+dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a101*k1[i]+a103*k3[i]+a104*k4[i]+a105*k5[i]+a106*k6[i]+a107*k7[i])
    end
    f(k10,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
  end
  if (allow_calc_end && length(k)< 16) || force_calc_end # Have not added the extra stages yet
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab
    uidx = eachindex(uprev)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1101*k[1][i]+a1104*k[4][i]+a1105*k[5][i]+a1106*k[6][i]+a1107*k[7][i]+a1108*k[8][i]+a1109*k[9][i])
    end
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1201*k[1][i]+a1204*k[4][i]+a1205*k[5][i]+a1206*k[6][i]+a1207*k[7][i]+a1208*k[8][i]+a1209*k[9][i]+a1211*k[11][i])
    end
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1301*k[1][i]+a1304*k[4][i]+a1305*k[5][i]+a1306*k[6][i]+a1307*k[7][i]+a1308*k[8][i]+a1309*k[9][i]+a1311*k[11][i]+a1312*k[12][i])
    end
    f(rtmp,tmp,p,t+c13*dt); copyat_or_push!(k,13,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1401*k[1][i]+a1404*k[4][i]+a1405*k[5][i]+a1406*k[6][i]+a1407*k[7][i]+a1408*k[8][i]+a1409*k[9][i]+a1411*k[11][i]+a1412*k[12][i]+a1413*k[13][i])
    end
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @tight_loop_macros for i in uidx
      tmp[i]=  uprev[i]+dt*(a1501*k[1][i]+a1504*k[4][i]+a1505*k[5][i]+a1506*k[6][i]+a1507*k[7][i]+a1508*k[8][i]+a1509*k[9][i]+a1511*k[11][i]+a1512*k[12][i]+a1513*k[13][i])
    end
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1601*k[1][i]+a1604*k[4][i]+a1605*k[5][i]+a1606*k[6][i]+a1607*k[7][i]+a1608*k[8][i]+a1609*k[9][i]+a1611*k[11][i]+a1612*k[12][i]+a1613*k[13][i])
    end
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
  end
  nothing
end
=#

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern8ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) <13 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310 = cache
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a0201*k[1]),p,t+c2*dt ))
    copyat_or_push!(k,3,f(uprev+dt*(a0301*k[1]+a0302*k[2]),p,t+c3*dt ))
    copyat_or_push!(k,4,f(uprev+dt*(a0401*k[1]       +a0403*k[3]),p,t+c4*dt ))
    copyat_or_push!(k,5,f(uprev+dt*(a0501*k[1]       +a0503*k[3]+a0504*k[4]),p,t+c5*dt ))
    copyat_or_push!(k,6,f(uprev+dt*(a0601*k[1]                +a0604*k[4]+a0605*k[5]),p,t+c6*dt ))
    copyat_or_push!(k,7,f(uprev+dt*(a0701*k[1]                +a0704*k[4]+a0705*k[5]+a0706*k[6]),p,t+c7*dt ))
    copyat_or_push!(k,8,f(uprev+dt*(a0801*k[1]                +a0804*k[4]+a0805*k[5]+a0806*k[6]+a0807*k[7]),p,t+c8*dt ))
    copyat_or_push!(k,9,f(uprev+dt*(a0901*k[1]                +a0904*k[4]+a0905*k[5]+a0906*k[6]+a0907*k[7]+a0908*k[8]),p,t+c9*dt ))
    copyat_or_push!(k,10,f(uprev+dt*(a1001*k[1]                +a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9]),p,t+c10*dt))
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]                +a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]                +a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11]),p,t+    dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1301*k[1]                +a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1310*k[10]),p,t+    dt))
  end
  if (allow_calc_end && length(k)< 21) || force_calc_end # Have not added the extra stages yet
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache
    copyat_or_push!(k,14,f(uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12]),p,t+c14*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14]),p,t+c15*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15]),p,t+c16*dt))
    copyat_or_push!(k,17,f(uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16]),p,t+c17*dt))
    copyat_or_push!(k,18,f(uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17]),p,t+c18*dt))
    copyat_or_push!(k,19,f(uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17]),p,t+c19*dt))
    copyat_or_push!(k,20,f(uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17]),p,t+c20*dt))
    copyat_or_push!(k,21,f(uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17]),p,t+c21*dt))
  end
  nothing
end

#=
"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern8Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) <13 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,tmp = cache
    uidx = eachindex(uprev)
    f(k1,uprev,p,t)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0201*k1[i])
    end
    f(k2,tmp,p,t+c2*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0301*k1[i]+a0302*k2[i])
    end
    f(k3,tmp,p,t+c3*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
    end
    f(k4,tmp,p,t+c4*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
    end
    f(k5,tmp,p,t+c5*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
    end
    f(k6,tmp,p,t+c6*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
    end
    f(k7,tmp,p,t+c7*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0801*k1[i]+a0804*k4[i]+a0805*k5[i]+a0806*k6[i]+a0807*k7[i])
    end
    f(k8,tmp,p,t+c8*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0901*k1[i]+a0904*k4[i]+a0905*k5[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
    end
    f(k9,tmp,p,t+c9*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1001*k1[i]+a1004*k4[i]+a1005*k5[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
    end
    f(k10,tmp,p,t+c10*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1101*k1[i]+a1104*k4[i]+a1105*k5[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
    end
    f(k11,tmp,p,t+c11*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1201*k1[i]+a1204*k4[i]+a1205*k5[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
    end
    f(k12,tmp,p,t+dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1301*k1[i]+a1304*k4[i]+a1305*k5[i]+a1306*k6[i]+a1307*k7[i]+a1308*k8[i]+a1309*k9[i]+a1310*k10[i])
    end
    f(k13,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
    copyat_or_push!(k,11,k11)
    copyat_or_push!(k,12,k12)
    copyat_or_push!(k,13,k13)
  end
  if (allow_calc_end && length(k)< 21) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab
    @unpack tmp = cache
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1401*k[1][i]+a1406*k[6][i]+a1407*k[7][i]+a1408*k[8][i]+a1409*k[9][i]+a1410*k[10][i]+a1411*k[11][i]+a1412*k[12][i])
    end
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1501*k[1][i]+a1506*k[6][i]+a1507*k[7][i]+a1508*k[8][i]+a1509*k[9][i]+a1510*k[10][i]+a1511*k[11][i]+a1512*k[12][i]+a1514*k[14][i])
    end
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1601*k[1][i]+a1606*k[6][i]+a1607*k[7][i]+a1608*k[8][i]+a1609*k[9][i]+a1610*k[10][i]+a1611*k[11][i]+a1612*k[12][i]+a1614*k[14][i]+a1615*k[15][i])
    end
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1701*k[1][i]+a1706*k[6][i]+a1707*k[7][i]+a1708*k[8][i]+a1709*k[9][i]+a1710*k[10][i]+a1711*k[11][i]+a1712*k[12][i]+a1714*k[14][i]+a1715*k[15][i]+a1716*k[16][i])
    end
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,17,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1801*k[1][i]+a1806*k[6][i]+a1807*k[7][i]+a1808*k[8][i]+a1809*k[9][i]+a1810*k[10][i]+a1811*k[11][i]+a1812*k[12][i]+a1814*k[14][i]+a1815*k[15][i]+a1816*k[16][i]+a1817*k[17][i])
    end
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,18,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1901*k[1][i]+a1906*k[6][i]+a1907*k[7][i]+a1908*k[8][i]+a1909*k[9][i]+a1910*k[10][i]+a1911*k[11][i]+a1912*k[12][i]+a1914*k[14][i]+a1915*k[15][i]+a1916*k[16][i]+a1917*k[17][i])
    end
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,19,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2001*k[1][i]+a2006*k[6][i]+a2007*k[7][i]+a2008*k[8][i]+a2009*k[9][i]+a2010*k[10][i]+a2011*k[11][i]+a2012*k[12][i]+a2014*k[14][i]+a2015*k[15][i]+a2016*k[16][i]+a2017*k[17][i])
    end
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,20,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2101*k[1][i]+a2106*k[6][i]+a2107*k[7][i]+a2108*k[8][i]+a2109*k[9][i]+a2110*k[10][i]+a2111*k[11][i]+a2112*k[12][i]+a2114*k[14][i]+a2115*k[15][i]+a2116*k[16][i]+a2117*k[17][i])
    end
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,21,rtmp)
  end
  nothing
end
=#

"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern9ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if (length(k) < 10 || always_calc_begin) || ((allow_calc_end && length(k)< 20) || force_calc_end)
  end
  if length(k) < 10 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613 = cache
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a0201*k[1]),p,t+c1*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a0301*k[1]+a0302*k[2]),p,t+c2*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a0401*k[1]       +a0403*k[3]),p,t+c3*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a0501*k[1]       +a0503*k[3]+a0504*k[4]),p,t+c4*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a0601*k[1]                +a0604*k[4]+a0605*k[5]),p,t+c5*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a0701*k[1]                +a0704*k[4]+a0705*k[5]+a0706*k[6]),p,t+c6*dt))
    copyat_or_push!(k,2,f(uprev+dt*(a0801*k[1]                                  +a0806*k[6]+a0807*k[7]),p,t+c7*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a0901*k[1]                                  +a0906*k[6]+a0907*k[7]+a0908*k[2]),p,t+c8*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a1001*k[1]                                  +a1006*k[6]+a1007*k[7]+a1008*k[2]+a1009*k[3]),p,t+c9*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a1101*k[1]                                  +a1106*k[6]+a1107*k[7]+a1108*k[2]+a1109*k[3]+a1110*k[4]),p,t+c10*dt))
    temp6 = recursivecopy(k[6])
    temp7 = recursivecopy(k[7])
    copyat_or_push!(k,6,f(uprev+dt*(a1201*k[1]                                  +a1206*temp6+a1207*temp7+a1208*k[2]+a1209*k[3]+a1210*k[4]+a1211*k[5]),p,t+c11*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a1301*k[1]                                  +a1306*temp6+a1307*temp7+a1308*k[2]+a1309*k[3]+a1310*k[4]+a1311*k[5]+a1312*k[6]),p,t+c12*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a1401*k[1]                                  +a1406*temp6+a1407*temp7+a1408*k[2]+a1409*k[3]+a1410*k[4]+a1411*k[5]+a1412*k[6]+a1413*k[7]),p,t+c13*dt))
    copyat_or_push!(k,9,f(uprev+dt*(a1501*k[1]                                  +a1506*temp6+a1507*temp7+a1508*k[2]+a1509*k[3]+a1510*k[4]+a1511*k[5]+a1512*k[6]+a1513*k[7]+a1514*k[8]),p,t+dt))
    copyat_or_push!(k,10,f(uprev+dt*(a1601*k[1]                                  +a1606*temp6+a1607*temp7+a1608*k[2]+a1609*k[3]+a1610*k[4]+a1611*k[5]+a1612*k[6]+a1613*k[7]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 20) || force_calc_end # Have not added the extra stages yet
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache
    copyat_or_push!(k,11,f(uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9]),p,t+c17*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11]),p,t+c18*dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12]),p,t+c19*dt))
    copyat_or_push!(k,14,f(uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13]),p,t+c20*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14]),p,t+c21*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15]),p,t+c22*dt))
    copyat_or_push!(k,17,f(uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15]),p,t+c23*dt))
    copyat_or_push!(k,18,f(uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15]),p,t+c24*dt))
    copyat_or_push!(k,19,f(uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15]),p,t+c25*dt))
    copyat_or_push!(k,20,f(uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15]),p,t+c26*dt))
  end
  nothing
end

#=
"""

"""
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern9Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,tmp = cache
    uidx = eachindex(uprev)
    f(k1,uprev,p,t)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0201*k1[i])
    end
    f(k2,tmp,p,t+c1*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0301*k1[i]+a0302*k2[i])
    end
    f(k3,tmp,p,t+c2*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
    end
    f(k4,tmp,p,t+c3*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
    end
    f(k5,tmp,p,t+c4*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
    end
    f(k6,tmp,p,t+c5*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
    end
    f(k7,tmp,p,t+c6*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0801*k1[i]+a0806*k6[i]+a0807*k7[i])
    end
    f(k8,tmp,p,t+c7*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a0901*k1[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
    end
    f(k9,tmp,p,t+c8*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1001*k1[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
    end
    f(k10,tmp,p,t+c9*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1101*k1[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
    end
    f(k11,tmp,p,t+c10*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1201*k1[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
    end
    f(k12,tmp,p,t+c11*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1301*k1[i]+a1306*k6[i]+a1307*k7[i]+a1308*k8[i]+a1309*k9[i]+a1310*k10[i]+a1311*k11[i]+a1312*k12[i])
    end
    f(k13,tmp,p,t+c12*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1401*k1[i]+a1406*k6[i]+a1407*k7[i]+a1408*k8[i]+a1409*k9[i]+a1410*k10[i]+a1411*k11[i]+a1412*k12[i]+a1413*k13[i])
    end
    f(k14,tmp,p,t+c13*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1501*k1[i]+a1506*k6[i]+a1507*k7[i]+a1508*k8[i]+a1509*k9[i]+a1510*k10[i]+a1511*k11[i]+a1512*k12[i]+a1513*k13[i]+a1514*k14[i])
    end
    f(k15,tmp,p,t+dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1601*k1[i]+a1606*k6[i]+a1607*k7[i]+a1608*k8[i]+a1609*k9[i]+a1610*k10[i]+a1611*k11[i]+a1612*k12[i]+a1613*k13[i])
    end
    f(k16,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k8)
    copyat_or_push!(k,3,k9)
    copyat_or_push!(k,4,k10)
    copyat_or_push!(k,5,k11)
    copyat_or_push!(k,6,k12)
    copyat_or_push!(k,7,k13)
    copyat_or_push!(k,8,k14)
    copyat_or_push!(k,9,k15)
    copyat_or_push!(k,10,k16)
  end
  if (allow_calc_end && length(k)< 20) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @unpack tmp = cache
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1701*k[1][i]+a1708*k[2][i]+a1709*k[3][i]+a1710*k[4][i]+a1711*k[5][i]+a1712*k[6][i]+a1713*k[7][i]+a1714*k[8][i]+a1715*k[9][i])
    end
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,11,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1801*k[1][i]+a1808*k[2][i]+a1809*k[3][i]+a1810*k[4][i]+a1811*k[5][i]+a1812*k[6][i]+a1813*k[7][i]+a1814*k[8][i]+a1815*k[9][i]+a1817*k[11][i])
    end
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,12,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a1901*k[1][i]+a1908*k[2][i]+a1909*k[3][i]+a1910*k[4][i]+a1911*k[5][i]+a1912*k[6][i]+a1913*k[7][i]+a1914*k[8][i]+a1915*k[9][i]+a1917*k[11][i]+a1918*k[12][i])
    end
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,13,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2001*k[1][i]+a2008*k[2][i]+a2009*k[3][i]+a2010*k[4][i]+a2011*k[5][i]+a2012*k[6][i]+a2013*k[7][i]+a2014*k[8][i]+a2015*k[9][i]+a2017*k[11][i]+a2018*k[12][i]+a2019*k[13][i])
    end
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,14,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2101*k[1][i]+a2108*k[2][i]+a2109*k[3][i]+a2110*k[4][i]+a2111*k[5][i]+a2112*k[6][i]+a2113*k[7][i]+a2114*k[8][i]+a2115*k[9][i]+a2117*k[11][i]+a2118*k[12][i]+a2119*k[13][i]+a2120*k[14][i])
    end
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,15,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2201*k[1][i]+a2208*k[2][i]+a2209*k[3][i]+a2210*k[4][i]+a2211*k[5][i]+a2212*k[6][i]+a2213*k[7][i]+a2214*k[8][i]+a2215*k[9][i]+a2217*k[11][i]+a2218*k[12][i]+a2219*k[13][i]+a2220*k[14][i]+a2221*k[15][i])
    end
    f(rtmp,tmp,p,t+c22*dt); copyat_or_push!(k,16,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2301*k[1][i]+a2308*k[2][i]+a2309*k[3][i]+a2310*k[4][i]+a2311*k[5][i]+a2312*k[6][i]+a2313*k[7][i]+a2314*k[8][i]+a2315*k[9][i]+a2317*k[11][i]+a2318*k[12][i]+a2319*k[13][i]+a2320*k[14][i]+a2321*k[15][i])
    end
    f(rtmp,tmp,p,t+c23*dt); copyat_or_push!(k,17,rtmp)
    @tight_loop_macros for i in uidx
      tmp[i]  = uprev[i]+dt*(a2401*k[1][i]+a2408*k[2][i]+a2409*k[3][i]+a2410*k[4][i]+a2411*k[5][i]+a2412*k[6][i]+a2413*k[7][i]+a2414*k[8][i]+a2415*k[9][i]+a2417*k[11][i]+a2418*k[12][i]+a2419*k[13][i]+a2420*k[14][i]+a2421*k[15][i])
    end
    f(rtmp,tmp,p,t+c24*dt); copyat_or_push!(k,18,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2501*k[1][i]+a2508*k[2][i]+a2509*k[3][i]+a2510*k[4][i]+a2511*k[5][i]+a2512*k[6][i]+a2513*k[7][i]+a2514*k[8][i]+a2515*k[9][i]+a2517*k[11][i]+a2518*k[12][i]+a2519*k[13][i]+a2520*k[14][i]+a2521*k[15][i])
    end
    f(rtmp,tmp,p,t+c25*dt); copyat_or_push!(k,19,rtmp)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a2601*k[1][i]+a2608*k[2][i]+a2609*k[3][i]+a2610*k[4][i]+a2611*k[5][i]+a2612*k[6][i]+a2613*k[7][i]+a2614*k[8][i]+a2615*k[9][i]+a2617*k[11][i]+a2618*k[12][i]+a2619*k[13][i]+a2620*k[14][i]+a2621*k[15][i])
    end
    f(rtmp,tmp,p,t+c26*dt); copyat_or_push!(k,20,rtmp)
 end
  nothing
end
=#