text:
# jni_md.h
typealias jint Cint
#ifdef _LP64 /* 64-bit Solaris */
# typedef long jlong;
typealias jlong Clonglong
typealias jbyte Cchar

# jni.h

typealias jboolean Cuchar
typealias jchar Cushort
typealias jshort Cshort
typealias jfloat Cfloat
typealias jdouble Cdouble
typealias jsize jint
jprimitive = @compat Union{jboolean, jchar, jshort, jfloat, jdouble, jint, jlong}

immutable JavaMetaClass{T}
    ptr::Ptr{Void}
end

#The metaclass, sort of equivalent to a the
JavaMetaClass(T, ptr) = JavaMetaClass{T}(ptr)

type JavaObject{T}
    ptr::Ptr{Void}

    function JavaObject(ptr)
        j=new(ptr)
        finalizer(j, deleteref)
        return j
    end

    JavaObject(argtypes::Tuple, args...) = jnew(T, argtypes, args...)

end

JavaObject(T, ptr) = JavaObject{T}(ptr)

function deleteref(x::JavaObject)
    if x.ptr == C_NULL; return; end
    if (penv==C_NULL); return; end
    #ccall(:jl_,Void,(Any,),x)
    ccall(jnifunc.DeleteLocalRef, Void, (Ptr{JNIEnv}, Ptr{Void}), penv, x.ptr)
    x.ptr=C_NULL #Safety in case this function is called direcly, rather than at finalize
    return
end

"""
```
isnull(obj::JavaObject)
```
Checks if the passed JavaObject is null or not

### Args
* obj: The object of type JavaObject

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaObject) = obj.ptr == C_NULL

"""
```
isnull(obj::JavaMetaClass)
```
Checks if the passed JavaMetaClass is null or not

### Args
* obj: The object of type JavaMetaClass

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaMetaClass) = obj.ptr == C_NULL

typealias JClass JavaObject{Symbol("java.lang.Class")}
typealias JObject JavaObject{Symbol("java.lang.Object")}
typealias JMethod JavaObject{Symbol("java.lang.reflect.Method")}
typealias JString JavaObject{Symbol("java.lang.String")}

function JString(str::AbstractString)
    jstring = ccall(jnifunc.NewStringUTF, Ptr{Void}, (Ptr{JNIEnv}, Ptr{UInt8}), penv, String(str))
    if jstring == C_NULL
        geterror()
    else
        return JString(jstring)
    end
end

# jvalue(v::Integer) = int64(v) << (64-8*sizeof(v))
jvalue(v::Integer) = @compat Int64(v)
jvalue(v::Float32) = jvalue(reinterpret(Int32, v))
jvalue(v::Float64) = jvalue(reinterpret(Int64, v))
jvalue(v::Ptr) = jvalue(@compat Int(v))

macro jimport(class)
    if isa(class, Expr)
        juliaclass=sprint(Base.show_unquoted, class)
    elseif isa(class, Symbol)
        juliaclass=string(class)
    elseif isa(class, AbstractString)
        juliaclass=class
    else
        error("Macro parameter is of type $(typeof(class))!\nShould be Expr, Symbol or String")
    end
    juliaclass = replace(juliaclass, " ", "") #handle $ for innerclass
    quote
        JavaObject{(Base.Symbol($juliaclass))}
    end
end

function jnew(T::Symbol, argtypes::Tuple, args...)
    sig = method_signature(Void, argtypes...)
    jmethodId = ccall(jnifunc.GetMethodID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T), String("<init>"), sig)
    if (jmethodId == C_NULL)
        error("No constructor for $T with signature $sig")
    end
    return  _jcall(metaclass(T), jmethodId, jnifunc.NewObjectA, JavaObject{T}, argtypes, args...)
end

# Call static methods
function jcall{T}(typ::Type{JavaObject{T}}, method::AbstractString, rettype::Type, argtypes::Tuple, args... )
    try
        gc_enable(false)
        sig = method_signature(rettype, argtypes...)

        jmethodId = ccall(jnifunc.GetStaticMethodID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T), String(method), sig)
        if jmethodId==C_NULL; geterror(true); end

        _jcall(metaclass(T), jmethodId, C_NULL, rettype, argtypes, args...)
    finally
        gc_enable(true)
    end

end

# Call instance methods
function jcall(obj::JavaObject, method::AbstractString, rettype::Type, argtypes::Tuple, args... )
    try
        gc_enable(false)
        sig = method_signature(rettype, argtypes...)
        jmethodId = ccall(jnifunc.GetMethodID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(obj), String(method), sig)
        if jmethodId==C_NULL; geterror(true); end
        _jcall(obj, jmethodId, C_NULL, rettype,  argtypes, args...)
    finally
        gc_enable(true)
    end
end

function jfield{T}(typ::Type{JavaObject{T}}, field::AbstractString, fieldType::Type)
    jfieldID  = ccall(jnifunc.GetStaticFieldID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T), String(field), signature(fieldType))
    if jfieldID==C_NULL; geterror(true); end
    _jfield(metaclass(T), jfieldID, fieldType)
end

function jfield(obj::JavaObject, field::AbstractString, fieldType::Type)
    jfieldID  = ccall(jnifunc.GetFieldID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(obj), String(field), signature(fieldType))
    if jfieldID==C_NULL; geterror(true); end
    _jfield(obj, jfieldID, fieldType)
end

for (x, y, z) in [ (:jboolean, :(jnifunc.GetBooleanField), :(jnifunc.GetStaticBooleanField)),
                  (:jchar, :(jnifunc.GetCharField), :(jnifunc.GetStaticCharField)),
                  (:jbyte, :(jnifunc.GetByteField), :(jnifunc.GetStaticBypeField)),
                  (:jshort, :(jnifunc.GetShortField), :(jnifunc.GetStaticShortField)),
                  (:jint, :(jnifunc.GetIntField), :(jnifunc.GetStaticIntField)),
                  (:jlong, :(jnifunc.GetLongField), :(jnifunc.GetStaticLongField)),
                  (:jfloat, :(jnifunc.GetFloatField), :(jnifunc.GetStaticFloatField)),
                  (:jdouble, :(jnifunc.GetDoubleField), :(jnifunc.GetStaticDoubleField)) ]

    m = quote
        function _jfield(obj, jfieldID::Ptr{Void}, fieldType::Type{$(x)})
            callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            result = ccall(callmethod, $x, (Ptr{JNIEnv}, Ptr{Void}, Ptr{Void}), penv, obj.ptr, jfieldID)
            if result==C_NULL; geterror(); end
            return convert_result(fieldType, result)
        end
    end
    eval(m)
end

function _jfield(obj, jfieldID::Ptr{Void}, fieldType::Type)
    callmethod = ifelse( typeof(obj)<:JavaObject, jnifunc.GetObjectField , jnifunc.GetStaticObjectField )
    result = ccall(callmethod, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{Void}), penv, obj.ptr, jfieldID)
    if result==C_NULL; geterror(); end
    return convert_result(fieldType, result)
end

#Generate these methods to satisfy ccall's compile time constant requirement
#_jcall for primitive and Void return types
for (x, y, z) in [ (:jboolean, :(jnifunc.CallBooleanMethodA), :(jnifunc.CallStaticBooleanMethodA)),
                  (:jchar, :(jnifunc.CallCharMethodA), :(jnifunc.CallStaticCharMethodA)),
                  (:jbyte, :(jnifunc.CallByteMethodA), :(jnifunc.CallStaticByteMethodA)),
                  (:jshort, :(jnifunc.CallShortMethodA), :(jnifunc.CallStaticShortMethodA)),
                  (:jint, :(jnifunc.CallIntMethodA), :(jnifunc.CallStaticIntMethodA)),
                  (:jlong, :(jnifunc.CallLongMethodA), :(jnifunc.CallStaticLongMethodA)),
                  (:jfloat, :(jnifunc.CallFloatMethodA), :(jnifunc.CallStaticFloatMethodA)),
                  (:jdouble, :(jnifunc.CallDoubleMethodA), :(jnifunc.CallStaticDoubleMethodA)),
                  (:Void, :(jnifunc.CallVoidMethodA), :(jnifunc.CallStaticVoidMethodA)) ]
    m = quote
        function _jcall(obj,  jmethodId::Ptr{Void}, callmethod::Ptr{Void}, rettype::Type{$(x)}, argtypes::Tuple, args... )
            if callmethod == C_NULL #!
                callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            end
            @assert callmethod != C_NULL
            @assert jmethodId != C_NULL
            if(isnull(obj)); error("Attempt to call method on Java NULL"); end
            savedArgs, convertedArgs = convert_args(argtypes, args...)
            result = ccall(callmethod, $x , (Ptr{JNIEnv}, Ptr{Void}, Ptr{Void}, Ptr{Void}), penv, obj.ptr, jmethodId, convertedArgs)
            if result==C_NULL; geterror(); end
            if result == nothing; return; end
            return convert_result(rettype, result)
        end
    end
    eval(m)
end

#_jcall for Object return types
#obj -- reciever - Class pointer or object prointer
#jmethodId -- Java method ID
#callmethod -- the C method pointer to call
function _jcall(obj,  jmethodId::Ptr{Void}, callmethod::Ptr{Void}, rettype::Type, argtypes::Tuple, args... )
    if callmethod == C_NULL
        callmethod = ifelse( typeof(obj)<:JavaObject, jnifunc.CallObjectMethodA , jnifunc.CallStaticObjectMethodA )
    end
    @assert callmethod != C_NULL
    @assert jmethodId != C_NULL
    if(isnull(obj)); error("Attempt to call method on Java NULL"); end
    savedArgs, convertedArgs = convert_args(argtypes, args...)
    result = ccall(callmethod, Ptr{Void} , (Ptr{JNIEnv}, Ptr{Void}, Ptr{Void}, Ptr{Void}), penv, obj.ptr, jmethodId, convertedArgs)
    if result==C_NULL; geterror(); end
    return convert_result(rettype, result)
end

global const _jmc_cache = Dict{Symbol, JavaMetaClass}()

function _metaclass(class::Symbol)
    jclass=javaclassname(class)
    jclassptr = ccall(jnifunc.FindClass, Ptr{Void}, (Ptr{JNIEnv}, Ptr{UInt8}), penv, jclass)
    if jclassptr == C_NULL; error("Class Not Found $jclass"); end
    return JavaMetaClass(class, jclassptr)
end

function metaclass(class::Symbol)
    if !haskey(_jmc_cache, class)
        _jmc_cache[class] = _metaclass(class)
    end
    return _jmc_cache[class]
end

metaclass{T}(::Type{JavaObject{T}}) = metaclass(T)
metaclass{T}(::JavaObject{T}) = metaclass(T)

javaclassname(class::Symbol) = replace(string(class), '.', '/')

