text:
# Interpolate data from input format to desired output format (vtu version)
function interpolate_data(::Val{:vtu}, input_data, mesh::Trixi.TreeMesh, n_visnodes, verbose)
  # Calculate equidistant output nodes (visualization nodes)
  dx = 2 / n_visnodes
  nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_visnodes))

  return raw2interpolated(input_data, nodes_out)
end

# Interpolate data from input format to desired output format (CurvedMesh or UnstructuredQuadMesh version)
function interpolate_data(::Val{:vtu}, input_data,
                          mesh::Union{Trixi.CurvedMesh, Trixi.UnstructuredQuadMesh},
                          n_visnodes, verbose)
  # Calculate equidistant output nodes
  nodes_out = collect(range(-1, 1, length=n_visnodes))

  return raw2interpolated(input_data, nodes_out)
end

# Interpolate data from input format to desired output format (vti version)
function interpolate_data(::Val{:vti}, input_data, mesh, n_visnodes, verbose)
  coordinates, levels, center_level_0, length_level_0 = extract_mesh_information(mesh)

  # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
  normalized_coordinates = similar(coordinates)
  for element_id in axes(coordinates, 2)
    @views normalized_coordinates[:, element_id] .= (
        (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 ))
  end

  # Determine level-wise resolution
  max_level = maximum(levels)
  resolution = n_visnodes * 2^max_level

  # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate
  # level-0-cell) that contains the number of visualization nodes for any
  # refinement level to visualize on an equidistant grid
  nvisnodes_per_level = [2^(max_level - level)*n_visnodes for level in 0:max_level]

  # Interpolate unstructured DG data to structured data
  structured_data = unstructured2structured(input_data, normalized_coordinates, levels,
                                            resolution, nvisnodes_per_level)

  return structured_data
end

# Interpolate unstructured DG data to structured data (cell-centered)
function unstructured2structured(unstructured_data::AbstractArray{Float64},
                                 normalized_coordinates::AbstractArray{Float64},
                                 levels::AbstractArray{Int}, resolution::Int,
                                 nvisnodes_per_level::AbstractArray{Int})
  # Extract number of spatial dimensions
  ndims_ = size(normalized_coordinates, 1)

  # Extract data shape information
  n_nodes_in, _, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Calculate interpolation vandermonde matrices for each level
  max_level = length(nvisnodes_per_level) - 1
  vandermonde_per_level = []
  for l in 0:max_level
    n_nodes_out = nvisnodes_per_level[l + 1]
    dx = 2 / n_nodes_out
    nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out))
    push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out))
  end

  # For each element, calculate index position at which to insert data in global data structure
  lower_left_index = element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level)

  # Create output data structure
  structured = Array{Float64}(undef, resolution, resolution, n_variables)

  # For each variable, interpolate element data and store to global data structure
  for v in 1:n_variables
    # Reshape data array for use in interpolate_nodes function
    reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements)

    for element_id in 1:n_elements
      # Extract level for convenience
      level = levels[element_id]

      # Determine target indices
      n_nodes_out = nvisnodes_per_level[level + 1]
      first = lower_left_index[:, element_id]
      last = first .+ (n_nodes_out - 1)

      # Interpolate data
      vandermonde = vandermonde_per_level[level + 1]
      structured[first[1]:last[1], first[2]:last[2], v] .= (
          reshape(interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1),
                  n_nodes_out, n_nodes_out))
    end
  end

  # Return as one 1D array for each variable
  return reshape(structured, resolution^ndims_, n_variables)
end

# For a given normalized element coordinate, return the index of its lower left
# contribution to the global data structure
function element2index(normalized_coordinates::AbstractArray{Float64}, levels::AbstractArray{Int},
                       resolution::Int, nvisnodes_per_level::AbstractArray{Int})
  # Extract number of spatial dimensions
  ndims_ = size(normalized_coordinates, 1)

  n_elements = length(levels)

  # First, determine lower left coordinate for all cells
  dx = 2 / resolution
  lower_left_coordinate = Array{Float64}(undef, ndims_, n_elements)
  for element_id in 1:n_elements
    nvisnodes = nvisnodes_per_level[levels[element_id] + 1]
    lower_left_coordinate[1, element_id] = (
        normalized_coordinates[1, element_id] - (nvisnodes - 1)/2 * dx)
    lower_left_coordinate[2, element_id] = (
        normalized_coordinates[2, element_id] - (nvisnodes - 1)/2 * dx)
  end

  # Then, convert coordinate to global index
  indices = coordinate2index(lower_left_coordinate, resolution)

  return indices
end

# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1])
function coordinate2index(coordinate, resolution::Integer)
  # Calculate 1D normalized coordinates
  dx = 2/resolution
  mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))

  # Find index
  id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], lt=(x,y)->x .< y .- dx/2)
  id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], lt=(x,y)->x .< y .- dx/2)
  return transpose(hcat(id_x, id_y))
end

# Interpolate to specified output nodes
function raw2interpolated(data_gl::AbstractArray{Float64}, nodes_out)
  # Extract number of spatial dimensions
  ndims_ = ndims(data_gl) - 2

  # Extract data shape information
  n_nodes_in = size(data_gl, 1)
  n_nodes_out = length(nodes_out)
  n_elements = size(data_gl, ndims_ + 1)
  n_variables = size(data_gl, ndims_ + 2)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Calculate Vandermonde matrix
  vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)

  if ndims_ == 2
    # Create output data structure
    data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_elements, n_variables)

    # For each variable, interpolate element data and store to global data structure
    for v in 1:n_variables
      # Reshape data array for use in interpolate_nodes function
      @views reshaped_data = reshape(data_gl[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements)

      # Interpolate data to visualization nodes
      for element_id in 1:n_elements
        @views data_vis[:, :, element_id, v] .= reshape(
            interpolate_nodes(reshaped_data[:, :, :, element_id], vandermonde, 1),
            n_nodes_out, n_nodes_out)
      end
    end
  elseif ndims_ == 3
    # Create output data structure
    data_vis = Array{Float64}(undef, n_nodes_out, n_nodes_out, n_nodes_out, n_elements, n_variables)

    # For each variable, interpolate element data and store to global data structure
    for v in 1:n_variables
      # Reshape data array for use in interpolate_nodes function
      @views reshaped_data = reshape(data_gl[:, :, :, :, v],
                                     1, n_nodes_in, n_nodes_in, n_nodes_in, n_elements)

      # Interpolate data to visualization nodes
      for element_id in 1:n_elements
        @views data_vis[:, :, :, element_id, v] .= reshape(
            interpolate_nodes(reshaped_data[:, :, :, :, element_id], vandermonde, 1),
            n_nodes_out, n_nodes_out, n_nodes_out)
      end
    end
  else
    error("Unsupported number of spatial dimensions: ", ndims_)
  end

  # Return as one 1D array for each variable
  return reshape(data_vis, n_nodes_out^ndims_ * n_elements, n_variables)
end

# Interpolate data using the given Vandermonde matrix and return interpolated values (2D version).
function interpolate_nodes(data_in::AbstractArray{T, 3},
                           vandermonde, n_vars) where T
  n_nodes_out = size(vandermonde, 1)
  data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out)
  interpolate_nodes!(data_out, data_in, vandermonde, n_vars)
end

function interpolate_nodes!(data_out::AbstractArray{T, 3}, data_in::AbstractArray{T, 3},
                            vandermonde, n_vars) where T
  n_nodes_out = size(vandermonde, 1)
  n_nodes_in  = size(vandermonde, 2)

  for j in 1:n_nodes_out
    for i in 1:n_nodes_out
      for v in 1:n_vars
        acc = zero(eltype(data_out))
        for jj in 1:n_nodes_in
          for ii in 1:n_nodes_in
            acc += vandermonde[i, ii] * data_in[v, ii, jj] * vandermonde[j, jj]
          end
        end
        data_out[v, i, j] = acc
      end
    end
  end

  return data_out
end

# Interpolate data using the given Vandermonde matrix and return interpolated values (3D version).
function interpolate_nodes(data_in::AbstractArray{T, 4},
                           vandermonde, n_vars) where T
  n_nodes_out = size(vandermonde, 1)
  data_out = zeros(eltype(data_in), n_vars, n_nodes_out, n_nodes_out, n_nodes_out)
  interpolate_nodes!(data_out, data_in, vandermonde, n_vars)
end

function interpolate_nodes!(data_out::AbstractArray{T, 4}, data_in::AbstractArray{T, 4},
                            vandermonde, n_vars) where T
  n_nodes_out = size(vandermonde, 1)
  n_nodes_in  = size(vandermonde, 2)

  for k in 1:n_nodes_out, j in 1:n_nodes_out, i in 1:n_nodes_out
    for v in 1:n_vars
      acc = zero(eltype(data_out))
      for kk in 1:n_nodes_in, jj in 1:n_nodes_in, ii in 1:n_nodes_in
        acc += vandermonde[i, ii] * vandermonde[j, jj] * vandermonde[k, kk] * data_in[v, ii, jj, kk]
      end
      data_out[v, i, j, k] = acc
    end
  end

  return data_out
end


text:
module magpoly
include("structures_and_functions.jl")

#############################################################################################################
#############################################################################################################

"""
    MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int)
Runs a MC simulation of the magnetic polymer by indexing the monte carlo
steps from 'start' to 'finish'. Modifications are done on the pol::Magnetic_polymer
and the data is stored in the traj::MC_data at the positions from 'start' to 'finish'
"""
function MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int)
    MC_run_base!(pol, traj, start, finish)
end

"""
    MC_run!(pol::Magnetic_polymer, traj::MC_data)
Same thins as MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int)
but start and finish are defaulted to 1 and n_steps
"""
function MC_run!(pol::Magnetic_polymer, traj::MC_data)
    MC_run_base!(pol, traj, 1, traj.n_steps)
end

function MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int, sample_lag::Int, n_samples::Int)
    MC_run_base!(pol, traj, start, min(finish, sample_lag*n_samples))
end

"""
    MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int, spins_configs::Matrix{Int}, sample_lag::Int)
Still performs the monte carlo simulation but one every 'sample_lag' steps the spins configuration is stored.
The configurations are stored in spins_configs::Matrix{Int}
"""
function MC_run!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int, spins_configs::Matrix{Int}, sample_lag::Int, n_samples::Int)
    MC_run_save!(pol, traj, start, finish, spins_configs, sample_lag, n_samples)
end

"""
    MC_run!(pol::Magnetic_polymer, traj::MC_data, spins_configs::Matrix{Int}, sample_lag::Int)
Monte Carlo that saves configurations as above but defaults starting and ending steps to 1 and n_steps
"""
function MC_run!(pol::Magnetic_polymer, traj::MC_data, spins_configs::Matrix{Int}, sample_lag::Int, n_samples::Int)
    MC_run_save!(pol, traj, 1, traj.n_steps, spins_configs, sample_lag, n_samples)
end

#############################################################################################################
#############################################################################################################

"""
    MMC_run!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
    n_strides::Int, stride::Int, inv_temps::Array{Float64})
Performs simulation with Multiple Markov Chains (MMC) attempting a swap between Chains
at two contiguous temperatures every 'stride' steps. The number of attempet
swaps is 'n_strides' and the values of the inverse temperatures of the
different Markov chains are in the array 'inv_temps'
"""
function MMC_run!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
    n_strides::Int, stride::Int, inv_temps::Array{Float64})
    MMC_run_base!(polymers, trajs, n_strides ,stride, inv_temps)
end

"""
    MMC_run!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
    n_strides::Int, stride::Int, inv_temps::Array{Float64}, spins_configs::Matrix{Int}, sample_lag::Int)
Same thing as above but the spins configurations are saved once every 'sample_lag' steps
"""
function MMC_run!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
    n_strides::Int, stride::Int, inv_temps::Array{Float64}, spins_configs::Matrix{Int}, sample_lag::Int,n_samples::Int)
    MMC_run_save!(polymers, trajs, n_strides ,stride, inv_temps,spins_configs,sample_lag,n_samples)
end

#############################################################################################################
#############################################################################################################

function write_results(pol::Magnetic_polymer, traj::MC_data)
    !isdir("simulation_data") && mkdir("simulation_data")

    open("simulation_data/final_config_$(pol.n_mono)_$(pol.inv_temp).txt", "w") do io
        writedlm(io, [transpose(pol.coord) pol.spins])
    end
    open("simulation_data/MC_data_$(pol.n_mono)_$(pol.inv_temp).txt", "w") do io
        writedlm(io, [traj.energies traj.magnetization traj.rg2])
    end
end

function write_results(pol::Magnetic_polymer)
    !isdir("simulation_data") && mkdir("simulation_data")

    open("simulation_data/final_config_$(pol.n_mono)_$(pol.inv_temp).txt", "w") do io
        writedlm(io, [transpose(pol.coord) pol.spins])
    end
end

#############################################################################################################
#############################################################################################################

function simulation(n_mono::Int, n_steps::Int ,beta_temp::Float64, spins_coupling::Float64, alpha::Float64)
    simulation_data = MC_data(n_steps)
    polymer = Magnetic_polymer(n_mono, beta_temp, spins_coupling, alpha)
    initialize_poly!(polymer)
    MC_run!(polymer, simulation_data)
    write_results(polymer, simulation_data)
end

function simulation(n_mono::Int, n_steps::Int ,beta_temp::Float64, spins_coupling::Float64, alpha::Float64, conf_file::String)
    simulation_data = MC_data(n_steps)
    polymer = Magnetic_polymer(n_mono, beta_temp, spins_coupling, alpha)
    initialize_poly!(polymer, conf_file)
    MC_run!(polymer, simulation_data)
    write_results(polymer, simulation_data)
end

