text:
# --------------------------------------------------------------------
# LabelEncoding aggregate logic

aggregate_fraction(numerator, denominator, labels, avgmode) =
    throw(ArgumentError("$avgmode is not supported"))

aggregate_fraction(numerator, denominator, labels, ::AvgMode.None) =
    Dict(Pair.(labels, numerator ./ denominator))

aggregate_fraction(numerator, denominator, labels, ::AvgMode.Macro) =
    mean(Base.Broadcast.broadcasted(/, numerator, denominator))

aggregate_fraction(numerator, denominator, labels, ::AvgMode.Micro) =
    sum(numerator) / sum(denominator)

# TODO: add weighted versions

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

"""
This will allow for readable function definition such as

    @reduce_fraction \"\"\"
    some documentation for the function
    \"\"\" ->
    precision := true_positives / predicted_positive

with the added bonus of the functions being implemented in such
a way as to avoid memory allocation and being able to compute
their results in just one pass.
"""
macro reduce_fraction(all)
    @assert all.head == :->
    docstr = all.args[1]
    expr = all.args[2].args[2]
    @assert expr.head == :(:=)
    fun = expr.args[1]
    @assert expr.args[2].args[1] == :/
    numer_fun = expr.args[2].args[2]
    denom_fun = expr.args[2].args[3]
    esc(quote

        # explicit binary case
        ($fun)(targets, outputs, avgmode::AverageMode) =
            reduce_fraction($numer_fun, $denom_fun,
                            targets, outputs, avgmode)

        ($fun)(targets, outputs; avgmode=AvgMode.None()) =
            reduce_fraction($numer_fun, $denom_fun,
                            targets, outputs,
                            convert(AverageMode, avgmode))

        ($fun)(targets, outputs, encoding; avgmode=AvgMode.None()) =
            reduce_fraction($numer_fun, $denom_fun,
                            targets, outputs, encoding,
                            convert(AverageMode, avgmode))

        ($fun)(targets, outputs, encoding, avgmode) =
            reduce_fraction($numer_fun, $denom_fun,
                            targets, outputs, encoding, avgmode)

        # for objects like BinaryMetricsCache
        ($fun)(object) = ($numer_fun)(object) / ($denom_fun)(object)

        # add documentation to function
        @doc """
            $($(string(fun)))(targets, outputs, [encoding], [avgmode = :none]) -> Union{Float64, Dict}

        $($docstr)
        """ ($fun)
    end)
end

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

# explicit binary case
function reduce_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        encoding::BinaryLabelEncoding,
        avgmode::AvgMode.None)
    @_dimcheck length(targets) == length(outputs)
    numer::Int = 0; denom::Int = 0
    @inbounds for I in eachindex(targets, outputs)
        target = targets[I]
        output = outputs[I]
        numer += numer_fun(target, output, encoding)
        denom += denom_fun(target, output, encoding)
    end
    (numer / denom)::Float64
end

# multiclass case. The function aggregate_fraction dispatches on the mode
function reduce_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        encoding::LabelEncoding,
        avgmode::AverageMode)
    @_dimcheck length(targets) == length(outputs)
    labels = label(encoding)
    n = length(labels)
    numer = zeros(n); denom = zeros(n)
    ovr = [LabelEnc.OneVsRest(l) for l in labels]
    @inbounds for I in eachindex(targets, outputs)
        target = targets[I]
        output = outputs[I]
        for j in 1:n
            numer[j] += numer_fun(target, output, ovr[j])
            denom[j] += denom_fun(target, output, ovr[j])
        end
    end
    aggregate_fraction(numer, denom, labels, avgmode)
end

# fallback for native labels as plain vector
function reduce_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        labels::AbstractVector,
        avgmode::AverageMode)
    reduce_fraction(numer_fun, denom_fun,
                    targets, outputs,
                    LabelEnc.NativeLabels(labels),
                    avgmode)
end

# fallbacks to try and choose compare mode automatically
function reduce_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        avgmode::AverageMode)
    reduce_fraction(numer_fun, denom_fun,
                    targets, outputs,
                    _labelenc(targets, outputs),
                    avgmode)
end

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

"""
This will allow for readable function definition such as

    @map_fraction \"\"\"
    some documentation for the function
    \"\"\" ->
    positive_likelihood_ratio := true_positive_rate / false_positive_rate

with the added bonus of the functions being implemented in such
a way as to avoid memory allocation and being able to compute
their results in just one pass.
"""
macro map_fraction(all)
    @assert all.head == :->
    docstr = all.args[1]
    expr = all.args[2].args[2]
    @assert expr.head == :(:=)
    fun = expr.args[1]
    @assert expr.args[2].args[1] == :/
    numer_fun = expr.args[2].args[2]
    denom_fun = expr.args[2].args[3]
    esc(quote

        # explicit binary case
        ($fun)(targets, outputs, avgmode::AverageMode) =
            map_fraction($numer_fun, $denom_fun,
                         targets, outputs, avgmode)

        ($fun)(targets, outputs; avgmode=AvgMode.None()) =
            map_fraction($numer_fun, $denom_fun,
                         targets, outputs,
                         convert(AverageMode, avgmode))

        ($fun)(targets, outputs, encoding; avgmode=AvgMode.None()) =
            map_fraction($numer_fun, $denom_fun,
                         targets, outputs, encoding,
                         convert(AverageMode, avgmode))

        ($fun)(targets, outputs, encoding, avgmode) =
            map_fraction($numer_fun, $denom_fun,
                         targets, outputs, encoding, avgmode)

        # for objects like BinaryMetricsCache
        ($fun)(object) = ($numer_fun)(object) / ($denom_fun)(object)

        # add documentation to function
        @doc """
            $($(string(fun)))(targets, outputs, [encoding], [avgmode = :none]) -> Union{Float64, Dict}

        $($docstr)
        """ ($fun)
    end)
end

# generic binary and multiclass case (except None). returns float
function map_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        encoding::LabelEncoding,
        avgmode::AverageMode)
    numer = numer_fun(targets, outputs, encoding, avgmode)
    denom = denom_fun(targets, outputs, encoding, avgmode)
    numer / denom
end

# binary case for AvgMode.None, because it also returns a float
function map_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        encoding::BinaryLabelEncoding,
        avgmode::AvgMode.None)
    numer = numer_fun(targets, outputs, encoding, avgmode)
    denom = denom_fun(targets, outputs, encoding, avgmode)
    numer / denom
end

# multiclass case for AvgMode.None should return a Dict again
function map_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        encoding::LabelEncoding,
        avgmode::AvgMode.None)
    numer = numer_fun(targets, outputs, encoding, avgmode)
    denom = denom_fun(targets, outputs, encoding, avgmode)
    Dict(label => numer[label] / denom[label] for label in keys(numer))
end

# fallback for native labels as plain vector
function map_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        labels::AbstractVector,
        avgmode::AverageMode)
    map_fraction(numer_fun, denom_fun,
                 targets, outputs,
                 LabelEnc.NativeLabels(labels),
                 avgmode)
end

# fallbacks to try and choose compare mode automatically
function map_fraction(
        numer_fun::Function,
        denom_fun::Function,
        targets::AbstractArray,
        outputs::AbstractArray,
        avgmode::AverageMode)
    map_fraction(numer_fun, denom_fun,
                 targets, outputs,
                 _labelenc(targets, outputs),
                 avgmode)
end


text:
function FieldEventState(arg0::FieldEventDetector)
    return FieldEventState((FieldEventDetector,), arg0)
end

function do_event(obj::FieldEventState, arg0::FieldSpacecraftState)
    return jcall(obj, "doEvent", FieldEventState_EventOccurrence, (FieldSpacecraftState,), arg0)
end

function evaluate_step(obj::FieldEventState, arg0::FieldOrekitStepInterpolator)
    return jcall(obj, "evaluateStep", jboolean, (FieldOrekitStepInterpolator,), arg0)
end

function get_event_date(obj::FieldEventState)
    return jcall(obj, "getEventDate", FieldAbsoluteDate, ())
end

function get_event_detector(obj::FieldEventState)
    return jcall(obj, "getEventDetector", FieldEventDetector, ())
end

function get_pending_event(obj::FieldEventState)
    return jcall(obj, "getPendingEvent", jboolean, ())
end

function init(obj::FieldEventState, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate), arg0, arg1)
end

function reinitialize_begin(obj::FieldEventState, arg0::FieldOrekitStepInterpolator)
    return jcall(obj, "reinitializeBegin", void, (FieldOrekitStepInterpolator,), arg0)
end

function try_advance(obj::FieldEventState, arg0::FieldSpacecraftState, arg1::FieldOrekitStepInterpolator)
    return jcall(obj, "tryAdvance", jboolean, (FieldSpacecraftState, FieldOrekitStepInterpolator), arg0, arg1)
end

ext: jl


text:
# upper bound tests
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();
#pmap(i -> importIpopt(),1:30);

dataList = Dict();
iterMax = 20;

caseList = [13,33,123];
T = 96;
τ = Int64(1/6*T);
Δt = 0.25;
N = 5;
pathTrain = load("pathHist_600.jld");

for ci in 1:length(caseList)
    for j in procs()
        remotecall_fetch(readInData,j,ci,caseList,T,τ);
    end
    pathDict = pathTrain["pathDict"][ci][T];

    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, 30, false,false, 2, 2);

    startPGT = time();
    cutDictPG = preGen(T, Δt, N, iterMax);
    preGenT = time() - startPGT;
    cutDictPGOri = deepcopy(cutDictPG);

    startT = time();
    # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false,
    #     Int64(round(100/N)),Int64(round(100/N)), cutDictPG, false, 200, [], 0, pathDict);
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
        max(Int64(round(500/N)),20), max(Int64(round(500/N)),20), cutDictPG, true, 500, [], 0, pathDict);
    elapsedT = time() - startT;
    dataList = [LBHist,UBHist,UBuHist,UBlHist,timeHist,elapsedT,preGenT];

    save("UBResults_$(ci).jld","UBOut",dataList);
end


text:
using Turing
using Test
using Random

Random.seed!(125)

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt.(s))
  x[1] ~ Normal(m, sqrt.(s))
  x[2] ~ Normal(m, sqrt.(s))
  return s, m
end

chain = sample(gdemo([1.5, 2.0]), SGLD(10000, 0.5))
# Note: samples are weigthed by step sizes cf 4.2 in paper
s_res1weightedMean = sum(chain[:lf_eps] .* chain[:s]) / sum(chain[:lf_eps])
m_res1weightedMean = sum(chain[:lf_eps] .* chain[:m]) / sum(chain[:lf_eps])
@test s_res1weightedMean ≈ 49/24 atol=0.2
@test m_res1weightedMean ≈ 7/6 atol=0.2


