text:
abstract type Smoother end
abstract type Sweep end
struct SymmetricSweep <: Sweep
end
struct ForwardSweep <: Sweep
end
struct BackwardSweep <: Sweep
end
struct GaussSeidel{S} <: Smoother
    sweep::S
    iter::Int
end
GaussSeidel(; iter = 1) = GaussSeidel(SymmetricSweep(), iter)
GaussSeidel(f::ForwardSweep) = GaussSeidel(f, 1)
GaussSeidel(b::BackwardSweep) = GaussSeidel(b, 1)
GaussSeidel(s::SymmetricSweep) = GaussSeidel(s, 1)

function (s::GaussSeidel{S})(A, x, b) where {S<:Sweep}
    for i in 1:s.iter
        if S === ForwardSweep || S === SymmetricSweep
            gs!(A, b, x, 1, 1, size(A, 1))
        end
        if S === BackwardSweep || S === SymmetricSweep
            gs!(A, b, x, size(A, 1), -1, 1)
        end
    end
end

function gs!(A, b, x, start, step, stop)
    n = size(A, 1)
    z = zero(eltype(A))
    @inbounds for col in 1:size(x, 2)
        for i in start:step:stop
            rsum = z
            d = z
            for j in nzrange(A, i)
                row = A.rowval[j]
                val = A.nzval[j]
                d = ifelse(i == row, val, d)
                rsum += ifelse(i == row, z, val * x[row, col])
            end
            x[i, col] = ifelse(d == 0, x[i, col], (b[i, col] - rsum) / d)
        end
    end
end

struct Jacobi{T,TX} <: Smoother
    ω::T
    temp::TX
    iter::Int
end
Jacobi(ω, x::TX; iter=1) where {T, TX<:AbstractArray{T}} = Jacobi{T,TX}(ω, similar(x), iter)
Jacobi(x::TX, ω=0.5; iter=1) where {T, TX<:AbstractArray{T}} = Jacobi{T,TX}(ω, similar(x), iter)

function (jacobi::Jacobi)(A, x, b)

    ω = jacobi.ω
    one = Base.one(eltype(A))
    temp = jacobi.temp
    z = zero(eltype(A))

    for i in 1:jacobi.iter
        @inbounds for col = 1:size(x, 2)
            for i = 1:size(A, 1)
                temp[i, col] = x[i, col]
            end

            for i = 1:size(A, 1)
                rsum = z
                diag = z

                for j in nzrange(A, i)
                    row = A.rowval[j]
                    val = A.nzval[j]

                    diag = ifelse(row == i, val, diag)
                    rsum += ifelse(row == i, z, val * temp[row, col])
                end

                xcand = (one - ω) * temp[i, col] + ω * ((b[i, col] - rsum) / diag)
                x[i, col] = ifelse(diag == 0, x[i, col], xcand)
            end
        end
    end
end

#=
using KissThreading: tmap!

struct ParallelJacobi{T,TX} <: Smoother
    ω::T
    temp::TX
end
ParallelJacobi(ω, x::TX) where {T, TX<:AbstractArray{T}} = ParallelJacobi{T,TX}(ω, similar(x))
ParallelJacobi(x::TX, ω = 0.5) where {T, TX<:AbstractArray{T}} = ParallelJacobi{T,TX}(ω, similar(x))

struct ParallelTemp{TX, TI}
    temp::TX
    col::TI
end
(ptemp::ParallelTemp)(i) = ptemp.temp[i, ptemp.col]

struct ParallelJacobiMapper{TA, TX, TB, TTemp, TF, TI}
    A::TA
    x::TX
    b::TB
    temp::TTemp
    ω::TF
    col::TI
end
function (pjacobmapper::ParallelJacobiMapper)(i)
    A = pjacobmapper.A
    x = pjacobmapper.x
    b = pjacobmapper.b
    temp = pjacobmapper.temp
    ω = pjacobmapper.ω
    col = pjacobmapper.col

    one = Base.one(eltype(A))
    z = zero(eltype(A))
    rsum = z
    diag = z

    for j in nzrange(A, i)
        row = A.rowval[j]
        val = A.nzval[j]

        diag = ifelse(row == i, val, diag)
        rsum += ifelse(row == i, z, val * temp[row, col])
    end
    xcand = (one - ω) * temp[i, col] + ω * ((b[i, col] - rsum) / diag)

    return ifelse(diag == 0, x[i, col], xcand)
end

function (jacobi::ParallelJacobi)(A, x, b)
    ω = jacobi.ω
    temp = jacobi.temp
    for col = 1:size(x, 2)
        @views tmap!(ParallelTemp(temp, col), x[1:size(A, 1), col], 1:size(A, 1))
        @views tmap!(ParallelJacobiMapper(A, x, b, temp, ω, col),
            x[1:size(A, 1), col], 1:size(A, 1))
    end
end
=#

struct JacobiProlongation{T}
    ω::T
end

struct DiagonalWeighting
end
struct LocalWeighting
end

function (j::JacobiProlongation)(A, T, S, B, degree = 1, weighting = LocalWeighting())
    D_inv_S = weight(weighting, A, j.ω)
    P = T
    for i = 1:degree
        P = P - (D_inv_S * P)
    end
    P
end

function weight(::DiagonalWeighting, S, ω)
    D_inv = 1 ./ diag(S)
    D_inv_S = scale_rows(S, D_inv)
    (eltype(S)(ω) / approximate_spectral_radius(D_inv_S)) * D_inv_S
    # (ω) * D_inv_S
end

function weight(::LocalWeighting, S, ω)
    #=D = abs.(S) * ones(eltype(S), size(S, 1))
    D_inv = 1 ./ D[find(D)]
    D_inv_S = scale_rows(S, D_inv)
    eltype(S)(ω) * D_inv_S=#
    D = zeros(eltype(S), size(S,1))
    for i = 1:size(S, 1)
        for j in nzrange(S, i)
            row = S.rowval[j]
            val = S.nzval[j]
            D[row] += abs(val)
        end
    end
    for i = 1:size(D, 1)
        if D[i] != 0
            D[i] = 1/D[i]
        end
    end
    D_inv_S = scale_rows(S, D)
    # eltype(S)(ω) * D_inv_S
    rmul!(D_inv_S, eltype(S)(ω))
end

#approximate_spectral_radius(A) =
#    eigs(A, maxiter = 15, tol = 0.01, ritzvec = false)[1][1] |> real

function scale_rows!(ret, S, v)
    n = size(S, 1)
    for i = 1:n
        for j in nzrange(S, i)
            row = S.rowval[j]
            ret.nzval[j] *= v[row]
        end
    end
    ret
end
scale_rows(S, v) = scale_rows!(deepcopy(S), S,  v)