function simulation(input::String)
    # First read parameters from input
    io = open("input.txt", "r")
    n_mono = 1
    n_steps = 1
    beta_temp = 1.0
    spins_coupling = 0.0
    alpha = 1.0
    conf_file = ""
    from_file = false
    while ! eof(io)
        line = split(readline(io))
        if line[1] == "n_mono"
            n_mono = parse(Int, line[2])
        elseif line[1] == "n_steps"
            n_steps = parse(Int, line[2])
        elseif line[1] == "spins_coupling"
            spins_coupling = parse(Float64, line[2])
        elseif line[1] == "alpha"
            alpha = parse(Float64, line[2])
        elseif line[1] == "inverse_temperature"
            beta_temp = parse(Float64, line[2])
        elseif line[1] == "initialize_from_files"
            from_file = parse(Bool, line[2])
        elseif line[1] == "conf_file"
            conf_file = String(line[2])
        end
    end

    close(io)

    simulation_data = MC_data(n_steps)
    polymer = Magnetic_polymer(n_mono, beta_temp, spins_coupling, alpha)

    if from_file
        initialize_poly!(polymer, conf_file)
    else
        initialize_poly!(polymer)
    end

    MC_run!(polymer, simulation_data)
    write_results(polymer, simulation_data)
end

#############################################################################################################
#############################################################################################################