text:
### ntt.jl
#
# Copyright (C) 2019 Jakub Wronowski.
#
# Maintainer: Jakub Wronowski <jakubwro@users.noreply.github.com>
# Keywords: number theoretic transform
#
# This file is a part of NumberTheoreticTransforms.jl.
#
# License is MIT.
#
### Commentary:
#
# This file contains implementation of general Number Theoretic Transform.
#
### Code:

export ntt, intt

"""
Constraints on NTT params to ensure that inverse can be computed
"""
function validate(N, g, q)
    @assert mod(q - 1, N) == 0
    @assert powermod(g, N, q) == 1
    @assert !(1 in powermod.(g, 2:N-1, q)) # this may be redundant
    @assert gcd(N, q) == 1
end

"""
    ntt(x::Array{T,1}, g, q) -> Array{T,1}

The [Number Theoretic Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general)#Number-theoretic_transform)
transforms data in a similar fashion to DFT, but instead complex roots of unity
it uses integer roots with all operation defined in a finite field (modulo an
integer number).

`ntt` function implements Number Theoretic Transform directly from the formula,
so it is flexible about choosing transformation params but lacks performance.

``\\bar{x}_k = \\sum_{n=1}^N{x_n g^{(n-1)(k-1)} } \\mod q``

There is also a few constraints on choosing parameters and input length to
ensure that inverse exists and equals to the original input:
-   ``g`` must ``N``-th root of one in modulo ``q`` arithmetic
-   ``q-1`` mod ``N`` must be equal zero
-   ``q`` must be grater than maximum element present in transformed array

To find parameter set you may use fint-ntt script.

The arguments of `ntt` function are

-   `x`: input data, its elements must be smaller than q
-   `g`: transform power base, must have inversion modulo q
-   `q`: defines modulo arithmetic
"""
function ntt(x::Array{T,1}, g::T, q::T) where {T<:Integer}
    N = length(x)
    #TODO: more validation of p,q, decompose it to struct
    #TODO: create transform object that validates input in the constructor

    validate(N, g, q)
    (lo, hi) = extrema(x)
    @assert lo >= 0
    @assert hi <= q-1

    t = [powermod(g, n * k, q) for n in 0:N-1, k in 0:N-1]
    #TODO: make result of ntt a struct that will hold infrmation about g
    return mod.(t * x, q)
end

# TODO: change implementation to handle any number of dimensions
# it can be done calling N-1 dimensional transform in a loop

function ntt(x::Array{T,2}, g::T, q::T) where {T<:Integer}
    N, M = size(x)
    @assert N == M #TODO: make it work for N != M (need different g for each dim)
    y = zeros(T, size(x))

    for n in 1:N
        y[n, :] = ntt(x[n, :], g, q)
    end

    for m in 1:M
        y[:, m] = ntt(y[:, m], g, q)
    end

    return y
end

"""
    intt(y::Array{T,1}, g, q) -> Array{T,1}

Inverse Number Theoretic Transform implementation directly from the formula.

``x_k = N^{-1} \\sum_{n=1}^N{\\bar{x}_n g^{-(n-1)(k-1)} } \\mod q``

The same input parameters constraints as for `ntt` function must be applied
"""
function intt(y::Array{T,1}, g::T, q::T) where {T<:Integer}
    N = length(y)

    validate(N, g, q)

    inv_g = invmod(g, q)
    inv_N = invmod(N, q)
    t = [powermod(inv_g, l * k, q) for k in 0:N-1, l in 0:N-1]
    return mod.(inv_N * t * y, q)
end

function intt(y::Array{T,2}, g::T, q::T) where {T<:Integer}
    N, M = size(y)
    x = zeros(T, size(y))

    for m in 1:M
        x[:, m] = intt(y[:, m], g, q)
    end

    for n in 1:N
        x[n, :] = intt(x[n, :], g, q)
    end

    return x
end

# those implementations are incorrect but transform results were interesing anyway
# TODO: check if those calculations are useful

# function ntt(g::T, q::T, x::Array{T,N}) where {T <: Integer,N}
#     return reshape(ntt(g, q, reshape(x, length(x))), size(x))
# end

# function intt(g::T, q::T, y::Array{T,N}) where {T <: Integer,N}
#     return reshape(intt(g, q, reshape(y, length(y))), size(y))
# end


text:
using HMMGradients, LinearAlgebra, Random, Statistics, Distributions, Zygote
use_log = true # do computation in log domain
plot_stuff = true

# define Gaussian probability distribution
logdist(μ::T,σ::T,x::T) where {T} = -( ((x-μ)/σ)^2/2 + log(σ*sqrt(T(2*π))))
dist(μ::T,σ::T,x::T) where {T} = exp(-((x-μ)/σ)^2 / 2 )/(σ*sqrt(T(2*π)))

# model parameters
mutable struct HMMGaussian{T<:AbstractFloat}
  a::Vector{T}
  A::Matrix{T}
  μs::Vector{T}
  σ::T
  use_log::Bool
  likelihood::Function
end

HMMGaussian(a,A,μs,σ;use_log=false) = use_log ?
  HMMGaussian(log.(a),log.(A),μs,σ,use_log,nlogMLlog) : HMMGaussian(a,A,μs,σ,use_log,nlogML)

# calculate observation probability
function predict(model::HMMGaussian{T},x::AbstractArray{T}) where {T}
  #hcat([dist.(model.μs,model.σs,x[t]) for t in 1:length(x)]...)
  X = x*ones(T,length(model.μs))'
  M = ones(T,length(x))*model.μs'
  if model.use_log
    return logdist.(M,model.σ,X)
  else
    return dist.(M,model.σ,X)
  end
end

# obtain a random signal from HMM model
function Base.rand(model::HMMGaussian{T},Nt::Int) where {T}
  if model.use_log # go back to probability
    model.a .= exp.(model.a)
    model.A .= exp.(model.A)
  end
  x = zeros(T,Nt)
  states = zeros(Int,Nt)
  s = rand(Categorical(a))
  for i = 1:Nt
    states[i] = s
    x[i] = model.μs[s] + model.σ * randn(T)
    s = rand(Categorical(A[s,:])) # new state
  end
  if model.use_log # go back to log
    model.a .= log.(model.a)
    model.A .= log.(model.A)
  end
  return x,states
end

# define maximum likelihood function
loss(model::HMMGaussian,x) = model.likelihood(length(x),model.a,model.A,predict(model,x))

T=Float32
Random.seed!(123)
μs = T[0.0;0.1;0.7]
σ = T(1e-1)         # standard deviation
Ns = 3              # number of states
a = T[1.,0.,0.] # initial probability distribution
A = T[0.5 0.9 0.0; 0.0 0.2 0.8; 0.5 0.0 0.5] # transition matrix
#a = rand(T,Ns)
#A = rand(T,Ns,Ns)
for s in 1:Ns A[s,:] ./= norm(A[s,:],1) end # normalize
a ./= norm(a,1) # normalize
@assert HMMGradients.is_transition_matrix(A)
@assert sum(a)-1.0 < 1e-5

model_gt = HMMGaussian(a,A,μs,σ;use_log=use_log)       # ground truth model
model = HMMGaussian(a,A,zeros(T,Ns),one(T);use_log=use_log) # initial model

Nt = 10 # number of samples
x, s = rand(model_gt,Nt)
fy_gt = loss(model_gt,x)
fy = loss(model,x)
println("cost ground truth model: $fy_gt")
println("cost initial model     : $fy   ")

# gradient descent
function update!(model::HMMGaussian{T},grads,γ=T(1e-3)) where {T}
  grads = grads[1][]
  model.μs .-= γ .* grads.μs
  model.σ -= γ * grads.σ
end

# generrate data
N = 200 # number od samples
X = [rand(model_gt,Nt)[1] for n=1:N]

epochs=10
for e in 1:epochs
  cost = 0
  for i = 1:N
    grads = gradient(model) do m
      f = loss(m,X[i])
      cost += f
      return f
    end
    update!(model,grads)
  end
  println("epoch: $e cost: $(cost/N)")
end

println("ground truth μs: $(model_gt.μs) ")
println("learned μs     : $(model.μs) ")
println("error μs       : $(norm(model.μs-model_gt.μs)) ")
println("ground truth σ : $(model_gt.σ) ")
println("learned σ      : $(model.σ) ")
println("error σ        : $(norm(model.σ-model_gt.σ)) ")