function geterror(allow=false)
    isexception = ccall(jnifunc.ExceptionCheck, jboolean, (Ptr{JNIEnv},), penv )

    if isexception == JNI_TRUE
        jthrow = ccall(jnifunc.ExceptionOccurred, Ptr{Void}, (Ptr{JNIEnv},), penv)
        if jthrow==C_NULL ; error("Java Exception thrown, but no details could be retrieved from the JVM"); end
        ccall(jnifunc.ExceptionDescribe, Void, (Ptr{JNIEnv},), penv ) #Print java stackstrace to stdout
        ccall(jnifunc.ExceptionClear, Void, (Ptr{JNIEnv},), penv )
        jclass = ccall(jnifunc.FindClass, Ptr{Void}, (Ptr{JNIEnv},Ptr{UInt8}), penv, "java/lang/Throwable")
        if jclass==C_NULL; error("Java Exception thrown, but no details could be retrieved from the JVM"); end
        jmethodId=ccall(jnifunc.GetMethodID, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), penv, jclass, "toString", "()Ljava/lang/String;")
        if jmethodId==C_NULL; error("Java Exception thrown, but no details could be retrieved from the JVM"); end
        res = ccall(jnifunc.CallObjectMethodA, Ptr{Void}, (Ptr{JNIEnv}, Ptr{Void}, Ptr{Void}, Ptr{Void}), penv, jthrow, jmethodId,C_NULL)
        if res==C_NULL; error("Java Exception thrown, but no details could be retrieved from the JVM"); end
        msg = unsafe_string(JString(res))
        ccall(jnifunc.DeleteLocalRef, Void, (Ptr{JNIEnv}, Ptr{Void}), penv, jthrow)
        error(string("Error calling Java: ",msg))
    else
        if allow==false
            return #No exception pending, legitimate NULL returned from Java
        else
            error("Null from Java. Not known how")
        end
    end
end

#get the JNI signature string for a method, given its
#return type and argument types
function method_signature(rettype, argtypes...)
    s=IOBuffer()
    write(s, "(")
    for arg in argtypes
        write(s, signature(arg))
    end
    write(s, ")")
    write(s, signature(rettype))
    return takebuf_string(s)
end

#get the JNI signature string for a given type
function signature(arg::Type)
    if is(arg, jboolean)
        return "Z"
    elseif is(arg, jbyte)
        return "B"
    elseif is(arg, jchar)
        return "C"
    elseif is(arg, jshort)
        return "S"
    elseif is(arg, jint)
        return "I"
    elseif is(arg, jlong)
        return "J"
    elseif is(arg, jfloat)
        return "F"
    elseif is(arg, jdouble)
        return "D"
    elseif is(arg, Void)
        return "V"
    elseif issubtype(arg, Array)
        return string("[", signature(eltype(arg)))
    end
end

signature{T}(arg::Type{JavaObject{T}}) = string("L", javaclassname(T), ";")


text:
Silixa_fmt = FormatDesc(
  "Silixa TDMS",
  "\"silixa\"",
  "Silixa, Hertfordshire, UK",
  "adapted from https://silixa.com/resources/software-downloads/",
  "https://silixa.com/about-us/contact/",
  HistVec(),
  ["Silixa variant on NI LabVIEW TDMS file format"],
  ["Silixa (nodal array data)"],
  [""],
  0xfd
  )
Silixa_fmt.ver = [  FmtVer(1, "2018-06-28", false) ]
formats["Silixa TDMS"] = Silixa_fmt

ext: jl


text:
using LsqFit
using Statistics
using StatsBase

"""
fit_straight_line(tdata::Vector{Float64}, vdata::Vector{Float64}; pa0=[0.0, 0.5], i0=1)

# Fields

`tdata`: time-like (x) data
`vdata`` : voltage-like (y) data
`pa0`: initial values of parameters
`i0` : index of first point in array to be used in the fit

# Returns
- coefficient of the field (a vector)
- standard errors of coefficients (a vector)
- yfit = fit(tdata), that is the prediction of the values corresponding to each tdata
  for the fit parameters

"""
function fit_straight_line(tdata::Vector{Float64}, vdata::Vector{Float64}; pa0=[0.0, 0.5], i0=1)
	tfun(t, a, b) = a + b * t
	pfun(t, p) = p[1] .+ p[2] .* t
	il = length(tdata)
	fit = curve_fit(pfun, tdata[i0:il], vdata[i0:il], pa0)
	coef(fit), stderror(fit), tfun.(tdata, coef(fit)...)
end

"""
fit_expo(tdata::Vector{Float64}, vdata::Vector{Float64}; pa0=[0.0, 0.5], i0=1)

# Fields

`tdata`: time-like (x) data
`vdata`` : voltage-like (y) data
`pa0`: initial values of parameters
`i0` : index of first point in array to be used in the fit

# Returns
- coefficient of the field (a vector)
- standard errors of coefficients (a vector)
- yfit = fit(tdata), that is the prediction of the values corresponding to each tdata
  for the fit parameters

"""
function fit_expo(tdata::Vector{Float64}, vdata::Vector{Float64}; pa0=[0.0, 0.5], i0=1)
	tfun(t, N, λ) = N*exp(-t/λ)
	pfun(t, p) = p[1] * exp.(-t/p[2])
	il = length(tdata)
	fit = curve_fit(pfun, tdata[i0:il], vdata[i0:il], pa0)
	coef(fit), stderror(fit), tfun.(tdata, coef(fit)...)
end

"""
    fit_time_histo(htime::Histogram; pa0=[100.0, 0.5], i0=1)

# Fields

	`htime` : An histogram of times
	`pa0`   : initial value of fit parameters
	`i0`    : first bin to fit

# Returns

# Returns
- coefficient of the field (a vector)
- standard errors of coefficients (a vector)
- yfit = fit(tdata), that is the prediction of the values corresponding to each tdata
  for the fit parameters

"""
function fit_time_histo(htime::Histogram; pa0=[100.0, 0.5], i0=1)
	tdata = centers(htime)
	vdata = htime.weights
    fit_expo(tdata, vdata; pa0=pa0, i0=i0)
end

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
"""
Return an IO that will always read in the right endianness.
"""
function correct_endianess(io::IO)
    sfnt = Base.peek(io, UInt32)
    if sfnt == 0x00000100
        SwapStream(io)
    else
        io
    end
end

version_16_dot_16(version::UInt32) = VersionNumber(version >> 16 + (version & 0x0000ffff) >> 8)

function word_align(size)
    4 * cld(size, 4)
end

include("parsing/table_records.jl")
include("parsing/font_header.jl")
include("parsing/maximum_profile.jl")
include("parsing/char_to_glyph.jl")
include("parsing/metrics.jl")
include("parsing/loca.jl")


text:
# heatmaps

# Plots
import Pkg
Pkg.add("Plots")
Pkg.add("PyPlot")
using Plots

using Random
Random.seed!(2018)

a = rand(5,5)
xlabel = string.(collect('A':'E'))
ylabel = string.(collect('a':'e'))
Plots.heatmap(a, xticks=(1:5, xlabel), yticks=(1:5, ylabel),
    fill_z=a, aspect_ratio=:equal)

fontsize = 15
nrow, ncol = size(a)
ann = [(i,j, Plots.text(round(a[i,j], digits=2), fontsize, :white, :center))
            for i in 1:nrow for j in 1:ncol]
annotate!(ann, linecolor=:white)

# makie
using Random

let
    Random.seed!(123)
    fig = Figure(resolution = (600, 400))
    ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "y")
    hmap = heatmap!(2rand(20, 20) .- 1; colormap = :RdPu_9)
    Colorbar(fig[1, 2], hmap; label = "values", width = 15, ticksize = 15, tickalign = 1)
    colsize!(fig.layout, 1, Aspect(1, 1.0))
    colgap!(fig.layout, 7)
    display(fig)
end;

# gadfly
using Gadfly

parameters = ["ξ₁", "η₁",  "ξ₂", "η₂"]
palettef = Scale.lab_gradient("darkgreen", "orange", "blue")

Gadfly.spy(rand(4,4), Scale.y_discrete(labels = i->parameters[i]), Scale.x_discrete,
Guide.ylabel("Parameters"), Guide.xlabel("Mode"), Scale.color_continuous(colormap=palettef, minvalue=0, maxvalue=1.5))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
struct TimePolynomial{T, P<:AbstractPolynomial{T}} <: AbstractPolynomial{T}
    p::P
    function TimePolynomial{T, P}(p::P) where {T, P<:AbstractPolynomial{T}}
        @assert length(variables(p)) == 1
        new{T, P}(p)
    end
end

TimePolynomial(p::P) where {T, P<:AbstractPolynomial{T}} = TimePolynomial{T, P}(p)

(p::TimePolynomial)(x) = p.p(variable(p) => x)
variable(p::TimePolynomial) = first(variables(p))
differentiate(p::TimePolynomial) = TimePolynomial(differentiate(p, first(variables(p))))
polynomial(p::TimePolynomial) = p.p
nvariables(::TimePolynomial) = 1
variables(p::TimePolynomial) = variables(p.p)
Base.show(io::IO, p::TimePolynomial) = print(io, p.p)

function PiecewiseFunctions.PiecewiseFunction(ps::AbstractVector{<:AbstractPolynomial}, breaks::AbstractVector)
    if eltype(ps) <: TimePolynomial
        error("huh?")
    end
    PiecewiseFunction(TimePolynomial.(ps), breaks)
end

PiecewiseFunctions.PiecewiseFunction(ps::AbstractVector{<:TimePolynomial}, breaks::AbstractVector) =
    PiecewiseFunction{eltype(ps), eltype(breaks), typeof(ps), typeof(breaks)}(ps, breaks)

function cubic_spline(Δt, x0, xf, ẋ0, ẋf)
    @polyvar t
    # p(0) = x0
    # p'(0) = x'0
    # p = a + bt + ct^2 + dt^3
    # a = x0
    # b = x'0
    # b + 2ct + 3dt^2 = x'f
    # a + bt + ct^2 + dt^3 = xf
    # bt + 2ct^2 + 3dt^3 = tx'f
    # 3a + 3bt + 3ct^2 + 3dt^3 = 3xf
    # 3a + 2bt + ct^2 = 3xf - tx'f
    a = x0
    b = ẋ0
    c = (3*xf - Δt*ẋf - 3*a - 2*b*Δt) / (Δt^2)
    d = (ẋf - b - 2*c*Δt) / (3*Δt^2)
    TimePolynomial(a + b*t + c*t^2 + d*t^3)
end

function cubic_spline(times::AbstractVector, positions::AbstractVector, velocities::AbstractVector)
    @assert length(times) > 1
    @assert length(times) == length(positions) == length(velocities)
    PiecewiseFunction(
        [cubic_spline(times[i] - times[i-1], positions[i-1], positions[i], velocities[i-1], velocities[i]) for i in 2:length(times)],
        times
    )