function MMC_simulation(n_mono::Int, n_strides::Int, stride::Int)
    spins_coupling = 1.0
    alpha = 1.0

    #inv_temps = [1.1, 1.05, 1.0, 0.95, 0.9, 0.8, 0.74, 0.69, 0.67, 0.64, 0.63, 0.625, 0.617,
    #0.606, 0.595, 0.588, 0.581, 0.571, 0.555, 0.541, 0.531, 0.518, 0.5, 0.4]

    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.57, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)
    polymers = Array{Magnetic_polymer}(undef, n_temps)
    trajs = Array{MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = MC_data(n_strides*stride)
        initialize_poly!(polymers[i_temp])#,"simulation_data/final_config_$(n_mono)_$(inv_temps[i_temp]).txt")
    end

    #=for i_strides in 1:n_strides
        println("Istride = ", i_strides)
        for i_temp in 1:n_temps
            MC_run!(polymers[i_temp], trajs[i_temp],(i_strides-1)*stride+1,i_strides*stride)
        end
    end=#

    MMC_run!(polymers,trajs,n_strides,stride,inv_temps)

    for i_temp in 1:n_temps
        write_results(polymers[i_temp],trajs[i_temp])
    end
end

end # end of module

ext: jl


text:
using DifferentialEquations, Plots

#Solving the simple pendulum with a traditional ODE method
#==========================================================#
function pendulum(du, u, params, t)
#reference: http://www.pgccphy.net/ref/advmech.pdf page 6

    g = params[1] #gravitational acceleration
    m = params[2] #mass
    l = params[3] #length

    θ = u[1]
    ℒ = u[2]

    d_θ = ℒ/(m*l^2)
    d_ℒ = -m*g*l*sin(θ)

    du .= [d_θ, d_ℒ]
    return nothing
end

g = 9.81
m = 2.
l = 1.
params = [g,m,l]
u0 = [1.0,1.0]

prob = ODEProblem(pendulum, u0, (0., 100.), params)
sol1 = solve(prob, AutoVern7(Rodas5()), dt = .05)

#==========================================================#

#Solving the simple pendulum with the DiffEqPhysics.jl HamiltonianProblem()
#==========================================================#
function H(θ, ℒ, params)
    g = params[1] #gravitational acceleration
    m = params[2] #mass
    l = params[3] #length

    return ℒ^2/(2*m*l^2) + m*g*l*(1-cos(θ))
end

g = 9.81
m = 2.
l = 1.
params = [g,m,l]
θ₀ = 1.
ℒ₀ = 1.

prob = HamiltonianProblem(H, θ₀, ℒ₀, (0., 100.), p=params)
sol2 = solve(prob, SofSpa10(), dt = .05);
#==========================================================#

#Plotting
#==========================================================#
plot(sol1, vars=1, tspan=(0,20), label="ODE Method")
plot!(sol2.t, sol2[1,:], xlim=(0,20), label="Hamiltonian Method")
#They produce the same solution!


text:
include("../../../linalg/src/linalg_module.jl")
include("../../../../std/util/src/util_module.jl")

module spline
using ..linalg
using ..util

struct NaturalSpline{T}
    x::Vector{T}
    y::Vector{T}
    c::Vector{T}
    function NaturalSpline{T}(x_::Vector{T}, y_::Vector{T}, c_::Vector{T}) where {T}
        return new( copy(x_), copy(y_), copy(c_) )
    end
end
function NaturalSpline(x::Vector{T}, y::Vector{T}, vec::Vector{T}) where {T}
    return NaturalSpline{T}(x, y, vec)
end

function Base.length(natural_spline::NaturalSpline)
    return length(natural_spline.x)
end

function fit_natural(x::Vector{T}, y::Vector{T}) where {T}
    n_points        = length(x)
    dx              = diff(x)
    dy              = diff(y)
    v1              = view(dx, 2:(n_points-1))
    v2              = view(dx, 1:(n_points-2))
    left_diag       = view(dx, 2:(n_points-2))
    main_diag       = 2 .* (v1 .+ v2)
    dydx            = dy ./ dx
    diff_dydx       = diff(dydx)
    rhs             = 6 * diff_dydx
    c               = linalg.linsolve_tridiagonal(left_diag, main_diag, left_diag, rhs)
    natural_spline  = NaturalSpline(x, y, [0.0; c; 0.0])
    return natural_spline
end

# function interpolate(natural_spline::NaturalSpline, xq::T; init_index::Int = 1) where {T <:Number}
function interpolate(spline_obj, xq::T; init_index::Int = 1) where {T <:Number}
    (x, a0, a1, a2, a3) = __find_interp_values(spline_obj, xq, init_index)
    return __poly_eval(x, a0, a1, a2, a3)
end

function __find_interp_values(spline_obj, xq, init_index)

    # search for the coeffs
    # use binary search for it
    left_idx = util.binary_search(
        spline_obj.x,
        xq;
        left_init_index = init_index)

    n = length(spline_obj)

    if left_idx >= n
        return spline_obj.y[n]
    elseif left_idx < 1
        return spline_obj.y[1]
    end

    x1 = spline_obj.x[left_idx]
    x2 = spline_obj.x[left_idx+1]

    y1 = spline_obj.y[left_idx]
    y2 = spline_obj.y[left_idx+1]

    c1 = spline_obj.c[left_idx]
    c2 = spline_obj.c[left_idx+1]

    dx = x2 - x1
    dy = y2 - y1

    # def of the coeffs
    a0 = y1
    a1 = dy / dx - (c2 + 2 * c1) / 6 * dx
    a2 = c1 / 2
    a3 = (c2 - c1) / (6 * dx)

    # eval of the poly
    x = xq - x1

    return (x, a0, a1, a2, a3)
end

function __poly_eval(x, a0, a1, a2, a3)
    return ((a3 .* x .+ a2) .* x + a1) .* x .+ a0
end

function __diff_poly_eval(x, a0, a1, a2, a3)
    return (3 .* a3 .* x .+ 2 .* a2) .* x + a1
end

function __poly_eval(x, a0, a1, a2)
    return (a2 .* x .+ a1) .* x + a0
end

# function interpolate(natural_spline::NaturalSpline, xq::Vector{T}) where {T <:Number}
function interpolate(spline_obj, xq::Vector{T}) where {T <:Number}
    # todo: rewrite it in a more better way
    n  = length(xq)
    yq = zeros(n)
    for ii = 1:n
        yq[ii] = interpolate(spline_obj, xq[ii])
    end
    return yq
end

function diff_interpolate(spline_obj, xq::T; init_index::Int = 1) where {T <:Number}
    (x, a0, a1, a2, a3) = __find_interp_values(spline_obj, xq, init_index)
    return __diff_poly_eval(x, a0, a1, a2, a3)
end

function diff_interpolate(spline_obj, xq::Vector{T}) where {T <:Number}
    # todo: rewrite it in a more better way
    n  = length(xq)
    yq = zeros(n)
    for ii = 1:n
        yq[ii] = diff_interpolate(spline_obj, xq[ii])
    end
    return yq
end

struct ClampedSpline{T}
    x::Vector{T}
    y::Vector{T}
    c::Vector{T}
    function ClampedSpline{T}(x_, y_, c_) where {T}
        return new( copy(x_), copy(y_), copy(c_) )
    end
end
function ClampedSpline(x::Vector{T}, y::Vector{T}, vec::Vector{T}) where {T}
    return ClampedSpline{T}(x, y, vec)
end

function Base.length(clamped_spline::ClampedSpline)
    return length(clamped_spline.x)
end

function fit_clamped(x::Vector{T}, y::Vector{T}, g1, g2) where {T}
    n_points        = length(x)
    dx              = diff(x)
    dy              = diff(y)
    v1              = view(dx, 2:(n_points-1))
    v2              = view(dx, 1:(n_points-2))
    left_diag       = dx
    main_diag       = [2 * dx[1]; 2 .* (v1 .+ v2); 2 * dx[end]]
    dydx            = dy ./ dx
    diff_dydx       = diff(dydx)
    rhs             = 6 .* [-g1 + dydx[1]; diff_dydx; g2 - dydx[end]]
    c               = linalg.linsolve_tridiagonal(left_diag, main_diag, left_diag, rhs)
    natural_spline  = ClampedSpline(x, y, c)
    return natural_spline
end

mutable struct ModifiedAkimaSpline{T}
    x::Vector{T}
    y::Vector{T}
    function ModifiedAkimaSpline{T}(x_, y_) where {T}
        return new( copy(x_), copy(y_) )
    end
end
function ModifiedAkimaSpline(x::Vector{T}, y::Vector{T}) where {T}
    return ModifiedAkimaSpline{T}(x, y)
end

function Base.length(akima_spline::ModifiedAkimaSpline)
    return length(akima_spline.x)
end

function fit_modified_akima(x, y)
    return ModifiedAkimaSpline(x, y)
end

function interpolate(
    spline_obj::ModifiedAkimaSpline,
    xq::T;
    init_index::Int = 1
    ) where {T <:Number}

    (a0, a1, a2, a3, x3) = __mod_akima_interp(
        spline_obj,
        xq;
        init_index = init_index)

    return __poly_eval(xq - x3, a0, a1, a2, a3)
end

function diff_interpolate(
    spline_obj::ModifiedAkimaSpline,
    xq::T;
    init_index::Int = 1
    ) where {T <: Number}

    (a0, a1, a2, a3, x3) = __mod_akima_interp(
        spline_obj,
        xq;
        init_index = init_index)

    return __diff_poly_eval(xq - x3, a0, a1, a2, a3)
end

function __mod_akima_interp(
    spline_obj::ModifiedAkimaSpline,
    xq::T;
    init_index::Int = 1
    ) where {T <:Number}

    left_idx = util.binary_search(
        spline_obj.x,
        xq;
        left_init_index = init_index)

    n = length(spline_obj)

    if left_idx >= n
        return spline_obj.y[n]
    elseif left_idx < 1
        return spline_obj.y[1]
    end

    # see https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
    # first two segments
    # quadratic interpolation
    t1 = 0.0
    t2 = 0.0
    g3 = 0.0
    x3 = 0.0
    x4 = 0.0
    y3 = 0.0

    ii = left_idx
    if left_idx == 1

        x3 = spline_obj.x[ii]
        x4 = spline_obj.x[ii+1]
        x5 = spline_obj.x[ii+2]
        x6 = spline_obj.x[ii+3]

        y3 = spline_obj.y[ii]
        y4 = spline_obj.y[ii+1]
        y5 = spline_obj.y[ii+2]
        y6 = spline_obj.y[ii+3]

        g3 = (y4 - y3) / (x4 - x3)
        g4 = (y5 - y4) / (x5 - x4)
        g5 = (y6 - y5) / (x6 - x5)

        g2 = 2 * g3 - g4
        g1 = 2 * g2 - g3

    elseif left_idx == 2

        x2 = spline_obj.x[ii-1]
        x3 = spline_obj.x[ii]
        x4 = spline_obj.x[ii+1]
        x5 = spline_obj.x[ii+2]
        x6 = spline_obj.x[ii+3]

        y2 = spline_obj.y[ii-1]
        y3 = spline_obj.y[ii]
        y4 = spline_obj.y[ii+1]
        y5 = spline_obj.y[ii+2]
        y6 = spline_obj.y[ii+3]

        g2 = (y3 - y2) / (x3 - x2)
        g3 = (y4 - y3) / (x4 - x3)
        g4 = (y5 - y4) / (x5 - x4)
        g5 = (y6 - y5) / (x6 - x5)

        g1 = 2 * g2 - g3

    elseif left_idx == n-1

        x1 = spline_obj.x[ii-2]
        x2 = spline_obj.x[ii-1]
        x3 = spline_obj.x[ii]
        x4 = spline_obj.x[ii+1]

        y1 = spline_obj.y[ii-2]
        y2 = spline_obj.y[ii-1]
        y3 = spline_obj.y[ii]
        y4 = spline_obj.y[ii+1]

        g1 = (y2 - y1) / (x2 - x1)
        g2 = (y3 - y2) / (x3 - x2)
        g3 = (y4 - y3) / (x4 - x3)

        g4 = 2 * g3 - g2
        g5 = 2 * g4 - g3

    elseif left_idx == n-2

        x1 = spline_obj.x[ii-2]
        x2 = spline_obj.x[ii-1]
        x3 = spline_obj.x[ii]
        x4 = spline_obj.x[ii+1]
        x5 = spline_obj.x[ii+2]

        y1 = spline_obj.y[ii-2]
        y2 = spline_obj.y[ii-1]
        y3 = spline_obj.y[ii]
        y4 = spline_obj.y[ii+1]
        y5 = spline_obj.y[ii+2]

        g1 = (y2 - y1) / (x2 - x1)
        g2 = (y3 - y2) / (x3 - x2)
        g3 = (y4 - y3) / (x4 - x3)
        g4 = (y5 - y4) / (x5 - x4)

        g5 = 2 * g4 - g3

    else
        # bulk
        x1 = spline_obj.x[ii-2]
        x2 = spline_obj.x[ii-1]
        x3 = spline_obj.x[ii]
        x4 = spline_obj.x[ii+1]
        x5 = spline_obj.x[ii+2]
        x6 = spline_obj.x[ii+3]

        y1 = spline_obj.y[ii-2]
        y2 = spline_obj.y[ii-1]
        y3 = spline_obj.y[ii]
        y4 = spline_obj.y[ii+1]
        y5 = spline_obj.y[ii+2]
        y6 = spline_obj.y[ii+3]

        g1 = (y2 - y1) / (x2 - x1)
        g2 = (y3 - y2) / (x3 - x2)
        g3 = (y4 - y3) / (x4 - x3)
        g4 = (y5 - y4) / (x5 - x4)
        g5 = (y6 - y5) / (x6 - x5)

    end

    w1 = abs(g4 - g3) + 0.5 * abs(g4 + g3)
    w2 = abs(g2 - g1) + 0.5 * abs(g2 + g1)
    t1 = w1 / (w1 + w2) * g2 + w2 / (w1 + w2) * g3

    w1 = abs(g5 - g4) + 0.5 * abs(g5 + g4)
    w2 = abs(g3 - g2) + 0.5 * abs(g3 + g2)
    t2 = w1 / (w1 + w2) * g3 + w2 / (w1 + w2) * g4

    if isnan(t1)
        t1 = 0
    end
    if isnan(t2)
        t2 = 0
    end

    a0 = y3
    a1 = t1
    a2 = (3 * g3 - 2 * t1 - t2) / (x4 - x3)
    a3 = (t1 + t2 - 2 * g3) / (x4 - x3)^2

    return (a0, a1, a2, a3, x3)
end

struct AkimaSpline{T}
    x::Vector{T}
    y::Vector{T}
end
function AkimaSpline(x::Vector{T}, y::Vector{T}) where {T}
    return AkimaSpline{T}(x, y)
end

function Base.length(akima_spline::AkimaSpline)
    return length(akima_spline.x)
end

function fit_akima(x, y)
    return AkimaSpline(x, y)
end

function interpolate(spline_obj::AkimaSpline, xq::T; init_index::Int = 1) where {T <:Number}
    left_idx = util.binary_search(
        spline_obj.x,
        xq;
        left_init_index = init_index)

    n = length(spline_obj)

    if left_idx >= n
        return spline_obj.y[n]
    elseif left_idx < 1
        return spline_obj.y[1]
    end

    # first two segments
    # quadratic interpolation
    if left_idx < 3
        return spline_obj.y[1]
    elseif left_idx > n-3
        return spline_obj.y[n-3]
    end

    # bulk
    ii = left_idx
    x1 = spline_obj.x[ii-2]
    x2 = spline_obj.x[ii-1]
    x3 = spline_obj.x[ii]
    x4 = spline_obj.x[ii+1]
    x5 = spline_obj.x[ii+2]
    x6 = spline_obj.x[ii+3]

    y1 = spline_obj.y[ii-2]
    y2 = spline_obj.y[ii-1]
    y3 = spline_obj.y[ii]
    y4 = spline_obj.y[ii+1]
    y5 = spline_obj.y[ii+2]
    y6 = spline_obj.y[ii+3]

    g1 = (y2 - y1) / (x2 - x1)
    g2 = (y3 - y2) / (x3 - x2)
    g3 = (y4 - y3) / (x4 - x3)
    g4 = (y5 - y4) / (x5 - x4)
    g5 = (y6 - y5) / (x6 - x5)

    t1 = abs(g4 - g3) * g2 + abs(g2 - g1) * g3
    t1 = t1 / ( abs(g4 - g3) + abs(g2 - g1) )

    t2 = abs(g5 - g4) * g3 + abs(g3 - g2) * g4
    t2 = t2 / ( abs(g5 - g4) + abs(g3 - g2))

    if isapprox(g1, g2)
        t1 = g2
    end
    if isapprox(g3, g4)
        t1 = g3
    end
    if isapprox(g1, g2) && isapprox(g3, g4) && !isapprox(g1, g3)
        t1 = 0.5 * (g1 + g3)
    end

    if isapprox(g2, g3)
        t2 = g3
    end
    if isapprox(g4, g5)
        t2 = g4
    end
    if isapprox(g2, g3) && isapprox(g4, g5) && !isapprox(g2, g4)
        t2 = 0.5 * (g2 + g4)
    end

    a0 = y3
    a1 = t1
    a2 = (3 * g3 - 2 * t1 - t2) / (x4 - x3)
    a3 = (t1 + t2 - 2 * g3) / (x4 - x3)^2

    return __poly_eval(xq - x3, a0, a1, a2, a3)
end

#=
struct HermiteSpline{T}
    # lazy version
    x::Vector{T}
    y::Vector{T}
    g::Vector{T}
end
function HermiteSpline(x::Vector{T}, y::Vector{T}, g1, g2) where {T}
    g = [g1; diff(y) ./ diff(x); g2]
    return HermiteSpline{T}(x, y, g)
end
function Base.length(hermite_spline::HermiteSpline)
    return length(hermite_spline.x)
end

function fit_hermite(
    x::Vector{T},
    y::Vector{T},
    g1 = (y[2] - y[1]) / (x[2] - x[1]),
    g2 = (y[end] - y[end-1]) / (x[end] - x[end-1])
    ) where {T}

    return HermiteSpline(x, y, g1, g2)

end

function interpolate(spline_obj::HermiteSpline, xq::T; init_index::Int = 1) where {T <:Number}
    left_idx = util.binary_search(
        spline_obj.x,
        xq;
        left_init_index = init_index)

    n = length(spline_obj)

    if left_idx >= n
        return spline_obj.y[n]
    elseif left_idx < 1
        return spline_obj.y[1]
    end

    x1 = spline_obj.x[left_idx]
    x2 = spline_obj.x[left_idx+1]
    y1 = spline_obj.y[left_idx]
    y2 = spline_obj.y[left_idx+1]

    g1 = spline_obj.g[left_idx]
    g2 = spline_obj.g[left_idx+2]

    a0 = y1
    a1 = g1

    dx_ = x2 - x1
    dy_ = y2 - y1

    a3 = ((g2 - g1) / dx_ - (dy_ - g1 * dx_) / dx_^2) / dx_
    a2 = 0.5 * ((g2 - g1) / dx_ - 3 * a3 * dx_)

    return __poly_eval(xq - x1, a0, a1, a2, a3)
end
=#

#=
struct HermiteSpline{T}
    # eager version
    x::Vector{T}
    y::Vector{T}
    a0::Vector{T}
    a1::Vector{T}
    a2::Vector{T}
    a3::Vector{T}
end
function HermiteSpline(x::Vector{T}, y::Vector{T}, a0, a1, a2, a3) where {T}
    return HermiteSpline{T}(x, y, a0, a1, a2, a3)
end

function Base.length(hermite_spline::HermiteSpline)
    return length(hermite_spline.x)
end

function fit_hermite(
    x::Vector{T},
    y::Vector{T},
    g1 = (y[2] - y[1]) / (x[2] - x[1]),
    g2 = (y[end] - y[end-1]) / (x[end] - x[end-1])
    ) where {T}

    L = length(x)
    dx = diff(x) # L-1
    dy = diff(y) # L-1

    # vy1 = view(y, 1:L-2) # L-2
    # vy2 = view(y, 3:L) # L-2
    # vx1 = view(x, 1:L-2) # L-2
    # vx2 = view(x, 3:L) # L-2

    # g = [g1; (vy2 .- vy1) ./ (vx2 .- vx1); g2] # L - 2 + 2 = L

    g      = zeros(L)
    g[1]   = g1
    g[end] = g2
    for kk = 2:L-1
        m   = x[kk]   - x[kk-1]
        p   = x[kk+1] - x[kk]
        dym = y[kk]   - y[kk-1]
        dyp = y[kk+1] - y[kk]

        g[kk] = (dyp * m^2 + dym * p^2) / ( m*p * (m + p) )
    end

    a0 = copy(y)
    pop!(a0) # L-1

    a1 = g[1:end-1] # L - 1
    a2 = zeros(L-1)
    a3 = zeros(L-1)
    for kk = 1:L-1
        g1  = g[kk]
        g2  = g[kk+1]
        dx_ = dx[kk]
        dy_ = dy[kk]

        a2[kk] = - (2 * g2 + g1) / dx_ + 3 * dy_ / dx_^2
        a3[kk] = (g1 + g2) / dx_^2 - 2 * dy_ / dx_^3
    end

    return HermiteSpline(x, y, a0, a1, a2, a3)
end

function interpolate(spline_obj::HermiteSpline, xq::T; init_index::Int = 1) where {T <:Number}
    left_idx = util.binary_search(
        spline_obj.x,
        xq;
        left_init_index = init_index)

    n = length(spline_obj)

    if left_idx >= n
        return spline_obj.y[n]
    elseif left_idx < 1
        return spline_obj.y[1]
    end

    a0 = spline_obj.a0[left_idx]
    a1 = spline_obj.a1[left_idx]
    a2 = 0 * spline_obj.a2[left_idx]
    a3 = 0 * spline_obj.a3[left_idx]

    return __poly_eval(xq - spline_obj.x[left_idx], a0, a1, a2, a3)
end
=#

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Random
using Distributions
#using FunctionWrappers

const Age=UInt8
abstract type AbstractSimParams end

abstract type InfectionModulation end

abstract type AbstractOutsideCases end

include("params/age_coupling.jl")
include("params/households.jl")
include("params/friendship.jl")
include("params/modulations.jl")
include("params/progression.jl")
include("params/hospital.jl")
include("params/phonetracing.jl")
include("params/spreading.jl")
include("params/strains.jl")
include("params/outside_cases.jl")
include("params/screening.jl")

struct SimParams <: AbstractSimParams
  household_ptrs::Vector{Tuple{PersonIdx,PersonIdx}}  # (i1,i2) where i1 and i2 are the indices of first and last member of the household

  ages::Vector{Age}
  genders::BitVector

  progressions::Vector{Progression} # not sure if progressions should be there

  strain_table::StrainTable

  constant_kernel_param::Float32
  household_kernel_param::Float32
  age_coupling_params::Union{Nothing, AgeCouplingParams} # nothing if kernel not active
  hospital_kernel_params::Union{Nothing, HospitalInfectionParams}  # nothing if hospital kernel not active

  hospital_detections::Bool
  mild_detection_prob::Float64

  backward_tracing_prob::Float32
  backward_detection_delay::TimeDiff

  forward_tracing_prob::Float32
  forward_detection_delay::TimeDiff

  quarantine_length::Float32
  testing_time::TimeDiff

  phone_tracing_params::Union{Nothing, PhoneTracingParams}

  infection_modulation_function::Union{Nothing,InfectionModulation}
  screening_params::Union{Nothing, ScreeningParams}
  spreading_params::Union{Nothing, SpreadingParams}
end

numindividuals(params::SimParams) = length(params.household_ptrs)
straindata(params::SimParams, strain::StrainKind) = getdata(params.strain_table, strain)

progressionof(params::SimParams, person_id::Integer) = params.progressions[person_id]
severityof(params::SimParams, person_id::Integer) = progressionof(params, person_id).severity
householdof(params::SimParams, person_id::Integer) = UnitRange(params.household_ptrs[person_id]...)
age(params::SimParams, person_id::Integer) = params.ages[person_id]
gender(params::SimParams, person_id::Integer) = params.genders[person_id]

socialcompetence(params::SimParams, person_id::Integer) =
  nothing!=params.friendship_kernel_params && socialcompetence(params.friendship_kernel_params, person_id)

ishealthcare(params::SimParams, person_id::Integer) =
  nothing!=params.hospital_kernel_params && ishealthcare(params.hospital_kernel_params, person_id)
uses_phone_tracing(params::SimParams, person_id::Integer) =
  nothing!=params.phone_tracing_params && uses_phone_tracing(params.phone_tracing_params, person_id)

spreading(params::SimParams, person_id::Integer) = isnothing(params.spreading_params) ? 1.0 : spreading(params.spreading_params, person_id)

function load_params(rng=MersenneTwister(0);
        population::DataFrame,
        infection_modulation_name::Union{Nothing,AbstractString}=nothing,
        infection_modulation_params::NamedTuple=NamedTuple{}(),
        kwargs...
        )

  individuals_df::DataFrame = population

  num_individuals = individuals_df |> nrow

  dist_incubation_time = LogNormal(1.3669786931887833, 0.5045104580676582)
  dist_symptom_onset_time = Gamma(0.8738003969079596, 2.9148873266517685)
  dist_hospitalization_time = Gamma(1.1765988120148885, 2.6664347368236787)
  dist_mild_recovery_time = Uniform(11, 17)
  dist_severe_recovery_time = Uniform(4*7, 8*7)
  dist_death_time = LogNormal(2.610727719719777, 0.44476420066780653)

  progressions = Vector{Progression}(undef, num_individuals);
  resample!(rng, progressions, individuals_df.age,
    dist_incubation_time,
    dist_symptom_onset_time,
    dist_hospitalization_time,
    dist_mild_recovery_time,
    dist_severe_recovery_time,
    dist_death_time
  )

  infection_modulation_function = isnothing(infection_modulation_name) ? nothing : make_infection_modulation(infection_modulation_name; infection_modulation_params...)

  make_params(
    rng;
    individuals_df=individuals_df,
    progressions=progressions,
    infection_modulation_function=infection_modulation_function,
    kwargs...
  )
end

function make_params(
  rng::AbstractRNG=MersenneTwister(0);
  individuals_df::DataFrame,
  progressions::AbstractArray{Progression},

  infection_modulation_function=nothing,

  screening_params::Union{Nothing,ScreeningParams}=nothing,

  hospital_kernel_param::Float64=0.0,
  healthcare_detection_prob::Float64=0.8,
  healthcare_detection_delay::Float64=1.0,

  constant_kernel_param::Float64=1.0,
  household_kernel_param::Float64=1.0,

  hospital_detections::Bool=true,
  mild_detection_prob::Float64=0.0,

  backward_tracing_prob::Float64=0.0,
  backward_detection_delay::Float64=1.0,

  forward_tracing_prob::Float64=0.0,
  forward_detection_delay::Float64=1.0,

  quarantine_length::Float64=14.0,
  testing_time::Float64=1.0,

  phone_tracing_usage::Real=0.0,
  phone_detection_delay::Real=0.25,
  phone_tracing_usage_by_household::Bool=false,

  spreading_alpha::Union{Nothing,Real}=nothing,
  spreading_x0::Real=1,
  spreading_truncation::Real=Inf,

  british_strain_multiplier::Real=1.70,
  delta_strain_multiplier::Real=1.7*1.5,

  age_coupling_thresholds::Union{Nothing, AbstractArray{T} where T<:Real}=nothing,
  age_coupling_weights::Union{Nothing, AbstractMatrix{T} where T<:Real}=nothing,
  age_coupling_use_genders::Bool=false,
  age_coupling_param::Union{Nothing, Real}=nothing,
)
  sort!(individuals_df, :household_index)

  num_individuals = individuals_df |> nrow

  @assert num_individuals == length(progressions)

  household_ptrs = make_household_ptrs(individuals_df.household_index)

  strain_table = make_strains(british_multiplier=british_strain_multiplier, delta_multiplier=delta_strain_multiplier)

  age_coupling_kernel_params =
    if nothing === age_coupling_weights && nothing === age_coupling_thresholds && nothing === age_coupling_param; nothing
    elseif age_coupling_weights !== nothing && age_coupling_thresholds !== nothing
      @assert minimum(individuals_df.age) >= 0
      @assert maximum(individuals_df.age) * (age_coupling_use_genders+1) < typemax(GroupIdx)
      AgeCouplingParams(
        individuals_df.age,
        age_coupling_use_genders === nothing ? individuals_df.gender : nothing,
        age_coupling_thresholds, age_coupling_weights,
        age_coupling_param
        )
    else error("age couplig params not fully given")
    end

  hospital_kernel_params =
    if 0 == hospital_kernel_param; nothing
    elseif 0.0 < hospital_kernel_param;
      HospitalInfectionParams(
        individuals_df.ishealthcare,
        (UInt32(1):UInt32(num_individuals))[individuals_df.ishealthcare],
        hospital_kernel_param,
        healthcare_detection_prob,
        healthcare_detection_delay
      )
    else error("hospital_kernel_param must be postive or 0, got $hospital_kernel_param")
    end

  phone_tracing_params =
    if 0 == phone_tracing_usage; nothing
    elseif 0.0 < phone_tracing_usage <= 1.0
      phone_tracing_usage_by_household ?
        PhoneTracingParams(rng, num_individuals, phone_tracing_usage, phone_detection_delay, 1, household_ptrs) :
        PhoneTracingParams(rng, num_individuals, phone_tracing_usage, phone_detection_delay, 1)
    else error("tracing_app_usage must be nonnegative, got $phone_tracing_usage")
    end

  spreading_params =
    if nothing === spreading_alpha; nothing
    elseif 0.0 < spreading_alpha
      SpreadingParams(rng, num_individuals, alpha=spreading_alpha, x0=spreading_x0, truncation=spreading_truncation)
    else error("spreading_alpha must be larger than 0, got $spreading_alpha")
    end

  params = SimParams(
    household_ptrs,
    individuals_df.age,
    individuals_df.gender,

    progressions,

    strain_table,

    constant_kernel_param,
    household_kernel_param,
    age_coupling_kernel_params,
    hospital_kernel_params,

    hospital_detections,
    mild_detection_prob,

    backward_tracing_prob,
    backward_detection_delay,

    forward_tracing_prob,
    forward_detection_delay,

    quarantine_length, # quarantine length
    testing_time, # testing time

    phone_tracing_params,
    infection_modulation_function,
    screening_params,
    spreading_params
  )
  params
end

function saveparams(dict, p::SimParams)
  dict["num_individuals"] = numindividuals(p)
  dict["constant/kernel_param"] = p.constant_kernel_param
  dict["household/kernel_param"] = p.household_kernel_param
  dict["quarantine/duration"] = p.quarantine_length
  dict["detections/hospital"] = p.hospital_detections
  dict["detections/mild"] = p.mild_detection_prob

  dict["infection_modulation"] = p.infection_modulation_function

  dict["tracing/backward_prob"] = p.backward_tracing_prob
  dict["tracing/backward_detection_delay"] = p.backward_detection_delay
  dict["tracing/forward_tracing_prob"] = p.forward_tracing_prob
  dict["tracing/forward_detection_delay"] = p.forward_detection_delay
  dict["tracing/testing_time"] = p.testing_time

  saveparams(dict, p.progressions, "progressions/")
  nothing===p.hospital_kernel_params || saveparams(dict, p.hospital_kernel_params, "hospital/")
  nothing===p.phone_tracing_params || saveparams(dict, p.phone_tracing_params, "phone_tracing/")
  nothing===p.spreading_params || saveparams(dict, p.spreading_params, "spreading")
  nothing
end


text:
"""
This submodule contains ready-to-use Scatterer models of zooplankton.
"""
module Models

using SDWBA
using DelimitedFiles

export  krill_mcgeehee,
		krill_conti,
		calanoid_copepod,
		sandeel,
		daphnia,
		nauplius

DATA_DIR = joinpath(dirname(@__FILE__), "data/")

"""
Krill shape from McGeehee et al. 1998, 'Effects of orientation on acoustic scattering
from Antarctic krill at 120 kHz.' Deep-Sea Research II 45(7), 1273-1294.
"""
krill_mcgeehee = from_csv(joinpath(DATA_DIR, "generic_krill_McGeehee1998.csv"))

"""
Krill shape from Conti and Demer, 2006, 'Improved parameterization of the SDWBA for
estimating krill target strength.' ICES Journal of Marine Science 63(5), 928-935.
"""
krill_conti = from_csv(joinpath(DATA_DIR, "generic_krill_Conti2006.csv"))

"""
Generic sand eel/sand lance shape (Ammodytes spp.).  Shape is taken from illustration in
Bigelow and Schroeder, 'Fishes of the Gulf of Maine.'  Material properties from
Yasuma et al. (2009), 'Density and sound-speed contrasts, and target strength of
Japanese sandeel Ammodytes personatus.' Fisheries Science 75 (3) 545-552.
"""
sandeel = from_csv(joinpath(DATA_DIR, "sand_eel.csv"))

"""
Generic calanoid copepod shape, 1 mm long, with g and h for "typical" marine zooplankton
(both 1.04, as per Stanton and Chu 2000, 'Review and recommendations for the modelling
of acoustic scattering by fluid-like elongated zooplankton: euphausiids and copepods.'
ICES Journal of Marine Science 57, 793-807).
"""
calanoid_copepod = from_csv(joinpath(DATA_DIR, "generic_acartia.csv"))

"""
Generic Daphnia shape used in Warren et al. 2016, 'Measuring the distribution,
abundance, and biovolume of zooplankton in an oligotrophic freshwater lake
with a 710 kHz scientific echosounder.' Limnology and Oceanography: Methods
"""
daphnia = from_csv(joinpath(DATA_DIR, "daphnia.csv"))

"""
Alternative (more realistic?) daphnia shape
"""
daphnia2 = from_csv(joinpath(DATA_DIR, "daphnia2.csv"))

"""
Generic nauplius larva shape, 0.5 mm long, with g and h for "typical" marine zooplankton
(both 1.04, as per Stanton and Chu 2000, 'Review and recommendations for the modelling
of acoustic scattering by fluid-like elongated zooplankton: euphausiids and copepods.'
ICES Journal of Marine Science 57, 793-807).
"""
nauplius = from_csv(joinpath(DATA_DIR, "nauplius.csv"))

end


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, InteractiveUtils
using LinearMaps: FiveArg, ThreeArg

@testset "block maps" begin
    @testset "hcat" begin
        for elty in (Float32, ComplexF64), n2 = (0, 20)
            A11 = rand(elty, 10, 10)
            A12 = rand(elty, 10, n2)
            v = rand(elty, 10)
            L = @inferred hcat(LinearMap(A11), LinearMap(A12))
            @test occursin("10×$(10+n2) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            @test L isa LinearMaps.BlockMap{elty}
            if elty <: Complex
                @test_throws ErrorException LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,))
            end
            A = [A11 A12]
            x = rand(10+n2)
            @test size(L) == size(A)
            @test Matrix(L) ≈ A
            @test L * x ≈ A * x
            L = @inferred hcat(LinearMap(A11), LinearMap(A12), LinearMap(A11))
            A = [A11 A12 A11]
            @test Matrix(L) ≈ A
            A = [I I I A11 A11 A11 v]
            @test (@which [A11 A11 A11]).module != LinearMaps
            @test (@which [I I I A11 A11 A11]).module != LinearMaps
            @test (@which hcat(I, I, I)).module != LinearMaps
            @test (@which hcat(I, I, I, LinearMap(A11), A11, A11)).module == LinearMaps
            maps = @inferred LinearMaps.promote_to_lmaps(ntuple(i->10, 7), 1, 1, I, I, I, LinearMap(A11), A11, A11, v)
            @inferred LinearMaps.rowcolranges(maps, (7,))
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v)
            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(v)]
            x = rand(elty, 61)
            @test L isa LinearMaps.BlockMap{elty}
            @test L * x ≈ A * x
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v, v, v, v)
            @test occursin("10×64 LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v, v, v, v, v, v, v)
            @test occursin("10×67 LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            A11 = rand(elty, 11, 10)
            A12 = rand(elty, 10, n2)
            @test_throws DimensionMismatch hcat(LinearMap(A11), LinearMap(A12))
        end
    end

    @testset "vcat" begin
        for elty in (Float32, ComplexF64)
            A11 = rand(elty, 10, 10)
            v = rand(elty, 10)
            L = @inferred vcat(LinearMap(A11))
            @test L == [LinearMap(A11);]
            @test Matrix(L) ≈ A11
            A21 = rand(elty, 20, 10)
            L = @inferred vcat(LinearMap(A11), LinearMap(A21))
            @test occursin("30×10 LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            @test L isa LinearMaps.BlockMap{elty}
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            @test (@which [A11; A21]).module != LinearMaps
            A = [A11; A21]
            x = rand(10)
            @test size(L) == size(A)
            @test Matrix(L) ≈ A
            @test L * x ≈ A * x
            A = [I; I; I; A11; A11; A11; v v v v v v v v v v]
            @test (@which [I; I; I; A11; A11; A11; v v v v v v v v v v]).module != LinearMaps
            L = @inferred vcat(I, I, I, LinearMap(A11), LinearMap(A11), LinearMap(A11), reduce(hcat, fill(v, 10)))
            @test L == [I; I; I; LinearMap(A11); LinearMap(A11); LinearMap(A11); reduce(hcat, fill(v, 10))]
            x = rand(elty, 10)
            @test L isa LinearMaps.BlockMap{elty}
            @test L * x ≈ A * x
            A11 = rand(elty, 10, 11)
            A21 = rand(elty, 20, 10)
            @test_throws DimensionMismatch vcat(LinearMap(A11), LinearMap(A21))
        end
    end

    @testset "hvcat" begin
        for elty in (Float32, ComplexF64)
            A11 = rand(elty, 10, 10)
            A12 = rand(elty, 10, 20)
            A21 = rand(elty, 20, 10)
            A22 = rand(elty, 20, 20)
            A = [A11 A12; A21 A22]
            @test (@which [A11 A12; A21 A22]).module != LinearMaps
            @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22))
            L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)]
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            @test @inferred !issymmetric(L)
            @test @inferred !ishermitian(L)
            x = rand(30)
            @test L isa LinearMaps.BlockMap{elty}
            @test size(L) == size(A)
            @test L * x ≈ A * x
            @test Matrix(L) == A
            @test convert(AbstractMatrix, L) == A
            A = [I A12; A21 I]
            @test (@which [I A12; A21 I]).module != LinearMaps
            @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
            L = @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
            @test L isa LinearMaps.BlockMap{elty}
            @test size(L) == (30, 30)
            @test Matrix(L) ≈ A
            @test L * x ≈ A * x
            y = randn(elty, size(L, 1))
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(y), L, x, α, β) ≈ y*β .+ A*x*α
            end
            X = rand(elty, 30, 10)
            Y = randn(elty, size(L, 1), 10)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), L, X, α, β) ≈ Y*β .+ A*X*α
            end
            A = rand(elty, 10,10); LA = LinearMap(A)
            B = rand(elty, 20,30); LB = LinearMap(B)
            @test [LA LA LA; LB] isa LinearMaps.BlockMap{elty}
            @test Matrix([LA LA LA; LB]) ≈ [A A A; B]
            @test [LB; LA LA LA] isa LinearMaps.BlockMap{elty}
            @test Matrix([LB; LA LA LA]) ≈ [B; A A A]
            @test [I; LA LA LA] isa LinearMaps.BlockMap{elty}
            @test Matrix([I; LA LA LA]) ≈ [I; A A A]
            A12 = LinearMap(rand(elty, 10, 21))
            A21 = LinearMap(rand(elty, 20, 10))
            @test_throws DimensionMismatch A = [I A12; A21 I]
            @test_throws DimensionMismatch A = [I A21; A12 I]
            @test_throws DimensionMismatch A = [A12 A12; A21 A21]
            @test_throws DimensionMismatch A = [A12 A21; A12 A21]

            # basic test of "misaligned" blocks
            M = ones(elty, 3, 2) # non-square
            A = LinearMap(M)
            B = [I A; A I]
            C = [I M; M I]
            @test B isa LinearMaps.BlockMap{elty}
            @test Matrix(B) == C
            @test Matrix(transpose(B)) == transpose(C)
            @test Matrix(adjoint(B)) == C'
        end
    end

    @testset "adjoint/transpose" begin
        for elty in (Float32, ComplexF64), transform in (transpose, adjoint)
            A12 = rand(elty, 10, 10)
            A = [I A12; transform(A12) I]
            L = [I LinearMap(A12); transform(LinearMap(A12)) I]
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            if elty <: Complex
                if transform == transpose
                    @test @inferred issymmetric(L)
                else
                    @test @inferred ishermitian(L)
                end
            end
            if elty <: Real
                @test @inferred ishermitian(L)
                @test @inferred issymmetric(L)
            end
            x = rand(elty, 20)
            @test L isa LinearMaps.LinearMap{elty}
            @test size(L) == size(A)
            @test L * x ≈ A * x
            @test Matrix(L) == A
            @test convert(AbstractMatrix, L) == A
            @test sparse(L) == sparse(A)
            Lt = @inferred transform(L)
            @test Lt isa LinearMaps.LinearMap{elty}
            @test Lt * x ≈ transform(A) * x
            @test convert(AbstractMatrix, Lt) == transform(A)
            @test sparse(transform(L)) == transform(A)
            Lt = @inferred transform(LinearMap(L))
            @test Lt * x ≈ transform(A) * x
            @test Matrix(Lt) == Matrix(transform(A))
            A21 = rand(elty, 10, 10)
            A = [I A12; A21 I]
            L = [I LinearMap(A12); LinearMap(A21) I]
            Lt = @inferred transform(L)
            @test Lt isa LinearMaps.LinearMap{elty}
            @test Lt * x ≈ transform(A) * x
            @test Matrix(Lt) ≈ Matrix(transform(LinearMap(L))) ≈ Matrix(transform(A))
            @test Matrix(transform(LinearMap(L))+transform(LinearMap(L))) ≈ 2Matrix(transform(A))
            X = rand(elty, size(L, 1), 10)
            Y = randn(elty, size(L, 2), 10)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), Lt, X, α, β) ≈ Y*β .+ transform(A)*X*α
            end
        end
    end

    @testset "block diagonal maps" begin
        for elty in (Float32, ComplexF64)
            m = 5; n = 6
            M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1)
            M2 = randn(elty, m, n+2); L2 = LinearMap(M2)
            M3 = randn(elty, m, n+3); L3 = LinearMap(M3)

            # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse:
            if elty <: Complex
                @test_throws ErrorException LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1))
            end
            Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
            @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps
            @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps
            x = randn(elty, size(Md, 2))
            Bd = @inferred blockdiag(L1, L2, L3, L2, L1)
            @test @inferred(LinearMaps.MulStyle(Bd)) === FiveArg()
            @test occursin("25×39 LinearMaps.BlockDiagonalMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), Bd))
            @test Matrix(Bd) == Md
            @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC
            @test sparse(Bd) == Md
            @test Matrix(@inferred blockdiag(L1)) == M1
            @test Matrix(@inferred blockdiag(L1, L2)) == blockdiag(sparse.((M1, M2))...)
            Bd2 = @inferred cat(L1, L2, L3, L2, L1; dims=(1,2))
            @test_throws ArgumentError cat(L1, L2, L3, L2, L1; dims=(2,2))
            @test Bd == Bd2
            @test Bd == blockdiag(L1, M2, M3, M2, M1)
            @test size(Bd) == (25, 39)
            @test !issymmetric(Bd)
            @test !ishermitian(Bd)
            @test @inferred Bd * x ≈ Md * x
            for transform in (identity, adjoint, transpose)
                @test Matrix(@inferred transform(Bd)) == transform(Md)
                @test Matrix(@inferred transform(LinearMap(Bd))) == transform(Md)
            end
            y = randn(elty, size(Md, 1))
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(y), Bd, x, α, β) ≈ y*β .+ Md*x*α
            end
            X = randn(elty, size(Md, 2), 10)
            Y = randn(elty, size(Md, 1), 10)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), Bd, X, α, β) ≈ Y*β .+ Md*X*α
            end
        end
    end

    @testset "function block map" begin
        N = 100
        T = ComplexF64
        CS! = LinearMap{T}(cumsum!,
            (y, x) -> (copyto!(y, x); reverse!(cumsum!(y, reverse!(y)))), N;
            ismutating=true)
        A = rand(T, N, N)
        B = rand(T, N, N)
        LT = LowerTriangular(ones(T, N, N))
        L1 = [CS! CS! CS!; CS! CS! CS!; CS! CS! CS!]
        M1 = [LT LT LT; LT LT LT; LT LT LT]
        L2 = [CS! LinearMap(A) CS!; LinearMap(B) CS! CS!; CS! CS! CS!]
        M2 = [LT A LT
              B LT LT
              LT LT LT]
        u = rand(T, 3N)
        v = rand(T, 3N)
        for α in (false, true, rand(T)), β in (false, true, rand(T))
            for transform in (identity, adjoint), (L, M) in ((L1, M1), (L2, M2))
                # @show α, β, transform
                @test mul!(copy(v), transform(L), u, α, β) ≈ transform(M)*u*α + v*β
                @test mul!(copy(v), transform(LinearMap(L)), u, α, β) ≈ transform(M)*u*α + v*β
                @test mul!(copy(v), LinearMap(transform(L)), u, α, β) ≈ transform(M)*u*α + v*β
                bmap = @benchmarkable mul!($(copy(v)), $(transform(L)), $u, $α, $β)
                transform != adjoint && @test run(bmap, samples=3).memory < 2sizeof(u)
            end
        end
    end