if plot_stuff
  # generate random samples from ground truth HMM
  # plots a random signal x, the state sequence and observation probability y
  Nt = 50
  x_gt,s_gt = rand(model_gt,Nt)
  y_gt = predict(model_gt,x_gt)
  x,s = rand(model,Nt)
  y = predict(model,x)
  using Plots
  p1 = plot(x_gt, title="ground truth HMM", label="x")
  q1 = plot(s_gt, marker=(:star), title="states", label="s")
  r1 = heatmap(y_gt', title="y")
  p2 = plot(x, title="learned HMM", label="x")
  q2 = plot(s, marker=(:star), label="s")
  r2 = heatmap(y')
  plot(p1,q1,r1,p2,q2,r2)
end


text:
# # A Stochastic Splicing Simulator

# ## Overview
# Slicing is a process by which subsequences of preRNA (called introns, by definition) are removed.
# The remaining subsequences (called exons) are rejoined together.
# The resulting RNA sequence with introns removed is called mature RNA (mRNA).

# The splicing process is often stochastic.  Given a preRNA sequence, there is a
# distribution over corresponding mRNA sequences.
# As is almost always the case, the stochasticity results from incomplete information.
# Sources of uncertainty:
# - Structure: RNA modifications / bulged nucleotides make pairing of rNA and snRNPs variable.
# - Mutation: Single nucleotide variation (Q: At what point do mutations happen?)
# - Spliceosome composition: Variability in RNA-Binding Protein concentrations of cell

# ## Model

# At a high level, this is an abstract biophysical model that simulates splicing using a discrete-time simulation.

# __Trans-elements__ are small nucleus RNA, proteins and their combination called small- nuclear
# ribonuclear proteins (snRNP).  Trans elements perform splicing by binding to the RNA
# in a sytematic fashion.

using BioSequences, Random

"A Trans-Element / Factor: Proten or snRNP that can bind to primanry RNA transcript"
abstract type TransElement end

# There are many properties of a trans-element that may affect its role in splicing.
# The simplest representation associates a transelement only with a name:

"Symbolic TransElement -- represents only its name"
struct OpaqueTransElement <: TransElement
  name::Symbol
end

# A __spliceosome__ is a complex of transelements.
# We model a spliceosome as a partial function from trans-elements to real values representing their
# relative concentrations:

"Collection of transelements at different concentrations"
struct Spliceosome{T <: TransElement, R <: Real}
  elems::Dict{T, R}   # Maps from transelement to its concentration in spliceosome
  # FIXME: Ensure concentratios are normalized
end

"Is transelement `te` in the splicesosome `sos`?"
insos(te, sos::Spliceosome) = te in keys(sos.elems)
# Spliceosome(elems::Dict{T, R}) where {T, R} = Spliceosome{T, R}(elems)

"concentration of transelement `te` in spliceosome `so`"
concentration(sos, te) = sos.elems[te]

"Number of trans-elements in splicesosome"
ntranselems(sos::Spliceosome) = length(sos.elems)

# The __primary transcript__ is a sequence of RNA to be spliced by the spliceosome
const RNASequence = BioSequence{RNAAlphabet{4}}

# Trans-elements can bind to sequences in primary transcript.
# A __bound sequence__ or transcript is a primary transcript with any number of trans-elements bound to it.

"Position -- Natural number 1:length(seq)"
const Pos = Int

"Primary transcript with tran-elements bound to it"
struct BoundSeq{T <: TransElement}
  seq::RNASequence             # Primary transcriipt
  binds::Vector{Tuple{Pos, T}} # position -> transelement bound at that position
end

BoundSeq(seq::RNASequence, ::Type{T} = TransElement) where T =
  BoundSeq(seq, Tuple{Pos, T}[])

# ### Interactions
# Splicing literature makes a distinction between the core splicing machinery
# and __splicing regulary elements__.
# The core splicing machinery is carried out by the __major splicesosome__, which
# consists of exactly 5 rnSNAs and around 150 proteins, which combine to make
# 5 rnSNPs: U1, U2, U4, U5, U6.
# The core splicing mechanism is an ordered, stepwise assembly of discrete snRNP particles
# on the pre-mRNA primary transcript.  First U1 binds to the 5' splice site, U2 binds
# to a branch site, and so on until a split occurs.

# Splicing regulartory elements (SREs) are patterns on RNA that inhibit or promote splicing machinery.
# They are believed to be the main mechanism for alternative splcing.

# We will model both the core machinery and SREs within the same framework.
# We model Interactions as tuples of the form $(p, e)$, where $p$ is a __pattern__ (described below)
# and $e$ is a set of __effects__:

struct Interaction{P, TE, F}
  pattern::P
  te::TE
  loc::Pos
  factor::F
end

# #### Patterns

# A pattern (or motif) is a property of an RNA sequence.
# Formally a pattern is a set of strings, i.e., a language.
# For instance, a g-run is the set of strings $\{G, GG, GGG, GGGG, GGGGG, \dots\}$.
# A pattern $p$ __matches__ a string $s$ if any element of $p$ is a subsequence of $s$.

# We assume that the extent to which a pattern has occured is graded.
# A __pattern weight function__ maps patterns to a non-negative real number, denoting the degree to which it has matched.
# p(seq) = 1 indicates a perfect match and p(seq) = 0 indicates no match at all.

# #### EFfects

# An effect $E \subseteq \textrm{TransElements} \times \textrm{Location} \times \mathbb{R}$ is a ternary relation:
# a set of tuples of the form $(t, l, f)$, where $t$ is a transelement, $l$ is an integer valued location, and $f$ is a factor.
# Informally, the semantics of an interaction $(p, \{(t_1, l_1, f_1), (t_2, l_2, f_2), \dots\})$.
# are that if a pattern $p$ is found, then the probability that transelement $t_i$ will bind to location $l_i$ is increated (or decreased) multiplicatively by factor $f_i$.

"Returns a protein extractor that returns the singleton `te` iff its in the splicesosome"
only(te::T) where T = sos -> insos(te, sos) ? [te] : T[]

"Initial probabilities for (TransElement, Position)"
initT(sos::Spliceosome{T, R}, seq) where {T, R} = Dict{Tuple{T, Pos}, R}()

"Absolute position where pattern matched"
matchpos(m::BioSequences.RE.RegexMatch) = m.captured[1]

"Simulate one binding step"
function step(rng, bseq, sos, intrs)
  T = initT(sos, bseq.seq)

  # For each interaction
  for intr in intrs
    # and for each match of the interaction's pattern
    for x in eachmatch(intr.pattern, bseq.seq)
      # For each
      abspos = matchpos(x) + intr.loc
      for protein in intr.te(sos)
        if (protein, abspos) ∉ keys(T)
          T[protein, abspos] = 1.0
        end
        T[protein, abspos] *= intr.factor
      end
    end
  end
  bseq = bind(rng, T, bseq, sos)
  maybesplice(bseq)
end

"Simulate a binding process"
function bind(rng, T, bseq, sos)
  @show T
  # Select transelement with probability proportional to concentration
  # bseq

  # How
end

"Simulate core splicing machinery -- checks for presence of u1 and u2"
maybesplice(bseq) = bseq

# The splicosome is then simulated by simply repeating this process until some criterion of convergence is me
"Simulate the splicing"
function splice(rng, bseq::BoundSeq, sos, intrs; n = 1000)
  for i = 1:n
    bseq = step(rng, bseq, sos, intrs)
    bseq = maybesplice(bseq)
  end
  bseq
end

splice(bseq, sos, intrs) = splice(Random.GLOBAL_RNG, bseq, sos, intrs)
splice(rng, rnaseq::RNASequence, sos, intrs) = splice(rng, BoundSeq(rnaseq), sos, intrs)

# Model
# A splicing round occurs when U1 and U2AF proteins reach a global minimum

# # Measurement (Sequencing)

# Unfortunately, we do not have exact data, and hence we must model the stochasticity in the sequencing process.
# The following model of sequencing assumes that from a given `rna`, `n` (possibly differing, due to alternative splicing) mRNAs
# are produced, and the sequencing measures a single sunsequence (of random length) from each

"Sequence RNA process"
function sequence(rng, seq, n)
  sequences = RNASequence[]

  splice_start(rng) = uniform(rng, 1:length(seq))

  splce_length(rng) = uniform(rng, 1:100)

  for i = 1:n
    mRNA = simsplice(rng, seq, n)
    lb = splice_start(rng)
    ub = lb + splice_length(rng)
    push!(sequences, mRNA[lb:ub]) # FIXME: Might go out of bounds
  end
  sequences
end


text:
Pkg.update()
Pkg.add("StatsBase")
using StatsBase
using RDatasets
iris_dataframe = dataset("datasets", "iris")
sample(iris_dataframe[:SepalLength], 5)

max_stars_repo_path: Chapter 3/1.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 3/1.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 3/1.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
# function BB_BB_intersect(tt::TT_Cache, a::AABB, b::AABB)
#     t = tt.R_a_b * b.c + tt.t_a_b - a.c
#     return BB_BB_intersect(a.e, b.e, t, tt.R_a_b, tt.abs_R_a_b)
# end

function BB_BB_intersect(tt::TT_Cache, a::OBB, b::OBB)
    i_dh_a = basic_dh(a.R', -a.R' * a.c)   # transform from frame A to   OBB A
    dh_a_b = basic_dh(tt.R_a_b, tt.t_a_b)  # transform from frame B to frame A
    dh_b   = basic_dh(b.R, b.c)            # transform from   OBB B to frame B

    dh_final = i_dh_a * dh_a_b * dh_b
    R_tot, t = dh_R_t(dh_final)

    abs_R_tot = abs.(R_tot) .+ 1.0e-14
    return BB_BB_intersect(a.e, b.e, t, R_tot, abs_R_tot)
end

s221(r::SVector{3, Float64}) = SVector{3, Float64}(r[3], r[3], r[2])
s100(r::SVector{3, Float64}) = SVector{3, Float64}(r[2], r[1], r[1])

function BB_BB_intersect(e_a::SVector{3, Float64}, e_b::SVector{3, Float64}, t::SVector{3, Float64}, R::SMatrix{3,3,Float64,9}, abs_R::SMatrix{3,3,Float64,9})
    # This implementation of the OBB_OBB_SAT test is based off of the test described
    # in table 4.1 on page 103 of Ericson's Real-Time Collision Detection.

    R0_012 = SVector{3, Float64}(R[1], R[4], R[7])
    R1_012 = SVector{3, Float64}(R[2], R[5], R[8])
    R2_012 = SVector{3, Float64}(R[3], R[6], R[9])
    abs_R0_012 = SVector{3, Float64}(abs_R[1], abs_R[4], abs_R[7])
    abs_R1_012 = SVector{3, Float64}(abs_R[2], abs_R[5], abs_R[8])
    abs_R2_012 = SVector{3, Float64}(abs_R[3], abs_R[6], abs_R[9])

    # face test 1/2
    T_dot_L = abs.(t)
    r_a = e_a
    r_b = abs_R * e_b
    any((r_a + r_b) .< T_dot_L) && (return false)

    # face test 2/2
    T_dot_L = abs.(transpose(R) * t)  # these are transposes
    r_a = transpose(abs_R) * e_a
    r_b = e_b
    any((r_a + r_b) .< T_dot_L) && (return false)

    eb_100     = s100(e_b)
    eb_221     = s221(e_b)
    abs_R0_221 = s221(abs_R0_012)
    abs_R0_100 = s100(abs_R0_012)
    abs_R1_221 = s221(abs_R1_012)
    abs_R1_100 = s100(abs_R1_012)
    abs_R2_221 = s221(abs_R2_012)
    abs_R2_100 = s100(abs_R2_012)

    t0 = t[1]
    t1 = t[2]
    t2 = t[3]
    ea0 = e_a[1]
    ea1 = e_a[2]
    ea2 = e_a[3]

    # cross test 1/3
    T_dot_L = abs.(t2      * R1_012     - t1      * R2_012)
    r_a     =      ea1     * abs_R2_012 + ea2     * abs_R1_012
    r_b     =      eb_100 .* abs_R0_221 + eb_221 .* abs_R0_100
    any((r_a + r_b) .< T_dot_L) && (return false)

    # cross test 2/3
    T_dot_L = abs.(t0  * R2_012     -     t2  * R0_012)
    r_a     =      ea0 * abs_R2_012 +    ea2  * abs_R0_012
    r_b     =  eb_100 .* abs_R1_221 + eb_221 .* abs_R1_100
    any((r_a + r_b) .< T_dot_L) && (return false)

    # cross test 3/3
    T_dot_L =  abs.(t1 * R0_012     -     t0  * R1_012)
    r_a     =      ea0 * abs_R1_012 +    ea1  * abs_R0_012
    r_b     =  eb_100 .* abs_R2_221 + eb_221 .* abs_R2_100
    any((r_a + r_b) .< T_dot_L) && (return false)
    return true
end


text:
function initialize!(integrator, cache::BS3ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

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

@muladd function perform_step!(integrator, cache::BS3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4 = cache
  k1 = integrator.fsalfirst
  a1 = dt*a21
  k2 = f(uprev+a1*k1, p, t+c1*dt)
  a2 = dt*a32
  k3 = f(uprev+a2*k2, p, t+c2*dt)
  u = uprev+dt*(a41*k1+a42*k2+a43*k3)
  k4 = f(u, p, t+dt); integrator.fsallast = k4
  integrator.destats.nf += 3
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::BS3Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k4
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::BS3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3,k4,utilde,tmp,atmp = cache
  @unpack a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4 = cache.tab
  # k1 = cache.fsalfirst
  k1 = integrator.fsalfirst
  a1 = dt*a21
  @.. tmp = uprev+a1*k1
  f(k2, tmp, p, t+c1*dt)
  a2 = dt*a32
  @.. tmp = uprev+a2*k2
  f(k3, tmp, p, t+c2*dt)
  @.. u = uprev+dt*(a41*k1+a42*k2+a43*k3)
  f(k4, u, p, t+dt)
  integrator.destats.nf += 3
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

function initialize!(integrator, cache::OwrenZen3ConstantCache)
    integrator.kshortsize = 4
    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
    integrator.destats.nf += 1

    # Avoid undefined entries if k is an array of arrays
    integrator.fsallast = zero(integrator.fsalfirst)
    integrator.k[1] = integrator.fsalfirst
    @inbounds for i in 2:integrator.kshortsize-1
      integrator.k[i] = zero(integrator.fsalfirst)
    end
    integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::OwrenZen3ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3 = cache
  k1 = integrator.fsalfirst
  a1 = dt*a21
  k2 = f(uprev+a1*k1, p, t+c1*dt)
  tmp =  uprev+ dt*(a31*k1 + a32*k2)
  k3 = f(tmp, p, t+c2*dt)
  u =  uprev+dt*(a41*k1+a42*k2+a43*k3)
  k4 = f(u, p, t+dt); integrator.fsallast = k4
  integrator.destats.nf += 3
  if integrator.opts.adaptive
    utilde =  dt*(btilde1*k1 + btilde2*k2 + btilde3*k3)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4
  integrator.u = u
end

function initialize!(integrator, cache::OwrenZen3Cache)
    integrator.kshortsize = 4
    resize!(integrator.k, integrator.kshortsize)
    integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
    integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
    integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k4  # setup pointers
    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
    integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::OwrenZen3Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k1,k2,k3,k4,utilde,tmp,atmp = cache
  @unpack a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3 = cache.tab
  a1 = dt*a21
  @.. tmp = uprev+a1*k1
  f(k2, tmp, p, t+c1*dt)
  @.. tmp = uprev+dt*(a31*k1+a32*k2)
  f(k3, tmp, p, t+c2*dt)
  @.. u = uprev+dt*(a41*k1+a42*k2+a43*k3)
  f(k4, u, p, t+dt)
  integrator.destats.nf += 3
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

function initialize!(integrator, cache::OwrenZen4ConstantCache)
  integrator.kshortsize = 6
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::OwrenZen4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t+c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
  k4 = f(uprev+dt*(a41*k1+a42*k2+a43*k3), p, t+c3*dt)
  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
  u = uprev+dt*(a61*k1+a63*k3+a64*k4+a65*k5)
  k6 = f(u, p, t+dt); integrator.fsallast = k6
  integrator.destats.nf += 5
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4;
  integrator.k[5]=k5; integrator.k[6]=k6
  integrator.u = u
end

function initialize!(integrator,cache::OwrenZen4Cache)
  integrator.kshortsize = 6
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k6  # setup pointers
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::OwrenZen4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k1,k2,k3,k4,k5,k6,utilde,tmp,atmp = cache
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5 = cache.tab
  a = dt*a21
  @.. tmp = uprev+a*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+a42*k2+a43*k3)
  f(k4, tmp, p, t+c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
  f(k5, tmp, p, t+c4*dt)
  @.. u = uprev+dt*(a61*k1+a63*k3+a64*k4+a65*k5)
  f(k6, u, p, t+dt)
  integrator.destats.nf += 5
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  return nothing
end

#=
@muladd function perform_step!(integrator, cache::OwrenZen4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k1,k2,k3,k4,k5,k6,utilde,tmp,atmp = cache
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5 = cache.tab
  a = dt*a21
  @tight_loop_macros for i in uidx
    @inbounds tmp[i] = uprev[i]+a*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]+a42*k2[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]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
  f(k5, tmp, p, t+c4*dt)
  @tight_loop_macros for i in uidx
    @inbounds u[i] = uprev[i]+dt*(a61*k1[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(k6, u, p, t+dt)
  integrator.destats.nf += 5
  if integrator.opts.adaptive
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end
=#

function initialize!(integrator, cache::OwrenZen5ConstantCache)
  integrator.kshortsize = 8
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::OwrenZen5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a31,a32,a41,a42,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t+c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
  k4 = f(uprev+dt*(a41*k1+a42*k2+k3), p, t+c3*dt)
  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
  k6 = f(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5), p, t+c5*dt)
  k7 = f(uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6), p, t+c6*dt)
  u = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  k8 = f(u, p, t+dt); integrator.fsallast = k8
  integrator.destats.nf += 7
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4
  integrator.k[5]=k5; integrator.k[6]=k6; integrator.k[7]=k7; integrator.k[8]=k8
  integrator.u = u
end

function initialize!(integrator, cache::OwrenZen5Cache)
  integrator.kshortsize = 8
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;
  integrator.k[7]=cache.k7; integrator.k[8]=cache.k8;
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k8  # setup pointers
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::OwrenZen5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp = cache
  @unpack a21,a31,a32,a41,a42,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  a = dt*a21
  @.. tmp = uprev+a*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+a42*k2+k3)
  f(k4, tmp, p, t+c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
  f(k5, tmp, p, t+c4*dt)
  @.. tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  f(k6, tmp, p, t+c5*dt)
  @.. tmp = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
  f(k7, tmp, p, t+c6*dt)
  @.. u = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  f(k8, u, p, t+dt)
  integrator.destats.nf += 7
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  return nothing
end

#=
@muladd function perform_step!(integrator, cache::OwrenZen5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp = cache
  @unpack a21,a31,a32,a41,a42,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  a = dt*a21
  @tight_loop_macros for i in uidx
    @inbounds tmp[i] = uprev[i]+a*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]+a42*k2[i]+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]+a52*k2[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]+a62*k2[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]+a72*k2[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 u[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, u, p, t+dt)
  integrator.destats.nf += 7
  if integrator.opts.adaptive
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end
=#

function initialize!(integrator, cache::BS5ConstantCache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 8) : (integrator.kshortsize = 11)
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:7
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast

  if !alg.lazy
    @inbounds for i in 9:11
      integrator.k[i] = zero(integrator.fsalfirst)
    end
  end
end

@muladd function perform_step!(integrator, cache::BS5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t+c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
  k4 = f(uprev+dt*(a41*k1+a42*k2+a43*k3), p, t+c3*dt)
  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
  k6 = f(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5), p, t+c5*dt)
  k7 = f(uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6), p, t+dt)
  integrator.destats.nf += 6
  u = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  integrator.fsallast = f(u, p, t+dt); k8 = integrator.fsallast
  integrator.destats.nf += 1
  if integrator.opts.adaptive
    uhat   = dt*(bhat1*k1 + bhat3*k3 + bhat4*k4 + bhat5*k5 + bhat6*k6)
    utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8)
    atmp = calculate_residuals(uhat, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst1 = integrator.opts.internalnorm(atmp,t)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst2 = integrator.opts.internalnorm(atmp,t)
    integrator.EEst = max(EEst1,EEst2)
  end
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3;integrator.k[4]=k4;integrator.k[5]=k5;integrator.k[6]=k6;integrator.k[7]=k7;integrator.k[8]=k8
  integrator.u = u

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || accept_step_controller(integrator, integrator.opts.controller))
    @unpack c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = cache
    k = integrator.k
    k[9] = f(uprev+dt*(a91*k[1]+a92*k[2]+a93*k[3]+a94*k[4]+a95*k[5]+a96*k[6]+a97*k[7]+a98*k[8]),p,t+c6*dt)
    k[10] = f(uprev+dt*(a101*k[1]+a102*k[2]+a103*k[3]+a104*k[4]+a105*k[5]+a106*k[6]+a107*k[7]+a108*k[8]+a109*k[9]),p,t+c7*dt)
    k[11] = f(uprev+dt*(a111*k[1]+a112*k[2]+a113*k[3]+a114*k[4]+a115*k[5]+a116*k[6]+a117*k[7]+a118*k[8]+a119*k[9]+a1110*k[10]),p,t+c8*dt)
    integrator.destats.nf += 3
  end
end

function initialize!(integrator, cache::BS5Cache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 8) : (integrator.kshortsize = 11)
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;
  integrator.k[7]=cache.k7; integrator.k[8]=cache.k8

  if !alg.lazy
    integrator.k[9]= similar(cache.k1)
    integrator.k[10]= similar(cache.k1)
    integrator.k[11]= similar(cache.k1)
  end

  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k8  # setup pointers
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::BS5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp = cache
  @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
  a = dt*a21
  @.. tmp = uprev+a*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+a42*k2+a43*k3)
  f(k4, tmp, p, t+c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
  f(k5, tmp, p, t+c4*dt)
  @.. tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  f(k6, tmp, p, t+c5*dt)
  @.. tmp = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
  f(k7, tmp, p, t+dt)
  @.. u = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  f(k8, u, p, t+dt)
  integrator.destats.nf += 7
  if integrator.opts.adaptive
    @.. utilde = dt*(bhat1*k1 + bhat3*k3 + bhat4*k4 + bhat5*k5 + bhat6*k6)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst1 = integrator.opts.internalnorm(atmp,t)
    @.. utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst2 = integrator.opts.internalnorm(atmp,t)
    integrator.EEst = max(EEst1,EEst2)
  end
  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || accept_step_controller(integrator, integrator.opts.controller))
    k = integrator.k
    @unpack c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = cache.tab
    @.. tmp = uprev+dt*(a91*k[1]+a92*k[2]+a93*k[3]+a94*k[4]+a95*k[5]+a96*k[6]+a97*k[7]+a98*k[8])
    f(k[9],tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a101*k[1]+a102*k[2]+a103*k[3]+a104*k[4]+a105*k[5]+a106*k[6]+a107*k[7]+a108*k[8]+a109*k[9])
    f(k[10],tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a111*k[1]+a112*k[2]+a113*k[3]+a114*k[4]+a115*k[5]+a116*k[6]+a117*k[7]+a118*k[8]+a119*k[9]+a1110*k[10])
    f(k[11],tmp,p,t+c8*dt)
    integrator.destats.nf += 3
  end
  return nothing
end

#=
@muladd function perform_step!(integrator, cache::BS5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp = cache
  @unpack c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8 = cache.tab
  a = dt*a21
  @tight_loop_macros for i in uidx
    @inbounds tmp[i] = uprev[i]+a*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]+a42*k2[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]+a52*k2[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]+a62*k2[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]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
  end
  f(k7, tmp, p, t+dt)
  @tight_loop_macros for i in uidx
    @inbounds u[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, u, p, t+dt)
  integrator.destats.nf += 7
  if integrator.opts.adaptive
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(bhat1*k1[i] + bhat3*k3[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst1 = integrator.opts.internalnorm(atmp,t)
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i] + btilde8*k8[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    EEst2 = integrator.opts.internalnorm(atmp,t)
    integrator.EEst = max(EEst1,EEst2)
  end

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || accept_step_controller(integrator, integrator.opts.controller))
    k = integrator.k
    @unpack c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = cache.tab
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a91*k[1][i]+a92*k[2][i]+a93*k[3][i]+a94*k[4][i]+a95*k[5][i]+a96*k[6][i]+a97*k[7][i]+a98*k[8][i])
    end
    f(k[9],tmp,p,t+c6*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a101*k[1][i]+a102*k[2][i]+a103*k[3][i]+a104*k[4][i]+a105*k[5][i]+a106*k[6][i]+a107*k[7][i]+a108*k[8][i]+a109*k[9][i])
    end
    f(k[10],tmp,p,t+c7*dt)
    @tight_loop_macros for i in uidx
      @inbounds tmp[i] = uprev[i]+dt*(a111*k[1][i]+a112*k[2][i]+a113*k[3][i]+a114*k[4][i]+a115*k[5][i]+a116*k[6][i]+a117*k[7][i]+a118*k[8][i]+a119*k[9][i]+a1110*k[10][i])
    end
    f(k[11],tmp,p,t+c8*dt)
    integrator.destats.nf += 3
  end
end
=#

function initialize!(integrator, cache::Tsit5ConstantCache)
  integrator.kshortsize = 7
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::Tsit5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t+c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
  k4 = f(uprev+dt*(a41*k1+a42*k2+a43*k3), p, t+c3*dt)
  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
  g6 = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  k6 = f(g6, p, t+dt)
  u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
  integrator.fsallast = f(u, p, t+dt); k7 = integrator.fsallast
  integrator.destats.nf += 6
  if typeof(integrator.alg) <: CompositeAlgorithm
    g7 = u
    # Hairer II, page 22
    integrator.eigen_est = integrator.opts.internalnorm(k7 - k6,t)/integrator.opts.internalnorm(g7 - g6,t)
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = k1
  integrator.k[2] = k2
  integrator.k[3] = k3
  integrator.k[4] = k4
  integrator.k[5] = k5
  integrator.k[6] = k6
  integrator.k[7] = k7
  integrator.u = u
end

function initialize!(integrator, cache::Tsit5Cache)
  integrator.kshortsize = 7
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7 # setup pointers
  resize!(integrator.k, integrator.kshortsize)
  # Setup k pointers
  integrator.k[1] = cache.k1
  integrator.k[2] = cache.k2
  integrator.k[3] = cache.k3
  integrator.k[4] = cache.k4
  integrator.k[5] = cache.k5
  integrator.k[6] = cache.k6
  integrator.k[7] = cache.k7
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::Tsit5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp = cache
  a = dt*a21
  @.. tmp = uprev+a*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+a42*k2+a43*k3)
  f(k4, tmp, p, t+c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
  f(k5, tmp, p, t+c4*dt)
  @.. tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  f(k6, tmp, p, t+dt)
  @.. u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
  f(k7, u, p, t+dt)
  integrator.destats.nf += 6
  if integrator.alg isa CompositeAlgorithm
    g7 = u
    g6 = tmp
    # Hairer II, page 22
    @.. utilde = k7 - k6
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g7 - g6
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  return nothing
end

#=
@muladd function perform_step!(integrator, cache::Tsit5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp = cache
  a = dt*a21
  @tight_loop_macros for i in uidx
    @inbounds tmp[i] = uprev[i]+a*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]+a42*k2[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]+a52*k2[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]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(k6, tmp, p, t+dt)
  @tight_loop_macros for i in uidx
    @inbounds u[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
  end
  f(k7, u, p, t+dt)
  integrator.destats.nf += 6
  if typeof(integrator.alg) <: CompositeAlgorithm
    g7 = u
    g6 = tmp
    # Hairer II, page 22
    @.. utilde = k7 - k6
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g7 - g6
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end
=#

function initialize!(integrator, cache::DP5ConstantCache)
  integrator.kshortsize = 4
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  @inbounds for i in eachindex(integrator.k)
    integrator.k[i] = zero(integrator.fsalfirst)
  end
end

@muladd function perform_step!(integrator, cache::DP5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6 = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t+c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
  k4 = f(uprev+dt*(a41*k1+a42*k2+a43*k3), p, t+c3*dt)
  k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
  g6 = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  k6 = f(g6, p, t+dt)
  update = a71*k1+a73*k3+a74*k4+a75*k5+a76*k6
  u = uprev+dt*update
  integrator.fsallast = f(u, p, t+dt); k7 = integrator.fsallast
  integrator.destats.nf += 6
  if typeof(integrator.alg) <: CompositeAlgorithm
    g7 = u
    # Hairer II, page 22
    integrator.eigen_est = integrator.opts.internalnorm(k7 - k6,t)/integrator.opts.internalnorm(g7 - g6,t)
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = update
  bspl = k1 - update
  integrator.k[2] = bspl
  integrator.k[3] = update - k7 - bspl
  integrator.k[4] = d1*k1+d3*k3+d4*k4+d5*k5+d6*k6+d7*k7
  integrator.u = u
end

function initialize!(integrator, cache::DP5Cache)
  integrator.kshortsize = 4
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = cache.update
  integrator.k[2] = cache.bspl
  integrator.k[3] = cache.dense_tmp3
  integrator.k[4] = cache.dense_tmp4
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::DP5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache.tab
  a = dt*a21
  @.. tmp = uprev+a*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+a42*k2+a43*k3)
  f(k4, tmp, p, t+c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
  f(k5, tmp, p, t+c4*dt)
  @.. tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
  f(k6, tmp, p, t+dt)
  @.. update = a71*k1+a73*k3+a74*k4+a75*k5+a76*k6
  @.. u = uprev+dt*update
  f(k7, u, p, t+dt)
  integrator.destats.nf += 6
  if integrator.alg isa CompositeAlgorithm
    g6 = tmp
    g7 = u
    # Hairer II, page 22
    ϱu, ϱd = zero(eltype(k7))^2, zero(eltype(g7))^2
    @inbounds for i in eachindex(k7)
      ϱu += (k7[i] - k6[i])^2
      ϱd += (g7[i] - g6[i])^2
    end
    integrator.eigen_est = sqrt(ϱu/ϱd)*oneunit(t)
  end
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  if integrator.opts.calck
    #integrator.k[4] == k5
    @.. integrator.k[4] = d1*k1+d3*k3+d4*k4+d5*k5+d6*k6+d7*k7
    #bspl == k3
    @.. bspl = k1 - update
    # k6 === integrator.k[3] === k2
    @.. integrator.k[3] = update - k7 - bspl
  end
  return nothing
end

#=
@muladd function perform_step!(integrator, cache::DP5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache.tab
  a = dt*a21
  @tight_loop_macros for i in uidx
    @inbounds tmp[i] = uprev[i]+a*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]+a42*k2[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]+a52*k2[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]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
  f(k6, tmp, p, t+dt)
  @tight_loop_macros for i in uidx
    @inbounds update[i] = a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
    @inbounds u[i] = uprev[i]+dt*update[i]
  end
  f(k7, u, p, t+dt)
  integrator.destats.nf += 6
  if typeof(integrator.alg) <: CompositeAlgorithm
    g6 = tmp
    g7 = u
    # Hairer II, page 22
    ϱu, ϱd = zero(eltype(k7))^2, zero(eltype(g7))^2
    @inbounds for i in eachindex(k7)
      ϱu += (k7[i] - k6[i])^2
      ϱd += (g7[i] - g6[i])^2
    end
    integrator.eigen_est = sqrt(ϱu/ϱd)*oneunit(t)
  end
  if integrator.opts.adaptive
    @tight_loop_macros for i in uidx
      @inbounds utilde[i] = dt*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i])
    end
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  if integrator.opts.calck
    @tight_loop_macros for i in uidx
      #integrator.k[4] == k5
      @inbounds integrator.k[4][i] = d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i]
      #bspl == k3
      @inbounds bspl[i] = k1[i] - update[i]
      # k6 === integrator.k[3] === k2
      @inbounds integrator.k[3][i] = update[i] - k7[i] - bspl[i]
    end
  end
end
=#

function initialize!(integrator,cache::KYK2014DGSSPRK_3S2_ConstantCache)
 integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
 integrator.destats.nf += 1
 integrator.kshortsize = 2
 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

 # Avoid undefined entries if k is an array of arrays
 integrator.fsallast = zero(integrator.fsalfirst)
 return nothing
end

@muladd function perform_step!(integrator,cache::KYK2014DGSSPRK_3S2_ConstantCache,repeat_step=false)
 @unpack t,dt,uprev,u,f,p = integrator
 @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache
 u_1 =  α_10 * uprev + dt*β_10*integrator.fsalfirst
 u_2 = (
  α_20 * uprev +
  α_21 * u_1 + dt*β_21*f(u_1, p, t + c_1*dt)
 )
 integrator.u = (
  α_30 * uprev + dt*β_30*integrator.fsalfirst +
  α_32 * u_2 + dt*β_32*f(u_2, p, t + c_2*dt)
 )
 integrator.k[1] = integrator.fsalfirst
 integrator.k[2] = f(integrator.u, p, t+dt) # For interpolation, then FSAL'd
 integrator.destats.nf += 3
 integrator.fsallast = integrator.k[2]
 return nothing
end

function initialize!(integrator,cache::KYK2014DGSSPRK_3S2_Cache)
 @unpack k,fsalfirst = cache
 integrator.fsalfirst = fsalfirst
 integrator.fsallast = k
 integrator.kshortsize = 2
 resize!(integrator.k, integrator.kshortsize)
 integrator.k[1] = integrator.fsalfirst
 integrator.k[2] = integrator.fsallast
 integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
 integrator.destats.nf += 1
 return nothing
end

@muladd function perform_step!(integrator,cache::KYK2014DGSSPRK_3S2_Cache,repeat_step=false)
 @unpack t,dt,uprev,u,f,p = integrator
 @unpack k,fsalfirst,u_1,u_2, kk_1, kk_2 = cache
 @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache.tab

 @.. u_1 =  α_10 * uprev + dt*β_10*integrator.fsalfirst
 f(kk_1, u_1, p, t + c_1*dt)
 @.. u_2 = (
  α_20 * uprev +
  α_21 * u_1 + dt*β_21*kk_1
 )
 f(kk_2, u_2, p, t + c_2*dt)
 @.. integrator.u = (
  α_30 * uprev + dt*β_30*integrator.fsalfirst +
  α_32 * u_2 + dt*β_32*kk_2
 )
 f(integrator.k[2], integrator.u, p, t+dt) # For interpolation, then FSAL'd
 integrator.destats.nf += 3
 return nothing
end

function initialize!(integrator, cache::RKO65ConstantCache)
  integrator.kshortsize = 6
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::RKO65ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p = integrator
  @unpack α21, α31, α41, α51, α32, α42, α52, α62, α43, α53, α63, α54, α64, α65, β2, β3, β4, β5, β6, c1, c2, c3, c4, c5, c6 = cache

  #k1=integrator.fsalfirst #f(uprev,p,t)
  k1 = f(uprev, p, t+c1*dt)
  k2 = f(uprev+α21*dt*k1, p, t+c2*dt)
  k3 = f(uprev+α31*dt*k1+α32*dt*k2, p, t+c3*dt)
  k4 = f(uprev+α41*dt*k1+α42*dt*k2+α43*dt*k3, p, t+c4*dt)
  k5 = f(uprev+α51*dt*k1+α52*dt*k2+α53*dt*k3+α54*dt*k4, p, t+c5*dt)
  k6 = f(uprev+α62*dt*k2+α63*dt*k3+α64*dt*k4+α65*dt*k5, p, t+c6*dt)
  u = uprev+dt*(β2*k2+β3*k3+β4*k4+β5*k5+β6*k6)

  integrator.fsallast = f(u, p, t+dt)  # For interpolation, then FSAL'd

  integrator.destats.nf += 6
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4; integrator.k[5]=k5; integrator.k[6]=k6
  integrator.u = u
end

function initialize!(integrator, cache::RKO65Cache)
  @unpack k,fsalfirst = cache
  integrator.kshortsize = 6
  resize!(integrator.k, integrator.kshortsize)

  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k6  # setup pointers

  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;

  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k6  # setup pointers

  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # Pre-start fsal

  integrator.destats.nf += 1

end

@muladd function perform_step!(integrator, cache::RKO65Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p = integrator
  @unpack tmp, k, k1, k2, k3, k4, k5, k6  = cache
  @unpack α21, α31, α41, α51, α32, α42, α52, α62, α43, α53, α63, α54, α64, α65, β2, β3, β4, β5, β6, c1, c2, c3, c4, c5, c6 = cache.tab
  #println("L221: tmp", tmp)
  f(k1, uprev, p, t+c1*dt)
  @.. tmp = uprev+α21*dt*k1
  #println("L224: tmp/k", tmp, k1)
  f(k2, tmp, p, t+c2*dt)
  @.. tmp = uprev+α31*dt*k1+α32*dt*k2
  f(k3, tmp, p, t+c3*dt)
  @.. tmp = uprev+α41*dt*k1+α42*dt*k2+α43*dt*k3
  f(k4, tmp, p, t+c4*dt)
  @.. tmp = uprev+α51*dt*k1+α52*dt*k2+α53*dt*k3+α54*dt*k4
  f(k5, tmp, p, t+c5*dt)
  @.. tmp = uprev+α62*dt*k2+α63*dt*k3+α64*dt*k4+α65*dt*k5
  f(k6, tmp, p, t+c6*dt)

  @.. u = uprev+dt*(β2*k2+β3*k3+β4*k4+β5*k5+β6*k6)
  #println("L238: tmp/u", tmp, u)
  integrator.destats.nf += 6

  #return nothing
end

function initialize!(integrator, cache::FRK65ConstantCache)
  integrator.kshortsize = 9
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::FRK65ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p = integrator
  @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11 = cache
  alg = unwrap_alg(integrator, true)
  ν = alg.omega*dt
  νsq = ν^2
  β4 = (d1 + νsq*(d2 + νsq*(d3 + νsq*(d4 + νsq*(d5 + νsq*(d6 + +νsq*d7))))))/(1 + νsq*(d8 + νsq*(d9 + νsq*(d10 + νsq*(d11 + νsq*(d12 + +νsq*d13))))))
  β5 = (e1 + νsq*(e2 + νsq*(e3 + νsq*(e4 + νsq*(e5 + νsq*e6)))))/(1 + νsq*(e8 + νsq*(e9 + νsq*(e10 + νsq*e11))))
  β6 = (f1 + νsq*(f2 + νsq*(f3 + νsq*(f4 + νsq*(f5 + νsq*f6)))))/(1 + νsq*(f8 + νsq*(f9 + νsq*(f10 + νsq*f11))))

  k1 = integrator.fsalfirst
  k2 = f(uprev+α21*dt*k1, p, t+c2*dt)
  k3 = f(uprev+α31*dt*k1+α32*dt*k2, p, t+c3*dt)
  k4 = f(uprev+α41*dt*k1+α43*dt*k3, p, t+c4*dt)
  k5 = f(uprev+α51*dt*k1+α53*dt*k3+α54*dt*k4, p, t+c5*dt)
  k6 = f(uprev+α61*dt*k1+α63*dt*k3+α64*dt*k4+α65*dt*k5, p, t+c6*dt)
  k7 = f(uprev+α71*dt*k1+α73*dt*k3+α74*dt*k4+α75*dt*k5+α76*dt*k6, p, t+c7*dt)
  k8 = f(uprev+α81*dt*k1+α83*dt*k3+α84*dt*k4+α85*dt*k5+α86*dt*k6+α87*dt*k7, p, t+c8*dt)
  u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8)
  integrator.fsallast = f(u, p, t+dt)
  k9 = integrator.fsallast
  if integrator.opts.adaptive
    utilde = dt*(β1tilde*k1 + β4tilde*k4 + β5tilde*k5 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.destats.nf += 8
  integrator.k[1]=k1
  integrator.k[2]=k2
  integrator.k[3]=k3
  integrator.k[4]=k4
  integrator.k[5]=k5
  integrator.k[6]=k6
  integrator.k[7]=k7
  integrator.k[8]=k8
  integrator.k[9]=k9
  integrator.u = u
end

function initialize!(integrator, cache::FRK65Cache)
  integrator.kshortsize = 9
  integrator.fsalfirst = cache.k1
  integrator.fsallast = cache.k9

  resize!(integrator.k, integrator.kshortsize)

  integrator.k[1]=cache.k1
  integrator.k[2]=cache.k2
  integrator.k[3]=cache.k3
  integrator.k[4]=cache.k4
  integrator.k[5]=cache.k5
  integrator.k[6]=cache.k6
  integrator.k[7]=cache.k7
  integrator.k[8]=cache.k8
  integrator.k[9]=cache.k9

  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::FRK65Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p = integrator
  @unpack tmp, k1, k2, k3, k4, k5, k6, k7, k8, k9, utilde, atmp  = cache
  @unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11 = cache.tab
  alg = unwrap_alg(integrator, true)

  ν = alg.omega*dt
  νsq = ν^2
  β4 = (d1 + νsq*(d2 + νsq*(d3 + νsq*(d4 + νsq*(d5 + νsq*(d6 + +νsq*d7))))))/(1 + νsq*(d8 + νsq*(d9 + νsq*(d10 + νsq*(d11 + νsq*(d12 + +νsq*d13))))))
  β5 = (e1 + νsq*(e2 + νsq*(e3 + νsq*(e4 + νsq*(e5 + νsq*e6)))))/(1 + νsq*(e8 + νsq*(e9 + νsq*(e10 + νsq*e11))))
  β6 = (f1 + νsq*(f2 + νsq*(f3 + νsq*(f4 + νsq*(f5 + νsq*f6)))))/(1 + νsq*(f8 + νsq*(f9 + νsq*(f10 + νsq*f11))))

  @.. tmp = uprev+α21*dt*k1
  f(k2, tmp, p, t+c2*dt)
  @.. tmp = uprev+α31*dt*k1+α32*dt*k2
  f(k3, tmp, p, t+c3*dt)
  @.. tmp = uprev+α41*dt*k1+α43*dt*k3
  f(k4, tmp, p, t+c4*dt)
  @.. tmp = uprev+α51*dt*k1+α53*dt*k3+α54*dt*k4
  f(k5, tmp, p, t+c5*dt)
  @.. tmp = uprev+α61*dt*k1+α63*dt*k3+α64*dt*k4+α65*dt*k5
  f(k6, tmp, p, t+c6*dt)
  @.. tmp = uprev+α71*dt*k1+α73*dt*k3+α74*dt*k4+α75*dt*k5+α76*dt*k6
  f(k7, tmp, p, t+c7*dt)
  @.. tmp = uprev+α81*dt*k1+α83*dt*k3+α84*dt*k4+α85*dt*k5+α86*dt*k6+α87*dt*k7
  f(k8, tmp, p, t+c8*dt)
  @.. u = uprev+dt*(β1*k1+β4*k4+β5*k5+β6*k6+β7*k7+β8*k8)
  f(k9, u, p, t+dt)

  integrator.destats.nf += 8
  if integrator.opts.adaptive
    @.. utilde = dt*(β1tilde*k1 + β4tilde*k4 + β5tilde*k5 + β6tilde*k6 + β7tilde*k7 + β8tilde*k8 + β9tilde*k9)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  return nothing
end

function initialize!(integrator,cache::RKMConstantCache)
  integrator.kshortsize = 6
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
	integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

	# Avoid undefined entries if k is an array of arrays
	integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::RKMConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α2, α3, α4, α5, α6, β1, β2, β3, β4, β6, c2, c3, c4, c5, c6 = cache

  #k1 = f(uprev, p, t)
  k1 = integrator.fsalfirst
  k2 = f(uprev+α2*dt*k1, p, t+c2*dt)
  k3 = f(uprev+α3*dt*k2, p, t+c3*dt)
  k4 = f(uprev+α4*dt*k3, p, t+c4*dt)
  k5 = f(uprev+α5*dt*k4, p, t+c5*dt)
  k6 = f(uprev+α6*dt*k5, p, t+c6*dt)
  u = uprev+dt*(β1*k1+β2*k2+β3*k3+β4*k4+β6*k6)

  integrator.fsallast = f(u, p, t+dt) #interpolation then FSAL'd
  integrator.destats.nf += 6
  integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3; integrator.k[4]=k4; integrator.k[5]=k5; integrator.k[6]=k6
  # integrator.k[1] = integrator.fsalfirst
  # integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator,cache::RKMCache)
  @unpack k,fsalfirst = cache
  integrator.kshortsize = 6
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1]=cache.k1; integrator.k[2]=cache.k2;
  integrator.k[3]=cache.k3; integrator.k[4]=cache.k4;
  integrator.k[5]=cache.k5; integrator.k[6]=cache.k6;

  integrator.fsalfirst = cache.k1; integrator.fsallast = zero(integrator.fsalfirst)  # setup pointers

  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # Pre-start fsal

  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::RKMCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k,k1,k2,k3,k4,k5,k6 = cache
  @unpack α2,α3,α4,α5,α6,β1,β2,β3,β4,β6,c2,c3,c4,c5,c6 = cache.tab

  @.. tmp = uprev+α2*dt*k1
  f(k2, tmp, p, t+c2*dt)
  @.. tmp = uprev+α3*dt*k2
  f(k3, tmp, p, t+c3*dt)
  @.. tmp = uprev+α4*dt*k3
  f(k4, tmp, p, t+c4*dt)
  @.. tmp = uprev+α5*dt*k4
  f(k5, tmp, p, t+c5*dt)
  @.. tmp = uprev+α6*dt*k5
  f(k6, tmp, p, t+c6*dt)
  @.. u = uprev+dt*(β1*k1+β2*k2+β3*k3+β4*k4+β6*k6)
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 6
  return nothing
end


text:
using Base.Test

# write your own tests here
include("mod1ch2.jl")
include("julia_doc.jl")
include("builtin_functions.jl")
include("data_structures.jl")
include("blas.jl")

@test 2 == 2


text:
using LinearAlgebra: I, nullspace
using UniqueVectors: UniqueVector
using StaticArrays: MMatrix

"""
    parametrize_subspace(M::AbstractMatrix{<:Real}, v::AbstractVector{<:Real})

`M` is a `K-n × K` full rank matrix and `v` is a `K-n` length vector. The equation `Mx = v`
describes an `n`-dimensional affine subspace. Compute `a` and `b` so that this space is
parametrized as `x = a + b * x[1:n]`. The columns of `b` are a basis for the null space of
`M` and `a` is a solution of `Mx=v`.
"""
function parametrize_subspace(M::AbstractMatrix{<:Real}, v::AbstractVector{<:Real})
    n = size(M, 2) - size(M, 1)
    @assert n >= 0
    b = nullspace(M)
    N = inv([transpose(b); M]) # this is square because M is full rank
    a = N[:, (n+1):end] * v
    return a, b
end

"""
    intersect_linear_quadratic(A::AbstractMatrix{<:Real},
    B::AbstractVector{<:Real}, C::Real, a::AbstractVector{<:Real},
    b::AbstractVector{<:Real})

Intersect the line `x = a + b * t` with the quadratic `xᵀAx + Bᵀx + C = 0` and return the
two intersection points.
"""
function intersect_linear_quadratic(A::AbstractMatrix{<:Real},
    B::AbstractVector{<:Real}, C::Real, a::AbstractVector{<:Real},
    b::AbstractVector{<:Real})
    quadratic = transpose(b) * A * b
    linear = transpose(a) * A * b + transpose(b) * A * a + transpose(B) * b
    constant = transpose(a) * A * a + transpose(B) * a + C
    t1, t2 = solve_quadratic(quadratic, linear, constant)
    return a + b * t1, a + b * t2
end

"""
    solve_quadratic(a::Real, b::Real, c::Real)

Find the two solutions of `ax² + bx + c = 0`.
"""
function solve_quadratic(a::Real, b::Real, c::Real)
    t1 = -b/(2 * a)
    t2 = sqrt(b^2 - 4 * a * c)/(2 * a)
    return t1 + t2, t1 - t2
end

"""
    sign_of_permutation(a1::AbstractVector{T}, a2::AbstractVector{T}) where T

Given two arrays which are permutations of each other, find the sign of the permutation.
"""
function sign_of_permutation(a1::AbstractVector{T}, a2::AbstractVector{T}) where T
    n = length(a1)
    @assert length(a2) == n
    inds = [findfirst(isequal(x), a2) for x in a1]
    m = MMatrix{n,n,Int}(zeros(n,n)) # det is faster for small matrices using StaticArrays
    for i in 1:n
        m[i, inds[i]] = 1
    end
    return sign(det(m))
end

"""
    first_setdiff_index(l1, l2)

Given two collections `l1` and `l2` find the index of the first element of `l1` that is
not in `l2`. If no such element exists return `0`.
"""
function first_setdiff_index(l1, l2)
    for (i, e) in enumerate(l1)
        if !(e in l2)
            return i
        end
    end
    return 0
end

"""
    first_setdiff(l1, l2)

Given two collections `l1` and `l2` find the first element of `l1` that is not in `l2`.
If no such element exists return `nothing`.
"""
function first_setdiff(l1, l2)
    i = first_setdiff_index(l1, l2)
    return i > 0 ? l1[i] : nothing
end

"""
    triangulated_lattice(v1::Vector{<:Real}, v2::Vector{<:Real}, n::Int, m::Int)

Subdivide a parallelogram with sides `v1` and `v2` into `n*m` similar parallelograms.
Divide each in two pieces to find a triangulation. Return the points and a
TriangulatedComplex.
"""
function triangulated_lattice(v1::Vector{<:Real}, v2::Vector{<:Real}, n::Int, m::Int)
    @assert length(v1) == length(v2) == 2
    inds = UniqueVector([(i, j) for i in 1:(n+1) for j in 1:(m+1)])
    points = [Point(v1 * (i-1)/n + v2 * (j-1)/m) for (i, j) in inds]
    simplices = Simplex{2, 3}[]
    for i in 1:n
        for j in 1:m
            p11 = points[findfirst(isequal((i,j)), inds)]
            p12 = points[findfirst(isequal((i,j+1)), inds)]
            p21 = points[findfirst(isequal((i+1,j)), inds)]
            p22 = points[findfirst(isequal((i+1,j+1)), inds)]
            push!(simplices, Simplex([p11, p12, p21]))
            push!(simplices, Simplex([p12, p22, p21]))
        end
    end
    return points, TriangulatedComplex(simplices)
end


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

baremodule Base

using Core.Intrinsics, Core.IR

# to start, we're going to use a very simple definition of `include`
# that doesn't require any function (except what we can get from the `Core` top-module)
const _included_files = Array{Tuple{Module,String},1}()
function include(mod::Module, path::String)
    ccall(:jl_array_grow_end, Cvoid, (Any, UInt), _included_files, UInt(1))
    Core.arrayset(true, _included_files, (mod, ccall(:jl_prepend_cwd, Any, (Any,), path)), arraylen(_included_files))
    Core.println(path)
    ccall(:jl_uv_flush, Nothing, (Ptr{Nothing},), Core.io_pointer(Core.stdout))
    Core.include(mod, path)
end
include(path::String) = include(Base, path)

# from now on, this is now a top-module for resolving syntax
const is_primary_base_module = ccall(:jl_module_parent, Ref{Module}, (Any,), Base) === Core.Main
ccall(:jl_set_istopmod, Cvoid, (Any, Bool), Base, is_primary_base_module)

# The @inline/@noinline macros that can be applied to a function declaration are not available
# until after array.jl, and so we will mark them within a function body instead.
macro inline()   Expr(:meta, :inline)   end
macro noinline() Expr(:meta, :noinline) end

# Try to help prevent users from shooting them-selves in the foot
# with ambiguities by defining a few common and critical operations
# (and these don't need the extra convert code)
getproperty(x::Module, f::Symbol) = (@inline; getfield(x, f))
setproperty!(x::Module, f::Symbol, v) = setfield!(x, f, v) # to get a decent error
getproperty(x::Type, f::Symbol) = (@inline; getfield(x, f))
setproperty!(x::Type, f::Symbol, v) = error("setfield! fields of Types should not be changed")
getproperty(x::Tuple, f::Int) = (@inline; getfield(x, f))
setproperty!(x::Tuple, f::Int, v) = setfield!(x, f, v) # to get a decent error

getproperty(x, f::Symbol) = (@inline; getfield(x, f))
setproperty!(x, f::Symbol, v) = setfield!(x, f, convert(fieldtype(typeof(x), f), v))

dotgetproperty(x, f) = getproperty(x, f)

getproperty(x::Module, f::Symbol, order::Symbol) = (@inline; getfield(x, f, order))
setproperty!(x::Module, f::Symbol, v, order::Symbol) = setfield!(x, f, v, order) # to get a decent error
getproperty(x::Type, f::Symbol, order::Symbol) = (@inline; getfield(x, f, order))
setproperty!(x::Type, f::Symbol, v, order::Symbol) = error("setfield! fields of Types should not be changed")
getproperty(x::Tuple, f::Int, order::Symbol) = (@inline; getfield(x, f, order))
setproperty!(x::Tuple, f::Int, v, order::Symbol) = setfield!(x, f, v, order) # to get a decent error

getproperty(x, f::Symbol, order::Symbol) = (@inline; getfield(x, f, order))
setproperty!(x, f::Symbol, v, order::Symbol) = (@inline; setfield!(x, f, convert(fieldtype(typeof(x), f), v), order))

swapproperty!(x, f::Symbol, v, order::Symbol=:notatomic) =
    (@inline; Core.swapfield!(x, f, convert(fieldtype(typeof(x), f), v), order))
modifyproperty!(x, f::Symbol, op, v, order::Symbol=:notatomic) =
    (@inline; Core.modifyfield!(x, f, op, v, order))
replaceproperty!(x, f::Symbol, expected, desired, success_order::Symbol=:notatomic, fail_order::Symbol=success_order) =
    (@inline; Core.replacefield!(x, f, expected, convert(fieldtype(typeof(x), f), desired), success_order, fail_order))

convert(::Type{Any}, Core.@nospecialize x) = x
convert(::Type{T}, x::T) where {T} = x
include("coreio.jl")

eval(x) = Core.eval(Base, x)
eval(m::Module, x) = Core.eval(m, x)

# init core docsystem
import Core: @doc, @__doc__, WrappedException, @int128_str, @uint128_str, @big_str, @cmd
if isdefined(Core, :Compiler)
    import Core.Compiler.CoreDocs
    Core.atdoc!(CoreDocs.docm)
end

include("exports.jl")

if false
    # simple print definitions for debugging. enable these if something
    # goes wrong during bootstrap before printing code is available.
    # otherwise, they just just eventually get (noisily) overwritten later
    global show, print, println
    show(io::IO, x) = Core.show(io, x)
    print(io::IO, a...) = Core.print(io, a...)
    println(io::IO, x...) = Core.println(io, x...)
end

"""
    time_ns()

Get the time in nanoseconds. The time corresponding to 0 is undefined, and wraps every 5.8 years.
"""
time_ns() = ccall(:jl_hrtime, UInt64, ())

start_base_include = time_ns()

## Load essential files and libraries
include("essentials.jl")
include("ctypes.jl")
include("gcutils.jl")
include("generator.jl")
include("reflection.jl")
include("options.jl")

# core operations & types
include("promotion.jl")
include("tuple.jl")
include("expr.jl")
Pair{A, B}(@nospecialize(a), @nospecialize(b)) where {A, B} = (@inline; Pair{A, B}(convert(A, a)::A, convert(B, b)::B))
#Pair{Any, B}(@nospecialize(a::Any), b) where {B} = (@inline; Pair{Any, B}(a, Base.convert(B, b)::B))
#Pair{A, Any}(a, @nospecialize(b::Any)) where {A} = (@inline; Pair{A, Any}(Base.convert(A, a)::A, b))
include("pair.jl")
include("traits.jl")
include("range.jl")
include("error.jl")

# core numeric operations & types
==(x, y) = x === y
include("bool.jl")
include("number.jl")
include("int.jl")
include("operators.jl")
include("pointer.jl")
include("refvalue.jl")
include("refpointer.jl")
include("checked.jl")
using .Checked

# Lazy strings
include("strings/lazy.jl")

# array structures
include("indices.jl")
include("array.jl")
include("abstractarray.jl")
include("subarray.jl")
include("views.jl")
include("baseext.jl")

include("ntuple.jl")

include("abstractdict.jl")
include("iddict.jl")
include("idset.jl")

include("iterators.jl")
using .Iterators: zip, enumerate, only
using .Iterators: Flatten, Filter, product  # for generators

include("namedtuple.jl")

# For OS specific stuff
# We need to strcat things here, before strings are really defined
function strcat(x::String, y::String)
    out = ccall(:jl_alloc_string, Ref{String}, (Csize_t,), Core.sizeof(x) + Core.sizeof(y))
    GC.@preserve x y out begin
        out_ptr = unsafe_convert(Ptr{UInt8}, out)
        unsafe_copyto!(out_ptr, unsafe_convert(Ptr{UInt8}, x), Core.sizeof(x))
        unsafe_copyto!(out_ptr + Core.sizeof(x), unsafe_convert(Ptr{UInt8}, y), Core.sizeof(y))
    end
    return out
end
include(strcat((length(Core.ARGS)>=2 ? Core.ARGS[2] : ""), "build_h.jl"))     # include($BUILDROOT/base/build_h.jl)
include(strcat((length(Core.ARGS)>=2 ? Core.ARGS[2] : ""), "version_git.jl")) # include($BUILDROOT/base/version_git.jl)

# numeric operations
include("hashing.jl")
include("rounding.jl")
using .Rounding
include("div.jl")
include("float.jl")
include("twiceprecision.jl")
include("complex.jl")
include("rational.jl")
include("multinverses.jl")
using .MultiplicativeInverses
include("abstractarraymath.jl")
include("arraymath.jl")

# SIMD loops
@pure sizeof(s::String) = Core.sizeof(s)  # needed by gensym as called from simdloop
include("simdloop.jl")
using .SimdLoop

# map-reduce operators
include("reduce.jl")

## core structures
include("reshapedarray.jl")
include("reinterpretarray.jl")
include("bitarray.jl")
include("bitset.jl")

if !isdefined(Core, :Compiler)
    include("docs/core.jl")
    Core.atdoc!(CoreDocs.docm)
end

include("multimedia.jl")
using .Multimedia

# Some type
include("some.jl")

include("dict.jl")
include("abstractset.jl")
include("set.jl")

# Strings
include("char.jl")
include("strings/basic.jl")
include("strings/string.jl")
include("strings/substring.jl")

# Initialize DL_LOAD_PATH as early as possible.  We are defining things here in
# a slightly more verbose fashion than usual, because we're running so early.
const DL_LOAD_PATH = String[]
let os = ccall(:jl_get_UNAME, Any, ())
    if os === :Darwin || os === :Apple
        if Base.DARWIN_FRAMEWORK
            push!(DL_LOAD_PATH, "@loader_path/Frameworks")
        end
        push!(DL_LOAD_PATH, "@loader_path")
    end
end

include("osutils.jl")
include("c.jl")

# Core I/O
include("io.jl")
include("iobuffer.jl")

# strings & printing
include("intfuncs.jl")
include("strings/strings.jl")
include("regex.jl")
include("parse.jl")
include("shell.jl")
include("show.jl")
include("arrayshow.jl")
include("methodshow.jl")

# multidimensional arrays
include("cartesian.jl")
using .Cartesian
include("multidimensional.jl")

include("broadcast.jl")
using .Broadcast
using .Broadcast: broadcasted, broadcasted_kwsyntax, materialize, materialize!,
                  broadcast_preserving_zero_d, andand, oror

# missing values
include("missing.jl")

# version
include("version.jl")

# system & environment
include("sysinfo.jl")
include("libc.jl")
using .Libc: getpid, gethostname, time

# Logging
include("logging.jl")
using .CoreLogging

# Concurrency
include("linked_list.jl")
include("condition.jl")
include("threads.jl")
include("lock.jl")
include("channels.jl")
include("task.jl")
include("threads_overloads.jl")
include("weakkeydict.jl")

include("env.jl")

# functions defined in Random
function rand end
function randn end

# I/O
include("libuv.jl")
include("asyncevent.jl")
include("iostream.jl")
include("stream.jl")
include("filesystem.jl")
using .Filesystem
include("cmd.jl")
include("process.jl")
include("ttyhascolor.jl")
include("secretbuffer.jl")

# RandomDevice support
include("randomdevice.jl")

# core math functions
include("floatfuncs.jl")
include("math.jl")
using .Math
const (√)=sqrt
const (∛)=cbrt

# now switch to a simple, race-y TLS, relative include for the rest of Base
delete_method(which(include, (Module, String)))
let SOURCE_PATH = ""
    global function include(mod::Module, path::String)
        prev = SOURCE_PATH::String
        path = normpath(joinpath(dirname(prev), path))
        Core.println(path)
        ccall(:jl_uv_flush, Nothing, (Ptr{Nothing},), Core.io_pointer(Core.stdout))
        push!(_included_files, (mod, abspath(path)))
        SOURCE_PATH = path
        result = Core.include(mod, path)
        SOURCE_PATH = prev
        return result
    end
end

# reduction along dims
include("reducedim.jl")  # macros in this file rely on string.jl
include("accumulate.jl")

include("permuteddimsarray.jl")
using .PermutedDimsArrays

# basic data structures
include("ordering.jl")
using .Order

# Combinatorics
include("sort.jl")
using .Sort

# BinaryPlatforms, used by Artifacts.  Needs `Sort`.
include("binaryplatforms.jl")

# Fast math
include("fastmath.jl")
using .FastMath

function deepcopy_internal end

# enums
include("Enums.jl")
using .Enums

# BigInts
include("gmp.jl")
using .GMP

# float printing: requires BigInt
include("ryu/Ryu.jl")
using .Ryu

# BigFloats
include("mpfr.jl")
using .MPFR

include("combinatorics.jl")

# irrational mathematical constants
include("irrationals.jl")
include("mathconstants.jl")
using .MathConstants: ℯ, π, pi

# metaprogramming
include("meta.jl")

# Stack frames and traces
include("stacktraces.jl")
using .StackTraces

# experimental API's
include("experimental.jl")

# utilities
include("deepcopy.jl")
include("download.jl")
include("summarysize.jl")
include("errorshow.jl")

include("initdefs.jl")

# worker threads
include("threadcall.jl")

# code loading
include("uuid.jl")
include("pkgid.jl")
include("toml_parser.jl")
include("loading.jl")

# misc useful functions & macros
include("timing.jl")
include("util.jl")

include("asyncmap.jl")

# deprecated functions
include("deprecated.jl")

# Some basic documentation
include("docs/basedocs.jl")

include("client.jl")

# Documentation -- should always be included last in sysimg.
include("docs/Docs.jl")
using .Docs
if isdefined(Core, :Compiler) && is_primary_base_module
    Docs.loaddocs(Core.Compiler.CoreDocs.DOCS)
end

# finally, now make `include` point to the full version
for m in methods(include)
    delete_method(m)
end
# These functions are duplicated in client.jl/include(::String) for
# nicer stacktraces. Modifications here have to be backported there
include(mod::Module, _path::AbstractString) = _include(identity, mod, _path)
include(mapexpr::Function, mod::Module, _path::AbstractString) = _include(mapexpr, mod, _path)

end_base_include = time_ns()

const _sysimage_modules = PkgId[]
in_sysimage(pkgid::PkgId) = pkgid in _sysimage_modules

# Precompiles for Revise
# TODO: move these to contrib/generate_precompile.jl
# The problem is they don't work there
for match = _methods(+, (Int, Int), -1, get_world_counter())
    m = match.method
    delete!(push!(Set{Method}(), m), m)
    copy(Core.Compiler.retrieve_code_info(Core.Compiler.specialize_method(match)))

    empty!(Set())
    push!(push!(Set{Union{GlobalRef,Symbol}}(), :two), GlobalRef(Base, :two))
    (setindex!(Dict{String,Base.PkgId}(), Base.PkgId(Base), "file.jl"))["file.jl"]
    (setindex!(Dict{Symbol,Vector{Int}}(), [1], :two))[:two]
    (setindex!(Dict{Base.PkgId,String}(), "file.jl", Base.PkgId(Base)))[Base.PkgId(Base)]
    (setindex!(Dict{Union{GlobalRef,Symbol}, Vector{Int}}(), [1], :two))[:two]
    (setindex!(IdDict{Type, Union{Missing, Vector{Tuple{LineNumberNode, Expr}}}}(), missing, Int))[Int]
    Dict{Symbol, Union{Nothing, Bool, Symbol}}(:one => false)[:one]
    Dict(Base => [:(1+1)])[Base]
    Dict(:one => [1])[:one]
    Dict("abc" => Set())["abc"]
    pushfirst!([], sum)
    get(Base.pkgorigins, Base.PkgId(Base), nothing)
    sort!([1,2,3])
    unique!([1,2,3])
    cumsum([1,2,3])
    append!(Int[], BitSet())
    isempty(BitSet())
    delete!(BitSet([1,2]), 3)
    deleteat!(Int32[1,2,3], [1,3])
    deleteat!(Any[1,2,3], [1,3])
    Core.svec(1, 2) == Core.svec(3, 4)
    any(t->t[1].line > 1, [(LineNumberNode(2,:none), :(1+1))])

    # Code loading uses this
    sortperm(mtime.(readdir(".")), rev=true)

    break   # only actually need to do this once
end

if is_primary_base_module
function __init__()
    # for the few uses of Libc.rand in Base:
    Libc.srand()
    # Base library init
    reinit_stdio()
    Multimedia.reinit_displays() # since Multimedia.displays uses stdout as fallback
    # initialize loading
    init_depot_path()
    init_load_path()
    init_active_project()
    append!(empty!(_sysimage_modules), keys(loaded_modules))
    if haskey(ENV, "JULIA_MAX_NUM_PRECOMPILE_FILES")
        MAX_NUM_PRECOMPILE_FILES[] = parse(Int, ENV["JULIA_MAX_NUM_PRECOMPILE_FILES"])
    end
    nothing
end

end

end # baremodule Base

max_stars_repo_name: fkastner/julia

max_issues_repo_name: fkastner/julia

max_forks_repo_name: fkastner/julia