text:
using Base64

const PROGRESS_ID = "weave_progress"

function run_doc(
    doc::WeaveDoc;
    doctype::Union{Nothing,AbstractString} = nothing,
    out_path::Union{Symbol,AbstractString} = :doc,
    args::Any = Dict(),
    mod::Union{Module,Nothing} = nothing,
    fig_path::Union{Nothing,AbstractString} = nothing,
    fig_ext::Union{Nothing,AbstractString} = nothing,
    cache_path::AbstractString = "cache",
    cache::Symbol = :off,
)
    # cache :all, :user, :off, :refresh

    doc.doctype = isnothing(doctype) ? (doctype = detect_doctype(doc.source)) : doctype
    doc.format = deepcopy(get_format(doctype))

    cwd = doc.cwd = get_cwd(doc, out_path)
    mkpath(cwd)

    # TODO: provide a way not to create `fig_path` ?
    if isnothing(fig_path)
        fig_path = if (endswith(doctype, "2pdf") && cache === :off) || endswith(doctype, "2html")
            basename(mktempdir(abspath(cwd)))
        else
            DEFAULT_FIG_PATH
        end
    end
    mkpath(normpath(cwd, fig_path))
    # This is needed for latex and should work on all output formats
    @static Sys.iswindows() && (fig_path = replace(fig_path, "\\" => "/"))
    set_rc_params(doc, fig_path, fig_ext)

    cache === :off || @eval import Serialization # XXX: evaluate in a more sensible module

    # New sandbox for each document with args exposed
    isnothing(mod) && (mod = sandbox = Core.eval(Main, :(module $(gensym(:WeaveSandBox)) end))::Module)
    Core.eval(mod, :(WEAVE_ARGS = $(args)))

    mimetypes = doc.format.mimetypes

    report = Report(cwd, doc.basename, doc.format, mimetypes)
    cd_back = let d = pwd(); () -> cd(d); end
    cd(cwd)
    pushdisplay(report)
    try
        if cache !== :off && cache !== :refresh
            cached = read_cache(doc, cache_path)
            isnothing(cached) && @info "No cached results found, running code"
        else
            cached = nothing
        end

        executed = []
        n = length(filter(chunk->isa(chunk,CodeChunk), doc.chunks))
        i = 0
        for chunk in doc.chunks
            if chunk isa CodeChunk
                options = merge(doc.chunk_defaults, chunk.options)
                merge!(chunk.options, options)

                @info "Weaving chunk $(chunk.number) from line $(chunk.start_line)" progress=(i)/n _id=PROGRESS_ID
                i+=1
            end

            restore = (cache === :user && chunk isa CodeChunk && chunk.options[:cache])
            result_chunks = if cached ≠ nothing && (cache === :all || restore)
                restore_chunk(chunk, cached)
            else
                run_chunk(chunk, doc, report, mod)
            end
            executed = [executed; result_chunks]
        end

        replace_header_inline!(doc, report, mod) # evaluate and replace inline code in header

        doc.header_script = report.header_script
        doc.chunks = executed

        cache !== :off && write_cache(doc, cache_path)

        @isdefined(sandbox) && clear_module!(sandbox)
    catch err
        rethrow(err)
    finally
        @info "Weaved all chunks" progress=1 _id=PROGRESS_ID
        cd_back()
        popdisplay(report) # ensure display pops out even if internal error occurs
    end

    return doc
end

run_doc(doc::WeaveDoc, doctype::Union{Nothing,AbstractString}; kwargs...) =
    run_doc(doc; doctype = doctype, kwargs...)

"""
    detect_doctype(path)

Detect the output format based on file extension.
"""
function detect_doctype(path)
    _, ext = lowercase.(splitext(path))

    match(r"^\.(jl|.?md|ipynb)", ext) !== nothing && return "md2html"
    ext == ".rst" && return "rst"
    ext == ".tex" && return "texminted"
    ext == ".txt" && return "asciidoc"

    return "pandoc"
end

function get_cwd(doc, out_path)
    return if out_path === :doc
        dirname(doc.path)
    elseif out_path === :pwd
        pwd()
    else
        path, ext = splitext(out_path)
        if isempty(ext) # directory given
            path
        else # file given
            dirname(path)
        end
    end |> abspath
end

function run_chunk(chunk::CodeChunk, doc, report, mod)
    result = eval_chunk(doc, chunk, report, mod)
    occursin("2html", doc.doctype) && (embed_figures!(result, report.cwd))
    return result
end

function embed_figures!(chunk::CodeChunk, cwd)
    for (i, fig) in enumerate(chunk.figures)
        chunk.figures[i] = img2base64(fig, cwd)
    end
end
embed_figures!(chunks, cwd) = embed_figures!.(chunks, Ref(cwd))

function img2base64(fig, cwd)
    ext = splitext(fig)[2]
    f = open(joinpath(cwd, fig), "r")
    raw = read(f)
    close(f)
    return ext == ".png" ? "data:image/png;base64," * stringmime(MIME("image/png"), raw) :
           ext == ".svg" ? "data:image/svg+xml;base64," * stringmime(MIME("image/svg"), raw) :
           ext == ".gif" ? "data:image/gif;base64," * stringmime(MIME("image/gif"), raw) :
           fig
end

function run_chunk(chunk::DocChunk, doc, report, mod)
    chunk.content = [run_inline(c, doc, report, mod) for c in chunk.content]
    return chunk
end

run_inline(inline::InlineText, ::WeaveDoc, ::Report, ::Module) = inline

const INLINE_OPTIONS = Dict(
    :term => false,
    :hold => true,
    :wrap => false
)

function run_inline(inline::InlineCode, doc::WeaveDoc, report::Report, mod::Module)
    # Make a temporary CodeChunk for running code. Collect results and don't wrap
    chunk = CodeChunk(inline.content, 0, 0, "", INLINE_OPTIONS)
    options = merge(doc.chunk_defaults, chunk.options)
    merge!(chunk.options, options)

    chunks = eval_chunk(doc, chunk, report, mod)
    occursin("2html", doc.doctype) && (embed_figures!(chunks, report.cwd))

    output = chunks[1].output
    endswith(output, "\n") && (output = output[1:end-1])
    inline.output = output
    inline.rich_output = chunks[1].rich_output
    inline.figures = chunks[1].figures
    return inline
end

function run_code(doc::WeaveDoc, chunk::CodeChunk, report::Report, mod::Module)
    code = chunk.content
    path = doc.path
    error = chunk.options[:error]
    codes = chunk.options[:term] ? split_code(code) : [code]
    capture = code -> capture_output(code, mod, path, error, report)
    return capture.(codes)
end

function split_code(code)
    res = String[]
    e = 1
    ex = :init
    while true
        s = e
        ex, e = Meta.parse(code, s)
        isnothing(ex) && break
        push!(res, strip(code[s:e-1]))
    end
    return res
end

function capture_output(code, mod, path, error, report)
    reset_report!(report)

    old = stdout
    rw, wr = redirect_stdout()
    reader = @async read(rw, String)

    local out = nothing
    task_local_storage(:SOURCE_PATH, path) do
        try
            obj = include_string(mod, code, path) # TODO: fix line number
            !isnothing(obj) && !REPL.ends_with_semicolon(code) && display(obj)
        catch _err
            err = unwrap_load_err(_err)
            error || throw(err)
            display(err)
            @warn "ERROR: $(typeof(err)) occurred, including output in Weaved document"
        finally
            redirect_stdout(old)
            close(wr)
            out = fetch(reader)
            close(rw)
        end
    end

    return ChunkOutput(code, remove_ansi_control_chars(out), report.rich_output, report.figures)
end

function reset_report!(report)
    report.rich_output = ""
    report.figures = String[]
end

unwrap_load_err(err) = return err
unwrap_load_err(err::LoadError) = return err.error

# https://stackoverflow.com/a/33925425/12113178
remove_ansi_control_chars(s) = replace(s, r"(\x9B|\x1B\[)[0-?]*[ -\/]*[@-~]" => "")

function eval_chunk(doc::WeaveDoc, chunk::CodeChunk, report::Report, mod::Module)
    if !chunk.options[:eval]
        chunk.output = ""
        chunk.options[:fig] = false
        return chunk
    end

    execute_prehooks!(chunk)

    report.fignum = 1
    report.cur_chunk = chunk

    if hasproperty(report.format, :out_width) && isnothing(chunk.options[:out_width])
        chunk.options[:out_width] = report.format.out_width
    end

    chunk.result = run_code(doc, chunk, report, mod)

    execute_posthooks!(chunk)

    return chunk.options[:term] ? collect_term_results(chunk) :
           chunk.options[:hold] ? collect_hold_results(chunk) :
           collect_results(chunk)
end

# Hooks to run before and after chunks, this is form IJulia,
const preexecution_hooks = Function[]
push_preexecution_hook!(f::Function) = push!(preexecution_hooks, f)
function pop_preexecution_hook!(f::Function)
    i = findfirst(x -> x == f, preexecution_hooks)
    isnothing(i) && error("this function has not been registered in the pre-execution hook yet")
    return splice!(preexecution_hooks, i)
end
function execute_prehooks!(chunk::CodeChunk)
    for prehook in preexecution_hooks
        Base.invokelatest(prehook, chunk)
    end
end

const postexecution_hooks = Function[]
push_postexecution_hook!(f::Function) = push!(postexecution_hooks, f)
function pop_postexecution_hook!(f::Function)
    i = findfirst(x -> x == f, postexecution_hooks)
    isnothing(i) && error("this function has not been registered in the post-execution hook yet")
    return splice!(postexecution_hooks, i)
end
function execute_posthooks!(chunk::CodeChunk)
    for posthook in postexecution_hooks
        Base.invokelatest(posthook, chunk)
    end
end

"""
    clear_module!(mod::Module)

Recursively sets variables in `mod` to `nothing` so that they're GCed.

!!! warning
    `const` variables can't be reassigned, as such they can't be cleared.
"""
function clear_module!(mod::Module)
    for name in names(mod; all = true)
        name === :eval && continue
        try
            v = getfield(mod, name)
            if v isa Module && v != mod
                clear_module!(v)
                continue
            end
            isconst(mod, name) && continue # can't clear constant
            Core.eval(mod, :($name = nothing))
        catch err
            @debug err
        end
    end
end

function get_figname(report::Report, chunk; fignum = nothing, ext = nothing)
    isnothing(ext) && (ext = chunk.options[:fig_ext])
    isnothing(fignum) && (fignum = report.fignum)

    chunkid = isnothing(chunk.options[:label]) ? chunk.number : chunk.options[:label]
    basename = string(report.basename, '_', chunkid, '_', fignum, ext)
    full_name = normpath(report.cwd, chunk.options[:fig_path], basename)
    rel_name = string(chunk.options[:fig_path], '/', basename) # Relative path is used in output
    return full_name, rel_name
end

function set_rc_params(doc::WeaveDoc, fig_path, fig_ext)
    doc.chunk_defaults[:fig_ext] = isnothing(fig_ext) ? doc.format.fig_ext : fig_ext
    doc.chunk_defaults[:fig_path] = fig_path
end

function collect_results(chunk::CodeChunk)
    content = ""
    result_chunks = CodeChunk[]
    for r in chunk.result
        content *= r.code
        # Check if there is any output from chunk
        if any(!isempty ∘ strip, (r.stdout, r.rich_output)) || !isempty(r.figures)
            rchunk = CodeChunk(
                content,
                chunk.number,
                chunk.start_line,
                chunk.optionstring,
                copy(chunk.options),
            )
            rchunk.output = r.stdout
            rchunk.rich_output = r.rich_output
            rchunk.figures = r.figures
            push!(result_chunks, rchunk)
            content = ""
        end
    end
    if !isempty(content)
        rchunk = CodeChunk(
            content,
            chunk.number,
            chunk.start_line,
            chunk.optionstring,
            copy(chunk.options),
        )
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function collect_term_results(chunk::CodeChunk)
    output = ""
    prompt = chunk.options[:prompt]
    result_chunks = CodeChunk[]
    for r in chunk.result
        output *= string('\n', indent_term_code(prompt, r.code), '\n', r.stdout)
        if !isempty(r.figures)
            rchunk = CodeChunk(
                "",
                chunk.number,
                chunk.start_line,
                chunk.optionstring,
                copy(chunk.options),
            )
            rchunk.output = output
            output = ""
            rchunk.figures = r.figures
            push!(result_chunks, rchunk)
        end
    end
    if !isempty(output)
        rchunk = CodeChunk(
            "",
            chunk.number,
            chunk.start_line,
            chunk.optionstring,
            copy(chunk.options),
        )
        rchunk.output = output
        push!(result_chunks, rchunk)
    end

    return result_chunks
end

function indent_term_code(prompt, code)
    prompt_with_space = string(prompt, ' ')
    n = length(prompt_with_space)
    pads = ' ' ^ n
    return map(enumerate(split(code, '\n'))) do (i,line)
        isone(i) ? string(prompt_with_space, line) : string(pads, line)
    end |> joinlines
end

function collect_hold_results(chunk::CodeChunk)
    for r in chunk.result
        chunk.output *= r.stdout
        chunk.rich_output *= r.rich_output
        chunk.figures = [chunk.figures; r.figures]
    end
    return [chunk]
end

const HEADER_INLINE = Regex("$(HEADER_INLINE_START)(?<code>.+)$(HEADER_INLINE_END)")

replace_header_inline!(doc, report, mod) = _replace_header_inline!(doc, doc.header, report, mod)

function _replace_header_inline!(doc, header, report, mod)
    replace!(header) do (k,v)
        return k =>
            v isa Dict ? _replace_header_inline!(doc, v, report, mod) :
            !isa(v, AbstractString) ? v :
            replace(v, HEADER_INLINE => s -> begin
                code = replace(s, HEADER_INLINE => s"\g<code>")
                run_inline_code(code, doc, report, mod)
            end)
    end
    return header
end

function run_inline_code(code, doc, report, mod)
    inline = InlineCode(code, 1, :inline)
    inline = run_inline(inline, doc, report, mod)
    return strip(inline.output, '"')
end


text:
# save symbols
compress_then_write(b::Vector{Char}, io) = begin
    meta = compress_then_write(Int.(b), io)
    (type = Char, len = meta.len, metadata = meta)
end

# load a Symbol column
column_loader(::Type{Char}, io, metadata) = begin
    strs = column_loader(Int, io, metadata.metadata)
    Char.(strs)
end

column_loader!(_, ::Type{Char}, io, metadata) = begin
    column_loader(Char, io, metadata)
end

some_elm(::Type{Char}) = 'J'

ext: jl


text:
# # Example 001: Periodic Advection
#
#md # !!! jupyter
#md #     This example is also available as a Jupyter notebook:
#md #     [`ex_001_periodic_advection.ipynb`](@__NBVIEWER_ROOT_URL__examples/DGmethods_old/generated/ex_001_periodic_advection.html)
#
# Key ideas of this tutorial:
#   - Setting up PDE
#   - Defining a numerical flux
#   - Defining finite element mesh
#   - Using the ODE solving framework
#   - Using VTK visualization
#
# ## Introduction
#
# In this example we will solve the constant coefficient advection equation on a
# periodic domain; the domain is taken to be the unit square or cube depending on
# whether the problem is two- or three-dimensional.
#
# The partial differential equation we wish to solve is
#
# ```math
#  \frac{\partial q}{\partial t} + \nabla \cdot (\vec{u} q) = 0,
# ```
#
# where $q$ is the advected scalar quantity and $\vec{u}$ is the constant
# velocity field. The quantity $\vec{u} q$ is more generally called the flux and
# denoted in the tutorial below as $\boldsymbol{F}(q) = \vec{u} q$.
#
# Below is a program interspersed with comments.
#md # The full program, without comments, can be found in the next
#md # [section](@ref ex_001_periodic_advection-plain-program).
#
# ## Commented Program

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

# ### Preliminaries
# Load in modules needed for solving the problem
using MPI
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VTK
using LinearAlgebra
using Logging
using Dates
using Printf
using StaticArrays

# Start up MPI if this has not already been done
MPI.Initialized() || MPI.Init()
#md nothing # hide

# define the velocity field for advection
const uvec = (1, 2, 3)
#md nothing # hide

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

# ### Physical Flux
#md # Now we define a function which given a value for $q$ computes the physical
#md # flux $\boldsymbol{F} = \vec{u} q$. Since we only have a single state, $q$,
#md # the state will come in as an `MVector` of length `1` and the flux to fill
#md # as an `MMatrix` of size `(3, 1)`.
#
#md # Note: for two-dimensional simulations the flux `MMatrix` will also be of
#md # size `(3, 1)` but the function need only fill the first two rows
function advectionflux!(F, state, _...)
  FT = eltype(state) # get the floating point type we are using
  @inbounds begin
    q = state[1]
    F[:, 1] = SVector{3, FT}(uvec) * q
  end
end
#md nothing # hide

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

# ### Numerical Flux
# In the discontinuous Galerkin method the continuity of the solution across
# element interfaces is imposed weakly through the use of a numerical flux. The
# numerical flux is a function that given the solution state on either side of
# an interface returns a unique value approximating $\boldsymbol{F}\cdot
# \vec{n}$ on the interface:
# ```math
# f^{*} = f^{*}(q^{-}, q^{+}; \vec{n}).
# ```
# We call the two sides of the interface the "minus side" and "plus side" hence
# the $\pm$ superscripts in $q$ in the notation above; the choice of which
# element is on the minus side and plus side is arbitrary. Here $\vec{n}$ is a
# unit normal to the interface and is typically taken to point from the minus
# side to the plus side. The numerical flux is required to be symmetric with
# respect to $q^{-}$ and $q^{+}$:
# ```math
# f^{*}(q^{-}, q^{+}; \vec{n}) = f^{*}(q^{+}, q^{-}; \vec{n}).
# ```
# and consistent with the physical flux in the sense that
# ```math
# \boldsymbol{F}(q) \cdot \vec{n} = f^{*}(q, q; \vec{n}).
# ```
# Taken together these two conditions also imply that
# ```math
# f^{*}(q^{-}, q^{+}; \vec{n}^{-}) = -f^{*}(q^{+}, q^{-}; \vec{n}^{+}),
# ```
# that is $f^{*}$ is skew-symmetry with respect to the unit normals $\vec{n}^{+}
# = -\vec{n}^{-}$.  The choice of numerical flux has important implications for
# the stability of the method. Though it is beyond the scope of these tutorials
# to dive into this in detail, often the "best" numerical fluxes are the ones
# that are constructed specifically for given set of equations by solving the
# Riemann problem; the Riemann problem is an initial value problem where the
# initial data are piecewise constant with a single discontinuity.
#
# In the CLIMA balance law solver the numerical flux function is a user-defined
# function that fills in an `MVector` for the numerical flux given two states,
# the "viscous state", a unit normal to the face, the simulation time, and a
# user-defined auxiliary state; the viscous and auxiliary states will be
# discussed in a subsequent examples. In the function below `F` is the numerical
# flux to fill, `nM` is the unit normal pointing away from the minus side and
# toward the plus side, `QM` and `QP` are `MVector`s of the solution state on
# the minus and plus sides of the interface, `viscM` and `viscP` are the viscous
# states on the minus and plus sides, `auxM` and `auxP` are the user-defined
# auxiliary state values on the minus and plus sides, and `t` is the simulation
# time.
#
# For linear advection the solution to the Riemann problem is trivial, since
# if $\vec{n}^{-} \cdot \vec{u} ≥ 0$ the state $q$ is being advected from
# the minus side to the plus side otherwise the reverse occurs. Thus the upwind
# numerical flux for advection is
# ```math
# f^{*}(q^{-}, q^{+}; \vec{n}^{-}) =
# \begin{cases}
#   \vec{n}^{-}\cdot\vec{u} \; q^{-}, &\text{ if } \vec{n} \cdot \vec{u} ≥ 0,\\
#   \vec{n}^{-}\cdot\vec{u} \; q^{+}, &\text{ otherwise.}
# \end{cases}
# ```
# This is done in the following function
function upwindflux!(fs, nM, stateM, viscM, auxM, stateP, viscP, auxP, t)
  FT = eltype(fs)
  @inbounds begin
    ## determine the advection speed and direction
    un = dot(nM, FT.(uvec))
    qM = stateM[1]
    qP = stateP[1]
    ## Determine which state is "upwind" of the minus side
    fs[1] = un ≥ 0 ? un * qM : un * qP
  end
end
#md nothing # hide

# In later examples we will demonstrate how to use the Rusanov flux which is
# included in the `CLIMA.DGBalanceLawDiscretizations.NumericalFluxes` submodule.
# This is a more general-purpose flux which approximates the solution Riemann
# problem by using an average of the flux on either side of the interface with
# additional dissipation added based on the local wave speed.

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

# ### Initial Condition
# In this example we take the initial condition to be
# ```math
# q(\vec{x}, t=0) = \prod_{i=1}^{d} \exp(\sin(2\pi x_{i})),
# ```
# where $d$ is the dimensionality of the problem. To use this initial condition
# we need a function that given $\vec{x}$ sets $q$.
#
# The initial condition is set by the solver through a function which takes $q$
# as an `MVector` to initialize based on the pointwise coordinate location
# `(x_1, x_2, x_3)`.
#
# Note: The initial condition will always be called as though the dimensionality
# of the problem is 3. For the domain used below `x_3 = 0` when the problem is
# actually two-dimensional; since when $x_3 = 0$ the function $\exp(\sin(2\pi
# x_{3})) = 1$ we can safely assume the dimensionality is always $3$ in our
# implication of the initial condition.
#
# Note: The last argument needs to be caught but not used for this example
function initialcondition!(Q, x_1, x_2, x_3, _...)
  @inbounds Q[1] = exp(sin(2π * x_1)) * exp(sin(2π * x_2)) * exp(sin(2π * x_3))
end
#md nothing # hide

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

# ### Exact Solution
# For periodic constant-velocity advection the exact solution is trivial to
# compute. Assuming that $\phi(x)$ is the periodically replicated initial
# condition, the analytic solution is
# ```math
# q(\vec{x}, t) = \phi(\vec{x} - \vec{u} t).
# ```
# This will be useful later since it will allow us to check our work by
# computing the error in our solution and estimating the convergence rate.
#
# For a general initial condition on the unit domain the following function can
# be used:
function exactsolution!(dim, Q, t, x_1, x_2, x_3, _...)
  @inbounds begin
    FT = eltype(Q)

    ## trace back the point (x_1, x_2, x_3) in the velocity field and
    ## determine where in our "original" [0, L_1] X [0, L_2] X [0, L_3] domain
    ## this point is located
    y_1 = mod(x_1 - FT(uvec[1]) * t, 1)
    y_2 = mod(x_2 - FT(uvec[2]) * t, 1)

    ## if we are really just 2-D we do not want to change the x_3 coordinate
    y_3 = dim == 3 ? mod(x_3 - FT(uvec[3]) * t, 1) : x_3

    initialcondition!(Q, y_1, y_2, y_3)
  end
end
#md nothing # hide
# The input argument `dim` is the "real" dimensionality of the problem and is
# needed in case `uvec[3] != 0`.

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

# ### Initialize the DG Method
# We are now at the point that we can initialize the structure for the DG
# method.  For convenience we define a function that initializes the DG method
# over a given MPI communicator `mpicomm`, for a given `polynomialorder`, using
# `dim` dimensions. The mesh used will be `Ne` x `Ne` elements when `dim == 2`
# and `Ne` x `Ne` x `Ne` elements when `dim == 3`. The floating point type of
# the computation is `FT` and whether the CPU or GPU is used is determined
# by `ArrayType`; `ArrayType === Array` is for the CPU and `ArrayType ===
# CuArray` is for NVIDIA GPUs.
#
# Note: This whole code chunk is in a function block
function setupDG(mpicomm, dim, Ne, polynomialorder, FT=Float64,
                 ArrayType=Array)

  @assert ArrayType === Array

  # We will use the `BrickTopology` from `CLIMA.Mesh.Topologies` to define the mesh.
  # The "topology" in CLIMA is the element connectivity information (e.g.,
  # neighbouring elements and interface data) along with coordinate locations
  # for corners of the elements. The `BrickTopology` creates a regular mesh of a
  # rectangular or regular hexahedral domain. This is done by specifying the
  # coordinate points of the element corners along each dimension. Here, we want
  # to mesh the unit square or cube with `Ne` elements in each dimension, thus
  # we specify the following `Tuple` of `range`s:
  brickrange = (range(FT(0); length=Ne+1, stop=1), # x_1 corner locations
                range(FT(0); length=Ne+1, stop=1), # x_2 corner locations
                range(FT(0); length=Ne+1, stop=1)) # x_3 corner locations
  # these coordinates will be combined in a tensor product fashion to define the
  # element corners.
  #
  # By default the `BrickTopology` is not periodic, so we need to define a
  # `Tuple` of boolean defining which dimensions are periodic.
  periodicity = (true, true, true)
  # Note: We have defined both the `brickrange` and `periodicity` as though we
  # are working in three-dimensions, and in when in two-dimensions we will
  # discard the third element of the `Tuple`; this could also be done from the
  # start by using the `ntuple` function.

  # Using `brickrange` and `periodicity` we can now initialize the topology
  # using the `BrickTopology` constructor. This will both create the topology as
  # well as do the partitioning of the elements across the MPI ranks available
  # in the mpi communicator `mpicomm`
  topology = BrickTopology(mpicomm, brickrange[1:dim];
                           periodicity=periodicity[1:dim])

  # The topology only has element connectivity and corner information, thus we
  # still need to create a grid (or mesh) of degrees of freedom. In CLIMA the
  # so-called discontinuous spectral element grid is used (aka tensor-product
  # quadrilateral and hexahedral elements with Legendre-Gauss-Lobatto
  # interpolation and quadrature weights). Given a topology and polynomial order
  # we can create the grid of degrees of free using
  grid = DiscontinuousSpectralElementGrid(topology; polynomialorder =
                                          polynomialorder, FloatType = FT,
                                          DeviceArray = ArrayType,)
  # Note: This constructor also takes in a `FloatType` which specifies the
  # floating point type to be used (e.g., for the coordinate points and geometry
  # metric terms). The argument `ArrayType` is used to determine the compute
  # device to use (i.e., `Array` will signify the CPU is being used and
  # `CuArray` will signal that an NVIDIA GPU is being used).

  # We can now define the discretization from the `grid` using the physical flux
  # `advectionflux!` and numerical flux `upwindflux!` defined above; we only
  # have a single state variable $q$ hence `length_state_vector = 1`
  spatialdiscretization = DGBalanceLaw(grid = grid, length_state_vector = 1,
                                       flux! = advectionflux!,
                                       numerical_flux! = upwindflux!)

  # (end of function)
end
#md nothing # hide

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

# ### Initializing and run the DG method
# Note: This whole code chunk is in a `let` block
let
  # We will just use the whole MPI communicator
  mpicomm = MPI.COMM_WORLD

  # Since this is an MPI enabled code, we use the Julia logging functionality to
  # ensure that only one MPI rank prints to the screen. Namely, MPI rank 0 does
  # all the logging and all other ranks dump their ouput to `devnull`.
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)
  # Note: The `NullLogger` should not be used because if any MPI functions are
  # called in a logging block deadlock will occur since `NullLogger` code is
  # not executed.

  # Dimensionality to run
  dim = 2

  # Mesh size along each dimension
  Ne = 20

  # order of polynomials to use
  polynomialorder = 4

  # Setup the DG discretization
  spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder)

  # Given the `spatialdiscretization` and the `initialcondition!` function we
  # can create and initialize storage for the solution. This is an MPI-aware
  # array
  Q = MPIStateArray(spatialdiscretization, initialcondition!)

  # A VTK file, which can be viewed in [ParaView](https://www.paraview.org/) or
  # [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit), can be
  # generated with the following command (the last `Tuple` gives strings for the
  # names of the state fields)
  filename = @sprintf("initialcondition_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))
  # Note: Currently the `writevtk` function writes one file for each MPI rank,
  # and the user is responsible for opening each of the files to "stitch"
  # together the image.

  # In order to run the simulation we need to use an ODE solver. In this example
  # we will use a low storage Runge-Kutta method which can be initialized with a
  # spatial discretization, solution vector (not stored but used to define
  # needed auxiliary arrays), initial solution time, and a time step size; this
  # particular Runge-Kutta method uses a fixed time step.
  #
  # Since we are using a regular mesh, with a constant wave speed, a "CFL"
  # restriction for the mesh is
  h = 1 / Ne                           # element size
  CFL = h / maximum(abs.(uvec[1:dim])) # time to cross the element once
  dt = CFL / polynomialorder^2         # DG time step scaling (for this
                                       ## particular RK scheme could go with a
                                       ## factor of ~2 larger time step)
  lsrk = LSRK54CarpenterKennedy(spatialdiscretization, Q; dt = dt, t0 = 0)

  # Here we run the ODE solver until the final time `timeend` using `Q` as the
  # initial condition `Q`. The solution will be updated in place so that the
  # final solution will also be stored in `Q`.
  finaltime = 1.0
  if (parse(Bool, lowercase(get(ENV,"TRAVIS","false")))       #src
      && "Test" == get(ENV,"TRAVIS_BUILD_STAGE_NAME","")) ||  #src
    parse(Bool, lowercase(get(ENV,"APPVEYOR","false")))       #src
    finaltime = 2dt                                           #src
  end                                                         #src
  solve!(Q, lsrk; timeend = finaltime)

  # The final solution can be visualized in a similar manner to the initial
  # condition
  filename = @sprintf("finalsolution_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))

  # Using the `finaltime` and `exactsolution!` we can calculate the exact
  # solution
  Qe = MPIStateArray(spatialdiscretization) do Qin, x, y, z, aux
    exactsolution!(dim, Qin, finaltime, x, y, z)
  end
  # and then compute the error by evaluating the Euclidean distance between the
  # computed solution `Q` and the exact solution `Qe`
  error = euclidean_distance(Q, Qe)
  with_logger(mpi_logger) do
    @info @sprintf("""Run with
                   dim              = %d
                   Ne               = %d
                   polynomial order = %d
                   error            = %e
                   """, dim, Ne, polynomialorder, error)
  end
end
#md nothing # hide

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

# ### Using ODE solver callback functions
# The above simulation run with `solve!` runs from the initial time to the final
# time. The ODE solver framework in CLIMA gives functionality that allows the
# user to *inject* code into the solver during execution. Here we show how to
# use some of the generic callback functions provided to
#
#  - Save diagnostic information
#  - display runtime simulation information
#  - save VTK files during the simulation
#
# Note: This whole code chunk is in a `let` block
let
  # code is the same as above until the `solve!` call
  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)
  dim = 2
  Ne = 20
  polynomialorder = 4
  spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder)
  Q = MPIStateArray(spatialdiscretization, initialcondition!)
  filename = @sprintf("initialcondition_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))
  h = 1 / Ne
  CFL = h / maximum(abs.(uvec[1:dim]))
  dt = CFL / polynomialorder^2
  lsrk = LSRK54CarpenterKennedy(spatialdiscretization, Q; dt = dt, t0 = 0)
  finaltime = 1.0
  if (parse(Bool, lowercase(get(ENV,"TRAVIS","false")))       #src
      && "Test" == get(ENV,"TRAVIS_BUILD_STAGE_NAME","")) ||  #src
    parse(Bool, lowercase(get(ENV,"APPVEYOR","false")))       #src
    finaltime = 2dt                                           #src
  end                                                         #src

  # The ODE solver callback functions are called both before the ODE solver
  # begins and then after each time step.
  #
  # For instance if a user wanted to store the norm of the solution every time
  # step the following callback could be used
  store_norm_index = 0
  normQ = Array{Float64}(undef, ceil(Int, finaltime / dt))
  function cb_store_norm()
    store_norm_index += 1
    normQ[store_norm_index] = norm(Q)
    nothing
  end
  # Note: that callbacks must return either `nothing` or `0` if the ODE solver
  # should continue, `1` is the ODE solver should stop after all the callbacks
  # have executed, or `2` is the time stepping should immediately stop with no
  # further callbacks executed.

  # Several generic callbacks are provided in the `CLIMA.GenericCallbacks`
  # submodule. For instance, the `EveryXSimulationSteps` callbacks will execute
  # every `X` time steps. This could be used to say write VTK output every `20`
  # time steps
  vtk_step = 0
  mkpath("vtk")
  cb_vtk = GenericCallbacks.EveryXSimulationSteps(20) do
    vtk_step += 1
    filename = @sprintf("vtk/advection_mpirank%04d_step%04d",
                         MPI.Comm_rank(mpicomm), vtk_step)
    writevtk(filename, Q, spatialdiscretization,
                                         ("q",))
    nothing
  end

  # Another provided generic callback is `EveryXWallTimeSeconds` which will be
  # called every `X` seconds of wall clock time (as opposed to simulation time).
  # This could be used to dump diagnostic information about the simulation. In
  # this case we display the norm of the simulation time, the run time, and the
  # norm of the solution.
  #
  # One unique feature of this call back is that it takes in a single optional
  # argument `init` which allows the ODE solver to call the callback for
  # initialization; all callbacks get called for initialization with a single
  # boolean argument set to `true`, but this occurs in a `try/catch` statement
  # in case the callback does not require initialization (such as the two
  # above).
  starttime = Ref(now())
  cb_info = GenericCallbacks.EveryXWallTimeSeconds(1, mpicomm) do (init=false)
    if init
      starttime[] = now()
    else
      with_logger(mpi_logger) do
        @info @sprintf("""Update
                       simtime = %.16e
                       runtime = %s
                       norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                       Dates.format(convert(Dates.DateTime,
                                            Dates.now()-starttime[]),
                                    Dates.dateformat"HH:MM:SS"),
                       norm(Q))
      end
    end
  end
  # Note that this callback also takes in the MPI communicator. This is
  # necessary because the callback needs to execute an `MPI.Allreduce` to ensure
  # that all the MPI ranks are using the same global run time.

  # the defined callbacks are based to the ODE `solve!` function through the
  # keyword argument `callbacks` as a tuple:
  solve!(Q, lsrk; timeend = finaltime,
         callbacks = (cb_store_norm, cb_vtk, cb_info))

  # The remainder of the function is the same as above
  filename = @sprintf("finalsolution_mpirank%04d", MPI.Comm_rank(mpicomm))
  writevtk(filename, Q, spatialdiscretization,
                                       ("q",))

  Qe = MPIStateArray(spatialdiscretization) do Qin, x, y, z, aux
    exactsolution!(dim, Qin, finaltime, x, y, z)
  end
  error = euclidean_distance(Q, Qe)
  with_logger(mpi_logger) do
    @info @sprintf("""Run with
                   dim              = %d
                   Ne               = %d
                   polynomial order = %d
                   error            = %e
                   """, dim, Ne, polynomialorder, error)
  end
end
#md nothing # hide

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

# ### Computing rates and errors
# If the above code is put in a loop over increasing `Ne` then a rate of
# convergence for the scheme can be established which we expect to be
# on the order of the polynomial order plus $\sim 1/2$.
let
  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)

  dim = 2
  polynomialorder = 4
  finaltime = 1.0

  with_logger(mpi_logger) do
    @info @sprintf("""Running with
                   dim              = %d
                   polynomial order = %d
                   """, dim, polynomialorder)
  end

  base_Ne = 5
  lvl_error = zeros(4) # number of levels to compute is length(lvl_error)
  for lvl = 1:length(lvl_error)
    ## `Ne` for this mesh level
    Ne = base_Ne * 2^(lvl-1)
    spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder)

    Q = MPIStateArray(spatialdiscretization, initialcondition!)
    h = 1 / Ne
    CFL = h / maximum(abs.(uvec[1:dim]))
    dt = CFL / polynomialorder^2
    if (parse(Bool, lowercase(get(ENV,"TRAVIS","false")))       #src
        && "Test" == get(ENV,"TRAVIS_BUILD_STAGE_NAME","")) ||  #src
      parse(Bool, lowercase(get(ENV,"APPVEYOR","false")))       #src
      finaltime = 2dt                                           #src
    end                                                         #src
    lsrk = LSRK54CarpenterKennedy(spatialdiscretization, Q; dt = dt, t0 = 0)

    solve!(Q, lsrk; timeend = finaltime)

    Qe = MPIStateArray(spatialdiscretization) do Qin, x, y, z, aux
      exactsolution!(dim, Qin, finaltime, x, y, z)
    end

    lvl_error[lvl] = euclidean_distance(Q, Qe)
    msg =  @sprintf   "Level      = %d" lvl
    msg *= @sprintf "\nNe               = %d" Ne
    msg *= @sprintf "\nerror            = %.4e" lvl_error[lvl]
    if lvl > 1
      rate = log2(lvl_error[lvl-1]) - log2(lvl_error[lvl])
      msg *= @sprintf "\nconvergence rate = %.4e" rate
    end
    with_logger(mpi_logger) do
      @info msg
    end
  end
end
#md nothing # hide

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

#md # ## [Plain Program](@id ex_001_periodic_advection-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [ex\_001\_periodic\_advection.jl](ex_001_periodic_advection.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
struct Image
    id::Integer
    name::String
    imagetype::String
    distribution::String
    slug::Union{Nothing, String}
    public::Bool
    regions::Array{String, 1}
    mindisksize::Integer
    sizegigabytes::Real
    created_at::DateTime

    function Image(data::Dict{String})
        # we assume all DO datetimes are in UTC
        data["created_at"] = DateTime(data["created_at"][1:end-1])

        new(
            data["id"],
            data["name"],
            data["type"],
            data["distribution"],
            data["slug"],
            data["public"],
            data["regions"],
            data["min_disk_size"],
            data["size_gigabytes"],
            data["created_at"]
        )
    end
end

function show(io::IO, i::Image)
    print(io, "Image ($(i.distribution) $(i.name))")
end

function getallimages!(client::AbstractClient)
    uri = joinpath(ENDPOINT, "images")
    data = getdata!(client, uri)
    images = Array{Image, 1}(UndefInitializer(), length(data))

    for (i, image) in enumerate(data)
        images[i] = Image(image)
    end
end


text:
using Parameters
using LinearAlgebra, StatsBase

sigmoid(x::Matrix, β::Vector) = @. 1 / (1 + exp(-$*(x,β)))

"""
    QDA(μ = μ, Σ = Σ)
Quadratic discrimination analysis

"""
@with_kw struct QDA
    μ::Array # 平均
    Σ::Matrix # 分散共分散行列
    invΣ::Matrix = inv(Σ) # 分散共分散行列の逆行列
    detΣ = det(Σ) # 分散共分散行列の行列式
end

"""
    function-like object of QDA

returns log likelihood of multivariate normal distribution

Example
```julia

```
"""
function (qda::QDA)(x...;prior=1.0)
    data = collect(x)
    @unpack μ,invΣ,detΣ = qda
    a = -0.5*(data-μ)' * invΣ * (data-μ)
    a[1] - log(detΣ) + log(prior)
end

function table_count(test,pred)
    @assert length(test) == length(pred)
    m = unique(test) |> length
    count = zeros(Int,m,m)
    for i in 1:length(test)
        count[test[i],pred[i]] += 1
    end
    return count
end

function knn(X_train::Matrix,y_train::Vector,X_test::Vector, k)
    n = size(X_train)[1]
    distance = [norm(X_train[i,:]-X_test) for i in 1:n]
    S = sortperm(distance)[1:k]
    u = counts(y_train[S], 1:k)
    u_max = maximum(u)
    m = findall(c->c==u_max,u)
    while length(m)!==1
        k -=1
        S = S[1:k]
        u = counts(S, 1:k)
        u_max = maximum(u)
        m = findall(c->c==u_max,u)
    end
    return m[1]
end

function knn(X_train::Matrix,y_train::Vector,X_test::Matrix, k)
    l = size(X_test)[1]
    w = Array{Int}(undef,l)
    for i in 1:l
        w[i] = knn(X_train,y_train,X_test[i,:], k)
    end
    return w
end


text:
# SimLynx.jl Event-Based Discrete-Event Simulation Example

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

struct Customer
    id::Int64
    Customer(id::Int64) = new(id)
end

Base.show(io::IO, customer::Customer) =
    print(io, "Customer($(customer.id))")

mutable struct Teller
    id::Int64
    serving::Union{Customer, Nothing}
    Teller(id::Int64) = new(id, nothing)
end

Base.show(io::IO, teller::Teller) =
    print(io, "Teller($(teller.id))")

tellers = nothing
teller_queue = nothing

function available_teller()::Union{Teller, Nothing}
    for teller in tellers
        if isnothing(teller.serving)
            return teller
        end
    end
    return nothing
end

generator_dist = Exponential(MEAN_INTERARRIVAL_TIME)

@event generate(i::Integer, n::Integer) begin
    @schedule now arrival(Customer(i))
    if i < n
        @schedule in rand(generator_dist) generate(i + 1, n)
    end
end

@event arrival(customer::Customer) begin
    teller = available_teller()
    if isnothing(teller)
        enqueue!(teller_queue, customer)
    else
        @schedule now service(teller, customer)
    end
end

@event service(teller::Teller, customer::Customer) begin
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    teller.serving = customer
    @schedule in rand(dist) departure(teller, customer)
end

@event departure(teller::Teller, customer::Customer) begin
    if isempty(teller_queue)
        teller.serving = nothing
    else
        next_customer = dequeue!(teller_queue)
        teller.serving = next_customer
        @schedule now service(teller, next_customer)
    end
end

function run_simulation(n::Integer)
    @assert n > 0
    println("SimLynx.jl Event-Based Discrete-Event Simulation Example")
    println("$N_TELLERS Teller, Single Queue Bank Model")
    println("  Inter-arrival time = Exponential($MEAN_INTERARRIVAL_TIME)")
    println("  Service time = Uniform($MIN_SERVICE_TIME, $MAX_SERVICE_TIME)")
    println("  Number of customers = $n")
    @simulation begin
        global tellers = [Teller(i) for i = 1:N_TELLERS]
        global teller_queue = FifoQueue{Customer}()
        current_trace!(true)
        @schedule at 0.0 generate(1, n)
        start_simulation()
        print_stats(teller_queue.n, title="Teller Queue Length Statistics")
        plot_history(teller_queue.n, title="Teller Queue Length History")
    end
end

run_simulation(10)


text:
@enum FPtyp  IEEE=1 IBM DEC
@enum ENDIAN BIG=1 LITTLE

expectedblksz(typ::Int32) = (typ in [4, 16]) ? 8 : ((typ in [3, 11]) ? 4 : 1)

function readinfo(io::IO, infodict)
    subtyp = read(io, Int32)
    blksz = read(io, Int32)
    blksz == expectedblksz(subtyp) || error("subtyp $subtyp gave blksz $blksz")
    contents = read(io, blksz * read(io, Int32))
    if subtyp == 3
        v = reinterpret(Int32, contents)
        infodict[:version] = VersionNumber(v[1], v[2], v[3])
        infodict[:FPrep] = FPtyp(v[5])
        infodict[:endian] = ENDIAN(v[7])
        infodict[:charactercode] = v[8]
    elseif subtyp == 4
        v = reinterpret(Float64, contents)
        infodict[:sysmiss] = v[1]
        infodict[:highest] = v[2]
        infodict[:lowest] = v[3]
    elseif subtyp == 10
        infodict[:product] = String(contents)
    elseif subtyp == 11
        infodict[:display] = reinterpret(Int32, contents)
    elseif subtyp == 13
        infodict[:longvarnames] =
            Dict{Symbol,Symbol}(map(x -> Symbol.(x),
                               filter(x -> length(x) == 2 && x[1] ≠ x[2],
                               split.(split(String(contents), '\t'), '='))))
    elseif subtyp == 14
        infodict[:stringlengths] = contents
    elseif subtyp == 16
        infodict[:ncases64] = reinterpret(Int64, contents)
    elseif subtyp == 17
        infodict[:datafileattributes] = String(contents)
    elseif subtyp == 18
        infodict[:variableattributes] = Dict{Symbol,String}(
            map(x -> (Symbol(x[1]), x[2]), split.(split(String(contents), '/'), ':')))
    elseif subtyp == 20
        infodict[:encoding] = String(contents)
    elseif subtyp == 21
        infodict[:longstringvalue] = contents
    elseif subtyp == 22
        infodict[:longstringmissing] = contents
    else
        error("unknown info subtype: ", subtyp)
    end
end

function readextensionrecords(io::IO, vdict::VariableDictionary)
    rectyp = Int32(7)
    infodict = Dict{Symbol,Any}()
    while rectyp == 7
        readinfo(io, infodict)
        rectyp = read(io, Int32)
    end
    return rectyp, infodict
end


text:
import Dates: TimeType
# Functions to be implemented, so that an array type
# can be used as an YAXArray. Here fakkback methods are implemented
# for the AbstractArray's
"""
    dimvals(x, i)

Returns the values associated to the i-th dimensions of the
array `x`
"""
dimvals(x::AbstractArray, i) = valfromaxis(axes(x,i))

"""
    dimname(x,i)

Returns the name of the i-th dimension of the array x as a symbol.
"""
dimname(::AbstractArray,i) = Symbol("Dim_",i)
#Mandatory interface ends here

# Optional methods
"""
    iscontdim(x, i)

Returns a boolean indicating if the i-th dimension of the array
should be interpreted as a continouus range (true) or as categorical (false).
"""
function iscontdim(x, i)
  iscontdimval(dimvals(x,i))
end
iscontdimval(v) = (eltype(v) <: Union{Number,TimeType}) && (issorted(v) || issorted(v,rev=true))

"""
    iscompressed(x)

Returns a boolean indicating if the underlying data is compressed
and therefore random access is slow.
"""
iscompressed(x) = false

"""
    dimnames(x)

Returns a tuple of dimension names for al dimensions of x.
"""
dimnames(x) = ntuple(i->dimname(x,i), ndims(getdata(x)))

"""
    getattributes(x)

If the data type supports storing additional attributes,
this returns a Dict with key-value pairs where the attribute names
(keys) are encoded as Strings.
"""
getattributes(x) = Dict{String,Any}()

"""
    getdata(x)

Returns a handle to the raw data of the array, which can be indexed
as a normal Julia Array. This might return the object itself, but ideally
if returns only the AbstractArray/DiskArray wrapped by the type.
"""
getdata(x) = x

"""
    valfromaxis(x)

If for a data type the dimension values are encoded in an object returned by `Base.axes(x)`
one can define a method for this function that returns the dimension values from that axis.
"""
valfromaxis(x) = x

"""
    yaxconvert(T::Type,x)

Converts an AbstractArray x that implements the interface to type T.
"""
function yaxconvert(T::Type{<:Any},x)
  data = getdata(x)
  yaxcreate(T, data, dimnames(x),dimvals.(Ref(x),1:ndims(data)),getattributes(x))
end

"""
    yaxcreate(T::Type,data,dnames,dvals,attributes)

Creates a new array with the given dimension names, values and attributes.
"""
function yaxcreate(T,data,dname,dvals,attributes) end


text:
@testset "inflate" begin
    input = Int8[-3, 3, 2, 4, 5]
    output = fill(0x00, 5)

    inflate!(Balloons.Packbits(), input, output)

    @test all(output .== [0x03, 0x03, 0x03, 0x04, 0x05])
end


text:
# Using PyCall is rather slow. Please compare to https://github.com/JuliaReinforcementLearning/ReinforcementLearningEnvironmentClassicControl.jl/blob/master/examples/cartpole.jl
using ReinforcementLearningEnvironmentGym, ReinforcementLearning

env = GymEnv("CartPole-v0")
rlsetup = RLSetup(ActorCriticPolicyGradient(ns = 4, na = 2, α = .02,
                                            nsteps = 25),
                  env, ConstantNumberSteps(400),
                  callbacks = [Visualize(wait = 0.)])
@info("Before learning.")
run!(rlsetup)
rlsetup.callbacks = [EvaluationPerEpisode(TotalReward())]
rlsetup.stoppingcriterion = ConstantNumberSteps(10^5)
@time learn!(rlsetup)
getvalue(rlsetup.callbacks[1])
@info("After learning.")
rlsetup.callbacks = [Visualize(wait = 0.)]
rlsetup.stoppingcriterion = ConstantNumberSteps(400)
run!(rlsetup)


text:
using LowRank, Test
using Random, Distributions

N = 50
d = 5
n = 10
k = 2

# non-zero error
A = zeros(d, N)
U = rand(Uniform(-1.0, 1.0), d-1, k)
V = rand(Uniform(-1.0, 1.0), k, n)
A[1:d-1, 1:n] = U*V
A[d,1:n] = rand(Uniform(-1.0, 1.0), 1, n)
A[:, n+1:N] = rand(Uniform(-1.0, 1.0), d , N-n)
error, _ = LowRank.truncated_svd(A[:,1:n], k)
data_set = collect(eachcol(A))

@test find_low_rank_subset_sample_rep(data_set,n,k)[1] <= 4*error
@test FastLowRank(data_set, n, k)[1] <= error + 0.0001


text:
# A table presents itself as a vector of named tuple

struct Table{names, T <: NamedTuple{names}, Vs <: Tuple{Vararg{AbstractVector}}, As <: Tuple{Vararg{AbstractIndex}}} <: AbstractVector{T}
    data::Vs
    accelerators::As
end

Table(;kwargs...) = Table(kwargs.data) # TODO `index` keyword argument?
Table(nt::NamedTuple{names}, index::Tuple{Vararg{AbstractIndex}} = ()) where {names} = Table{names}(Tuple(nt), index)

Table{names}(;kwargs...) where {names} = Table{names}(kwargs.data) # TODO `index` keyword argument?

function Table{names}(data::Tuple{Vararg{AbstractVector}}, index::Tuple{Vararg{AbstractIndex}} = ()) where {names}
    T = _eltypes(data)
    return Table{names, NamedTuple{names, T}, typeof(data), typeof(index)}(data, index)
end

function Table{names}(nt::NamedTuple{names, <:Tuple{Vararg{AbstractVector}}}, index::Tuple{Vararg{AbstractIndex}} = ()) where {names}
    Table{names}(Tuple(nt), index)
end

function Table{names}(nt::NamedTuple{names2, <:Tuple{Vararg{AbstractVector}}}, index::Tuple{Vararg{AbstractIndex}} = ()) where {names, names2}
    if _issetequal(names, names2)
        data = getindices(nt, names)
        return Table{names}(data, index)
    else
        error("Table column names $names do not match input data with names $names2")
    end
end

# Helpers to get the data directly from the Table struct
accelerators(::AbstractVector{<:NamedTuple}) = ()
accelerators(t::Table) = Core.getfield(t, :accelerators)
getdata(t::Table) = Core.getfield(t, :data)

# Simple column access via `table.columnname`
@inline Base.getproperty(t::Table, name::Symbol) = getdata(t)[_find(colnames(t), name)]

colnames(::AbstractArray{<:NamedTuple{names}}) where {names} = names
columns(t::Table{names}) where {names} = NamedTuple{names}(getdata(t))

@inline size(t::Table) = size(first(getdata(t)))
@inline axes(t::Table) = axes(first(getdata(t)))

@generated function getindex(t::Table{names}, i::Int) where {names}
    exprs = [:($(names[j]) = getdata(t)[$j][i]) for j = 1:length(names)]
    return quote
        @_propagate_inbounds_meta
        $(Expr(:tuple, exprs...))
    end
end

@generated function getindex(t::Table{names}, i::AbstractVector{Int}) where {names}
    exprs = [:(getdata(t)[$j][i]) for j = 1:length(names)]
    return quote
        @_propagate_inbounds_meta
        return Table{names}($(Expr(:tuple, exprs...)), ())  # Indexing always removes acceleration indexes, for now
    end
end

function getindex(t::Table{names}, ::Colon) where {names}
    return Table{names}(copy.(getdata(t)), ()) # Indexing always removes acceleration indexes, for now (compare with `copy`)
end

@generated function setindex!(t::Table{names}, v::NamedTuple{names2}, i) where {names, names2}
    if !issetequal(names, names2)
        return quote
            error("Attempted to assign named tuple with names $names2 to table with names $names")
        end
    end

    # TODO - rebuild indices? or make them immutable? or let the user "use with care"?

    exprs = [:(getdata(t)[$j][i] = getproperty(v, $(Expr(:quote, names[j])))) for j = 1:length(names)]
    return quote
        @_propagate_inbounds_meta
        if accelerators(t) !== ()
            error("Mutating tables with acceleration indexes currently unsupported")
        end
        $(Expr(:block, exprs...))
        return v
    end
end

function similar(t::Table, ::Type{NamedTuple{names, Ts}}, dims::Tuple{Int}) where {names, Ts}
    data = _makevectors(Ts, dims)
    return Table{names}(data)
end

function copy(t::Table{names}) where {names}
    return Table{names}(copy.(getdata(t)), copy.(accelerators(t)))
end

@inline function (p::Project{names})(t::Table) where {names}
    data = p(columns(t))
    indexes = p(accelerators(t))
    return Table{names}(Tuple(data), indexes)
end

@inline function (r::Rename{oldnames, newnames})(t::Table{names}) where {oldnames, newnames, names}
    names2 = _rename(Val(oldnames), Val(newnames), Val(names))
    return Table{names2}(getdata(t), r(accelerators(t)))
end


text:
function initialize!(integrator,cache::SSPRK22ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK22ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator

  # u1 -> stored as u
  u = uprev + dt*integrator.fsalfirst
  k = f(u, p, t+dt)
  # u
  u = (uprev + u + dt*k) / 2

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK22Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK22Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,stage_limiter!,step_limiter! = cache

  # u1 -> stored as u
  @. u = uprev + dt*fsalfirst
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  # u
  @. u = (uprev + u + dt*k) / 2
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::LDDRK64ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

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

  # u1
  tmp = dt*integrator.fsalfirst
  u   = uprev + β1*tmp
  # u2
  tmp = α2*tmp + dt*f(u, p, t+c2*dt)
  u   = u + β2*tmp
  # u3
  tmp = α3*tmp + dt*f(u, p, t+c3*dt)
  u   = u + β3*tmp
  # u4
  tmp = α4*tmp + dt*f(u, p, t+c4*dt)
  u   = u + β4*tmp
  # u5 = u
  tmp = α5*tmp + dt*f(u, p, t+c5*dt)
  u   = u + β5*tmp
  # u6
  tmp = α6*tmp + dt*f(u, p, t+c6*dt)
  u = u + β6*tmp

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::LDDRK64Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

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

  # u1
  @. tmp = dt*fsalfirst
  @. u   = uprev + β1*tmp
  # u2
  f( k,  u, p, t+c2*dt)
  @. tmp = α2*tmp + dt*k
  @. u   = u + β2*tmp
  # u3
  f( k,  u, p, t+c3*dt)
  @. tmp = α3*tmp + dt*k
  @. u   = u + β3*tmp
  # u4
  f( k,  u, p, t+c4*dt)
  @. tmp = α4*tmp + dt*k
  @. u   = u + β4*tmp
  # u5 = u
  f( k,  u, p, t+c5*dt)
  @. tmp = α5*tmp + dt*k
  @. u   = u + β5*tmp

  f( k,  u, p, t+c6*dt)
  @. tmp = α6*tmp + dt*k
  @. u   = u + β6*tmp

  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK33ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK33ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator

  # u1
  u = uprev + dt*integrator.fsalfirst
  k = f(u, p, t+dt)
  # u2
  u = (3*uprev + u + dt*k) / 4
  k = f(u,p,t+dt/2)
  # u
  u = (uprev + 2*u + 2*dt*k) / 3

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK33Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK33Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,stage_limiter!,step_limiter! = cache

  # u1
  @. u = uprev + dt*fsalfirst
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  # u2
  @. u = (3*uprev + u + dt*k) / 4
  stage_limiter!(u, f, t+dt/2)
  f(k,u,p,t+dt/2)
  # u
  @. u = (uprev + 2*u + 2*dt*k) / 3
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK53ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK53ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache

  # u1
  tmp = uprev + β10 * dt * integrator.fsalfirst
  k = f(tmp, p, t+c1*dt)
  # u2 -> stored as u
  u = tmp + β21 * dt * k
  k = f(u, p, t+c2*dt)
  # u3
  tmp = α30 * uprev + α32 * u + β32 * dt * k
  k = f(tmp, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α43 * tmp + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u
  u = α52 * u + α54 * tmp + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK53Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK53Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache.tab

  # u1
  @. tmp = uprev + β10 * dt * fsalfirst
  stage_limiter!(tmp, f, t+c1*dt)
  f( k,  tmp, p, t+c1*dt)
  # u2 -> stored as u
  @. u = tmp + β21 * dt * k
  stage_limiter!(u, f, t+c2*dt)
  f( k,  u, p, t+c2*dt)
  # u3
  @. tmp = α30 * uprev + α32 * u + β32 * dt * k
  stage_limiter!(tmp, f, t+c3*dt)
  f( k,  tmp, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α43 * tmp + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u
  @. u = α52 * u + α54 * tmp + β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK53_2N1ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK53_2N1ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α40,α43,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache
  #stores in u for all intermediate stages
  # u1
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2
  u = u + β21 * dt * k
  k = f(u, p, t+c2*dt)
  # u3
  u = u + β32 * dt * k
  k = f(u, p, t+c3*dt)
  # u4
  u= α40 * uprev + α43 * u + β43 * dt * k
  k = f(u, p, t+c4*dt)
  # u
  u =  u + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK53_2N1Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK53_2N1Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α40,α43,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache.tab
  #stores in u for all intermediate stages
  # u1
  @. u = uprev + β10 * dt * fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2
  @. u = u + β21 * dt * k
  stage_limiter!(u, f, t+c2*dt)
  f( k,  u, p, t+c2*dt)
  # u3
  @. u = u + β32 * dt * k
  stage_limiter!(u, f, t+c3*dt)
  f( k,  u, p, t+c3*dt)
  # u4
  @. u = α40 * uprev + α43 * u + β43 * dt * k
  stage_limiter!(u, f, t+c4*dt)
  f( k,  u, p, t+c4*dt)
  # u
  @. u = u + β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK53_2N2ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK53_2N2ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α30,α32,α50,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache
  #stores in u for all intermediate stages
  # u1
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2 -> stored as u
  u = u + β21 * dt * k
  k = f(u, p, t+c2*dt)
  # u3
  u = α30 * uprev + α32 * u + β32 * dt * k
  k = f(u, p, t+c3*dt)
  # u4
  u = u + β43 * dt * k
  k = f(u, p, t+c4*dt)
  # u
  u = α50 * uprev + α54 * u + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK53_2N2Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK53_2N2Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α30,α32,α50,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache.tab

  # u1
  @. u = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2 -> stored as u
  @. u = u + β21 * dt * k
  stage_limiter!(u, f, t+c2*dt)
  f( k,  u, p, t+c2*dt)
  # u3
  @. u = α30 * uprev + α32 * u + β32 * dt * k
  stage_limiter!(u, f, t+c3*dt)
  f( k,  u, p, t+c3*dt)
  # u4
  @. u = u + β43 * dt * k
  stage_limiter!(u, f, t+c4*dt)
  f( k,  u, p, t+c4*dt)
  # u
  @. u = α50* uprev + α54 * u+ β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK63ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK63ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5 = cache

  # u1 -> stored as u
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2
  u₂ = u + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u3
  tmp = u₂ + β32 * dt * k
  k = f(tmp, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α41 * u + α43 * tmp + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u
  u = α62 * u₂ + α65 * tmp + β65 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK63Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK63Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₂,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5 = cache.tab

  # u1 -> stored as u
  @. u = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2
  @. u₂ = u + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u3
  @. tmp = u₂ + β32 * dt * k
  stage_limiter!(tmp, f, t+c3*dt)
  f( k,  tmp, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α41 * u + α43 * tmp + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u
  @. u = α62 * u₂ + α65 * tmp + β65 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK73ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK73ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6 = cache

  # u1
  u₁ = uprev + β10 * dt * integrator.fsalfirst
  k = f(u₁,p,t+c1*dt)
  # u2
  tmp = u₁ + β21 * dt * k
  k = f(tmp, p, t+c2*dt)
  # u3 -> stored as u
  u = tmp + β32 * dt * k
  k = f(u, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α43 * u + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = α50 * uprev + α51 * u₁ + α54 * tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u6
  tmp = tmp + β65 * dt * k
  k = f(tmp, p, t+c6*dt)
  # u
  u = α73 * u + α76 * tmp + β76 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK73Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK73Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₁,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6 = cache.tab

  # u1
  @. u₁ = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u₁, f, t+c1*dt)
  f(k,u₁,p,t+c1*dt)
  # u2
  @. tmp = u₁ + β21 * dt * k
  stage_limiter!(tmp, f, t+c2*dt)
  f( k,  tmp, p, t+c2*dt)
  # u3 -> stored as u
  @. u = tmp + β32 * dt * k
  stage_limiter!(u, f, t+c3*dt)
  f( k,  u, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α43 * u + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = α50 * uprev + α51 * u₁ + α54 * tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u6
  @. tmp = tmp + β65 * dt * k
  stage_limiter!(tmp, f, t+c6*dt)
  f( k,  tmp, p, t+c6*dt)
  # u
  @. u = α73 * u + α76 * tmp + β76 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK83ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

@muladd function perform_step!(integrator,cache::SSPRK83ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7 = cache

  # u1 -> save as u
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2
  u₂ = u + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u3
  u₃ = u₂ + β32 * dt * k
  k = f(u₃,p,t+c3*dt)
  # u4
  tmp = u₃ + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = α50 * uprev + α51 * u + α54 * tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u6
  tmp = α61 * u + α65 * tmp + β65 * dt * k
  k = f(tmp, p, t+c6*dt)
  # u7
  tmp = α72 * u₂ + α73 * u₃ + α76 * tmp + β76 * dt * k
  k = f(tmp, p, t+c7*dt)
  # u
  u = tmp + β87 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK83Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK83Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₂,u₃,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7 = cache.tab

  # u1 -> save as u
  @. u = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2
  @. u₂ = u + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u3
  @. u₃ = u₂ + β32 * dt * k
  stage_limiter!(u₃, f, t+c3*dt)
  f(k,u₃,p,t+c3*dt)
  # u4
  @. tmp = u₃ + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = α50 * uprev + α51 * u + α54 * tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u6
  @. tmp = α61 * u + α65 * tmp + β65 * dt * k
  stage_limiter!(tmp, f, t+c6*dt)
  f( k,  tmp, p, t+c6*dt)
  # u7
  @. tmp = α72 * u₂ + α73 * u₃ + α76 * tmp + β76 * dt * k
  stage_limiter!(tmp, f, t+c7*dt)
  f( k,  tmp, p, t+c7*dt)
  # u
  @. u = tmp + β87 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

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

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

@muladd function perform_step!(integrator,cache::SSPRK432ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_2 = dt / 2

  # u1
  u = uprev + dt_2*integrator.fsalfirst
  k = f(u, p, t+dt_2)
  # u2
  u = u + dt_2*k
  k = f(u, p, t+dt)
  u = u + dt_2*k
  if integrator.opts.adaptive
    utilde = (uprev + 2*u) / 3
  end
  # u3
  u = (2*uprev + u) / 3
  k = f(u, p, t+dt_2)
  # u
  u = u + dt_2*k

  integrator.fsallast = f(u, p, t+dt)
  if integrator.opts.adaptive
    utilde = utilde - u
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

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

@muladd function perform_step!(integrator,cache::SSPRK432Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,utilde,atmp,stage_limiter!,step_limiter! = cache
  dt_2 = dt / 2

  # u1
  @. u = uprev + dt_2*fsalfirst
  stage_limiter!(u, f, t+dt_2)
  f( k,  u, p, t+dt_2)
  # u2
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  #
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt+dt_2)
  if integrator.opts.adaptive
    @. utilde = (uprev + 2*u) / 3
  end
  # u3
  @. u = (2*uprev + u) / 3
  f( k,  u, p, t+dt_2)
  #
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)

  if integrator.opts.adaptive
    @. utilde = utilde - u
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f( k,  u, p, t+dt)
end

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

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

@muladd function perform_step!(integrator,cache::SSPRKMSVS32ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack u_1,u_2,dts,dtf,μ,v_n = cache

  if integrator.iter == 1
    cache.dts[1] = dt
    cache.dts[2] = dt
    cache.dtf[1] = dt
  end
  accpt = true
  dt = dts[1]

  if cache.step < 3 #starting Procedure
    k = f(u,p,t+dt)
    u = uprev + dt*k
    k = f(u, p, t+dt)
    u = (uprev + u + dt*k) / 2
      if cache.step == 1
        u_2 = uprev
      else
        u_1 = uprev
      end
    if integrator.opts.adaptive
      v_n = dt/dts[2]*0.5
      cache.dtf[2] = dtf[1]
      cache.dtf[1] = dt/v_n*0.5
      if v_n > 0.5
        cache.step -= 1
        accpt = false
      end
      cache.dts[3] = dts[2]
      cache.dts[2] = dt
      dt = 0.9*dtf[1]
      μ = min(dtf[1],dtf[2])

    end
  else
    if integrator.opts.adaptive
      Ω = (dts[2] + dts[3])/dt
    else
      Ω = 2
    end
    u = (Ω*Ω - 1)/(Ω*Ω)*(uprev + Ω/(Ω-1)*dt*integrator.fsalfirst) + 1/(Ω*Ω)*u_2
    u_2 = u_1
    u_1 = uprev
    if integrator.opts.adaptive
      v_n = (dts[2]+dts[3]-dt)/(dts[2]+dts[3])*0.5
      dt = (dts[2] + dts[3])/(dts[2]+dts[3]+μ)*μ
      cache.dtf[2] = dtf[1]
      dtf[1]  = dt/v_n*0.5
      cache.dts[3] = dts[2]
      cache.dts[2] = dt
      μ = min(dtf[1],dtf[2])
    end
  end
  if accpt == true
    integrator.fsallast = f(u, p, t+dt)
    integrator.k[1] = integrator.fsalfirst
    integrator.u = u
  else
    integrator.fsallast = f(uprev, p, t+dt)
    integrator.k[1] = integrator.fsalfirst
    integrator.u = uprev
  end
  cache.dts[1] = dt
  cache.step += 1
  cache.u_1 = u_1
  cache.u_2 = u_2
  cache.μ = μ

end

function initialize!(integrator,cache::SSPRKMSVS32Cache)
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
end

@muladd function perform_step!(integrator,cache::SSPRKMSVS32Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,u_1,u_2,stage_limiter!,step_limiter! = cache

  println(cache.fsalfirst)
  if cache.step < 3
    @. u = uprev + dt*fsalfirst
    stage_limiter!(u, f, t+dt)
    f(k,u, p, t+dt)
    @. u = (uprev + u + dt*k) / 2
    stage_limiter!(u, f, t+dt)
    step_limiter!(u, f, t+dt)

    if cache.step == 1
      cache.u_2 .= uprev
    else
      cache.u_1 .= uprev
    end
  else
      Ω = 2
    @. u = ((Ω*Ω - 1)/(Ω*Ω))*(uprev + (Ω/(Ω-1))*dt*fsalfirst) + (1/(Ω*Ω))*cache.u_2
    cache.u_2 .= u_1
    cache.u_1 .= uprev
    print
    stage_limiter!(u, f, t+dt)
    step_limiter!(u, f, t+dt)
  end
  cache.step += 1
  f(k,u, p, t+dt)

end

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

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

@muladd function perform_step!(integrator,cache::SSPRKMSVS43ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack u_1,u_2,u_3,k1,k2,k3 = cache

  if cache.step < 4
    u = uprev + dt*integrator.fsalfirst
    k = f(u, p, t+dt)
    u = (uprev + u + dt*k) / 2
      if cache.step == 1
        u_3 = uprev
        cache.k3 = f(u_3,p,t+dt)
      end
      if cache.step == 2
        u_2 = uprev
        cache.k2 = f(u_2,p,t+dt)
      end
      if cache.step == 3
        u_1 = uprev
        cache.k1 = f(u_1,p,t+dt)
      end
  # u
  else
    u = (16/27)*(uprev + 3*dt*integrator.fsalfirst) + (11/27)*(u_3 + (12/11)*dt*k3)
    cache.k3 = k2
    cache.k2 = k1
    cache.k1 = integrator.fsalfirst
    u_3 = u_2
    u_2 = u_1
    u_1 = uprev
  end
  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
  cache.step += 1
  cache.u_1 = u_1
  cache.u_2 = u_2
  cache.u_3 = u_3
end

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

@muladd function perform_step!(integrator,cache::SSPRKMSVS43Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,u_1,u_2,u_3,stage_limiter!,step_limiter!,k1,k2,k3 = cache
  if cache.step < 4
    @. u = uprev + dt*fsalfirst
    stage_limiter!(u, f, t+dt)
    f(k,u, p, t+dt)
    @. u = (uprev + u + dt*k) / 2
    stage_limiter!(u, f, t+dt)
    step_limiter!(u, f, t+dt)
    if cache.step == 1
      cache.u_3 .= uprev
      f(k3,u_3,p,t+dt)
    end
    if cache.step == 2
      cache.u_2 .= uprev
      f(k2,u_2,p,t+dt)
    end
    if cache.step == 3
      cache.u_1 .= uprev
      f(k1,u_1,p,t+dt)
    end
  # u
  else
    @. u = (16/27)*(uprev + 3*dt*fsalfirst) + (11/27)*(u_3 + (12/11)*dt*k3)
    stage_limiter!(u, f, t+dt)
    step_limiter!(u, f, t+dt)
    cache.k3 .= k2
    cache.k2 .= k1
    cache.k1 .= fsalfirst
    cache.u_3 .= u_2
    cache.u_2 .= u_1
    cache.u_1 .= uprev
  end
  cache.step += 1
  f( k, u, p, t+dt)
end

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

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

@muladd function perform_step!(integrator,cache::SSPRK932ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_6 = dt / 6
  dt_3 = dt / 3
  dt_2 = dt / 2

  # u1
  u = uprev + dt_6*integrator.fsalfirst
  k = f(u, p, t+dt_6)
  # u2
  u = u + dt_6*k
  k = f(u, p, t+dt_3)
  # u3
  u = u + dt_6*k
  k = f(u, p, t+dt_2)
  # u4
  u = u + dt_6*k
  k = f(u, p, t+2*dt_3)
  # u5
  u = u + dt_6*k
  k = f(u, p, t+5*dt_6)
  # u6
  u = u + dt_6*k
  if integrator.opts.adaptive
    k = f(u, p, t+dt)
    utilde = (uprev + 6*u + 6*dt*k) / 7
  end
  # u6*
  u = (3*uprev + dt_2*integrator.fsalfirst + 2*u) / 5
  k = f(u, p, t+dt_2)
  # u7*
  u = u + dt_6*k
  k = f(u, p, t+2*dt_3)
  # u8*
  u = u + dt_6*k
  k = f(u, p, t+5*dt_6)
  # u
  u = u + dt_6*k

  integrator.fsallast = f(u, p, t+dt)
  if integrator.opts.adaptive
    utilde = utilde - u
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

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

@muladd function perform_step!(integrator,cache::SSPRK932Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,utilde,atmp,stage_limiter!,step_limiter! = cache
  dt_6 = dt / 6
  dt_3 = dt / 3
  dt_2 = dt / 2

  # u1
  @. u = uprev + dt_6*fsalfirst
  stage_limiter!(u, f, t+dt_6)
  f( k,  u, p, t+dt_6)
  # u2
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt_3)
  f( k,  u, p, t+dt_3)
  # u3
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt_2)
  f( k,  u, p, t+dt_2)
  # u4
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+2*dt_3)
  f( k,  u, p, t+2*dt_3)
  # u5
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+5*dt_6)
  f( k,  u, p, t+5*dt_6)
  # u6
  @. u = u + dt_6*k
  if integrator.opts.adaptive
    stage_limiter!(u, f, t+dt)
    f( k,  u, p, t+dt)
    @. utilde = (uprev + 6*u + 6*dt*k) / 7
    stage_limiter!(utilde, f, t+dt)
  end
  # u6*
  @. u = (3*uprev + dt_2*fsalfirst + 2*u) / 5
  stage_limiter!(u, f, t+dt_6)
  f( k,  u, p, t+dt_2)
  # u7*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+2*dt_3)
  f( k,  u, p, t+2*dt_3)
  # u8*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+5*dt_6)
  f( k,  u, p, t+5*dt_6)
  # u9*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)

  if integrator.opts.adaptive
    @. utilde = utilde - u
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK54ConstantCache)
  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::SSPRK54ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4 = cache

  # u₁
  u₂ = uprev + β10 * dt * integrator.fsalfirst
  k = f(u₂,p,t+c1*dt)
  # u₂
  u₂ = α20 * uprev + α21 * u₂ + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u₃
  u₃ = α30 * uprev + α32 * u₂ + β32 * dt * k
  k₃ = f(u₃,p,t+c3*dt)
  # u₄ -> stored as tmp
  tmp = α40 * uprev + α43 * u₃ + β43 * dt * k₃
  k = f(tmp, p, t+c4*dt)
  # u
  u = α52 * u₂ + α53 * u₃ + β53 * dt * k₃ + α54 * tmp + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

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

@muladd function perform_step!(integrator,cache::SSPRK54Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,k₃,u₂,u₃,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4 = cache.tab

  # u₁
  @. u₂ = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u₂, f, t+c1*dt)
  f(k,u₂,p,t+c1*dt)
  # u₂
  @. u₂ = α20 * uprev + α21 * u₂ + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u₃
  @. u₃ = α30 * uprev + α32 * u₂ + β32 * dt * k
  stage_limiter!(u₃, f, t+c3*dt)
  f(k₃,u₃,p,t+c3*dt)
  # u₄ -> stored as tmp
  @. tmp = α40 * uprev + α43 * u₃ + β43 * dt * k₃
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u
  @. u = α52 * u₂ + α53 * u₃ + β53 * dt * k₃ + α54 * tmp + β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK104ConstantCache)
  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::SSPRK104ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_6 = dt/6
  dt_3 = dt/3
  dt_2 = dt/2

  tmp = uprev + dt_6 * integrator.fsalfirst # u₁
  k = f(tmp, p, t+dt_6)
  tmp = tmp + dt_6 * k # u₂
  k = f(tmp, p, t+dt_3)
  tmp = tmp + dt_6 * k # u₃
  k = f(tmp, p, t+dt_2)
  u = tmp + dt_6 * k # u₄
  k₄ = f(u, p, t+2*dt_3)
  tmp = (3*uprev + 2*u + dt_3 * k₄) / 5 # u₅
  k = f(tmp, p, t+dt_3)
  tmp = tmp + dt_6 * k # u₆
  k = f(tmp, p, t+dt_2)
  tmp = tmp + dt_6 * k # u₇
  k = f(tmp, p, t+2*dt_3)
  tmp = tmp + dt_6 * k # u₈
  k = f(tmp, p, t+5*dt_6)
  tmp = tmp + dt_6 * k # u₉
  k = f(tmp, p, t+dt)
  u = (uprev + 9*(u + dt_6*k₄) + 15*(tmp + dt_6*k)) / 25

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

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

@muladd function perform_step!(integrator,cache::SSPRK104Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,k₄,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  dt_6 = dt/6
  dt_3 = dt/3
  dt_2 = dt/2

  @. tmp = uprev + dt_6 * integrator.fsalfirst
  stage_limiter!(tmp, f, t+dt_6)
  f( k,  tmp, p, t+dt_6)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_3)
  f( k,  tmp, p, t+dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_2)
  f( k,  tmp, p, t+dt_2)
  @. u = tmp + dt_6 * k
  stage_limiter!(u, f, t+2*dt_3)
  f(k₄,u,p,t+2*dt_3)
  @. tmp = (3*uprev + 2*u + 2*dt_6 * k₄) / 5
  stage_limiter!(tmp, f, t+dt_3)
  f( k,  tmp, p, t+dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_2)
  f( k,  tmp, p, t+dt_2)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+2*dt_3)
  f( k,  tmp, p, t+2*dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+5*dt_6)
  f( k,  tmp, p, t+5*dt_6)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt)
  f( k,  tmp, p, t+dt)
  @. u = (uprev + 9*(u + dt_6*k₄) + 15*(tmp + dt_6*k)) / 25
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f(k, u, p, t+dt)
end

function initialize!(integrator,cache::ORK256ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  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
end

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

  # u1
  tmp = dt*integrator.fsalfirst
  u   = uprev + β1*tmp
  # u2
  tmp = α2*tmp + dt*f(u, p, t+c2*dt)
  u   = u + β2*tmp
  # u3
  tmp = α3*tmp + dt*f(u, p, t+c3*dt)
  u   = u + β3*tmp
  # u4
  tmp = α4*tmp + dt*f(u, p, t+c4*dt)
  u   = u + β4*tmp
  # u
  tmp = α5*tmp + dt*f(u, p, t+c5*dt)
  u   = u + β5*tmp

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::ORK256Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::ORK256Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,tmp = cache
  @unpack α2,α3,α4,α5,β1,β2,β3,β4,β5,c2,c3,c4,c5 = cache.tab

  # u1
  @. tmp = dt*fsalfirst
  @. u   = uprev + β1*tmp
  # u2
  f( k,  u, p, t+c2*dt)
  @. tmp = α2*tmp + dt*k
  @. u   = u + β2*tmp
  # u3
  f( k,  u, p, t+c3*dt)
  @. tmp = α3*tmp + dt*k
  @. u   = u + β3*tmp
  # u4
  f( k,  u, p, t+c4*dt)
  @. tmp = α4*tmp + dt*k
  @. u   = u + β4*tmp
  # u
  f( k,  u, p, t+c5*dt)
  @. tmp = α5*tmp + dt*k
  @. u   = u + β5*tmp

  f( k,  u, p, t+dt)
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e6d6879e-048f-11eb-1328-2d8ae06f1f1f
# Include package.
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.envdir(), "dev-env"))
	using FractalTools
	using Plots
	plotly()
	default(:label, "")
	default(:lw, 2)
end

# ╔═╡ 2df6b090-0490-11eb-1093-ff3605ecef25
# Define function
# f(x) = sin(2π*1/2*x)
# f(x) = sin(2π*x + π/6)
f(x) = 4 - (x - 2)^2 + 1

# ╔═╡ 3f42bc90-0490-11eb-3570-3bafd44e6d9d
# Interpolation points
begin
	xi, dx, xf = 0, 1,  4
	x = collect(xi : dx : xf)
	y = f.(x)
	d = 0.5 * ones(length(x) - 1)
end

# ╔═╡ 4abce91a-0490-11eb-0841-9d7d87776e2c
# Construct interpolation
begin
	interp = interpolate(x, y, d)
	ws = interp.ifs.ws
end

# ╔═╡ 595c870a-0490-11eb-12fa-09f822a8d78d

function transform_plot(dval, niters)
	d = dval * ones(length(x) - 1)
	xc = [0, 4, 4, 0, 0]
	yc = [0, 0, 4, 4, 0]
	pc = [xc yc]
	xt = collect(xi : 0.1dx : xf)
	plt = plot()
	plot!(pc[:, 1], pc[:, 2], label="initial_domain")
	plot!(xt, f.(xt), label="true_function")
	for niter in 1 : niters
		pc = vcat([mapslices(wt, pc, dims=2) for wt in ws]...)
		plot!(pc[:, 1], pc[:, 2], label="tranformed_domain")
	end
	interp = interpolate(x, y, d, niter=niters)
	plot!(xt, interp.(xt))
	scatter!(x, y, marker=(:circle, 3))
	plt
end

# ╔═╡ 5e0c510c-04a1-11eb-3ba6-37f073caacfa
@bind dval html"<input type=range min=0 max=1 step=0.1>"

# ╔═╡ cb8afd34-0490-11eb-2d48-b3718511cf5d
@bind niters html"<input type=range min=1 max=5 step=1>"

# ╔═╡ f7275e38-0490-11eb-333d-5be95205a2f4
niters, dval

# ╔═╡ bddf0a98-04b0-11eb-0f06-37f3b02e5f1d
transform_plot(dval, niters)

# ╔═╡ Cell order:
# ╠═e6d6879e-048f-11eb-1328-2d8ae06f1f1f
# ╠═2df6b090-0490-11eb-1093-ff3605ecef25
# ╠═3f42bc90-0490-11eb-3570-3bafd44e6d9d
# ╠═4abce91a-0490-11eb-0841-9d7d87776e2c
# ╠═595c870a-0490-11eb-12fa-09f822a8d78d
# ╠═5e0c510c-04a1-11eb-3ba6-37f073caacfa
# ╠═cb8afd34-0490-11eb-2d48-b3718511cf5d
# ╠═f7275e38-0490-11eb-333d-5be95205a2f4
# ╠═bddf0a98-04b0-11eb-0f06-37f3b02e5f1d