end


text:
function buildconnectedgraph!(v::Array{Int,1}, G::DictGraph; weightfunction::Function = (a,b) -> 1)
    length(v) <= 1 && return
    x = pop!(v)
    for u in v
        addedge!(G, x, u; weight = weightfunction(x,u))
    end
    buildconnectedgraph!(v,G; weightfunction = weightfunction)
end

function buildconnectedgraph(v::Array{Int,1}, G::DictGraph = DictGraph(); weightfunction::Function = (a,b) -> 1)
    length(v) <= 1 && return G
    x = pop!(v)
    for u in v
        addedge!(G, x, u; weight = weightfunction(x,u))
    end
    buildconnectedgraph(v,G; weightfunction = weightfunction)
end


text:
# Overall sound pressure level
function f_oaspl(spl)

    # Compute
    oaspl = (10 * log10.( sum(10 .^(spl./10.), dims=2) ))

    return oaspl
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
struct TRBDF2Tableau{T,T2}
  γ::T2
  d::T
  ω::T
  btilde1::T
  btilde2::T
  btilde3::T
  α1::T
  α2::T
end

#=
Tableau:

c = [0; γ; 1]
A = [0  0  0
     d  d  0
     ω  ω  d]
b = [ω; ω; d]
bhat = [(1-ω)/3; (3ω+1)/3; d/3]