text:
using Base.Test, DataStreams, DataFrames, NullableArrays

if !isdefined(Core, :String)
    typealias String UTF8String
end

ROWS = 2
COLS = 3

sch = Data.Schema([Int, Int, Int], ROWS)
data = [ 1 2 ; 10 20 ]

@test Data.header(sch) == String["Column1", "Column2", "Column3"]
@test Data.types(sch) == [Int, Int, Int]
@test size(sch) == (ROWS, COLS)

"""
A generic `Sink` type that fulfills the DataStreams interface
wraps any kind of Julia structure `T`; by default `T` = Vector{NullableVector}
"""
type TableSink{T} <: Data.Sink
    schema::Data.Schema
    data::T
    other::Any # for other metadata, references, etc.
end

function TableSink(sch::Data.Schema)
    rows, cols = size(sch)
    TableSink(sch, NullableVector[NullableArray(T, rows) for T in Data.types(sch)], 0)
end

TableSink(s::Data.Source) = TableSink(schema(s))

import DataStreams.Data.stream!

function Data.stream!(src::DataFrame, snk::TableSink)
    # TODO: this could be improved considering different source Schema.
    snk.data = src.columns
end

src_tb = DataFrame()
snk = TableSink(sch)

snk_tb = Data.stream!(src_tb, snk)

ext: jl


text:
import Base.show

export NodeData

mutable struct NodeData
    # Properties
    id::UInt32
    name::String
    visible::Bool
    dirty::Bool
    world::Ranger.World

    parent::Ranger.AbstractNode

    function NodeData()
        NodeData(UInt32(0), "_NoName_")
    end

    function NodeData(id::UInt32, name::String)
        obj = new()

        obj.id = id
        obj.name = name
        obj.visible = true
        obj.dirty = true

        obj.parent = NodeNil()    # default to self

        obj
    end

    function NodeData(name::String, parent::Ranger.AbstractNode, world::Ranger.World)
        obj = NodeData(Ranger.gen_id(world), name)

        obj.parent = parent
        obj.world = world
        obj
    end
end

function Base.show(io::IO, node::NodeData)
    print(io, "|'", node.name, "' (", node.id, ")|");
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module TimeProbeSeismic

using PyCall
using Distributed, Reexport

@reexport using LinearAlgebra, Images, HDF5, SegyIO, JOLI, JUDI, Serialization
@reexport using JLD2
@reexport using SlimOptim, PyPlot, JOLI, Printf, Random, Statistics

export h5read, forward, backprop, smooth, combine_probes, simil
export typedict, set_r!, set_ps!

# python imports
const dv = PyNULL()
const wu = PyNULL()
const ker = PyNULL()
const geom = PyNULL()
const si = PyNULL()
const ut = PyNULL()

function __init__()
    pushfirst!(PyVector(pyimport("sys")."path"), joinpath(JUDIPATH, "pysource"))
    copy!(dv, pyimport("devito"))
    copy!(wu, pyimport("wave_utils"))
    copy!(ker,  pyimport("kernels"))
    copy!(geom, pyimport("geom_utils"))
    copy!(si, pyimport("sensitivity"))
    copy!(ut, pyimport("utils"))
end

# Propagators
include("propagators.jl")
include("interface.jl")

# JUDI functions
include("judi.jl")

# Utility functions
include("utils.jl")

end

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
import HTTP
import Pkg

const GIT_LFS_REPOS_TOML_FILE_URL = "https://raw.githubusercontent.com/UnofficialJuliaMirror/MirrorUpdater.jl/master/config/repositories/git-lfs-repos-src-dest-pairs.toml"

const GIT_LFS_REPOS = Pkg.TOML.parse(
    String(HTTP.request("GET", GIT_LFS_REPOS_TOML_FILE_URL).body)
    )


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

module TaylorIntegration

using Reexport
@reexport using TaylorSeries
using LinearAlgebra
using Markdown
using Requires

export taylorinteg, lyap_taylorinteg, @taylorize

include("explicitode.jl")
include("lyapunovspectrum.jl")
include("rootfinding.jl")
include("parse_eqs.jl")

function __init__()
    @require DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" include("common.jl")
end

end #module


text:
function initialize!(integrator, cache::Union{ImplicitEulerConstantCache,
                                              ImplicitMidpointConstantCache,
                                              TrapezoidConstantCache,
                                              TRBDF2ConstantCache,
                                              SDIRK2ConstantCache,
                                              SDIRK22ConstantCache,
                                              SSPSDIRK2ConstantCache,
                                              Cash4ConstantCache,
                                              Hairer4ConstantCache,
                                              ESDIRK54I8L2SAConstantCache,
                                              SFSDIRK4ConstantCache,
                                              SFSDIRK5ConstantCache,
                                              SFSDIRK6ConstantCache,
                                              SFSDIRK7ConstantCache,
                                              SFSDIRK8ConstantCache})
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

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

function initialize!(integrator, cache::Union{ImplicitEulerCache,
                                              ImplicitMidpointCache,
                                              TrapezoidCache,
                                              TRBDF2Cache,
                                              SDIRK2Cache,
                                              SDIRK22Cache,
                                              SSPSDIRK2Cache,
                                              Cash4Cache,
                                              Hairer4Cache,
                                              ESDIRK54I8L2SACache,
                                              SFSDIRK4Cache,
                                              SFSDIRK5Cache,
                                              SFSDIRK6Cache,
                                              SFSDIRK7Cache,
                                              SFSDIRK8Cache})
  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) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ImplicitEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  # initial guess
  if alg.extrapolant == :linear
    nlsolver.z = dt*integrator.fsalfirst
  else # :constant
    nlsolver.z = zero(u)
  end

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

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    tmp = r*integrator.opts.internalnorm.((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  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::ImplicitEulerCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack atmp,nlsolver = cache
  @unpack z,tmp = nlsolver
  mass_matrix = integrator.f.mass_matrix
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  # initial guess
  if alg.extrapolant == :linear
    @.. z = dt*integrator.fsalfirst
  else # :constant
    z .= zero(eltype(u))
  end

  nlsolver.tmp .= uprev
  nlsolver.γ = 1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = uprev + z

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    @.. tmp = r*integrator.opts.internalnorm((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end
  integrator.destats.nf += 1
  f(integrator.fsallast,u,p,t+dt)
end

@muladd function perform_step!(integrator, cache::ImplicitMidpointConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  γ = 1//2
  markfirststage!(nlsolver)

  # initial guess
  if alg.extrapolant == :linear
    nlsolver.z = dt*integrator.fsalfirst
  else # :constant
    nlsolver.z = zero(u)
  end

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

  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::ImplicitMidpointCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack nlsolver = cache
  @unpack z,tmp = nlsolver
  mass_matrix = integrator.f.mass_matrix
  alg = unwrap_alg(integrator, true)
  γ = 1//2
  markfirststage!(nlsolver)

  # initial guess
  if alg.extrapolant == :linear
    @.. z = dt*integrator.fsalfirst
  else # :constant
    z .= zero(eltype(u))
  end

  nlsolver.tmp = uprev
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = nlsolver.tmp + z

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

@muladd function perform_step!(integrator, cache::TrapezoidConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  # precalculations
  γ = 1//2
  γdt = γ*dt
  markfirststage!(nlsolver)

  # initial guess: constant extrapolation
  nlsolver.z = uprev

  if f.mass_matrix === I
    nlsolver.tmp = @.. uprev * inv(γdt) + integrator.fsalfirst
  else
    nlsolver.tmp = (f.mass_matrix * uprev) .* inv(γdt) .+ integrator.fsalfirst
  end
  nlsolver.α = 1
  nlsolver.γ = γ
  nlsolver.method = COEFFICIENT_MULTISTEP
  u = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      # tmp = r*abs(((u - uprev)/dt1 - (uprev - uprev2)/dt2) - ((uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4)/dt5)
      DD31 = (u - uprev)/dt1 - (uprev - uprev2)/dt2
      DD30 = (uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4
      tmp = r*integrator.opts.internalnorm((DD31 - DD30)/dt5,t)
      atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        cache.uprev3 = uprev2
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      cache.uprev3 = integrator.uprev2
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    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

@muladd function perform_step!(integrator, cache::TrapezoidCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack atmp,nlsolver = cache
  @unpack z,tmp = nlsolver
  alg = unwrap_alg(integrator, true)
  mass_matrix = integrator.f.mass_matrix

  # precalculations
  γ = 1//2
  γdt = γ*dt
  markfirststage!(nlsolver)

  # initial guess: constant extrapolation
  @.. z = uprev
  invγdt = inv(γdt)
  if mass_matrix === I
    @.. tmp = uprev * invγdt + integrator.fsalfirst
  else
    mul!(u, mass_matrix, uprev)
    @.. tmp = u * invγdt + integrator.fsalfirst
  end
  nlsolver.α = 1
  nlsolver.γ = γ
  nlsolver.method = COEFFICIENT_MULTISTEP
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = z

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      # @.. tmp = r*abs(((u - uprev)/dt1 - (uprev - uprev2)/dt2) - ((uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4)/dt5)
      @.. tmp = r*integrator.opts.internalnorm((
          ((u - uprev)/dt1 - (uprev- uprev2)/dt2) #DD31
          - ((uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4) #DD30
      )/dt5,t)
      calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        copyto!(cache.uprev3,uprev2)
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      copyto!(cache.uprev3,integrator.uprev2)
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    end
  end

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

@muladd function perform_step!(integrator, cache::TRBDF2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,d,ω,btilde1,btilde2,btilde3,α1,α2 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  # FSAL
  zprev = dt*integrator.fsalfirst

  ##### Solve Trapezoid Step

  # TODO: Add extrapolation
  zᵧ = zprev
  nlsolver.z = zᵧ
  nlsolver.c = γ

  nlsolver.tmp = uprev + d*zprev
  zᵧ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve BDF2 Step

  ### Initial Guess From Shampine
  z = α1*zprev + α2*zᵧ
  nlsolver.z = z
  nlsolver.c = 1

  nlsolver.tmp = uprev + ω*zprev + ω*zᵧ
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + d*z

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

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

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

@muladd function perform_step!(integrator, cache::TRBDF2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack zprev,zᵧ,atmp,nlsolver = cache
  @unpack z,tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  b = nlsolver.ztmp
  @unpack γ,d,ω,btilde1,btilde2,btilde3,α1,α2 = cache.tab
  alg = unwrap_alg(integrator, true)

  # FSAL
  @.. zprev = dt*integrator.fsalfirst
  markfirststage!(nlsolver)

  ##### Solve Trapezoid Step

  # TODO: Add extrapolation
  @.. zᵧ = zprev
  z .= zᵧ
  @.. tmp = uprev + d*zprev
  nlsolver.c = γ
  zᵧ .= nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve BDF2 Step

  ### Initial Guess From Shampine
  @.. z = α1*zprev + α2*zᵧ
  @.. tmp = uprev + ω*zprev + ω*zᵧ
  nlsolver.c = 1
  isnewton(nlsolver) && set_new_W!(nlsolver, false)
  nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + d*z

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

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

  @.. integrator.fsallast = z/dt
end

@muladd function perform_step!(integrator, cache::TRBDF2Cache{<:Array}, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack zprev,zᵧ,atmp,nlsolver = cache
  @unpack z,tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  b = nlsolver.ztmp
  @unpack γ,d,ω,btilde1,btilde2,btilde3,α1,α2 = cache.tab
  alg = unwrap_alg(integrator, true)

  # FSAL
  @inbounds @simd ivdep for i in eachindex(u)
    zprev[i] = dt*integrator.fsalfirst[i]
  end
  markfirststage!(nlsolver)

  ##### Solve Trapezoid Step

  # TODO: Add extrapolation
  copyto!(zᵧ,zprev)
  copyto!(z,zᵧ)
  @inbounds @simd ivdep for i in eachindex(u)
    tmp[i] = uprev[i] + d*zprev[i]
  end
  nlsolver.c = γ
  zᵧ .= nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve BDF2 Step

  ### Initial Guess From Shampine
  @inbounds @simd ivdep for i in eachindex(u)
    z[i] = α1*zprev[i] + α2*zᵧ[i]
  end
  @inbounds @simd ivdep for i in eachindex(u)
    tmp[i] = uprev[i] + ω*zprev[i] + ω*zᵧ[i]
  end
  nlsolver.c = 1
  isnewton(nlsolver) && set_new_W!(nlsolver, false)
  nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @inbounds @simd ivdep for i in eachindex(u)
    u[i] = tmp[i] + d*z[i]
  end

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

  if integrator.opts.adaptive

    @inbounds @simd ivdep for i in eachindex(u)
      tmp[i] = btilde1*zprev[i] + btilde2*zᵧ[i] + btilde3*z[i]
    end
    if alg.smooth_est && isnewton(nlsolver) # From Shampine
      est = nlsolver.cache.dz
      nlsolver.cache.linsolve(vec(est),W,vec(tmp),false)
      integrator.destats.nsolve += 1
    else
      est = tmp
    end
    calculate_residuals!(atmp, est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @inbounds @simd ivdep for i in eachindex(u)
    integrator.fsallast[i] = z[i]/dt
  end
end

@muladd function perform_step!(integrator, cache::SDIRK2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  # initial guess
  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
    z₁ = u - uprev
  elseif alg.extrapolant == :linear
    z₁ = dt*integrator.fsalfirst
  else
    z₁ = zero(u)
  end

  nlsolver.tmp = uprev
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ### Initial Guess Is α₁ = c₂/γ, c₂ = 0 => z₂ = α₁z₁ = 0
  z₂ = zero(u)
  nlsolver.z = z₂
  nlsolver.tmp = uprev - z₁
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  u = uprev + z₁/2 + z₂/2

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

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

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

@muladd function perform_step!(integrator, cache::SDIRK2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,atmp,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  # initial guess
  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
    @.. z₁ = u - uprev
  elseif alg.extrapolant == :linear
    @.. z₁ = dt*integrator.fsalfirst
  else
    z₁ .= zero(eltype(u))
  end
  nlsolver.z = z₁

  ##### Step 1
  nlsolver.tmp = uprev
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 2

  ### Initial Guess Is α₁ = c₂/γ, c₂ = 0 => z₂ = α₁z₁ = 0
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂
  isnewton(nlsolver) && set_new_W!(nlsolver, false)
  @.. tmp = uprev - z₁
  nlsolver.tmp = tmp
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = uprev + z₁/2 + z₂/2

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

  if integrator.opts.adaptive
    @.. tmp = z₁/2 - z₂/2
    if alg.smooth_est && isnewton(nlsolver) # From Shampine
      est = nlsolver.cache.dz
      nlsolver.cache.linsolve(vec(est),W,vec(tmp),false)
      integrator.destats.nsolve += 1
    else
      est = tmp
    end
    calculate_residuals!(atmp, est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

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

@muladd function perform_step!(integrator, cache::SDIRK22ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack a,α,β = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)

  # precalculations
  γ = a*dt
  γdt = γ*dt
  markfirststage!(nlsolver)

  # initial guess
  zprev = dt*integrator.fsalfirst
  nlsolver.z = zprev

  # first stage
  nlsolver.tmp = uprev + γdt*integrator.fsalfirst
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  uprev = α*nlsolver.tmp + β*z

  # final stage
  γ = dt
  γdt = γ*dt
  markfirststage!(nlsolver)
  nlsolver.tmp = uprev + γdt*integrator.fsalfirst
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      DD31 = (u - uprev)/dt1 - (uprev - uprev2)/dt2
      DD30 = (uprev - uprev2)/dt3 - (uprev2 - uprev3)/dt4
      tmp = r*abs((DD31 - DD30)/dt5)
      atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        cache.uprev3 = uprev2
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      cache.uprev3 = integrator.uprev2
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    end
  end

  integrator.destats.nf += 2
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

@muladd function perform_step!(integrator, cache::SDIRK22Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack atmp,nlsolver = cache
  @unpack z,tmp = nlsolver
  @unpack a,α,β = cache.tab
  alg = unwrap_alg(integrator, true)
  mass_matrix = integrator.f.mass_matrix

  # precalculations
  γ = a*dt
  γdt = γ*dt
  markfirststage!(nlsolver)

  # first stage
  @.. z = dt*integrator.fsalfirst
  @.. tmp = uprev + γdt*integrator.fsalfirst
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = α*tmp + β*z

  # final stage
  γ = dt
  γdt = γ*dt
  markfirststage!(nlsolver)
  @.. tmp = uprev + γdt*integrator.fsalfirst
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = nlsolver.tmp

  if integrator.opts.adaptive
    if integrator.iter > 2
      # local truncation error (LTE) bound by dt^3/12*max|y'''(t)|
      # use 3rd divided differences (DD) a la SPICE and Shampine

      # TODO: check numerical stability
      uprev2 = integrator.uprev2
      tprev = integrator.tprev
      uprev3 = cache.uprev3
      tprev2 = cache.tprev2

      dt1 = dt*(t+dt-tprev)
      dt2 = (t-tprev)*(t+dt-tprev)
      dt3 = (t-tprev)*(t-tprev2)
      dt4 = (tprev-tprev2)*(t-tprev2)
      dt5 = t+dt-tprev2
      c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
      r = c*dt^3/2 # by mean value theorem 3rd DD equals y'''(s)/6 for some s

      @inbounds for i in eachindex(u)
        DD31 = (u[i] - uprev[i])/dt1 - (uprev[i] - uprev2[i])/dt2
        DD30 = (uprev[i] - uprev2[i])/dt3 - (uprev2[i] - uprev3[i])/dt4
        tmp[i] = r*abs((DD31 - DD30)/dt5)
      end
      calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
      if integrator.EEst <= 1
        copyto!(cache.uprev3,uprev2)
        cache.tprev2 = tprev
      end
    elseif integrator.success_iter > 0
      integrator.EEst = 1
      copyto!(cache.uprev3,integrator.uprev2)
      cache.tprev2 = integrator.tprev
    else
      integrator.EEst = 1
    end
  end

  integrator.destats.nf += 2
  f(integrator.fsallast,u,p,t+dt)
end

@muladd function perform_step!(integrator, cache::SSPSDIRK2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)

  γ = eltype(u)(1//4)
  c2 = typeof(t)(3//4)

  markfirststage!(nlsolver)

  # initial guess
  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
    z₁ = u - uprev
  elseif alg.extrapolant == :linear
    z₁ = dt*integrator.fsalfirst
  else
    z₁ = zero(u)
  end
  nlsolver.z = z₁

  ##### Step 1

  tstep = t + dt
  u = uprev + γ*z₁

  nlsolver.c = 1
  nlsolver.tmp = uprev
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 2

  ### Initial Guess Is α₁ = c₂/γ
  z₂ = c2/γ
  nlsolver.z = z₂

  nlsolver.tmp = uprev + z₁/2
  nlsolver.c = 1
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + z₂/2

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

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

@muladd function perform_step!(integrator, cache::SSPSDIRK2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,nlsolver = cache
  @unpack tmp = nlsolver
  alg = unwrap_alg(integrator, true)

  γ = eltype(u)(1//4)
  c2 = typeof(t)(3//4)
  markfirststage!(nlsolver)

  # initial guess
  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
    @.. z₁ = u - uprev
  elseif alg.extrapolant == :linear
    @.. z₁ = dt*integrator.fsalfirst
  else
    z₁ .= zero(eltype(u))
  end
  nlsolver.z = z₁
  nlsolver.tmp = uprev

  ##### Step 1
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Solve Step 2

  ### Initial Guess Is α₁ = c₂/γ
  @.. z₂ = c2/γ
  nlsolver.z = z₂

  @.. tmp = uprev + z₁/2
  nlsolver.tmp = tmp
  isnewton(nlsolver) && set_new_W!(nlsolver, false)
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + z₂/2

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

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

@muladd function perform_step!(integrator, cache::Cash4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,c2,c3,c4 = cache.tab
  @unpack b1hat1,b2hat1,b3hat1,b4hat1,b1hat2,b2hat2,b3hat2,b4hat2 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use yhat2 for prediction
  z₅ = b1hat2*z₁ + b2hat2*z₂ + b3hat2*z₃ + b4hat2*z₄
  nlsolver.z = z₅

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

  u = nlsolver.tmp + γ*z₅

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

  if integrator.opts.adaptive
    if alg.embedding == 3
      btilde1 = b1hat2-a51; btilde2 = b2hat2-a52;
      btilde3 = b3hat2-a53; btilde4 = b4hat2-a54; btilde5 = -γ
    else
      btilde1 = b1hat1-a51; btilde2 = b2hat1-a52;
      btilde3 = b3hat1-a53; btilde4 = b4hat1-a54; btilde5 = -γ
    end

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

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

@muladd function perform_step!(integrator, cache::Cash4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,atmp,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,c2,c3,c4 = cache.tab
  @unpack b1hat1,b2hat1,b3hat1,b4hat1,b1hat2,b2hat2,b3hat2,b4hat2 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  @.. z₅ = b1hat2*z₁ + b2hat2*z₂ + b3hat2*z₃ + b4hat2*z₄
  nlsolver.z = z₅
  @.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = 1
  z₅ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₅

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

  if integrator.opts.adaptive
    if alg.embedding == 3
      btilde1 = b1hat2-a51; btilde2 = b2hat2-a52;
      btilde3 = b3hat2-a53; btilde4 = b4hat2-a54; btilde5 = -γ
    else
      btilde1 = b1hat1-a51; btilde2 = b2hat1-a52;
      btilde3 = b3hat1-a53; btilde4 = b4hat1-a54; btilde5 = -γ
    end

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

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

@muladd function perform_step!(integrator, cache::SFSDIRK4ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,c2,c3,c4 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

  ################################## Final Step

  u = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄

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

@muladd function perform_step!(integrator, cache::SFSDIRK4Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,c2,c3,c4 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)
  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

  @.. u = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄

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

@muladd function perform_step!(integrator, cache::SFSDIRK5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  z₅ = z₄
  nlsolver.z = z₅

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

  ################################## Final Step

  u = uprev + a61*z₁ + a62*z₂ + a63*z₃ + a64*z₄+ a65*z₅

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

@muladd function perform_step!(integrator, cache::SFSDIRK5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)
  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  @.. z₅ = z₄
  nlsolver.z = z₅

  @.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = c5
  z₅ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
################################### Finalize
  @.. u = uprev + a61*z₁ + a62*z₂ + a63*z₃ + a64*z₄ + a65*z₅

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

@muladd function perform_step!(integrator, cache::SFSDIRK6ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,c2,c3,c4,c5,c6 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  z₆ = z₅
  nlsolver.z = z₆

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

  ################################## Final Step

  u = uprev + a71*z₁ + a72*z₂ + a73*z₃ + a74*z₄+ a75*z₅ + a76*z₆

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

@muladd function perform_step!(integrator, cache::SFSDIRK6Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,z₆,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,c2,c3,c4,c5,c6 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)
  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  @.. z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  @.. z₆ = z₅
  nlsolver.z = z₆

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

  ################################### Finalize
  @.. u = uprev + a71*z₁ + a72*z₂ + a73*z₃ + a74*z₄ + a75*z₅ + a76*z₆

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

@muladd function perform_step!(integrator, cache::SFSDIRK7ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a82,a83,a84,a85,a86,a87,c2,c3,c4,c5,c6,c7 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  z₆ = z₅
  nlsolver.z = z₆

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

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

  # Use constant z prediction
  z₇ = z₆
  nlsolver.z = z₇

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

  ################################## Final Step

  u = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄+ a85*z₅ + a86*z₆ + a87*z₇

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

@muladd function perform_step!(integrator, cache::SFSDIRK7Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,z₆,z₇,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a82,a83,a84,a85,a86,a87,c2,c3,c4,c5,c6,c7 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)
  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  @.. z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  @.. z₆ = z₅
  nlsolver.z = z₆

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

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

  # Use constant z prediction
  @.. z₇ = z₆
  nlsolver.z = z₇

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

  ################################### Finalize
  @.. u = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇

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

@muladd function perform_step!(integrator, cache::SFSDIRK8ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a82,a83,a84,a85,a86,a87,a91,a92,a93,a94,a95,a96,a97,a98,c2,c3,c4,c5,c6,c7,c8 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)

  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z = z₁

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

  ##### Step 2

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

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

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

  # Guess starts from z₁
  z₃ = z₁
  nlsolver.z = z₃

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

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

  # Use constant z prediction
  z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  z₆ = z₅
  nlsolver.z = z₆

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

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

  # Use constant z prediction
  z₇ = z₆
  nlsolver.z = z₇

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

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

  # Use constant z prediction
  z₈ = z₇
  nlsolver.z = z₈

  nlsolver.tmp = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
  nlsolver.c = c8
  z₈ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################## Final Step

  u = uprev + a91*z₁ + a92*z₂ + a93*z₃ + a94*z₄+ a95*z₅ + a96*z₆ + a97*z₇ + a98*z₈

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

@muladd function perform_step!(integrator, cache::SFSDIRK8Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,nlsolver = cache
  @unpack tmp = nlsolver
  W = isnewton(nlsolver) ? get_W(nlsolver) : nothing
  @unpack γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a82,a83,a84,a85,a86,a87,a91,a92,a93,a94,a95,a96,a97,a98,c2,c3,c4,c5,c6,c7,c8 = cache.tab
  alg = unwrap_alg(integrator, true)
  markfirststage!(nlsolver)
  ##### Step 1

  # TODO: Add extrapolation for guess
  z₁ .= zero(eltype(z₁))
  nlsolver.z = z₁
  nlsolver.c = γ
  nlsolver.tmp = uprev

  # initial step of NLNewton iteration
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(z₂))
  nlsolver.z = z₂

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

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

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

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

  # Use constant z prediction
  @.. z₄ = z₃
  nlsolver.z = z₄

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

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

  # Use constant z prediction
  @.. z₅ = z₄
  nlsolver.z = z₅

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

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

  # Use constant z prediction
  @.. z₆ = z₅
  nlsolver.z = z₆

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

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

  # Use constant z prediction
  @.. z₇ = z₆
  nlsolver.z = z₇

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

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

  # Use constant z prediction
  @.. z₈ = z₇
  nlsolver.z = z₈

  @.. tmp = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
  nlsolver.c = c8
  z₈ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ################################### Finalize
  @.. u = uprev + a91*z₁ + a92*z₂ + a93*z₃ + a94*z₄ + a95*z₅ + a96*z₆ + a97*z₇ + a98*z₈

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

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

  # precalculations
  γdt = γ*dt
  markfirststage!(nlsolver)

  # TODO: Add extrapolation for guess
  z₁ = zero(u)
  nlsolver.z, nlsolver.tmp = z₁, uprev
  nlsolver.c = γ
  z₁ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  ##### Step 2

  z₂ = α21*z₁
  nlsolver.z = z₂
  nlsolver.tmp = uprev + a21*z₁
  nlsolver.c = c2
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

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

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

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

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

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

  # Use yhat2 for prediction
  z₅ = bhat1*z₁ + bhat2*z₂ + bhat3*z₃ + bhat4*z₄
  nlsolver.z = z₅
  nlsolver.tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = 1
  z₅ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₅

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

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

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

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

  # initial guess
  if integrator.success_iter > 0 && !integrator.reeval_fsal && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
    @.. z₁ = u - uprev
  elseif alg.extrapolant == :linear
    @.. z₁ = dt*integrator.fsalfirst
  else
    z₁ .= zero(eltype(z₁))
  end
  nlsolver.z = z₁
  nlsolver.tmp = uprev

  ##### Step 1

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

  ##### Step 2

  @.. z₂ = α21*z₁
  nlsolver.z = z₂
  @.. tmp = uprev + a21*z₁
  nlsolver.tmp = tmp
  nlsolver.c = c2
  isnewton(nlsolver) && set_new_W!(nlsolver, false)
  z₂ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

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

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

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

  # Use constant z prediction
  @.. z₄ = α41*z₁ + α43*z₃
  nlsolver.z = z₄
  @.. tmp = uprev + a41*z₁ + a42*z₂ + a43*z₃
  nlsolver.c = c4
  z₄ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

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

  # Use yhat prediction
  @.. z₅ = bhat1*z₁ + bhat2*z₂ + bhat3*z₃ + bhat4*z₄
  nlsolver.z = z₅
  @.. tmp = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄
  nlsolver.c = 1
  z₅ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₅

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

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

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

@muladd function perform_step!(integrator, cache::ESDIRK54I8L2SAConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack γ,
          a31, a32,
          a41, a42, a43,
          a51, a52, a53, a54,
          a61, a62, a63, a64, a65,
          a71, a72, a73, a74, a75, a76,
          a81, a82, a83, a84, a85, a86, a87,
                    c3,  c4,  c5,  c6,  c7,
          btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7, btilde8 = cache.tab
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt
  markfirststage!(nlsolver)

  # TODO: Add extrapolation for guess

  ##### Step 1

  z₁ = dt*integrator.fsalfirst

  ##### Step 2

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

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

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

  nlsolver.z = z₃ = zero(z₂)

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

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

  nlsolver.z = z₄ = zero(z₃)

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

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

  nlsolver.z = z₅ = zero(z₄)

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

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

  nlsolver.z = z₆ = zero(z₅)

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

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

  nlsolver.z = z₇ = zero(z₆)

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

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

  nlsolver.z = z₈ = zero(z₇)

  nlsolver.tmp = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
  nlsolver.c = 1
  z₈ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  u = nlsolver.tmp + γ*z₈

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

  if integrator.opts.adaptive
    est = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈
    atmp = calculate_residuals(est, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

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

@muladd function perform_step!(integrator, cache::ESDIRK54I8L2SACache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,atmp,nlsolver = cache
  @unpack tmp = nlsolver
  @unpack γ,
          a31, a32,
          a41, a42, a43,
          a51, a52, a53, a54,
          a61, a62, a63, a64, a65,
          a71, a72, a73, a74, a75, a76,
          a81, a82, a83, a84, a85, a86, a87,
                    c3,  c4,  c5,  c6,  c7,
          btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7, btilde8 = cache.tab
  alg = unwrap_alg(integrator, true)

  # precalculations
  γdt = γ*dt
  markfirststage!(nlsolver)

  ##### Step 1

  @.. z₁ = dt*integrator.fsalfirst

  ##### Step 2

  # TODO: Add extrapolation for guess
  z₂ .= zero(eltype(u))
  nlsolver.z = z₂

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

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

  nlsolver.z = fill!(z₃, zero(eltype(u)))

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

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

  # Use constant z prediction
  nlsolver.z = fill!(z₄, zero(eltype(u)))

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

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

  nlsolver.z = fill!(z₅, zero(eltype(u)))

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

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

  nlsolver.z = fill!(z₆, zero(eltype(u)))

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

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

  nlsolver.z = fill!(z₇, zero(eltype(u)))

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

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

  nlsolver.z = fill!(z₈, zero(eltype(u)))

  @.. nlsolver.tmp = uprev + a81*z₁ + a82*z₂ + a83*z₃ + a84*z₄ + a85*z₅ + a86*z₆ + a87*z₇
  nlsolver.c = oneunit(nlsolver.c)
  z₈ = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. u = tmp + γ*z₈

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

  if integrator.opts.adaptive
    @.. tmp = btilde1*z₁ + btilde2*z₂ + btilde3*z₃ + btilde4*z₄ + btilde5*z₅ + btilde6*z₆ + btilde7*z₇ + btilde8*z₈
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @.. integrator.fsallast = z₈/dt
  return
end


text:
"""
`UnivariateFeatureSelector` has the following fields:

* `method::Function` (required)- Method to calculate feature importance. The
  method chosen will determine the scoring. Below is the scring with available
  statistical method to obtain them.

    * Correlation - higher score means more important

        * `pearson_correlation`

    * P-value - lower score means more important

        * `f_test`
        * `chisq_test`

* `k::Union{Int64,Nothing}` - Select top `k` features with the highest correlation to target
  variable. You could ignore this by specifying k == nothing. This defaults to nothing.
* `threshold::Union{Float64,Nothing}` - Select features with correlation more than or equal
  to threshold. To ignore, simply set threshold to nothing (default behavior).
"""
mutable struct UnivariateFeatureSelector
    method::Function
    k::Union{Int64,Nothing}
    threshold::Union{Float64,Nothing}
end

function UnivariateFeatureSelector(;
    method::Function,
    k::Union{Int64,Nothing} = nothing,
    threshold::Union{Float64,Nothing} = nothing,
)
    UnivariateFeatureSelector(method, k, threshold)
end

"""
    select_features(selector,
                    X::DataFrame,
                    y::Vector;
                    verbose::Bool=false)

Select features based on the importance, which is defined by `selector.method`
to target `y`. if `verbose` is true, logs will be printed - this defaults to
false. This function will return only the feature names of selected features.

If you have feature `X_data` as matrix and feature names `X_features` as a
Vector, you can replace `X` with `X_data` and `X_features` (in this order).

# Example

```jldoctest
julia> using RDatasets, FeatureSelectors, DataFrames

julia> boston = dataset("MASS", "Boston");

julia> selector = UnivariateFeatureSelector(method=pearson_correlation, k=5)
UnivariateFeatureSelector(FeatureSelectors.pearson_correlation, 5, nothing)

julia> select_features(
           selector,
           boston[:, Not(:MedV)],
           boston.MedV
       )
5-element Vector{String}:
 "LStat"
 "Rm"
 "PTRatio"
 "Indus"
 "Tax"

```
"""
function select_features(
    selector::UnivariateFeatureSelector,
    X_data::Matrix,
    X_features::Vector,
    y::Vector;
    verbose::Bool = false,
)
    score_arr = selector.method(X_data, y)
    # Correlation should be sorted descending
    if selector.method == pearson_correlation
        rev = true
        filter_operator = >=
        if verbose
            @info "method is correlation score, result will be " * "descendingly sorted" rev
        end
    else
        rev = false
        filter_operator = <=
        if verbose
            @info "method is p-value, result will be " * "ascendingly sorted" rev
        end
    end
    sorted_tup =
        sort([x for x in zip(X_features, score_arr)], by = v -> abs(v[2]), rev = rev)
    if !isnothing(selector.k)
        if selector.k < 1
            @warn "k cannot be less than 1. Resetting value to 1"
            selector.k = 1
        end
        verbose && @info "Filtering top k features" selector.k
        sorted_tup = sorted_tup[1:selector.k]
    end
    if !isnothing(selector.threshold)
        verbose && @info "Filtering by threshold" selector.threshold
        sorted_tup = filter(v -> filter_operator(abs(v[2]), selector.threshold), sorted_tup)
    end
    first.(sorted_tup)
end

select_features(
    selector::UnivariateFeatureSelector,
    X::DataFrame,
    y::Vector;
    verbose::Bool = false,
) = select_features(
    selector,
    Matrix(X),
    names(X),
    y;
    verbose = verbose,
)

"""
    pearson_correlation(X_data::Matrix, y::Vector)

Calculate pearson's correlation on `X_data` to `y`.
"""
pearson_correlation(X_data::Matrix, y::Vector) = cor(X_data, y)

"""
    f_test(X_data::Matrix, y::Vector)

Calculate p-value using f-test method.
"""
function f_test(X_data::Matrix, y::Vector)
    pvals = [pvalue(VarianceFTest(X_col, y)) for X_col in eachcol(X_data)]
end

"""
    chisq_test(X_data::Matrix, y::Vector)

Calculate p-value using chi-square test.
"""
function chisq_test(X_data::Matrix, y::Vector)
    # DataFrame for easy aggregation
    data_df = DataFrame(hcat(X_data, y), :auto)
    # Rename cols to prevent name clash with count aggregation below
    rename!(data_df, [Symbol("input$i") for i = 1:size(data_df)[2]])
    y_name = names(data_df[!, [end]])[1]
    pvals = Vector{Float64}()
    for col_name in names(data_df[:, 1:end-1])
        # Hacky method to generate frequency count and put 0 for missing
        # combination. `by` function is not able to achieve this. See:
        # https://github.com/JuliaData/DataFrames.jl/issues/2136

        # Feature + target to process and put count = 1 for each row
        Xy_df = hcat(data_df[:, [col_name, y_name]], Int.(ones(size(data_df)[1])))
        rename!(Xy_df, :x1 => :Count)
        # This nested join will append missing in Count to missing combination
        # of value in col_name and y_name
        Xy_df = outerjoin(
            crossjoin(
                DataFrame(; Symbol(col_name) => unique(Xy_df[:, col_name])),
                DataFrame(; Symbol(y_name) => unique(Xy_df[:, y_name])),
            ),
            Xy_df,
            on = [col_name, y_name],
        )
        # sum count and recode missing to 0
        _y = replace(
            sort(combine(r -> sum(r.Count), groupby(Xy_df, [col_name, y_name]))).x1,
            missing => 0,
        )
        _y = reshape(
            _y,
            (length(unique(Xy_df[:, col_name])), length(unique(Xy_df[:, y_name]))),
        )
        push!(pvals, pvalue(ChisqTest(_y)))
    end
    pvals
end


text:
type CRT
    data::Ptr{Void}

    function CRT()
        c = new()
        c.data = Libc.malloc(1000)  # 552
        ccall((:mbedtls_x509_crt_init, MBED_X509), Void, (Ptr{Void},), c.data)
        finalizer(c, c->begin
            ccall((:mbedtls_x509_crt_free, MBED_X509), Void, (Ptr{Void},), c.data)
            Libc.free(c.data)
        end)
        c
    end
end

show(io::IO, crt::CRT) = print(io, crt_info(crt))

function crt_info(crt::CRT)
    buf = zeros(UInt8, 1000)
    ccall((:mbedtls_x509_crt_info, MBED_X509), Cint,
        (Ptr{Void}, Csize_t, Cstring, Ptr{Void}),
        buf, 1000, "", crt.data)
    String(pointer(buf))
end

function crt_parse!(chain, buf::String)
    ret = ccall((:mbedtls_x509_crt_parse, MBED_X509), Cint,
        (Ptr{Void}, Ptr{UInt8}, Csize_t),
        chain.data, buf, sizeof(buf)+1)
    ret == 0 || mbed_err(ret)
    chain
end

crt_parse!(chain, buf::IOStream) = crt_parse!(chain, readstring(buf))

crt_parse_file(path) = crt_parse(readstring(path))

function crt_parse(buf)
    crt = CRT()
    crt_parse!(crt, String(buf))
    crt
end


text:
function commandSet(ed::Editor, args::String)
    arg_arr = strip.(split(args, ' '))

    # Initial checks
    if length(arg_arr) != 2
        setStatusMessage(ed, "set: command requires two arguments")
        return
    elseif !Base.isidentifier(arg_arr[1])
        setStatusMessage(ed, "set: $(arg_arr[1]) is not a valid command name")
        return
    end

    sym = Symbol(arg_arr[1])

    # Check if it is a valid parameter
    if !configIsParam(sym)
        setStatusMessage(ed, "set: '$sym' is not a valid parameter name")
        return
    end

    # Attempt to assign the parameter
    try
        val = parse(arg_arr[2])
        configSet(sym, val)
    catch Exception
        setStatusMessage(ed, "set: invalid argument for $sym '$(arg_arr[2])'")
    end

    for row in ed.rows
        update!(row)
    end
    refreshScreen(ed)
end

addCommand(:set, commandSet,
           help="set <param> <value>: set the given parameter to a value")


text:
module FastIntegers

typealias IntTypes Union{Int8, Int16, Int32, Int64, Int128}
typealias UIntTypes Union{UInt8, UInt16, UInt32, UInt64, UInt128}
typealias IntegerTypes Union{IntTypes, UIntTypes}

# LLVM seems to ignore this
export fast_assert
function fast_assert(cond::Bool)
    Base.llvmcall("""
        br i1 %0, label %good, label %bad
    good:
        ret void
    bad:
        unreachable
    """, Void, Tuple{Bool}, cond)
end
@generated function fast_assert{T<:IntegerTypes}(x::T, cond::Bool)
    tp = "i$(8*sizeof(T))"
    quote
        Base.llvmcall("""
            br i1 %1, label %good, label %bad
        good:
            ret $($tp) %0
        bad:
            unreachable
        """, T, Tuple{T, Bool}, x, cond)
    end
end

# This leads to segfaults in LLVM
#=
export fast_assume
@generated function fast_assume{T<:IntegerTypes}(x::T, cond::Bool)
    tp = "i$(8*sizeof(T))"
    quote
        Base.llvmcall("""
            br i1 %1, label %good, label %bad
        good:
            ret $($tp) %0
        bad:
            ret $($tp) undef
        """, T, Tuple{T, Bool}, x, cond)
    end
end
=#

export fast_pos
fast_pos{T<:IntegerTypes}(x::T) = x

export fast_neg
fast_neg{T<:IntegerTypes}(x::T) = T(0) - x

export fast_abs
function fast_abs{T<:IntTypes}(x::T)
    # Note: -x = ~(x-1)
    m = -signbit(x)
    r = (x + m) $ m
    # fast_assert(r >= T(0))
    # r = fast_assert(r, r >= T(0))
    # r = fast_assume(r, r >= T(0))
    r
end
fast_abs{T<:UIntTypes}(x::T) = x

export fast_add
@generated function fast_add{T<:IntegerTypes}(x::T, y::T)
    tp = "i$(8*sizeof(T))"
    nw = T <: IntTypes ? "nsw" : "nuw"
    quote
        Base.llvmcall("""
            %res = add $($nw) $($tp) %0, %1
            ret $($tp) %res
        """, T, Tuple{T, T}, x, y)
    end
end

export fast_sub
@generated function fast_sub{T<:IntegerTypes}(x::T, y::T)
    tp = "i$(8*sizeof(T))"
    nw = T <: IntTypes ? "nsw" : "nuw"
    quote
        Base.llvmcall("""
            %res = sub $($nw) $($tp) %0, %1
            ret $($tp) %res
        """, T, Tuple{T, T}, x, y)
    end
end

export fast_mul
@generated function fast_mul{T<:IntegerTypes}(x::T, y::T)
    tp = "i$(8*sizeof(T))"
    nw = T <: IntTypes ? "nsw" : "nuw"
    quote
        Base.llvmcall("""
            %res = mul $($nw) $($tp) %0, %1
            ret $($tp) %res
        """, T, Tuple{T, T}, x, y)
    end
end

export fast_div
@generated function fast_div{T<:IntegerTypes}(x::T, y::T)
    tp = "i$(8*sizeof(T))"
    div = T <: IntTypes ? "sdiv" : "udiv"
    quote
        Base.llvmcall("""
            %res = $($div) $($tp) %0, %1
            ret $($tp) %res
        """, T, Tuple{T, T}, x, y)
    end
end

export fast_rem
@generated function fast_rem{T<:IntegerTypes}(x::T, y::T)
    tp = "i$(8*sizeof(T))"
    rem = T <: IntTypes ? "srem" : "urem"
    quote
        Base.llvmcall("""
            %res = $($rem) $($tp) %0, %1
            ret $($tp) %res
        """, T, Tuple{T, T}, x, y)
    end
end

export fast_fld
function fast_fld{T<:IntTypes}(x::T, y::T)
    # fld(x,y) = div(x,y) - ((x>=0) != (y>=0) && rem(x,y) != 0 ? 1 : 0)
    d = fast_div(x, y)
    d - (signbit(x $ y) & (d * y != x))
end
fast_fld{T<:UIntTypes}(x::T, y::T) = fast_div(x, y)

export fast_mod
function fast_mod{T<:IntTypes}(x::T, y::T)
    # x = fld(x,y)*y + mod(x,y)
    x - fast_fld(x, y) * y
end
fast_mod{T<:UIntTypes}(x::T, y::T) = fast_rem(x, y)

export fast_cld
function fast_cld{T<:IntTypes}(x::T, y::T)
    # cld(x,y) = div(x,y) + ((x>0) = (y>0) && rem(x,y) != 0 ? 1 : 0)
    d = fast_div(x, y)
    d + (((x > 0) == (y > 0)) & (d * y != x))
end
function fast_cld{T<:UIntTypes}(x::T, y::T)
    d = fast_div(x, y)
    d + (d * y != x)
end

end


text:
using Flux
using Flux: onehot, onehotbatch, crossentropy, reset!, throttle

corpora = Dict()

cd(@__DIR__)
for file in readdir("corpus")
  lang = Symbol(match(r"(.*)\.txt", file).captures[1])
  corpus = split(String(read("corpus/$file")), ".")
  corpus = strip.(normalize_string.(corpus, casefold=true, stripmark=true))
  corpus = filter(!isempty, corpus)
  corpora[lang] = corpus
end

langs = collect(keys(corpora))
alphabet = ['a':'z'; '0':'9'; ' '; '\n'; '_']

# See which chars will be represented as "unknown"
unique(filter(x -> x ∉ alphabet, join(vcat(values(corpora)...))))

dataset = [(onehotbatch(s, alphabet, '_'), onehot(l, langs))
           for l in langs for s in corpora[l]] |> shuffle

train, test = dataset[1:end-100], dataset[end-99:end]

N = 15

scanner = Chain(Dense(length(alphabet), N, σ), LSTM(N, N))
encoder = Dense(N, length(langs))

function model(x)
  state = scanner.(x.data)[end]
  reset!(scanner)
  softmax(encoder(state))
end

loss(x, y) = crossentropy(model(x), y)

testloss() = mean(loss(t...) for t in test)

opt = ADAM(params(scanner, encoder))
evalcb = () -> @show testloss()

Flux.train!(loss, train, opt, cb = throttle(evalcb, 10))

max_stars_repo_name: domluna/model-zoo

max_issues_repo_name: domluna/model-zoo

max_forks_repo_name: domluna/model-zoo