end


text:
import StatsFuns.poispdf

seed!(12345)

function testPoisson(μ, low, high, N)
  function pmf(x::Int64)
    return poispdf(μ, x)
  end
  function sampler()
    return samplePoisson(μ)
  end
  return testGOFDiscrete(pmf, sampler, low, high, N)
end

@test testPoisson(3.0, 0, 12, 2^21) > 0.01
@test testPoisson(15.0, 5, 30, 2^21) > 0.01
@test testPoisson(50.0, 25, 80, 2^21) > 0.01


text:
using seisimu

vp = [5000.0]
vs = [3500.0]
rho = 2.8
depth = 2000
Horizon = 2000
pkf = 30
T = 1
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

model = initmodel(vp,vs,rho,depth,Horizon,pkf,T);

model.wf.txz = rand(Float64,model.nwf.BDntxz[1],model.nwf.BDntxz[2]);

model.wf.txx = rand(Float64,model.nwf.BDntpp[1],model.nwf.BDntpp[2]);

model.wf.tzz = rand(Float64,model.nwf.BDntpp[1],model.nwf.BDntpp[2]);

@time begin
for i in 1:200
extrap2d!(model);
end
end


text:
# http://web.mit.edu/puzzle/www/2015/puzzle/school_of_fish_meta/
@testset "school of fish meta" begin
    f = best_feature(wordlist(IOBuffer("""
        cred
        ruble
        marble
        lodge
        regent
        barfed
        boone
        lunge
        toughen
        freeing
        ahead
        stolen
        mice
        elite
        rambles
        octal
        crook
        sedan
        amadeus
        dangle
        boucher
        clanlaird
        crasher
        eking
        dylan
        sturgeon
        wraps
        stooge
        rotorset
        ankles
        coloredice
        raptors
        spare net
        sick porno
        hawke
        limon
        taco
        tuba
        others
        triage
        heartpin
        rant
        blazer
        isabel
        baron
        paws
        opera
        internal
        pikers
        soaped
        draft
        papers
        serbia
        crowds
        cleaned
        piper
        """)), 10)
    @test f.description == "has a 1-letter transdeletion"
end


text:
#=
# Sample texts for usage
# * https://benlauwens.github.io/ThinkJulia.jl/latest/book.html
=#

SampleDocument = raw"""
In January 2018 I started the preparation of a programming course targeting students without programming experience. I wanted to use Julia, but I found that there existed no book with the purpose of learning to program with Julia as the first programming language. There are wonderful tutorials that explain Julia’s key concepts, but none of them pay sufficient attention to learning how to think like a programmer.

I knew the book Think Python by Allen Downey, which contains all the key ingredients to learn to program properly. However, this book was based on the Python programming language. My first draft of the course notes was a melting pot of all kinds of reference works, but the longer I worked on it, the more the content started to resemble the chapters of Think Python. Soon, the idea of developing my course notes as a port of that book to Julia came to fruition.

All the material was available as Jupyter notebooks in a GitHub repository. After I posted a message on the Julia Discourse site about the progress of my course, the feedback was overwhelming. A book about basic programming concepts with Julia as the first programming language was apparently a missing link in the Julia universe. I contacted Allen to ask if I could start an official port of Think Python to Julia, and his answer was immediate: “Go for it!” He put me in touch with his editor at O’Reilly Media, and a year later I was putting the finishing touches on this book.

It was a bumpy ride. In August 2018 Julia v1.0 was released, and like all my fellow Julia programmers I had to do a migration of the code. All the examples in the book were tested during the conversion of the source files to O’Reilly-compatible AsciiDoc files. Both the toolchain and the example code had to be made Julia v1.0–compliant. Luckily, there are no lectures to give in August….

I hope you enjoy working with this book, and that it helps you learn to program and think like a computer scientist, at least a little bit.

Ben Lauwens
"""


text:
export ChunkedTransferDecodeStream

mutable struct ChunkedTransferDecodeStream <: IO
    source::IO
    maxsize::UInt64
    nread::UInt64
    chunkstart::UInt64
    chunksize::UInt64   # of the current chunk
    ondone::Function
    onstart::Function
end

function ChunkedTransferDecodeStream(source::IO; maxsize=typemax(UInt64), ondone::Function=NOOP, onstart::Function=NOOP)
    ChunkedTransferDecodeStream(source, maxsize, 0, 0, 0, ondone, onstart)
end

function Base.read(s::ChunkedTransferDecodeStream, ::Type{UInt8})
    # has the side affect of parsing the next chunksize if necessary
    if eof(s)
        throw(EOFError())
    end
    b = read(s.source, UInt8)
    s.nread += 1
    if s.nread >= s.maxsize
        throw(HTTPError(413; message="Unwilling to process a request entity with total size greater than $(s.maxsize)", shouldclose=true))
    end
    b
end

# function Base.readbytes!(s::ChunkedTransferDecodeStream, dest::Vector{UInt8}, count::Int64)
#     for i in 1:count
#         push!(dest, read(s, UInt8))
#     end
#
#     count
# end

function Base.eof(s::ChunkedTransferDecodeStream)
    if s.nread == 0
        s.onstart()
    end
    if eof(s.source)
        true
    else
        chunkpos = s.nread - s.chunkstart
        if s.nread == 0 || chunkpos == s.chunksize
            readchunksize(s)
        end

        eof(s.source) || s.chunksize == 0
    end
end

function Base.isopen(s::ChunkedTransferDecodeStream)
    isopen(s.source)
end

function Base.close(s::ChunkedTransferDecodeStream)
    nothing
end

function readchunksize(s::ChunkedTransferDecodeStream)
    source = s.source
    complete = false
    hexbytes = Vector{UInt8}()
    if isopen(source)
        b = read(s.source, UInt8)
        s.nread += 1
        while in(b, HEX) && isopen(source)
            if s.nread >= s.maxsize
                throw(HTTPError(413; message="Unwilling to process a request entity with total size greater than $(s.maxsize)", shouldclose=true))
            end
            push!(hexbytes, b)
            b = read(source, UInt8)
            s.nread += 1
        end

        if length(hexbytes) == 0
            throw(HTTPError(400; message="Chunked transfer encoding is missing chunksize at position $(s.nread) (found $(b))", shouldclose=true))
        end

        extcount = 0
        # TODO: We really don't care about any chunked encoding extensions at this point, but
        #       we may want to in the future...also we should validate that at least the syntax
        #       is correct for the extension.  Currently limiting the total extensions size to
        #       1024 bytes, but we may want to make this configurable.
        while b != CR && b != LF && isopen(source)
            extcount += 1
            if extcount > 1024
                throw(HTTPError(400; message="Unwilling to parse a chunked encoding extension greater than 1024 bytes", shouldclose=true))
            end
            if s.nread >= s.maxsize
                throw(HTTPError(413; message="Unwilling to process a request entity with total size greater than $(s.maxsize)", shouldclose=true))
            end

            b = read(source, UInt8)
            s.nread += 1
        end

        if b == LF
            throw(HTTPError(400; message="Invalid line ending in chunked encoding after chunk size", shouldclose=true))
        end

        if s.nread >= s.maxsize
            throw(HTTPError(413; message="Unwilling to process a request entity with total size greater than $(s.maxsize)", shouldclose=true))
        end

        if b == CR && isopen(source)
            b = read(source, UInt8)
            s.nread += 1

            foundend = true
            if s.nread >= s.maxsize
                throw(HTTPError(413; message="Unwilling to process a request entity with total size greater than $(s.maxsize)", shouldclose=true))
            end
            if b != LF
                throw(HTTPError(400; message="Invalid line ending in chunked encoding after chunk size", shouldclose=true))
            end
        end
    end

    if !foundend
        # source must have closed unexpectedly
        throw(HTTPError(400; message="Unexpected end of entity data while decoding chunked transfer encoding", shouldclose=true))
    end

    try
        s.chunksize = parse(UInt64, String(hexbytes), base = 16)
    catch e
        if isa(e, OverflowError)
            throw(HTTPError(413; message="Chunk size too large to be parsed", shouldclose=true))
        else
            throw(HTTPError(500; message="An unexpected error occurred", shouldclose=true, data=Dict{AbstractString, Any}("error" => e)))
        end
    end

    s.chunkstart = s.nread

    # notify parser so trailers can be consumed
    if s.chunksize == 0
        s.ondone()
    end

    s.chunksize
end


text:
import Blosc

# port of https://github.com/Blosc/c-blosc/blob/3a668dcc9f61ad22b5c0a0ab45fe8dad387277fd/hdf5/blosc_filter.c (copyright 2010 Francesc Alted, license: MIT/expat)

const FILTER_BLOSC_VERSION = 2
const FILTER_BLOSC = 32001 # Filter ID registered with the HDF Group for Blosc
const blosc_name = "blosc"

const blosc_flags = Ref{Cuint}()
const blosc_values = Vector{Cuint}(undef,8)
const blosc_nelements = Ref{Csize_t}(length(blosc_values))
const blosc_chunkdims = Vector{API.hsize_t}(undef,32)

function blosc_set_local(dcpl::API.hid_t, htype::API.hid_t, space::API.hid_t)
    API.h5p_get_filter_by_id(dcpl, FILTER_BLOSC, blosc_flags, blosc_nelements, blosc_values, 0, C_NULL, C_NULL)
    flags = blosc_flags[]

    nelements = max(blosc_nelements[], 4) # First 4 slots reserved

    # Set Blosc info in first two slots
    blosc_values[1] = FILTER_BLOSC_VERSION
    blosc_values[2] = Blosc.VERSION_FORMAT

    ndims = API.h5p_get_chunk(dcpl, 32, blosc_chunkdims)
    chunksize = prod(resize!(blosc_chunkdims, ndims))
    if ndims < 0 || ndims > 32 || chunksize > Blosc.MAX_BUFFERSIZE
        return API.herr_t(-1)
    end

    htypesize = API.h5t_get_size(htype)
    if API.h5t_get_class(htype) == API.H5T_ARRAY
        hsuper = API.h5t_get_super(htype)
        basetypesize = API.h5t_get_size(hsuper)
        API.h5t_close(hsuper)
    else
        basetypesize = htypesize
    end

    # Limit large typesizes (they are pretty inefficient to shuffle
    # and, in addition, Blosc does not handle typesizes larger than
    # blocksizes).
    if basetypesize > Blosc.MAX_TYPESIZE
        basetypesize = 1
    end
    blosc_values[3] = basetypesize
    blosc_values[4] = chunksize * htypesize # size of the chunk

    API.h5p_modify_filter(dcpl, FILTER_BLOSC, flags, nelements, blosc_values)

    return API.herr_t(1)
end

function blosc_filter(flags::Cuint, cd_nelmts::Csize_t,
                      cd_values::Ptr{Cuint}, nbytes::Csize_t,
                      buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Cvoid}})
    typesize = unsafe_load(cd_values, 3) # The datatype size
    outbuf_size = unsafe_load(cd_values, 4)
    # Compression level:
    clevel = cd_nelmts >= 5 ? unsafe_load(cd_values, 5) : Cuint(5)
    # Do shuffle:
    doshuffle = cd_nelmts >= 6 ? unsafe_load(cd_values, 6) != 0 : true
    # to do: set compressor based on compcode in unsafe_load(cd_values, 7)?

    if (flags & API.H5Z_FLAG_REVERSE) == 0 # compressing
        # Allocate an output buffer exactly as long as the input data; if
        # the result is larger, we simply return 0. The filter is flagged
        # as optional, so HDF5 marks the chunk as uncompressed and proceeds.
        outbuf_size = unsafe_load(buf_size)
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return Csize_t(0)
        status = Blosc.blosc_compress(clevel, doshuffle, typesize, nbytes,
                                      unsafe_load(buf), outbuf, nbytes)
        status < 0 && (Libc.free(outbuf); return Csize_t(0))
    else # decompressing
        # Extract the exact outbuf_size from the buffer header.
        #
        # NOTE: the guess value got from "cd_values" corresponds to the
        # uncompressed chunk size but it should not be used in a general
        # cases since other filters in the pipeline can modify the buffer
        # size.
        outbuf_size, cbytes, blocksize = Blosc.cbuffer_sizes(unsafe_load(buf))
        outbuf = Libc.malloc(outbuf_size)
        outbuf == C_NULL && return Csize_t(0)
        status = Blosc.blosc_decompress(unsafe_load(buf), outbuf, outbuf_size)
        status <= 0 && (Libc.free(outbuf); return Csize_t(0))
    end

    if status != 0
        Libc.free(unsafe_load(buf))
        unsafe_store!(buf, outbuf)
        unsafe_store!(buf_size, outbuf_size)
        return Csize_t(status) # size of compressed/decompressed data
    end
    Libc.free(outbuf); return Csize_t(0)
end

# register the Blosc filter function with HDF5
function register_blosc()
    c_blosc_set_local = @cfunction(blosc_set_local, API.API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    c_blosc_filter = @cfunction(blosc_filter, Csize_t,
                                (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                 Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))
    API.h5z_register(API.H5Z_class_t(API.H5Z_CLASS_T_VERS, FILTER_BLOSC, 1, 1, pointer(blosc_name), C_NULL, c_blosc_set_local, c_blosc_filter))

    return nothing
end

const set_blosc_values = Cuint[0,0,0,0,5,1,0]
function set_blosc(p::Properties, level::Integer=5)
    0 <= level <= 9 || throw(ArgumentError("blosc compression $level not in [0,9]"))
    set_blosc_values[5] = level
    API.h5p_set_filter(p.id, FILTER_BLOSC, API.H5Z_FLAG_OPTIONAL, length(set_blosc_values), set_blosc_values)

    return nothing
end


text:
export hypercubicI, copyvertex, copytensor, copygate
export Cylinder, rand_maskedsquare

"""periodic in y direction."""
struct Cylinder <: Viznet.AbstractSquareLattice
    Nx::Int
    Ny::Int
end

Base.size(lt::Cylinder) = (lt.Nx, lt.Ny)
Base.size(lt::Cylinder, i::Int) = size(lt, i)
Viznet.unit(lt::Cylinder) = 1/(max(size(lt)...))
Viznet.vertices(lt::Cylinder) = 1:lt.Nx * lt.Ny
function Base.getindex(lt::Cylinder, i::Real, j::Real)
    step = unit(lt)
    (i-0.5)*step, (j-0.5)*step
end

function solve(::Type{TT}, sg::AbstractSpinglass{LT}; usecuda=false) where {LT<:Cylinder,TT}
    # Yao gates
    lt = sg.lattice
    Lx, Ly = lt.Nx, lt.Ny
    reg = _init_reg(TT, lt.Ny, Val(usecuda))
    k = 0

    for i=1:Lx
        println("Layer $i/$Lx")
        i!=1 && for j=1:Ly
            k += 1
            reg |> put(Ly, j=>Ghb(bondtensor(TT, sg, k)))
        end
        for j=1:Ly
            reg |> put(Ly, j=>Gh(vertextensor(TT, sg, (i-1)*Ly + j)))
        end
        for j=1:Ly
            k += 1
            reg |> put(Ly, (j,mod1(j+1, Ly))=>Gvb(bondtensor(TT, sg, k)))
        end
    end
    sum(state(reg))
end

Viznet.bonds(lt::Cylinder) = sgbonds(lt)

function sgvertices(lt::Cylinder)
    sgvertices(SquareLattice(size(lt)...))
end

function sgbonds(lt::Cylinder)
    edges = Tuple{Int, Int}[]
    LI = LinearIndices(size(lt))
    function trypush!(i, j)
        push!(edges, (LI[i...], LI[j...]))
    end
    Lx, Ly = size(lt)
    for i=1:Lx
        for j=1:Ly
            trypush!([i,j], [i,mod1(j+1,Ly)])
        end
        (i!=Lx) && for j=1:Ly
            trypush!([i,j], [i+1,j])
        end
    end
    edges
end

function assign_Js_hs(lt::Cylinder, grad_Js::AbstractVector{T}, grad_hs) where T
    grid = zeros(size(lt.mask))
    vorder = sgvertices(lt)
    for (i, hg) in enumerate(grad_hs)
        grid[vorder[i]] = hg
    end
    for ((i, j), Jg) in zip(sgbonds(lt), grad_Js)
        assign_one!(grid, i, j, Jg)
    end
    return grid
end

regsize(lt::Cylinder) = size(lt,2)
function _init_reg(::Type{T}, lt::Cylinder, ::Val{:false}) where T
    nbit = regsize(lt)
    state = ones(T, 1<<nbit)
    ArrayReg(state)
end


text:
###################################################################################################
#  Extremely manually keyed-in plaquette sign structure required for Hamiltonian
#
# (1) plaquette(uc,b)
#
#       WORKS ONLY FOR (10,3)a !!!
#
#
#    -> carries the signs for bonds between nearest neighbors
#    -> VERY manual/ Brute-force
#    -> The way the signs have been keyed in: Print the list of all bonds in a unitcell of (10,3)a
#    -> Refer to relevant paper for the signs a bond between each pair of nearest-sites carries
#    -> See which pair-of-sites it corresponds to in the bond-list printed
#    -> Make matrix of signs "plaquette_convention" accordingly
#
#
#
#  (2) plaquetteNNN
#
#      -> Sign structure for bonds between next nearest neighbors
#      -> Works for sign-structure of any unitcell once the plaquette(uc,b) is set-up properly
#      -> Not as brute-force as the plaquette(uc,b) function
#
#######################################################################################################

function plaquette(
        uc :: U,
        b  :: B
    ) where {L,N,S,B<:AbstractBond{L,N},U<:AbstractUnitcell{S,B}}

   bond_index =  findall(bd -> bd == b,bonds(uc))

  # sign-structure ONLY FOR (10,3)a!!!

  plaquette_convention =   -[-1 1 1 1 -1 -1 1 -1 1 -1 -1 1]

    return( plaquette_convention[bond_index])
end

#sign-structure for next-nearest neighbors
function plaquetteNNN(
        uc :: U,
        b  :: B
    ) where {L,N,S,B<:AbstractBond{L,N},U<:AbstractUnitcell{S,B}}

    #b is the bond between pair of next-nearest neighbors
    i_from = from(b)
    i_to = to(b)
    z_eff = 0

    # to find the site lying in between the pair of NNN neighbors i_from and i_to

    # list of all bonds in NN-list in which i_from is starting site and i_to is ending site
    from_bonds = findall(bd -> from(bd)==i_from, bonds(uc))
    to_bonds = findall(bd -> to(bd)==i_to, bonds(uc))

    for i in 1:length(from_bonds)
        current_frombond = bonds(uc)[from_bonds[i]]
        for j in 1:length(to_bonds)
            current_tobond = bonds(uc)[to_bonds[j]]

            #adding up position vectors
            uc_add = wrap(current_tobond) .+ wrap(current_frombond)
            uc_sub = wrap(current_tobond) .- wrap(current_frombond)

            if(to(current_frombond)==from(current_tobond))

              #check if uc_add matches position vector of NNN-pair that we have
              if(uc_add==wrap(b))

                  #extract signs of the two bonds involved
                  z1 = plaquette(uc,current_frombond)
                  z2 = plaquette(uc,current_tobond)

                  #to extract the type of bond (x,y or z)
                  dir_index = findall(bd -> from(bd)== i_from&&to(bd)==i_to, bonds(uc))
                   sign_bond = bonds(uc)[dir_index][1]

                   #enter the info of bond-labels to calculate handedness
                   perm = [label(current_frombond),label(current_tobond),label(sign_bond)]

                   #effective/resultant sign of NNN bond
                   z_eff = levicivita(perm).*z1.*z2

            end
        end
    end
 end
    return z_eff[1]
end


text:
export Pow

"""
`Pow([domainType=Float64::Type,] dim_in::Tuple)`

Elementwise power `p` non-linear operator with input dimensions `dim_in`.

"""
struct Pow{T,N,I<:Real} <: NonLinearOperator
	dim::NTuple{N,Int}
    p::I
end

function Pow(DomainType::Type, DomainDim::NTuple{N,Int}, p::I) where {N, I <: Real}
	Pow{DomainType, N, I}(DomainDim, p)
end

Pow(DomainDim::NTuple{N,Int}, p::I) where {N, I <: Real} = Pow{Float64,N,I}(DomainDim,p)

function mul!(y::AbstractArray{T,N}, L::Pow{T,N,I}, x::AbstractArray{T,N}) where {T,N,I}
	y .= x.^L.p
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{Pow{T, N, I},TT}},
              b::AbstractArray) where {T, N, I, TT <: AbstractArray{T,N}}
    L = J.A
    y .= conj.(L.A.p.*(L.x).^(L.A.p-1)).*b
end

fun_name(L::Pow) = "『"

size(L::Pow) = (L.dim, L.dim)

domainType(L::Pow{T,N}) where {T,N} = T
codomainType(L::Pow{T,N}) where {T,N} = T


text:
#!/usr/bin/julia

# This was written to merge the output of `samtools depth`. I had executed
# samtools depth on six BAM files, each the result of mapping a sample's reads
# to a common reference (as when comparing sample-specific transcription
# or transcription levels). Currently it goes through six files at the same time,
# averaging coverage per transcript and saving the results to another file.

using DataFrames

namen = ASCIIString[]

# Files to be read from
in85 = open("85_depth", "r")
in86 = open("86_depth", "r")
in87 = open("87_depth", "r")
in88 = open("88_depth", "r")
in89 = open("89_depth", "r")
in90 = open("90_depth", "r")

# One Float64 and one Int array for every file. Avgs are to contain the average
# coverage per transcript. tempcov collects per-base coverages (for one
# transcript at a time) that will be averaged once a whole transcript has been
# parsed.

avgs85 = Float64[]
tempcov85 = Int[]

avgs86 = Float64[]
tempcov86 = Int[]

avgs87 = Float64[]
tempcov87 = Int[]

avgs88 = Float64[]
tempcov88 = Int[]

avgs89 = Float64[]
tempcov89 = Int[]

avgs90 = Float64[]
tempcov90 = Int[]

# Read first line of every file
line85 = readline(in85)
line86 = readline(in86)
line87 = readline(in87)
line88 = readline(in88)
line89 = readline(in89)
line90 = readline(in90)

# currnam is the name of the transcript whose coverage is being calculated a
currnam = split(line85)[1]

# Push the first base's coverage into tempcov (for every file)
push!(tempcov85, parse(Int, split(line85)[3]))
push!(tempcov86, parse(Int, split(line86)[3]))
push!(tempcov87, parse(Int, split(line87)[3]))
push!(tempcov88, parse(Int, split(line88)[3]))
push!(tempcov89, parse(Int, split(line89)[3]))
push!(tempcov90, parse(Int, split(line90)[3]))

# count will count the lenth of each transcrip (necessary for calculation
# of the mean coverage).
count = 1

#dcount = 1 #for debug

# Main loop of this script

for line85 in eachline(in85)

# dcount and the if stuff is for debug only
#dcount += 1
#if dcount == 10000
#	break
#end

# newnam - transcript name of the next line. If it differs from the last line's
# transcript name (currnam), we have reached a new transcript and we can
# calculate last one's coverage and push it into avgs. Also count we be set back
# to 1 and tempcov will be emptied/
	newnam = split(line85)[1]

	if newnam != currnam

		push!(namen, currnam)
		push!(avgs85, sum(tempcov85)/count)
		push!(avgs86, sum(tempcov86)/count)
		push!(avgs87, sum(tempcov87)/count)
		push!(avgs88, sum(tempcov88)/count)
		push!(avgs89, sum(tempcov89)/count)
		push!(avgs90, sum(tempcov90)/count)

		println(currnam)
		currnam = newnam
		count = 1
		tempcov85 = Int[]
		tempcov86 = Int[]
		tempcov87 = Int[]
		tempcov88 = Int[]
		tempcov89 = Int[]
		tempcov90 = Int[]

	end

	line86 = readline(in86)
	line87 = readline(in87)
	line88 = readline(in88)
	line89 = readline(in89)
	line90 = readline(in90)

	push!(tempcov85, parse(Int, split(line85)[3]))
	push!(tempcov86, parse(Int, split(line86)[3]))
	push!(tempcov87, parse(Int, split(line87)[3]))
	push!(tempcov88, parse(Int, split(line88)[3]))
	push!(tempcov89, parse(Int, split(line89)[3]))
	push!(tempcov90, parse(Int, split(line90)[3]))

	count += 1

	push!(tempcov85, parse(Int, split(line85)[3]))

end

# Once we have read through the whole file, we calculate the last transcript's
# coverage and push it to avgs for every file. Then, we put all averages and
# names in a dataframe and dump to a file (very much like in R).
push!(namen, currnam)
push!(avgs85, sum(tempcov85)/count)
push!(avgs86, sum(tempcov86)/count)
push!(avgs87, sum(tempcov87)/count)
push!(avgs88, sum(tempcov88)/count)
push!(avgs89, sum(tempcov89)/count)
push!(avgs90, sum(tempcov90)/count)

df = DataFrame(names = namen, d85 = avgs85, d86 = avgs86, d87 = avgs87,
	d88 = avgs88, d89 = avgs89, d90 = avgs90)

writetable("coverages", df, separator = '\t')

# Close opened files.
close(in85)
close(in86)
close(in87)
close(in88)
close(in89)
close(in90)

max_stars_repo_name: hannesbecher/PodiSNP

max_issues_repo_name: hannesbecher/PodiSNP

max_forks_repo_name: hannesbecher/PodiSNP


text:
function generate_mvn_with_outliers(n_dim, n_observations,
    seed=123, normalized=true, incl_outliers=true)

    Random.seed!(seed)
    norm_distribution = MvNormal(zeros(n_dim), Matrix(1.0I, n_dim, n_dim))
    inliers = rand(norm_distribution, n_observations)
    tmp = [rand(MvNormal([x, y], Matrix(1.0I, 2, 2)), 2) for x in [4,-4] for y in [4,-4]]
    outliers = vcat(hcat(tmp...), zeros(n_dim - 2, 8))

    if incl_outliers
        x = hcat(inliers, outliers)
        labels = vcat(fill(:inlier, n_observations), fill(:outlier, 8))
    else
        x = inliers
        labels = vcat(fill(:inlier, n_observations))
    end
    if normalized
        x = mapslices(normalize, x, dims=2)
    end
    return (x, labels)
end


text:
# uniform random walk, with bounds check
function proposal1(current, tuning)
    lb_param_ub = [
        0.95    0.99   	0.999; 	    # beta
        0.0	    2      	5;		    # gam
        0    	0.9   	0.999;	    # rho1
        0.0001       0.02 	0.1;    # sigma1
        0    	0.7     0.999;      # rho2
        0.0001	    0.01  	0.1;    # sigma2
        6/24    8/24	9/24	    # nss
    ]
    lb = lb_param_ub[:,1]
    ub = lb_param_ub[:,3]
    trial = copy(current)
    if rand() > 0.0
        i = rand(1:size(current,1))
        trial[i] = current[i] + tuning[i].*randn()
    else
        trial = lb + (ub - lb).*rand(size(lb))
    end
    return trial
end

function proposal2(current, cholV)
    trial = copy(current)
    lb_param_ub = [
        0.95    0.99   	0.999; 	    # beta
        0.0	    2      	5;		    # gam
        0    	0.9   	0.999;	    # rho1
        0.0001       0.02 	0.1;    # sigma1
        0    	0.7     0.999;      # rho2
        0.0001	    0.01  	0.1;    # sigma2
        6/24    8/24	9/24	    # nss
    ]
    lb = lb_param_ub[:,1]
    ub = lb_param_ub[:,3]
    if rand() > 0.0
        trial += cholV'*randn(size(trial))
    else
        trial = lb + (ub - lb).*rand(size(lb))
    end

end

function prior(theta)
    lb_param_ub = [
        0.95    0.99   	0.999; 	# beta
        0.0	    2      	5;		    # gam
        0    	0.9   	0.999;	    # rho1
        0.0001       0.02 	0.1;    # sigma1
        0    	0.7     0.999;     # rho2
        0.0001	    0.01  	0.1;    # sigma2
        6/24    8/24	9/24	    # nss
    ]
    lb = lb_param_ub[:,1]
    ub = lb_param_ub[:,3]
    a = 0.0
    if(all((theta .>= lb) .& (theta .<= ub)))
        a = 1.0
    end
    return a
end

function logL(theta, data)
    n = size(data,1)
    gs = DSGEmoments(theta, data)
    sigma = cov(gs)
    siginv = inv(sigma)
    ghat = mean(gs,dims=1)
    #lnL = - 0.5*(ghat*siginv*ghat')[1,1]
    lnL = log(sqrt(det(siginv))) - 0.5*n*(ghat*siginv*ghat')[1,1]
end

function logL_ABC(theta, data)
    n = size(data,1)
    gs = DSGEmoments(theta, data)
    Σ = cov(gs)
    cholΣ = chol(Σ)
    Σinv = inv(Σ)
    #Σinv = 1.0  # using the full asymptotic likelihood contracts too much
                # the true parameters can be outside the posterior's support
                # due to the small sample bias. Using a simple sum of squares
                # criterion doesn't have the problem
    ghat = mean(gs,1)
    #lnL = - 0.5*(ghat*Σinv*ghat')[1,1]
    lnL = log(sqrt(det(Σinv))) - 0.5*n*(ghat*Σinv*ghat')[1,1]
    ghat += randn(size(ghat))*cholΣ/sqrt(n) # used for nonparametric fitting
    return lnL, ghat

end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
"""
    get_atoms(dimensions, N, rₘᵢₙ; kwargs...)

`dimensions` : 1D/2D/3D == `1`/`2`/`3`
`N` : number of atoms
`rₘᵢₙ` : radius of exclusion
"""
function get_atoms(dimensions, N, rₘᵢₙ; kwargs...)
    new_atom = zeros(dimensions)

    get_single_atom = kwargs[:createFunction]
    r = get_single_atom(;kwargs...)
    if N == 1
        return r
    end
    temp_atom = zeros(eltype(r), dimensions)
    nValid = 1
    for iLoop in 1:1_000_000
        temp_atom = get_single_atom(;kwargs...)

        if is_atom_valid(temp_atom, r, rₘᵢₙ)
            nValid = nValid + 1
            r = hcat(r, temp_atom)
        end
        if nValid == N
            break
        end
    end
    if nValid < N
        @error "Could not generate all data after 10^6 interactions"
    end

    return r
end

"""
    Checks if the `new atom` has a `minimum distance` for `all the other atoms`
"""
function is_atom_valid(new_atom, r, rₘᵢₙ)
    A = r
    b = new_atom
    allDistances = get_Distance_A_to_b(A, b)
    return all(allDistances .≥ rₘᵢₙ)
end

function get_Distance_A_to_b(A, b)
    N = size(A,2)
    distanceAb = zeros(N)
    i = 1
    for c = eachcol(A)
        distanceAb[i] = Distances.evaluate(Euclidean(), c, b)
        i += 1
    end
    return distanceAb
end

"""
    get_rₘᵢₙ(ρ) = (ρ^(-2 / 3)) / π^2
"""
get_rₘᵢₙ(ρ) = (ρ^(-2 / 3)) / π^2  # AD HOC

"""
    get_pairwise_matrix(r)
r = matrix with `N` atoms positions

returns a `NxN` Float64 matrix, with zeros at diagonal
"""
function get_pairwise_matrix(r)
    r_matrix = transpose(r)
    R_jk = Distances.pairwise(Euclidean(), r_matrix, r_matrix; dims=1)
    R_jk[diagind(R_jk)] .= 0
    return R_jk
end

### -----SHAPE CONSTRUCTOR GIVEN MATRIX -------
function Atom(geometry::T, r::Matrix, kR::Union{Real, Integer}) where T <: Dimension
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(geometry, r, N, Float64(kR))
end

### --------------- CONVERSIONS ---------------
"""
    cube_inputs(N::Integer, ρ::Real)
"""
function cube_inputs(N::Integer, ρ::Real)
    kL = (N / ρ)^(1 / 3)
    return (N, kL)
end
"""
    sphere_inputs(N::Integer, ρ::Real)
"""
function sphere_inputs(N::Integer, ρ::Real)
    kR = (N / ((4π / 3) * ρ))^(1 / 3)
    return (N, kR)
end

### --------------- IO ---------------
"""
If atomic position is `atomic.r` of `size(N, dimensions)`.

Select one column, effectively as : `atoms.r[:, direction]`
"""
function select_atoms_axes(atoms, direction)
    return select_matrix_axes(atoms.r, direction)
end


text:
using Hexagons
using Test

convert(HexagonOffsetOddR, HexagonAxial(2, 4))

x, y = center(HexagonAxial(2, 3))
@test isapprox(x,7.06217782649107)
@test isapprox(y,5.5)

h = cube_round(23.5, 4.67)

collect(vertices(HexagonAxial(2, 3)))


text:
using Nabla
using Test, LinearAlgebra, Statistics, Random
using Distributions, BenchmarkTools, SpecialFunctions, DualNumbers

using Nabla: unbox, pos, tape, oneslike, zeroslike

# Helper function for comparing `Ref`s, since they don't compare equal under `==`
ref_equal(a::Ref{T}, b::Ref{T}) where {T} = a[] == b[]
ref_equal(a::Ref, b::Ref) = false

@testset "Core" begin
    include("core.jl")
    include("code_transformation/util.jl")
    include("code_transformation/differentiable.jl")
    include("sensitivity.jl")
end

@testset "Sensitivities" begin
    include("finite_differencing.jl")

    # Test sensitivities for the basics.
    include("sensitivities/indexing.jl")
    include("sensitivities/scalar.jl")
    include("sensitivities/array.jl")

    # Test sensitivities for functionals.
    @testset "Functional" begin
        include("sensitivities/functional/functional.jl")
        include("sensitivities/functional/reduce.jl")
        include("sensitivities/functional/reducedim.jl")
    end

    # Test sensitivities for linear algebra optimisations.
    @testset "Linear algebra" begin
        include("sensitivities/linalg/generic.jl")
        include("sensitivities/linalg/symmetric.jl")
        include("sensitivities/linalg/uniformscaling.jl")
        include("sensitivities/linalg/diagonal.jl")
        include("sensitivities/linalg/triangular.jl")
        include("sensitivities/linalg/strided.jl")
        include("sensitivities/linalg/blas.jl")

        @testset "Factorisations" begin
            include("sensitivities/linalg/factorization/cholesky.jl")
            include("sensitivities/linalg/factorization/svd.jl")
        end
    end

    include("checkpointing.jl")
end


text:
function initialize!(integrator,cache::Union{AB3ConstantCache,
                                             AB4ConstantCache,
                                             AB5ConstantCache,
                                             ABM32ConstantCache,
                                             ABM43ConstantCache,
                                             ABM54ConstantCache})

  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function initialize!(integrator,cache::Union{AB3Cache,
                                             AB4Cache,
                                             AB5Cache,
                                             ABM32Cache,
                                             ABM43Cache,
                                             ABM54Cache})

  @unpack fsalfirst,k = 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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::AB3ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2, k3 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 2
    cache.step += 1
    ttmp = t + (2/3)*dt
    ralk2 = f(uprev + (2/3)*dt*k1, p, ttmp)       #Ralston Method
    integrator.destats.nf += 1
    u = uprev + (dt/4)*(k1 + 3*ralk2)
    if cnt == 1
      k3 = k1
    else
      k2 = k1
    end
  else
    u  = uprev + (dt/12)*(23*k1 - 16*k2 + 5*k3)
    k3 = k2
    k2 = k1
  end
  cache.k2 = k2
  cache.k3 = k3
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::AB3Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,ralk2,k = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 2
    cache.step += 1
    ttmp = t + (2/3)*dt
    @.. tmp = uprev + (2/3)*dt*k1
    f(ralk2, tmp, p, ttmp)
    integrator.destats.nf += 1
    @.. u = uprev + (dt/4)*(k1 + 3*ralk2)        #Ralston Method
    if cnt == 1
      cache.k3 .= k1
    else
      cache.k2 .= k1
    end
  else
    @.. u  = uprev + (dt/12)*(23*k1 - 16*k2 + 5*k3)
    cache.k2,cache.k3 = k3,k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::ABM32ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step == 1
    cache.step += 1
    ttmp = t + (2/3)*dt
    ralk2 = f(uprev + (2/3)*dt*k1, p, ttmp)     #Ralston Method
    integrator.destats.nf += 1
    u = uprev + (dt/4)*(k1 + 3*ralk2)
    k2 = k1
  else
    perform_step!(integrator, AB3ConstantCache(k2,k3,cnt))
    k = integrator.fsallast
    u = uprev + (dt/12)*(5*k + 8*k1 - k2)
    k3 = k2
    k2 = k1
  end
  cache.k2 = k2
  cache.k3 = k3
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::ABM32Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,ralk2,k = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step == 1
    cache.step += 1
    ttmp = t + (2/3)*dt
    @.. tmp = uprev + (2/3)*dt*k1
    f(ralk2, tmp, p, ttmp)
    integrator.destats.nf += 1
    @.. u = uprev + (dt/4)*(k1 + 3*ralk2)       #Ralston Method
    cache.k2 .= k1
  else
    if cnt == 2
      perform_step!(integrator, AB3Cache(u,uprev,fsalfirst,copy(k2),k3,ralk2,k,tmp,cnt))  #Here passing copy of k2, otherwise it will change in AB3()
    else
      perform_step!(integrator, AB3Cache(u,uprev,fsalfirst,k2,k3,ralk2,k,tmp,cnt))
    end
    k = integrator.fsallast
    @.. u = uprev + (dt/12)*(5*k + 8*k1 - k2)
    cache.k2,cache.k3 = k3,k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::AB4ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3,k4 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 3
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    k2 = f(uprev + halfdt*k1, p, ttmp)
    k3 = f(uprev + halfdt*k2, p, ttmp)
    k4 = f(uprev + dt*k3, p, t+dt)
    integrator.destats.nf += 3
    u = uprev + (dt/6)*(2*(k2 + k3) + (k1+k4))   #RK4
    if cnt == 1
      cache.k4 = k1
    elseif cnt == 2
      cache.k3 = k1
    else
      cache.k2 = k1
    end
  else
    u  = uprev + (dt/24)*(55*k1 - 59*k2 + 37*k3 - 9*k4)
    cache.k4 = k3
    cache.k3 = k2
    cache.k2 = k1
  end
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::AB4Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,k4,ralk2,k,t2,t3,t4 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 3
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    @.. tmp = uprev + halfdt*k1
    f(t2,tmp,p,ttmp)
    @.. tmp = uprev + halfdt*t2
    f(t3,tmp,p,ttmp)
    @.. tmp = uprev + dt*t3
    f(t4,tmp,p,t+dt)
    integrator.destats.nf += 3
    @.. u = uprev + (dt/6)*(2*(t2 + t3) + (k1 + t4))   #RK4
    if cnt == 1
      cache.k4 .= k1
    elseif cnt == 2
      cache.k3 .= k1
    else
      cache.k2 .= k1
    end
  else
    @.. u  = uprev + (dt/24)*(55*k1 - 59*k2 + 37*k3 - 9*k4)
    cache.k4, cache.k3 = k3, k4
    cache.k3 .= k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::ABM43ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3,k4 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 2
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    k2 = f(uprev + halfdt*k1, p, ttmp)
    k3 = f(uprev + halfdt*k2, p, ttmp)
    k4 = f(uprev + dt*k3, p, t+dt)
    integrator.destats.nf += 3
    u = uprev + (dt/6)*(2*(k2 + k3) + (k1+k4))   #RK4
    if cnt == 1
      cache.k3 = k1
    else
      cache.k2 = k1
    end
  else
    perform_step!(integrator, AB4ConstantCache(k2,k3,k4,cnt))
    k = integrator.fsallast
    u = uprev + (dt/24)*(9*k + 19*k1 - 5*k2 + k3)
    cache.k4 = k3
    cache.k3 = k2
    cache.k2 = k1
  end
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::ABM43Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,k4,ralk2,k,t2,t3,t4,t5,t6,t7 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 2
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    @.. tmp = uprev + halfdt*k1
    f(t2,tmp,p,ttmp)
    @.. tmp = uprev + halfdt*t2
    f(t3,tmp,p,ttmp)
    @.. tmp = uprev + dt*t3
    f(t4,tmp,p,t+dt)
    integrator.destats.nf += 3
    @.. u = uprev + (dt/6)*(2*(t2 + t3) + (k1 + t4))   #RK4
    if cnt == 1
      cache.k3 .= k1
    else
      cache.k2 .= k1
    end
  else
    t2 .= k2
    t3 .= k3
    t4 .= k4
    perform_step!(integrator, AB4Cache(u,uprev,fsalfirst,t2,t3,t4,ralk2,k,tmp,t5,t6,t7,cnt))
    k = integrator.fsallast
    @.. u = uprev + (dt/24)*(9*k + 19*k1 - 5*k2 + k3)
    cache.k4, cache.k3 = k3, k4
    cache.k3 .= k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::AB5ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3,k4,k5 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  if cache.step <= 4
    cnt = cache.step
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    k2 = f(uprev + halfdt*k1, p, ttmp)
    k3 = f(uprev + halfdt*k2, p, ttmp)
    k4 = f(uprev + dt*k3, p, t+dt)
    integrator.destats.nf += 3
    u = uprev + (dt/6)*(2*(k2 + k3) + (k1+k4))  #RK4
    if cnt == 1
      cache.k5 = k1
    elseif cnt == 2
      cache.k4 = k1
    elseif cnt == 3
      cache.k3 = k1
    else
      cache.k2 = k1
    end
  else
    u  = uprev + (dt/720)*(1901*k1 - 2774*k2 + 2616*k3 - 1274*k4 + 251*k5)
    cache.k5 = k4
    cache.k4 = k3
    cache.k3 = k2
    cache.k2 = k1
  end
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::AB5Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,k4,k5,k,t2,t3,t4 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 4
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    @.. tmp = uprev + halfdt*k1
    f(t2,tmp,p,ttmp)
    @.. tmp = uprev + halfdt*t2
    f(t3,tmp,p,ttmp)
    @.. tmp = uprev + dt*t3
    f(t4,tmp,p,t+dt)
    integrator.destats.nf += 1
    @.. u = uprev + (dt/6)*(2*(t2 + t3) + (k1 + t4))   #RK4
    if cnt == 1
      cache.k5 .= k1
    elseif cnt == 2
      cache.k4 .= k1
    elseif cnt == 3
      cache.k3 .= k1
    else
      cache.k2 .= k1
    end
  else
    @.. u  = uprev + (dt/720)*(1901*k1 - 2774*k2 + 2616*k3 - 1274*k4 + 251*k5)
    cache.k5, cache.k4 = k4, k5
    cache.k4 .= k3
    cache.k3 .= k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::ABM54ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,k3,k4,k5 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 3
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    k2 = f(uprev + halfdt*k1, p, ttmp)
    k3 = f(uprev + halfdt*k2, p, ttmp)
    k4 = f(uprev + dt*k3, p, t+dt)
    integrator.destats.nf += 3
    u = uprev + (dt/6)*(2*(k2 + k3) + (k1+k4))   #RK4
    if cnt == 1
      cache.k4 = k1
    elseif cnt == 2
      cache.k3 = k1
    else
      cache.k2 = k1
    end
  else
    perform_step!(integrator, AB5ConstantCache(k2,k3,k4,k5,cnt))
    k = integrator.fsallast
    u = uprev + (dt/720)*(251*k + 646*k1 - 264*k2 + 106*k3 - 19*k4)
    cache.k5 = k4
    cache.k4 = k3
    cache.k3 = k2
    cache.k2 = k1
  end
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::ABM54Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack tmp,fsalfirst,k2,k3,k4,k5,k,t2,t3,t4,t5,t6,t7,t8 = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  cnt = cache.step
  if cache.step <= 3
    cache.step += 1
    halfdt = dt/2
    ttmp = t+halfdt
    @.. tmp = uprev + halfdt*k1
    f(t2,tmp,p,ttmp)
    @.. tmp = uprev + halfdt*t2
    f(t3,tmp,p,ttmp)
    @.. tmp = uprev + dt*t3
    f(t4,tmp,p,t+dt)
    @.. u = uprev + (dt/6)*(2*(t2 + t3) + (k1 + t4))   #RK4
    integrator.destats.nf += 3
    if cnt == 1
      cache.k4 .= k1
    elseif cnt == 2
      cache.k3 .= k1
    else
      cache.k2 .= k1
    end
  else
    t2 .= k2
    t3 .= k3
    t4 .= k4
    t5 .= k5
    perform_step!(integrator, AB5Cache(u,uprev,fsalfirst,t2,t3,t4,t5,k,tmp,t6,t7,t8,cnt))
    k = integrator.fsallast
    @.. u = uprev + (dt/720)*(251*k + 646*k1 - 264*k2 + 106*k3 - 19*k4)
    cache.k5, cache.k4 = k4, k5
    cache.k4 .= k3
    cache.k3 .= k2
    cache.k2 .= k1
  end
  f(k, u, p, t+dt)
  integrator.destats.nf += 1
end

# Variable Step Size Multistep Methods

function initialize!(integrator,cache::VCAB3ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCAB3ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕstar_n,ϕstar_nm1,order,tab = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[3]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2
    perform_step!(integrator, tab)
  else
    g_coefs!(cache,k)
    u = uprev
    for i = 1:k
        u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      utilde = g[k] * ϕstar_n[k]      # Using lower order AB from subset of coefficients
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:2
          dts[i] = dts[i+1]
        end
        dts[3] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCAB3Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCAB3Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k4,dts,g,ϕstar_n,ϕstar_nm1,order,atmp,utilde,bs3cache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[3]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache, k1, k)
  if k == 1 || k == 2
    perform_step!(integrator, bs3cache)
    @unpack k4 = bs3cache
    integrator.fsallast .= k4
  else
    g_coefs!(cache, k)
    @.. u = uprev
    for i = 1:k
      @.. u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      @.. utilde = g[k] * ϕstar_n[k]    # Using lower order AB from subset of coefficients
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:2
          dts[i] = dts[i+1]
        end
        dts[3] = tmp
        return nothing
      end
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCAB4ConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCAB4ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕstar_n,ϕstar_nm1,order,rk4constcache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[4]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3
    perform_step!(integrator, rk4constcache)
  else
    g_coefs!(cache,k)
    u = uprev
    for i = 1:k
        u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      utilde = g[k] * ϕstar_n[k]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:3
          dts[i] = dts[i+1]
        end
        dts[4] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCAB4Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCAB4Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k4,dts,g,ϕ_n,ϕstar_n,ϕstar_nm1,order,atmp,utilde,rk4cache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[4]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3
    rk4cache.fsalfirst .= k1
    perform_step!(integrator, rk4cache)
    integrator.fsallast .= rk4cache.k
  else
    g_coefs!(cache,k)
    @.. u = uprev
    for i = 1:k
      @.. u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      @.. utilde = g[k] * ϕstar_n[k]
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:3
          dts[i] = dts[i+1]
        end
        dts[4] = tmp
        return nothing
      end
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

# VCAB5

function initialize!(integrator,cache::VCAB5ConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCAB5ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕstar_n,ϕstar_nm1,order,rk4constcache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[5]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 4
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[5] = dts[4]
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3 || k == 4
    perform_step!(integrator, rk4constcache)
  else
    g_coefs!(cache,k)
    u = uprev
    for i = 1:k
        u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      utilde = g[k] * ϕstar_n[k]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:4
          dts[i] = dts[i+1]
        end
        dts[5] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCAB5Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCAB5Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k4,dts,g,ϕ_n,ϕstar_n,ϕstar_nm1,order,atmp,utilde,rk4cache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[5]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 4
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[5] = dts[4]
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3 || k == 4
    rk4cache.fsalfirst .= k1
    perform_step!(integrator, rk4cache)
    integrator.fsallast .= rk4cache.k
  else
    g_coefs!(cache,k)
    @.. u = uprev
    for i = 1:k
      @.. u += g[i] * ϕstar_n[i]
    end
    if integrator.opts.adaptive
      @.. utilde = g[k] * ϕstar_n[k]
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:4
          dts[i] = dts[i+1]
        end
        dts[5] = tmp
        return nothing
      end
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

# VCABM3

function initialize!(integrator,cache::VCABM3ConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCABM3ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,tab = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[3]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2
    perform_step!(integrator, tab)
  else
    g_coefs!(cache,k+1)
    u = uprev
    for i = 1:(k-1)
        u += g[i] * ϕstar_n[i]
    end
    du_np1 = f(u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, du_np1, k+1)
    u += g[end-1] * ϕ_np1[end-1]
    if integrator.opts.adaptive
      utilde = (g[end] - g[end-1]) * ϕ_np1[end]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:2
          dts[i] = dts[i+1]
        end
        dts[3] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCABM3Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCABM3Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k4,dts,g,ϕstar_n,ϕ_np1,ϕstar_nm1,order,atmp,utilde,bs3cache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[3]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache, k1, k)
  if k == 1 || k == 2
    perform_step!(integrator, bs3cache)
    @unpack k4 = bs3cache
    integrator.fsallast .= k4
  else
    g_coefs!(cache, k+1)
    @.. u = uprev
    for i = 1:(k-1)
      @.. u += g[i] * ϕstar_n[i]
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, k4, k+1)
    @.. u += g[end-1] * ϕ_np1[end-1]
    if integrator.opts.adaptive
      @.. utilde = (g[end] - g[end-1]) * ϕ_np1[end]
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:2
          dts[i] = dts[i+1]
        end
        dts[3] = tmp
        return nothing
      end
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

# VCABM4

function initialize!(integrator,cache::VCABM4ConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCABM4ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,rk4constcache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[4]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3
    perform_step!(integrator, rk4constcache)
  else
    g_coefs!(cache,k+1)
    u = uprev
    for i = 1:(k-1)
        u += g[i] * ϕstar_n[i]
    end
    du_np1 = f(u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, du_np1, k+1)
    u += g[end-1] * ϕ_np1[end-1]
    if integrator.opts.adaptive
      utilde = (g[end] - g[end-1]) * ϕ_np1[end]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:3
          dts[i] = dts[i+1]
        end
        dts[4] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCABM4Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCABM4Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k4,dts,g,ϕstar_n,ϕ_np1,ϕstar_nm1,order,atmp,utilde,rk4cache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[4]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache, k1, k)
  if k == 1 || k == 2 || k == 3
    rk4cache.fsalfirst .= k1
    perform_step!(integrator, rk4cache)
    integrator.fsallast .= rk4cache.k
  else
    g_coefs!(cache, k+1)
    @.. u = uprev
    for i = 1:(k-1)
      @.. u += g[i] * ϕstar_n[i]
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, k4, k+1)
    @.. u += g[end-1] * ϕ_np1[end-1]
    if integrator.opts.adaptive
      @.. utilde = (g[end] - g[end-1]) * ϕ_np1[end]
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:3
          dts[i] = dts[i+1]
        end
        dts[4] = tmp
        return nothing
      end
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

# VCABM5

function initialize!(integrator,cache::VCABM5ConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCABM5ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,rk4constcache = cache
  k1 = integrator.fsalfirst
  if integrator.u_modified
    cache.step = 1
  end
  k = cache.step
  tmp = dts[5]
  if k == 1
    dts[1] = dt
    cache.step += 1
  elseif k == 2
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 3
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  elseif k == 4
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
    cache.step += 1
  else
    dts[5] = dts[4]
    dts[4] = dts[3]
    dts[3] = dts[2]
    dts[2] = dts[1]
    dts[1] = dt
  end
  ϕ_and_ϕstar!(cache,k1,k)
  if k == 1 || k == 2 || k == 3 || k == 4
    perform_step!(integrator, rk4constcache)
  else
    g_coefs!(cache,k+1)
    u = uprev
    for i = 1:(k-1)
        u += g[i] * ϕstar_n[i]
    end
    du_np1 = f(u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, du_np1, k+1)
    u += g[end-1] * ϕ_np1[end-1]
    if integrator.opts.adaptive
      utilde = (g[end] - g[end-1]) * ϕ_np1[end]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:4
          dts[i] = dts[i+1]
        end
        dts[5] = tmp
        return nothing
      end
    end
    integrator.fsallast = f(u, p, t+dt)
    integrator.destats.nf += 1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
  end
  cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
end

function initialize!(integrator,cache::VCABM5Cache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCABM5Cache,repeat_step=false)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack k4,dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,atmp,utilde,rk4cache = cache
    k1 = integrator.fsalfirst
    if integrator.u_modified
      cache.step = 1
    end
    k = cache.step
    tmp = dts[5]
    if k == 5
      dts[5] = dts[4]
      dts[4] = dts[3]
      dts[3] = dts[2]
      dts[2] = dts[1]
      dts[1] = dt
    elseif k == 1
      dts[1] = dt
      cache.step += 1
    elseif k == 2
      dts[2] = dts[1]
      dts[1] = dt
      cache.step += 1
    elseif k == 3
      dts[3] = dts[2]
      dts[2] = dts[1]
      dts[1] = dt
      cache.step += 1
    elseif k == 4
      dts[4] = dts[3]
      dts[3] = dts[2]
      dts[2] = dts[1]
      dts[1] = dt
      cache.step += 1
    end
    ϕ_and_ϕstar!(cache, k1, 5)
    if k <= 4
      rk4cache.fsalfirst .= k1
      perform_step!(integrator, rk4cache)
      integrator.fsallast .= rk4cache.k
    else
      g_coefs!(cache, 6)
      @.. u = muladd(g[1], ϕstar_n[1], uprev)
      @.. u = muladd(g[2], ϕstar_n[2], u)
      @.. u = muladd(g[3], ϕstar_n[3], u)
      @.. u = muladd(g[4], ϕstar_n[4], u)
      f(k4,u,p,t+dt)
      integrator.destats.nf += 1
      ϕ_np1!(cache, k4, 6)
      @.. u = muladd(g[6-1], ϕ_np1[6-1], u)
      if integrator.opts.adaptive
        @.. utilde = (g[6] - g[6-1]) * ϕ_np1[end]
        calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        integrator.EEst = integrator.opts.internalnorm(atmp,t)
        if integrator.EEst > one(integrator.EEst)
          dts[1] = dts[2]
          dts[2] = dts[3]
          dts[4] = dts[5]
          dts[5] = tmp
          return nothing
        end
      end
      cache.ϕstar_nm1[1] .= ϕstar_n[1]
      cache.ϕstar_nm1[2] .= ϕstar_n[2]
      cache.ϕstar_nm1[3] .= ϕstar_n[3]
      cache.ϕstar_nm1[4] .= ϕstar_n[4]
      f(k4,u,p,t+dt)
      integrator.destats.nf += 1
    end
    cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
    return nothing
  end # inbounds
end

# VCABM

function initialize!(integrator,cache::VCABMConstantCache)
  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)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VCABMConstantCache,repeat_step=false)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,max_order = cache
    k1 = integrator.fsalfirst
    step = integrator.iter
    k = order
    tmp = dts[13]
    for i = 12:-1:1
      dts[i+1] = dts[i]
    end
    dts[1] = dt
    ϕ_and_ϕstar!(cache,k1,k)
    g_coefs!(cache,k+1)
    u = muladd(g[1], ϕstar_n[1], uprev)
    for i = 2:k-1
      u = muladd(g[i], ϕstar_n[i], u)
    end
    du_np1 = f(u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, du_np1, k+1)
    u = muladd(g[k], ϕ_np1[k], u)
    if integrator.opts.adaptive
      utilde = (g[k+1]-g[k]) * ϕ_np1[k+1]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:12
          dts[i] = dts[i+1]
        end
        dts[13] = tmp
        return nothing
      end
      integrator.fsallast = f(u, p, t+dt)
      integrator.destats.nf += 1
      if step <= 4 || order < 3
        cache.order = min(order+1,3)
      else
        # utildem2 = dt * γstar[(k-2)+1] * ϕ_np1[k-1]
        utildem2 = (g[k-1]-g[k-2]) * ϕ_np1[k-1]
        # utildem1 = dt * γstar[(k-1)+1] * ϕ_np1[k]
        utildem1 = (g[k]-g[k-1]) * ϕ_np1[k]
        expand_ϕ_and_ϕstar!(cache, k+1)
        ϕ_np1!(cache, integrator.fsallast, k+2)
        utildep1 = dt * γstar[(k+1)+1] * ϕ_np1[k+2]
        atmpm2 = calculate_residuals(utildem2, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        atmpm1 = calculate_residuals(utildem1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        atmpp1 = calculate_residuals(utildep1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        errm2 = integrator.opts.internalnorm(atmpm2,t)
        errm1 = integrator.opts.internalnorm(atmpm1,t)
        errp1 = integrator.opts.internalnorm(atmpp1,t)
        if max(errm2,errm1) <= integrator.EEst
          cache.order = order - 1
        elseif errp1 < integrator.EEst
          cache.order = min(order+1,max_order)
          integrator.EEst = one(integrator.EEst)   # for keeping the stepsize constant in the next step
        end # if
      end # step <= 4
    end # integrator.opts.adaptive
    cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u
    return nothing
  end
end

function initialize!(integrator,cache::VCABMCache)
  @unpack fsalfirst,k4 = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k4
  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) # pre-start FSAL
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator,cache::VCABMCache,repeat_step=false)
  @inbounds begin
    @unpack t,dt,uprev,u,f,p = integrator
    @unpack k4,dts,g,ϕ_n,ϕ_np1,ϕstar_n,ϕstar_nm1,order,max_order,utilde,utildem2,utildem1,utildep1,atmp,atmpm1,atmpm2,atmpp1 = cache
    k1 = integrator.fsalfirst
    step = integrator.iter
    k = order
    tmp = dts[13]
    for i = 12:-1:1
      dts[i+1] = dts[i]
    end
    dts[1] = dt
    ϕ_and_ϕstar!(cache,k1,k)
    g_coefs!(cache,k+1)
    # unroll the predictor once
    @.. u = muladd(g[1], ϕstar_n[1], uprev)
    for i = 2:k-1
      @.. u = muladd(g[i], ϕstar_n[i], u)
    end
    f(k4,u,p,t+dt)
    integrator.destats.nf += 1
    ϕ_np1!(cache, k4, k+1)
    @.. u = muladd(g[k], ϕ_np1[k], u)
    if integrator.opts.adaptive
      @.. utilde = (g[k+1]-g[k]) * ϕ_np1[k+1]
      calculate_residuals!(atmp,utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst > one(integrator.EEst)
        for i = 1:12
          dts[i] = dts[i+1]
        end
        dts[13] = tmp
        return nothing
      end
      f(k4,u,p,t+dt)
      integrator.destats.nf += 1
      if step <= 4 || order < 3
        cache.order = min(order+1,3)
      else
        # @.. utildem2 = dt * γstar[(k-2)+1] * ϕ_np1[k-1]
        @.. utildem2 = (g[k-1]-g[k-2]) * ϕ_np1[k-1]
        # @.. utildem1 = dt * γstar[(k-1)+1] * ϕ_np1[k]
        @.. utildem1 = (g[k]-g[k-1]) * ϕ_np1[k]
        expand_ϕ_and_ϕstar!(cache, k+1)
        ϕ_np1!(cache, k4, k+2)
        @.. utildep1 = dt * γstar[(k+1)+1] * ϕ_np1[k+2]
        calculate_residuals!(atmpm2, utildem2, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        calculate_residuals!(atmpm1, utildem1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        calculate_residuals!(atmpp1, utildep1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        errm2 = integrator.opts.internalnorm(atmpm2,t)
        errm1 = integrator.opts.internalnorm(atmpm1,t)
        errp1 = integrator.opts.internalnorm(atmpp1,t)
        if max(errm2,errm1) <= integrator.EEst
          cache.order = order - 1
        elseif errp1 < integrator.EEst
          cache.order = min(order+1,max_order)
          integrator.EEst = one(integrator.EEst)    # for keeping the stepsize constant in the next step
        end
      end
    end
    cache.ϕstar_nm1, cache.ϕstar_n = ϕstar_n, ϕstar_nm1
    return nothing
  end # inbounds
end

# CNAB2

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

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

function perform_step!(integrator,cache::CNAB2ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,nlsolver = cache
  cnt = integrator.iter
  f1 = integrator.f.f1
  f2 = integrator.f.f2
  du₁ = f1(uprev, p, t)
  integrator.destats.nf += 1
  k1 = integrator.fsalfirst - du₁
  # Explicit part
  if cnt == 1
    tmp = uprev + dt * k1
  else
    tmp = uprev + dt * (3//2*k1 - 1//2*k2)
  end
  nlsolver.tmp = tmp
  # Implicit part
  # precalculations
  γ = 1//2
  γdt = γ*dt

  # initial guess
  zprev = dt*du₁
  nlsolver.z = z = zprev # Constant extrapolation

  nlsolver.tmp += γ*zprev
  markfirststage!(nlsolver)
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp + 1//2*z

  cache.k2 = k1
  integrator.fsallast = f1(u, p, t+dt) + f2(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::CNAB2Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
  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)
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::CNAB2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k1,k2,du₁,nlsolver = cache
  @unpack z,tmp = nlsolver
  @unpack f1 = f
  cnt = integrator.iter

  f1(du₁, uprev, p, t)
  integrator.destats.nf += 1
  @.. k1 = integrator.fsalfirst - du₁
  # Explicit part
  if cnt == 1
    @.. tmp = uprev + dt * k1
  else
    @.. tmp = uprev + dt * (3//2*k1 - 1//2*k2)
  end
  # Implicit part
  # precalculations
  γ = 1//2
  γdt = γ*dt

  # initial guess
  @.. z = dt*du₁
  @.. tmp += γ*z
  markfirststage!(nlsolver)
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = tmp + 1//2*z

  cache.k2 .= k1
  f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 1
end

# CNLF2

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

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

function perform_step!(integrator,cache::CNLF2ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k2,uprev2,nlsolver = cache
  cnt = integrator.iter
  f1 = integrator.f.f1
  f2 = integrator.f.f2
  du₁ = f1(uprev, p, t)
  integrator.destats.nf += 1
  # Explicit part
  if cnt == 1
    tmp = uprev + dt * (integrator.fsalfirst - du₁)
  else
    tmp = uprev2 + 2//1 * dt * (integrator.fsalfirst - du₁)
  end
  # Implicit part
  # precalculations
  γ = 1//1
  if cnt != 1
    tmp += γ*dt*k2
  end
  γdt = γ*dt
  nlsolver.tmp = tmp

  # initial guess
  zprev = dt*du₁
  nlsolver.z = z = zprev # Constant extrapolation

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

  cache.uprev2 = uprev
  cache.k2 = du₁
  integrator.fsallast = f1(u, p, t+dt) + f2(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::CNLF2Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
  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)
  integrator.destats.nf += 1
end

function perform_step!(integrator, cache::CNLF2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uprev2,k2,du₁,nlsolver = cache
  @unpack z,tmp = nlsolver
  @unpack f1 = f
  cnt = integrator.iter

  f1(du₁, uprev, p, t)
  integrator.destats.nf += 1
  # Explicit part
  if cnt == 1
    @.. tmp = uprev + dt * (integrator.fsalfirst - du₁)
  else
    @.. tmp = uprev2 + 2//1 * dt * (integrator.fsalfirst - du₁)
  end
  # Implicit part
  # precalculations
  γ = 1//1
  if cnt != 1
    @.. tmp += γ*dt*k2
  end
  γdt = γ*dt

  # initial guess
  @.. z = dt*du₁
  markfirststage!(nlsolver)
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = tmp + γ*z

  cache.uprev2 .= uprev
  cache.k2 .= du₁
  f(integrator.fsallast,u,p,t+dt)
  integrator.destats.nf += 1
end