where γ=2-√2, d=γ/2, and ω=(√2)/4. Hence

btilde = bhat-b = [(1-4ω)/3; 1/3; -2d/3] = [(1-√2)/3; 1/3; (√2-2)/3]

Let uᵧ and u be approximations to u(t+γ*dt) and u(t+dt). Then initial guesses of
approximations zᵧ and z to dt*f(t+γ*dt,uᵧ) and dt*f(t+dt,u) according to Shampine:

zᵧ = zprev

z = (1.5+√2)*zprev + (2.5+2√2)*zᵧ - (6+4.5√2)*(uᵧ-uprev) = (by def. of uᵧ)
  = (1.5+√2)*zprev + (2.5+2√2)*zᵧ - (6+4.5√2)*(d*zᵧ + d*zprev) =
  = (-√2)/2*zprev + (1+(√2)/2)*zᵧ
=#
function TRBDF2Tableau(T,T2)
  γ = T2(2-sqrt(2))
  d = T(1-sqrt(2)/2)
  ω = T(sqrt(2)/4)
  btilde1 = T((1-sqrt(2))/3)
  btilde2 = T(1//3)
  btilde3 = T((sqrt(2)-2)/3)
  α1 = T(-sqrt(2)/2)
  α2 = T(1+sqrt(2)/2)
  TRBDF2Tableau(γ,d,ω,btilde1,btilde2,btilde3,α1,α2)
end

struct ESDIRK4Tableau{T,T2}
    γ::T
    a31::T
    a32::T
    a41::T
    a42::T
    a43::T
    btilde1::T
    btilde2::T
    btilde3::T
    btilde4::T
    c3::T2
    α31::T
    α32::T
    α41::T
    α42::T
end

#=
Derivative of Hermite Polynomial
k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt

Extrapolation for ESDIRK interior step 3
dt = c2 since interval is [c1,c2] and c1 = 0
θ =  c3/c2 the extrapolation point
z = dt*k

z₁ + Θ*(-4dt*z₁ - 2dt*z₂ - 6y₀ + Θ*(3dt*z₁ + 3z₂ + 6y₀ - 6y₁ ) + 6y₁)/dt

Test Expression on TRBDF2:
c2 = 2 - sqrt(2)
c3 = 1
θ = c3/c2; dt = c2

Coefficient on z₁:
(1 + (-4θ + 3θ^2))*z₁
1 + (-4θ + 3θ^2) - (1.5 + sqrt(2)) # 1.5 + sqrt(2) given by Shampine

Coefficient on z₂:
(-2θ + 3θ^2)*z₂
(-2θ + 3θ^2) - (2.5 + 2sqrt(2)) # 2.5 + 2sqrt(2) given by Shampine

Coefficient on y₁-y₀:
θ*(θ*6(y₀-y₁)+6(y₁-y₀))/dt
θ*(-6θ(y₁-y₀)+6(y₁-y₀))/dt
(y₁-y₀)6θ*(1-θ)/dt

(6θ*(1-θ)/dt)*(y₁-y₀)

6θ*(1-θ)/dt - (- (6 + 4.5sqrt(2)))  # - (6 + 4.5sqrt(2)) given by Shampine

# Write only in terms of z primatives
y₀ = uprev
y₁ = uprev + γ*z₁ + γ*z₂
y₁-y₀ = γ*z₁ + γ*z₂

# Full Expression
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/dt)*γ)*z₁ + ((-2θ + 3θ^2) + (6θ*(1-θ)/dt)*γ)*z₂
=#

#=
# Kvaerno3
# Predict z4 from yhat

yhat = uprev + a31*z1 + a32*z2 + γ*z3
z₄ = yhat - uprev = a31*z1 + a32*z2 + γ*z3

# Note Hermite is too small of an interval for this one!!!
=#

function Kvaerno3Tableau(T,T2)
  γ   = T(0.4358665215)
  a31 = T(0.490563388419108)
  a32 = T(0.073570090080892)
  a41 = T(0.308809969973036)
  a42 = T(1.490563388254106)
  a43 = -T(1.235239879727145)
  # bhat1 = T(0.490563388419108)
  # bhat2 = T(0.073570090080892)
  # bhat3 = T(0.4358665215)
  # bhat4 = T(0.0)
  btilde1 = T(0.181753418446072) # bhat1-a41
  btilde2 = T(-1.416993298173214) # bhat2-a42
  btilde3 = T(1.671106401227145) # bhat3-a43
  btilde4 = -γ # bhat4-γ
  c3 = T2(1)
  c2 = 2γ
  θ = c3/c2
  α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
  α41 = T(0.0)
  α42 = T(0.0)
  ESDIRK4Tableau(γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,α32,α41,α42)
end

struct KenCarp3Tableau{T,T2}
    γ::T
    a31::T
    a32::T
    a41::T
    a42::T
    a43::T
    btilde1::T
    btilde2::T
    btilde3::T
    btilde4::T
    c3::T2
    α31::T
    α32::T
    α41::T
    α42::T
    ea21::T
    ea31::T
    ea32::T
    ea41::T
    ea42::T
    ea43::T
    eb1::T
    eb2::T
    eb3::T
    eb4::T
    ebtilde1::T
    ebtilde2::T
    ebtilde3::T
    ebtilde4 ::T
end

#=
# KenCarp3
# Predict z4 from Hermite z2 and z1
# Not z3 because c3 < c2 !

θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
θ = c4/c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
=#
function KenCarp3Tableau(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  γ  = T(0.435866521508459)
  a31 = T(0.2576482460664272)
  a32 = -T(0.09351476757488625)
  a41 = T(0.18764102434672383)
  a42 = -T(0.595297473576955)
  a43 = T(0.9717899277217721)
  # bhat1 = T(2756255671327//12835298489170)
  # bhat2 = -T(10771552573575//22201958757719)
  # bhat3 = T(9247589265047//10645013368117)
  # bhat4 = T(2193209047091//5459859503100)
  btilde1 = T(0.027099261876665316) # bhat1-a41
  btilde2 = T(0.11013520969201586) # bhat2-a42
  btilde3 = T(-0.10306492520138458) # bhat3-a43
  btilde4 = T(-0.0341695463672966) # bhat4-γ
  c3 = T2(0.6)
  c2 = 2γ
  θ = c3/c2
  α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
  θ = 1/c2
  α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

  ea21 = T(0.871733043016918)
  ea31 = T(0.5275890119763004)
  ea32 = T(0.0724109880236996)
  ea41 = T(0.3990960076760701)
  ea42 = -T(0.4375576546135194)
  ea43 = T(1.0384616469374492)
  eb1 = T(0.18764102434672383)
  eb2 = T(-0.595297473576955)
  eb3 = T(0.9717899277217721)
  eb4 = T(0.435866521508459)
  ebtilde1 = T(0.027099261876665316)
  ebtilde2 = T(0.11013520969201586)
  ebtilde3 = -T(0.10306492520138458)
  ebtilde4 = -T(0.0341695463672966)

  KenCarp3Tableau(γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,
                  α32,α41,α42,ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,
                  ebtilde1,ebtilde2,ebtilde3,ebtilde4)
end

function KenCarp3Tableau(T,T2)
  γ  = T(1767732205903//4055673282236)
  a31 = T(2746238789719//10658868560708)
  a32 = -T(640167445237//6845629431997)
  a41 = T(1471266399579//7840856788654)
  a42 = -T(4482444167858//7529755066697)
  a43 = T(11266239266428//11593286722821)
  # bhat1 = T(2756255671327//12835298489170)
  # bhat2 = -T(10771552573575//22201958757719)
  # bhat3 = T(9247589265047//10645013368117)
  # bhat4 = T(2193209047091//5459859503100)
  btilde1 = T(parse(BigInt,"681815649026867975666107")//parse(BigInt,"25159934323302256049469295")) # bhat1-a41
  btilde2 = T(parse(BigInt,"18411887981491912264464127")//parse(BigInt,"167175311446532472108584143")) # bhat2-a42
  btilde3 = T(parse(BigInt,"-12719313754959329011138489")//parse(BigInt,"123410692144842870217698057")) # bhat3-a43
  btilde4 = T(parse(BigInt,"-47289384293135913063989")//parse(BigInt,"1383962894467812063558225")) # bhat4-γ
  c3 = T2(3//5)
  c2 = 2γ
  θ = c3/c2
  α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
  θ = 1/c2
  α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

  # Explicit Tableau
  ea21 = T(1767732205903//2027836641118)
  ea31 = T(5535828885825//10492691773637)
  ea32 = T(788022342437//10882634858940)
  ea41 = T(6485989280629//16251701735622)
  ea42 = -T(4246266847089//9704473918619)
  ea43 = T(10755448449292//10357097424841)
  eb1 = T(1471266399579//7840856788654)
  eb2 = T(-4482444167858//7529755066697 )
  eb3 = T(11266239266428//11593286722821)
  eb4 = T(1767732205903//4055673282236)
  ebtilde1 = T(parse(BigInt,"681815649026867975666107")//parse(BigInt,"25159934323302256049469295"))
  ebtilde2 = T(parse(BigInt,"18411887981491912264464127")//parse(BigInt,"167175311446532472108584143"))
  ebtilde3 = -T(parse(BigInt,"12719313754959329011138489")//parse(BigInt,"123410692144842870217698057"))
  ebtilde4 = -T(parse(BigInt,"47289384293135913063989")//parse(BigInt,"1383962894467812063558225"))
  KenCarp3Tableau(γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,
                  α32,α41,α42,ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,
                  ebtilde1,ebtilde2,ebtilde3,ebtilde4)
end

# Flip them all!

# ebtilde1 = big(1471266399579)//7840856788654 - big(2756255671327)//12835298489170
# ebtilde2 = -big(4482444167858)//7529755066697 + big(10771552573575)//22201958757719
# ebtilde3 = big(11266239266428)//11593286722821 - big(9247589265047)//10645013368117
# ebtilde4 = big(1767732205903)//4055673282236 - big(2193209047091)//5459859503100

struct Cash4Tableau{T,T2}
  γ::T
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  b1hat1::T
  b2hat1::T
  b3hat1::T
  b4hat1::T
  b1hat2::T
  b2hat2::T
  b3hat2::T
  b4hat2::T
  c2::T2
  c3::T2
  c4::T2
end

#=
Extrapolation for Cash interior step 3
dt = c1-c2 since interval is [c2,c1] and c1 = 0
c2 < c1, so z₂ is the left
θ =  (c3-c1)/dt the extrapolation point
z = dt*k

z₂ + Θ*(-4dt*z₂ - 2dt*z₁ - 6y₀ + Θ*(3dt*z₂ + 3z₁ + 6y₀ - 6y₁ ) + 6y₁)/dt

Coefficient on z₁:
(-2θ + 3θ^2)

Coefficient on z₂:
(1 + (-4θ + 3θ^2))

Coefficient on y₁-y₀:
(6θ*(1-θ)/dt)

# Write only in terms of z primatives
y₁ = uprev + a21*z₁ + γ*z₂
y₀ = uprev + γ*z₁
y₁-y₀ = (a21-γ)*z₁ + γ*z₂

θ = 1.1
Full z₁ coefficient: (-2θ + 3θ^2) + (6θ*(1-θ)/dt)*(a21-γ)
Full z₂ coefficient: (1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/dt)*γ

f(θ)= (-2θ + 3θ^2) + (6θ*(1-θ)/dt)*(a21-γ)
g(θ) = (1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/dt)*γ
t = linspace(0,1.5,100)
y = f.(t)
z = g.(t)
plot(t,y)
plot!(t,z)

The extrapolation is really bad that far
Hairer's extrapolation is no better.
Using constant extrapolations
=#

function Cash4Tableau(T,T2)
  γ = T(0.435866521508)
  a21 = T(-1.13586652150)
  a31 = T(1.08543330679)
  a32 = -T(0.721299828287)
  a41 = T(0.416349501547)
  a42 = T(0.190984004184)
  a43 = T(-0.118643265417)
  a51 = T(0.896869652944)
  a52 = T(0.0182725272734)
  a53 = -T(0.0845900310706)
  a54 = -T(0.266418670647)
  b1hat1 = T(1.05646216107052)
  b2hat1 = -T(0.0564621610705236)
  b3hat1 = T(0)
  b4hat1 = T(0)
  b1hat2 = T(0.776691932910)
  b2hat2 = T(0.0297472791484)
  b3hat2 = -T(0.0267440239074)
  b4hat2 = T(0.220304811849)
  c2 = -T2(0.7)
  c3 = T2(0.8)
  c4 = T2(0.924556761814)
  Cash4Tableau(γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
               b1hat1,b2hat1,b3hat1,b4hat1,b1hat2,b2hat2,b3hat2,b4hat2,c2,c3,c4)
end

struct Hairer4Tableau{T,T2}
  γ::T
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  bhat1::T
  bhat2::T
  bhat3::T
  bhat4::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  btilde5::T
  c2::T2
  c3::T2
  c4::T2
  r11::T
  r12::T
  r13::T
  r14::T
  r21::T
  r22::T
  r23::T
  r24::T
  r31::T
  r32::T
  r33::T
  r34::T
  r41::T
  r42::T
  r43::T
  r44::T
  r51::T
  r52::T
  r53::T
  r54::T
  α21::T
  α31::T
  α32::T
  α41::T
  α43::T
end

function Hairer4Tableau(T,T2)
  γ = T(1//4)
  c2 = T(3//4)
  c3 = T(11//20)
  c4 = T(1//2)

  #=
  α21 = T(2)
  α31 = T(42//25)
  α32 = -T(4//25)
  α41 = T(89//68)
  α42 = -T(25//136)
  α43 = T(15//136)
  α51 = -T(37//12)
  α52 = -T(103//24)
  α53 = T(275//8)
  α54 = -T(85//3)

  alpha = -inv(A)*γ

  α = [-1 0 0 0 0
       α21 -1 0 0 0
       α31 α32 -1 0 0
       α41 α42 α43 -1 0
       α51 α52 α53 α54 -1]
  A α = -γ
  A = -γ*inv(α)

  Now zⱼ = fⱼ + ∑α_ⱼᵢ zᵢ
  =#

  a11 = T(1//4)
  a21 = T(1//2)
  a22 = T(1//4)
  a31 = T(17//50)
  a32 = T(-1//25)
  a33 = T(1//4)
  a41 = T(371//1360)
  a42 = T(-137//2720)
  a43 = T(15//544)
  a44 = T(1//4)
  a51 = T(25//24)
  a52 = T(-49//48)
  a53 = T(125//16)
  a54 = T(-85//12)

  #=
  e1 = -T(23//6)
  e2 = -T(17//12)
  e3 = T(125//4)
  e4 = -T(85//3)
  E = [e1 e2 e3 e4 0]

  bhat = [59//48,-17//96,225//32,-85//12,0]

  α = [-1 0 0 0 0
       α21 -1 0 0 0
       α31 α32 -1 0 0
       α41 α42 α43 -1 0
       e1 e2 e3 e4 -1]

  A = [γ 0 0 0 0
       a21 γ 0 0 0
       a31 a32 γ 0 0
       a41 a42 a43 γ 0
       a51 a52 a53 a54 γ]

  E = bhat'*inv(A)
  bhat = E*A
  =#

  bhat1 = T(59//48)
  bhat2 = T(-17//96)
  bhat3 = T(225//32)
  bhat4 = T(-85//12)

  btilde1 = T(3//16) # bhat1-a51
  btilde2 = T(27//32) # bhat2-a52
  btilde3 = T(-25//32) # bhat3-a53
  btilde4 = T(0) # bhat4-a54
  btilde5 = -γ

  #=
  d11 = T(61//27)
  d12 = T(-185//54)
  d13 = T(2525//18)
  d14 = T(-3740//27)
  d15 = T(-44//9)
  d21 = T(2315//81)
  d22 = T(1049//162)
  d23 = T(-27725//54)
  d24 = T(40460//81)
  d25 = T(557//27)
  d31 = T(-6178//81)
  d32 = T(-1607//81)
  d33 = T(20075//27)
  d34 = T(-56440//81)
  d35 = T(-718//27)
  d41 = T(3680//81)
  d42 = T(1360//81)
  d43 = T(-10000//27)
  d44 = T(27200//81)
  d45 = T(320//27)

  D = [d11 d12 d13 d14 d15
       d21 d22 d23 d24 d25
       d31 d32 d33 d34 d35
       d41 d42 d43 d44 d45]
  R = (D*A)'
  =#

  r11 = T(11//3)
  r12 = T(-463//72)
  r13 = T(217//36)
  r14 = T(-20//9)
  r21 = T(11//2)
  r22 = T(-385//16)
  r23 = T(661//24)
  r24 = T(-10//1)
  r31 = T(-125//18)
  r32 = T(20125//432)
  r33 = T(-8875//216)
  r34 = T(250//27)
  r41 = T(0)
  r42 = T(-85//4)
  r43 = T(85//6)
  r44 = T(0//1)
  r51 = T(-11//9)
  r52 = T(557//108)
  r53 = T(-359//54)
  r54 = T(80//27)

  # c2/γ
  α21 = T(3)
  #=
  # Prediction alphas from Hairer
  # Predict z3 from z1 and z2
  A = [c1 c2
  γ*c1 a21*c1+γ*c2]
  b = [c3,a31*c1+a32*c2+γ*c3]
  A\b
  =#
  α31 = T(88//100)
  α32 = T(44//100)
  #=
  # Predict z4 from z1 and z3
  A = [c1   c3
       γ*c1 a31*c1+a32*c2+γ*c3]
  b = [c4,a41*c1+a42*c2+a43*c3+γ*c4]
  A\b
  =#
  α41 = T(3//17)
  α43 = T(155//187)

 Hairer4Tableau(γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                bhat1,bhat2,bhat3,bhat4,btilde1,btilde2,btilde3,
                btilde4,btilde5,c2,c3,c4,r11,r12,r13,r14,
                r21,r22,r23,r24,r31,r32,r33,r34,r41,r42,r43,r44,r51,
                r52,r53,r54,α21,α31,α32,α41,α43)
end

function Hairer42Tableau(T,T2)
  γ  = T(4//15)
  c2 = T2(23//30)
  c3 = T2(17//30)
  c4 = T2(2881//28965)+γ
  #=
  α21 = T(15//8)
  α31 = T(1577061//922880)
  α32 = -T(23427//115360)
  α41 = T(647163682356923881//2414496535205978880)
  α42 = -T(593512117011179//3245291041943520)
  α43 = T(559907973726451//1886325418129671)
  α51 = T(724545451//796538880)
  α52 = -T(830832077//267298560)
  α53 = T(30957577//2509272)
  α54 = -T(69863904375173//6212571137048)

  α = [-1 0 0 0 0
       α21 -1 0 0 0
       α31 α32 -1 0 0
       α41 α42 α43 -1 0
       α51 α52 α53 α54 -1]
  A = -γ*inv(α)
  =#
  a21 = T(1//2)
  a31 = T(51069//144200)
  a32 = T(-7809//144200)
  a41 = T(12047244770625658//141474406359725325)
  a42 = T(-3057890203562191//47158135453241775)
  a43 = T(2239631894905804//28294881271945065)
  a51 = T(181513//86430)
  a52 = T(-89074//116015)
  a53 = T(83636//34851)
  a54 = T(-69863904375173//23297141763930)

  #=

  A = [γ 0 0 0 0
       a21 γ 0 0 0
       a31 a32 γ 0 0
       a41 a42 a43 γ 0
       a51 a52 a53 a54 γ]

  A = convert(Array{Rational{BigInt}},A)

  e1 = T(7752107607//11393456128)
  e2 = -T(17881415427//11470078208)
  e3 = T(2433277665//179459416)
  e4 = -T(96203066666797//6212571137048)
  E = [e1 e2 e3 e4 0]

  bhat = E*A
  =#

  bhat1 = T(33665407//11668050)
  bhat2 = T(-2284766//15662025)
  bhat3 = T(11244716//4704885)
  bhat4 = T(-96203066666797//23297141763930)

  btilde1 = T(4580576//5834025) # bhat1-a51
  btilde2 = T(9740224//15662025) # bhat2-a52
  btilde3 = T(-46144//4704885) # bhat3-a53
  btilde4 = T(-13169581145812//11648570881965) # bhat4-a54
  btilde5 = -γ

  #=
  d11 = T(24.74416644927758)
  d12 = -T(4.325375951824688)
  d13 = T(41.39683763286316)
  d14 = T(-61.04144619901784)
  d15 = T(-3.391332232917013)
  d21 = T(-51.98245719616925)
  d22 = T(10.52501981094525)
  d23 = T(-154.2067922191855)
  d24 = T(214.3082125319825)
  d25 = T(14.71166018088679)
  d31 = T(33.14347947522142)
  d32 = T(-19.72986789558523)
  d33 = T(230.4878502285804)
  d34 = T(-287.6629744338197)
  d35 = T(-18.99932366302254)
  d41 = T(-5.905188728329743)
  d42 = T(13.53022403646467)
  d43 = T(-117.6778956422581)
  d44 = T(134.3962081008550)
  d45 = T(8.678995715052762)
  =#

  r11 = T(6.776439256624082)
  r12 = T(-14.066831911883533)
  r13 = T(16.204808856162565)
  r14 = T(-6.8143005003361985)
  r21 = T(3.166691382949011)
  r22 = T(-14.034196189427504)
  r23 = T(15.497198116229603)
  r24 = T(-5.3974733381957005)
  r31 = T(-1.9310469085972866)
  r32 = T(11.146663701107887)
  r33 = T(-6.9009212321038405)
  r34 = T(0.085120800673252)
  r41 = T(-6.107728468864632)
  r42 = T(13.031255018633459)
  r43 = T(-19.734599430149146)
  r44 = T(9.812254180511282)
  r51 = T(-0.9043552621112034)
  r52 = T(3.9231093815698106)
  r53 = T(-5.066486310139344)
  r54 = T(2.3143988573474035)

  # c2/γ
  α21 = T(23//8)
  α31 = T(0.9838473040915402)
  α32 = T(0.3969226768377252)
  α41 = T(0.6563374010466914)
  α43 = T(0.3372498196189311)

  Hairer4Tableau(γ,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                 bhat1,bhat2,bhat3,bhat4,btilde1,btilde2,btilde3,btilde4,btilde5,
                 c2,c3,c4,r11,r12,r13,r14,
                 r21,r22,r23,r24,r31,r32,r33,r34,r41,r42,r43,r44,r51,
                 r52,r53,r54,α21,α31,α32,α41,α43)
end

struct Kvaerno4Tableau{T,T2}
  γ::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  btilde5::T
  c3::T2
  c4::T2
  α21::T
  α31::T
  α32::T
  α41::T
  α42::T
end

#=
# Kvaerno4
# Predict z3 from Hermite z2 and z1

c2 = 2γ
θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z4 from Hermite z2 and z1

θ = c4/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
=#

function Kvaerno4Tableau(T,T2)
  γ = T(0.4358665215)
  a31 = T(0.140737774731968)
  a32 = T(-0.108365551378832)
  a41 = T(0.102399400616089)
  a42 = T(-0.376878452267324)
  a43 = T(0.838612530151233)
  a51 = T(0.157024897860995)
  a52 = T(0.117330441357768)
  a53 = T(0.61667803039168)
  a54 = T(-0.326899891110444)
  btilde1 = T(-0.054625497244906) # a41 - a51
  btilde2 = T(-0.494208893625092) # a42 - a52
  btilde3 = T(0.221934499759553) # a43 - a53
  btilde4 = T(0.762766412610444) # γ - a54
  btilde5 = -γ
  c3 = T2(0.468238744853136)
  c4 = T2(1)
  α21 = T(2) # c2/γ
  α31 = T(0.462864521870446)
  α32 = T(0.537135478129554)
  α41 = T(-0.14714018016178376)
  α42 = T(1.1471401801617838)
  Kvaerno4Tableau(γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                  btilde1,btilde2,btilde3,btilde4,btilde5,
                  c3,c4,α21,α31,α32,α41,α42)
end

struct KenCarp4Tableau{T,T2}
  γ::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  c3::T2
  c4::T2
  c5::T2
  α21::T
  α31::T
  α32::T
  α41::T
  α42::T
  α51::T
  α52::T
  α53::T
  α54::T
  α61::T
  α62::T
  α63::T
  α64::T
  α65::T
  ea21::T
  ea31::T
  ea32::T
  ea41::T
  ea42::T
  ea43::T
  ea51::T
  ea52::T
  ea53::T
  ea54::T
  ea61::T
  ea62::T
  ea63::T
  ea64::T
  ea65::T
  eb1::T
  eb3::T
  eb4::T
  eb5::T
  eb6::T
  ebtilde1::T
  ebtilde3::T
  ebtilde4::T
  ebtilde5::T
  ebtilde6::T
end

#=
# KenCarp4
# Predict z3 from Hermite z2 and z1

c2 = 2γ
θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z4 from Hermite z2 and z1
θ = c4/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z5 from Hermite z2 and z1
θ = c5/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z5 from z1 and z4

θ = c5/c4
dt = c4

y₀ = uprev
y₁ = uprev + a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄
y₁-y₀ = a41*z₁ + a42*z₂ + a43*z₃ + γ*z₄

(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))*z₁ +
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))*z₅
+ (6θ*(1-θ)/dt)*(a52*z₂ + a53*z₃ + a54*z₄)

(1 + (-4θ + 3θ^2) + a41*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a42
(6θ*(1-θ)/dt)*a43
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

# Predict last stage from z1 and z5

θ = 1/c5
dt = c5

y₀ = uprev
y₁ = uprev + a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄ + γ*z₅
y₁-y₀ = a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄ + γ*z₅

(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))*z₁ +
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))*z₅
+ (6θ*(1-θ)/dt)*(a52*z₂ + a53*z₃ + a54*z₄)

(1 + (-4θ + 3θ^2) + a51*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a52
(6θ*(1-θ)/dt)*a53
(6θ*(1-θ)/dt)*a54
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

=#
function KenCarp4Tableau(T,T2)
  γ = T(1//4)
  a31 = T(8611//62500)
  a32 = -T(1743//31250)
  a41 = T(5012029//34652500)
  a42 = -T(654441//2922500)
  a43 = T(174375//388108)
  a51 = T(15267082809//155376265600)
  a52 = -T(71443401//120774400)
  a53 = T(730878875//902184768)
  a54 = T(2285395//8070912)
  a61 = T(82889//524892)
  a63 = T(15625//83664)
  a64 = T(69875//102672)
  a65 = -T(2260//8211)
  # bhat1 = T(4586570599//29645900160)
  # bhat3 = T(178811875//945068544)
  # bhat4 = T(814220225//1159782912)
  # bhat5 = -T(3700637//11593932)
  # bhat6 = T(61727//225920)
  btilde1 = T(-31666707//9881966720) # bhat1-a61
  btilde3 = T(256875//105007616) # bhat3-a63
  btilde4 = T(2768025//128864768) # bhat4-a64
  btilde5 = T(-169839//3864644) # bhat5-a65
  btilde6 = T(5247//225920) # bhat6-γ
  c3 = T2(83//250)
  c4 = T2(31//50)
  c5 = T2(17//20)
  α21 = T(2) # c2/γ
  α31 = T(42//125)
  α32 = T(83//125)
  α41 = T(-6//25)
  α42 = T(31//25)
  α51 = T(914470432//2064665255)
  α52 = T(798813//724780)
  α53 = T(-824765625//372971788)
  α54 = T(49640//29791)
  α61 = T(288521442795//954204491116)
  α62 = T(2224881//2566456)
  α63 = T(-1074821875//905317354)
  α64 = T(-3360875//8098936)
  α65 = T(7040//4913)

  ea21 = T(1//2)
  ea31 = T(13861//62500)
  ea32 = T(6889//62500)
  ea41 = -T(116923316275//2393684061468)
  ea42 = -T(2731218467317//15368042101831)
  ea43 = T(9408046702089//11113171139209)
  ea51 = -T(451086348788//2902428689909)
  ea52 = -T(2682348792572//7519795681897)
  ea53 = T(12662868775082//11960479115383)
  ea54 = T(3355817975965//11060851509271)
  ea61 = T(647845179188//3216320057751)
  ea62 = T(73281519250//8382639484533)
  ea63 = T(552539513391//3454668386233)
  ea64 = T(3354512671639//8306763924573)
  ea65 = T(4040//17871)

  eb1 = T(82889//524892)
  eb3 = T(15625//83664)
  eb4 = T(69875//102672)
  eb5 = -T(2260//8211)
  eb6 = T(1//4)

  ebtilde1 = -T(31666707//9881966720)
  ebtilde3 =  T(256875//105007616)
  ebtilde4 =  T(2768025//128864768)
  ebtilde5 = -T(169839//3864644)
  ebtilde6 =  T(5247//225920)

  KenCarp4Tableau(γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a63,a64,a65,
                  btilde1,btilde3,btilde4,btilde5,btilde6,
                  c3,c4,c5,
                  α21,α31,α32,α41,α42,α51,α52,α53,α54,α61,α62,α63,α64,α65,
                  ea21,ea31,ea32,ea41,ea42,ea43,ea51,ea52,ea53,ea54,ea61,ea62,
                  ea63,ea64,ea65,eb1,eb3,eb4,eb5,eb6,ebtilde1,ebtilde3,ebtilde4,
                  ebtilde5,ebtilde6)
end

# Flip them all!

# ebtilde1 = 82889//524892 - 4586570599//29645900160
# ebtilde3 = 15625//83664 - 178811875//945068544
# ebtilde4 = 69875//102672 - 814220225//1159782912
# ebtilde5 = -2260//8211 + 3700637//11593932
# ebtilde6 = 1//4 - 61727//225920

struct Kvaerno5Tableau{T,T2}
  γ::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  α31::T
  α32::T
  α41::T
  α42::T
  α43::T
  α51::T
  α52::T
  α53::T
  α61::T
  α62::T
  α63::T
end

#=
# Kvaerno5
# Predict z3 from Hermite z2 and z1

c2 = 2γ
θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict others from z1 and z3 since it covers [0,1.23]

dt = c3 since interval is [c1,c3] and c1 = 0
θ =  c4/c3, c5/c3, c6/c3, c7/c3
z = dt*k

z₁ + Θ*(-4dt*z₁ - 2dt*z₃ - 6y₀ + Θ*(3dt*z₁ + 3z₃ + 6y₀ - 6y₁ ) + 6y₁)/dt

(1 + (-4θ + 3θ^2))*z₁ + (-2θ + 3θ^2)*z₃ + (6θ*(1-θ)/dt)*(y₁-y₀)

y₀ = uprev
y₁ = uprev + a31*z₁ + a32*z₂ + γ*z₃
y₁-y₀ = a31*z₁ + a32*z₂ + γ*z₃

(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))*z₁ +
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))*z₃ + (6θ*(1-θ)/dt)*a32*z₂

dt = c3
θ = c4/c3
(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a32
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

θ = c5/c3
(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a32
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

θ = c6/c3
(1 + (-4θ + 3θ^2) + a31*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a32
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))
=#

function Kvaerno5Tableau(T,T2)
  γ = T(0.26)
  a31 = T(0.13)
  a32 = T(0.84033320996790809)
  a41 = T(0.22371961478320505)
  a42 = T(0.47675532319799699)
  a43 = -T(0.06470895363112615)
  a51 = T(0.16648564323248321)
  a52 = T(0.10450018841591720)
  a53 = T(0.03631482272098715)
  a54 = -T(0.13090704451073998)
  a61 = T(0.13855640231268224)
  a63 = -T(0.04245337201752043)
  a64 = T(0.02446657898003141)
  a65 = T(0.61943039072480676)
  a71 = T(0.13659751177640291)
  a73 = -T(0.05496908796538376)
  a74 = -T(0.04118626728321046)
  a75 = T(0.62993304899016403)
  a76 = T(0.06962479448202728)
  btilde1 = T(0.00195889053627933) # a61-a71
  btilde3 = T(0.01251571594786333) # a63-a73
  btilde4 = T(0.06565284626324187) # a64-a74
  btilde5 = -T(0.01050265826535727) # a65-a75
  btilde6 = T(0.19037520551797272) # γ-a76
  btilde7 = -γ
  α21 = T(2) # c2/γ
  α31 = T(-1.366025403784441)
  α32 = T(2.3660254037844357)
  α41 = T(-0.19650552613122207)
  α42 = T(0.8113579546496623)
  α43 = T(0.38514757148155954)
  α51 = T(0.10375304369958693)
  α52 = T(0.937994698066431)
  α53 = T(-0.04174774176601781)
  α61 = T(-0.17281112873898072)
  α62 = T(0.6235784481025847)
  α63 = T(0.5492326806363959)
  c3 = T(1.230333209967908)
  c4 = T(0.895765984350076)
  c5 = T(0.436393609858648)
  c6 = T(1)
  Kvaerno5Tableau(γ,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                  a61,a63,a64,a65,a71,a73,a74,a75,a76,
                  btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,
                  c3,c4,c5,c6,α31,α32,α41,α42,α43,α51,α52,α53,
                  α61,α62,α63)
end

struct KenCarp5Tableau{T,T2}
  γ::T
  a31::T
  a32::T
  a41::T
  a43::T
  a51::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a84::T
  a85::T
  a86::T
  a87::T
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  c7::T2
  α31::T
  α32::T
  α41::T
  α42::T
  α51::T
  α52::T
  α61::T
  α62::T
  α71::T
  α72::T
  α73::T
  α74::T
  α75::T
  α81::T
  α82::T
  α83::T
  α84::T
  α85::T
  btilde1::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  btilde8::T
  ea21::T
  ea31::T
  ea32::T
  ea41::T
  ea43::T
  ea51::T
  ea53::T
  ea54::T
  ea61::T
  ea63::T
  ea64::T
  ea65::T
  ea71::T
  ea73::T
  ea74::T
  ea75::T
  ea76::T
  ea81::T
  ea83::T
  ea84::T
  ea85::T
  ea86::T
  ea87::T
  eb1::T
  eb4::T
  eb5::T
  eb6::T
  eb7::T
  eb8::T
  ebtilde1::T
  ebtilde4::T
  ebtilde5::T
  ebtilde6::T
  ebtilde7::T
  ebtilde8::T
end

#=
# KenCarp5
# Predict z3 from Hermite z2 and z1

c2 = 2γ
θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z4 from z2 and z1
θ = c4/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z5 from z2 and z1
θ = c5/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z6 from z2 and z1
θ = c6/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

# Predict z7 from z5 and z1
θ = c7/c5
dt = c5

(1 + (-4θ + 3θ^2))*z₁ + (-2θ + 3θ^2)*z₃ + (6θ*(1-θ)/dt)*(y₁-y₀)
y₁-y₀ = a51*z₁ + a52*z₂ + a53*z₃ + a54*z₄ + γ*z₅

(1 + (-4θ + 3θ^2) + a51*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a52
(6θ*(1-θ)/dt)*a53
(6θ*(1-θ)/dt)*a54
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

# Predict z8 from z5 and z1
θ = 1/c5
dt = c5

(1 + (-4θ + 3θ^2) + a51*(6θ*(1-θ)/dt))
(6θ*(1-θ)/dt)*a52
(6θ*(1-θ)/dt)*a53
(6θ*(1-θ)/dt)*a54
(-2θ + 3θ^2 + γ*(6θ*(1-θ)/dt))

=#

function KenCarp5Tableau(T,T2)
  γ = T(41//200)
  a31 = T(41//400)
  a32 = -T(567603406766//11931857230679)
  a41 = T(683785636431//9252920307686)
  a43 = -T(110385047103//1367015193373)
  a51 = T(3016520224154//10081342136671)
  a53 = T(30586259806659//12414158314087)
  a54 = -T(22760509404356//11113319521817)
  a61 = T(218866479029//1489978393911)
  a63 = T(638256894668//5436446318841)
  a64 = -T(1179710474555//5321154724896)
  a65 = -T(60928119172//8023461067671)
  a71 = T(1020004230633//5715676835656)
  a73 = T(25762820946817//25263940353407)
  a74 = -T(2161375909145//9755907335909)
  a75 = -T(211217309593//5846859502534)
  a76 = -T(4269925059573//7827059040749)
  a81 = -T(872700587467//9133579230613)
  a84 = T(22348218063261//9555858737531)
  a85 = -T(1143369518992//8141816002931)
  a86 = -T(39379526789629//19018526304540)
  a87 = T(32727382324388//42900044865799)
  # bhat1 = -T(975461918565//9796059967033)
  # bhat4 = T(78070527104295//32432590147079)
  # bhat5 = -T(548382580838//3424219808633)
  # bhat6 = -T(33438840321285//15594753105479)
  # bhat7 = T(3629800801594//4656183773603)
  # bhat8 = T(4035322873751//18575991585200)
  btilde1 = -T(360431431567533808054934//89473089856732078284381229) # bhat1-a81
  btilde4 = T(21220331609936495351431026//309921249937726682547321949) # bhat4-a84
  btilde5 = -T(42283193605833819490634//2144566741190883522084871) # bhat5-a85
  btilde6 = -T(21843466548811234473856609//296589222149359214696574660) # bhat6-a86
  btilde7 = T(3333910710978735057753642//199750492790973993533703797) # bhat7-a87
  btilde8 = T(45448919757//3715198317040) # bhat8-γ
  c3 = T2(2935347310677//11292855782101)
  c4 = T2(1426016391358//7196633302097)
  c5 = T2(92//100)
  c6 = T2(24//100)
  c7 = T2(3//5)
  α31 = T(169472355998441//463007087066141)
  α32 = T(293534731067700//463007087066141)
  α41 = T(152460326250177//295061965385977)
  α42 = T(142601639135800//295061965385977)
  α51 = T(-51//41)
  α52 = T(92//41)
  α61 = T(17//41)
  α62 = T(24//41)
  α71 = T(13488091065527792//122659689776876057)
  α72 = T(-3214953045//3673655312)
  α73 = T(550552676519862000//151043064207496529)
  α74 = T(-409689169278408000//135215758621947439)
  α75 = T(3345//12167)
  α81 = T(1490668709762032//122659689776876057)
  α82 = T(5358255075//14694621248)
  α83 = T(-229396948549942500//151043064207496529)
  α84 = T(170703820532670000//135215758621947439)
  α85 = T(30275//24334)

  ea21 = T(41//100)
  ea31 = T(367902744464//2072280473677)
  ea32 = T(677623207551//8224143866563)
  ea41 = T(1268023523408//10340822734521)
  ea43 = T(1029933939417//13636558850479)
  ea51 = T(14463281900351//6315353703477)
  ea53 = T(66114435211212//5879490589093)
  ea54 = -T(54053170152839//4284798021562)
  ea61 = T(14090043504691//34967701212078)
  ea63 = T(15191511035443//11219624916014)
  ea64 = -T(18461159152457//12425892160975)
  ea65 = -T(281667163811//9011619295870)
  ea71 = T(19230459214898//13134317526959)
  ea73 = T(21275331358303//2942455364971)
  ea74 = -T(38145345988419//4862620318723)
  ea75 = -T(1//8)
  ea76 = -T(1//8)
  ea81 = -T(19977161125411//11928030595625)
  ea83 = -T(40795976796054//6384907823539)
  ea84 = T(177454434618887//12078138498510)
  ea85 = T(782672205425//8267701900261)
  ea86 = -T(69563011059811//9646580694205)
  ea87 = T(7356628210526//4942186776405)

  eb1 = -T(872700587467//9133579230613)
  eb4 = T(22348218063261//9555858737531)
  eb5 = -T(1143369518992//8141816002931)
  eb6 = -T(39379526789629//19018526304540)
  eb7 = T(32727382324388//42900044865799)
  eb8 = T(41//200)

  ebtilde1 = -T(360431431567533808054934//89473089856732078284381229)
  ebtilde4 = T(21220331609936495351431026//309921249937726682547321949)
  ebtilde5 = -T(42283193605833819490634//2144566741190883522084871)
  ebtilde6 = -T(21843466548811234473856609//296589222149359214696574660)
  ebtilde7 = T(3333910710978735057753642//199750492790973993533703797)
  ebtilde8 = T(45448919757//3715198317040)

  KenCarp5Tableau(γ,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,
                  a71,a73,a74,a75,a76,a81,a84,a85,a86,a87,
                  c3,c4,c5,c6,c7,α31,α32,α41,α42,α51,α52,
                  α61,α62,α71,α72,α73,α74,α75,α81,α82,α83,α84,α85,
                  btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,
                  ea21,ea31,ea32,ea41,ea43,ea51,ea53,ea54,ea61,ea63,
                  ea64,ea65,ea71,ea73,ea74,ea75,ea76,ea81,ea83,ea84,
                  ea85,ea86,ea87,eb1,eb4,eb5,eb6,eb7,eb8,ebtilde1,
                  ebtilde4,ebtilde5,ebtilde6,ebtilde7,ebtilde8)
end

# Flip them all!

# ebtilde1 = -Int128(872700587467)//9133579230613 + Int128(975461918565)//9796059967033
# ebtilde4 = Int128(22348218063261)//9555858737531 - Int128(78070527104295)//32432590147079
# ebtilde5 = -Int128(1143369518992)//8141816002931 + Int128(548382580838)//3424219808633
# ebtilde6 = -Int128(39379526789629)//19018526304540 + Int128(33438840321285)//15594753105479
# ebtilde7 = Int128(32727382324388)//42900044865799 - Int128(3629800801594)//4656183773603
# ebtilde8 = Int128(41)//200 - Int128(4035322873751)//18575991585200


text:
packages = [
    let (name, version) = split(line)
        (name, VersionNumber(version))
    end
    for line in split(strip(read("REQUIRE", String)), '\n')[2:end]
]

for (name, version) in packages
    info("Pkg.add($name)")
    Pkg.add(name)  # to make the following work
end
names = map(first, packages)
info("Pkg.free & .update: $(join(names, ", "))")
Pkg.free(names)            # required for Pkg.update
Pkg.update(names...)       # to update METADATA and package repository

for (name, version) in packages
    info("Pkg.pin($name, $version)")
    Pkg.pin(name, version)
end

println("Pkg.resolve()")
Pkg.resolve()


text:
"""
    LOCF(; limit=nothing)

Last observation carried forward (LOCF) iterates forwards through the `data` and fills
missing data with the last existing observation. The current implementation is univariate,
so each variable in a table or matrix will be handled independently.

See also:
- [`Impute.NOCB`](@ref): Next Observation Carried Backward

!!! Missing elements at the head of the array may not be imputed if there is no
existing observation to carry forward. As a result, this method does not guarantee
that all missing values will be imputed.

# Keyword Arguments
* `limit::Union{UInt, Nothing}`: Optionally limits the amount of consecutive missing values to replace.

# Example
```jldoctest
julia> using Impute: LOCF, impute

julia> M = [1.0 2.0 missing missing missing 6.0; 1.1 missing missing 4.4 5.5 6.6]
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0       missing   missing   missing  6.0
 1.1   missing  missing  4.4       5.5       6.6

julia> impute(M, LOCF(); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.0  2.0  2.0  6.0
 1.1  1.1  1.1  4.4  5.5  6.6

julia> impute(M,  LOCF(; limit=2); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.0  2.0   missing  6.0
 1.1  1.1  1.1  4.4  5.5       6.6
```
"""
struct LOCF <: Imputor
    limit::Union{UInt, Nothing}
end

LOCF(; limit=nothing) = LOCF(limit)

function _impute!(data::AbstractVector{Union{T, Missing}}, imp::LOCF) where T
    @assert !all(ismissing, data)
    start_idx = findfirst(!ismissing, data)
    count = 1

    for i in start_idx + 1:lastindex(data)
        if ismissing(data[i])
            if imp.limit === nothing
                data[i] = data[i-1]
            elseif count <= imp.limit
                data[i] = data[start_idx]
                count += 1
            end
        else
            start_idx = i
            count = 1
        end
    end

    return data
end


text:
using EventSimulation
using Distributions
using Statistics
using Random
using Test

mutable struct MMQueue <: AbstractState
    sp::Exponential{Float64} # service rate
    qlen::Int   # queue length
    busy::Bool  # is server busy?
    arrival_times::Vector{Float64} # when customers arrived
    departure_times::Vector{Float64} # when customers left
end

# customer arrival event
function arrival(s)
    push!(s.state.arrival_times, s.now)
    if s.state.busy # server is working?
        s.state.qlen += 1
    else
        s.state.busy = true
        register!(s, leave, rand(s.state.sp))
    end
end

# customer leave the system event
function leave(s)
    push!(s.state.departure_times, s.now)
    if s.state.qlen > 0 # any customers waiting?
        s.state.qlen -= 1
        register!(s, leave, rand(s.state.sp))
    else
        s.state.busy = false
    end
end

function run(ar, sr)
    q = MMQueue(Exponential(1 / sr), 0, false, Float64[], Float64[])
    s = Scheduler(q, Float64)
    repeat_register!(s, arrival, x -> rand(Exponential(1/ar)))
    go!(s, 1_000_000)
    mean(d-a for (a,d) in zip(s.state.arrival_times, s.state.departure_times))
end

Random.seed!(1)
for ar in 0.1:0.1:0.9
    @test abs(run(ar, 1.0)*(1-ar)-1) < 0.03
end


text:
"""
dₖ = solve_tCG(M, x, gradf_k, hessf_x_h, ν=1e-2, ϵ_residual = 1e-5)

Conjugate gradient method to solve the linear problem:
find d s.t. hessf_x_h(d) + gradf_k = 0

Note that we are looking for a direction d in the tangent space of M at x, and the kernel of
∇²f_M contains at least the orthogonal to this set. CG should naturally give a solution
living in the tangent space.
"""
function solve_tCG(M, x, gradfₖ, hessf_x_h; ν=1e-3, ϵ_residual = 1e-13, maxiter=1e5, safe_projection=false, printlev=0)
    norm_x(η) = norm(M, x, η)

    dₖ = zero_tangent_vector(M, x)
    rⱼ, rⱼ_prev = zero_tangent_vector(M, x), zero_tangent_vector(M, x)
    vⱼ, vⱼ_prev = zero_tangent_vector(M, x), zero_tangent_vector(M, x)

    d_type = :Unsolved
    j = 0

    (printlev>0) && @printf "\nj     norm(rⱼ)             norm(vⱼ)             ⟨vⱼ, hessf(x)[vⱼ]⟩   ν * norm(vⱼ)^2\n"
    while true
        # current residual, conjugated direction
        rⱼ = hessf_x_h(dₖ) + gradfₖ
        if safe_projection
            rⱼ = project(M, x, rⱼ)
        end

        vⱼ = - rⱼ
        if j ≥ 1
            βⱼ = norm_x(rⱼ)^2 / norm_x(rⱼ_prev)^2
            vⱼ += βⱼ * vⱼ_prev
        end

        hessf_x_vⱼ = hessf_x_h(vⱼ)
        (printlev>0) && @printf "%5i %.10e    %.10e    % .10e     %.10e\n" j norm_x(rⱼ) norm_x(vⱼ) inner(M, x, vⱼ, hessf_x_vⱼ) ν * norm_x(vⱼ)^2
        if norm_x(rⱼ) < ϵ_residual || inner(M, x, vⱼ, hessf_x_vⱼ) < ν * norm_x(vⱼ)^2 || j > maxiter
            ## Satisfying point obtained
            if j == 0
                dₖ = -gradfₖ
            end

            if norm_x(rⱼ) < ϵ_residual
                d_type = :Solved
                # printstyled("Exiting: ||rⱼ|| < ϵ : $(norm_x(rⱼ)) < $ϵ_residual\n", color=:red)
            elseif j > maxiter
                d_type = :MaxIter
                # printstyled("Exiting: j > maxiter : $j > $maxiter\n", color=:red)
            else
                d_type = :QuasiNegCurvature
                a = inner(M, x, vⱼ, hessf_x_vⱼ)
                b = ν * norm_x(vⱼ)^2
                # printstyled("Exiting: ⟨vⱼ, hessf_x(vⱼ)⟩ < ν * ||vⱼ||² : $a < $b\n", color=:red)
            end

            break
        end

        tⱼ = - inner(M, x, rⱼ, vⱼ) / inner(M, x, vⱼ, hessf_x_vⱼ)
        dₖ += tⱼ * vⱼ

        rⱼ_prev = deepcopy(rⱼ)
        vⱼ_prev = deepcopy(vⱼ)

        j += 1
    end

    return dₖ, j, d_type
end


text:
using LazySets, MathematicalSystems, Optim, Random, LinearAlgebra
import Distributions, Polyhedra, CDDLib, DifferentialEquations
using LazySets: dim, HalfSpace, sample, basetype

# use Reachability or ReachabilityAnalysis library depending on Julia version
const USE_REACHABILITY = (VERSION < v"1.3")
if USE_REACHABILITY
    using Reachability
else
    using ReachabilityAnalysis
end

include("AffDyn.jl")

include("julia_library.jl")

include("time_series_to_symbolic.jl")

include("analyze_sequence.jl")

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth