text:
export dcs_parallel
using DensestCommonSubgraph
using Devectorize

function reduce_dcs(sA::(Array{Int64, 1}, Float64), sB::(Array{Int64, 1}, Float64))
  if sA[2] > sB[2]
    return sA
  else
    return sB
  end
end

function dcs_parallel{T}(G::Array{SparseMatrixCSC{T, Int64}, 1}, n::Number)
  @parallel reduce_dcs for j=[1:n;]
    dcs_exhaustive_worker(G, n, j)
  end
end

function dcs_exhaustive_worker{T}(G::Array{SparseMatrixCSC{T, Int64}, 1}, num_parts::Number, curr_part::Number)
  m = length(G)
  n = size(G[1], 1)

  prev_state = zeros(Bool, n)
  prev_nodes = zeros(Int64, n)
  prev_num_nodes = 0

  binary_state = zeros(Bool, n)

  curr_state = zeros(Bool, n)
  curr_nodes = zeros(Int64, n)
  curr_num_nodes = 0

  opt_idx = 0
  opt_deg = 0.0

  curr_ne = zeros(m, 1)
  curr_deg = zeros(m, 1)

  # Which part of state space to look at
  if VERSION.minor == 3
    j_begin = ifloor(2^(n)/num_parts*(curr_part - 1)) + 1
    j_end = min(iceil(2^(n)/num_parts*(curr_part)), 2^n - 1)
  else # avoids deprecated warning!
    j_begin = floor(Integer, 2^(n)/num_parts*(curr_part - 1)) + 1
    j_end = min(ceil(Integer, 2^(n)/num_parts*(curr_part)), 2^n - 1)
  end
  # update previous state
  if j_begin > 1
    int_to_binary_bool!(binary_state, j_begin - 1, n)
    binary_to_gray!(prev_state, binary_state)
    prev_num_nodes = find_set_bits!(prev_nodes, prev_state)
    for k=1:m
      for i=1:prev_num_nodes
        for j=(i+1):prev_num_nodes
          ni = prev_nodes[i]
          nj = prev_nodes[j]
          curr_ne[k] = curr_ne[k] + G[k][ni, nj] + G[k][nj, ni]
        end
      end
      curr_deg[k] = curr_ne[k]/prev_num_nodes
    end
    opt_deg = minimum(curr_deg)
    opt_idx = j_begin - 1
    # println("ne: $(curr_ne') deg: $(curr_deg') nodes: $(prev_nodes[1:prev_num_nodes])")
  end

  # start the main loop
  for j=j_begin:j_end
    int_to_binary_bool!(binary_state, j, n)
    binary_to_gray!(curr_state, binary_state)
    curr_num_nodes = find_set_bits!(curr_nodes, curr_state)

    # Find which node has been added or deleted
    changed_idx = 1
    while (curr_state[changed_idx] == prev_state[changed_idx])
      changed_idx = changed_idx + 1
    end

    # Has the node been added or deleted?
    if curr_state[changed_idx]
      m_sign = 1
    else
      m_sign = -1
    end
    @inbounds begin
      for k=1:m
        delta_deg = 0
        if prev_num_nodes > 0
          for ll=1:prev_num_nodes
            other_node = prev_nodes[ll]
            delta_deg = delta_deg + G[k][other_node, changed_idx] + G[k][changed_idx, other_node]
          end
        end
        curr_ne[k] = curr_ne[k] + m_sign  * delta_deg
        curr_deg[k] = curr_ne[k] / curr_num_nodes
        #      println("G[$k] j: $j: $(curr_ne[k]), $(curr_deg[k])")
      end
    end
    @devec curr_min_deg = minimum(curr_deg)
    if curr_min_deg > opt_deg
      opt_deg = curr_min_deg
      opt_idx = j
    end

    for k=[1:n;]
      prev_state[k] = curr_state[k]
      prev_nodes[k] = curr_nodes[k]
    end
    prev_num_nodes = curr_num_nodes
  end
  int_to_binary_bool!(binary_state, opt_idx, n)
  binary_to_gray!(curr_state, binary_state)
  curr_num_nodes = find_set_bits!(curr_nodes, curr_state)
  opt_S = curr_nodes[1:curr_num_nodes]
  # println("$(lpad(curr_part, 3))/$(rpad(num_parts, 3)) min_deg: $(rpad(opt_deg, 25)) S: $opt_S")

  return opt_S , opt_deg
end


text:
using Documenter
using fug

push!(LOAD_PATH,"../src/")
makedocs(
    sitename = "Fenske-Underwood-Gilliand Documation",
    pages = ["Index" => "index.md"],
    format = Documenter.HTML(prettyurls = false),
    modules = [fug]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "https://github.com/jmox0351/fug.git",
    devbranch = "gh-pages"
)

ext: jl


text:
# PolyJuMP Data
mutable struct Data
    # Default set for NonNegPoly
    nonnegpoly_default
    # Default set for PosDefPolyMatrix
    posdefpolymatrix_default
    function Data()
        new(nothing, nothing)
    end
end

function getpolydata(m::JuMP.Model)
    if !haskey(m.ext, :Poly)
        m.ext[:Poly] = Data()
    end
    m.ext[:Poly]
end


text:
# # Optimal transport

#=
In this example, we show how to differentiate through the solution of the entropy-regularized optimal transport problem.
=#

using Distances
using FiniteDiff
using ImplicitDifferentiation
using OptimalTransport
using Zygote
using Test, LinearAlgebra #src

#=
## Introduction

Here we give a brief introduction to optimal transport, see the [book by Gabriel Peyré and Marco Cuturi](https://optimaltransport.github.io/book/) for more details.

### Problem description

Suppose we have a distribution of mass ``\mu \in \Delta^{n}`` over points ``x_1, ..., x_{n} \in \mathbb{R}^d`` (where ``\Delta`` denotes the probability simplex).
We want to transport it to a distribution ``\nu \in \Delta^{m}`` over points ``y_1, ..., y_{m} \in \mathbb{R}^d``.
The unit moving cost from ``x`` to ``y`` is proportional to the squared Euclidean distance ``c(x, y) = \lVert x - y \rVert_2^2``.

A transportation plan can be described by a coupling ``p = \Pi(\mu, \nu)``, i.e. a probability distribution on the product space with the right marginals:
```math
\Pi(\mu, \nu) = \{p \in \Delta^{n \times m}: p \mathbf{1} = \mu, p^\top \mathbf{1} = \nu\}
```
Let ``C \in \mathbb{R}^{n \times m}`` be the moving cost matrix, with ``C_{ij} = c(x_i, y_j)``.
The basic optimization problem we want to solve is a linear program:
```math
\hat{p}(C) = \min_{p \in \Pi(\mu, \nu)} \sum_{i,j} p_{ij} C_{ij}
```
In order to make it smoother, we add an entropic regularization term:
 ```math
\hat{p}_{\varepsilon}(C) = \min_{p \in \Pi(\mu, \nu)} \sum_{i,j} \left(p_{ij} C_{ij} + \varepsilon p_{ij} \log \frac{p_{ij}}{\mu_i \nu_j} \right)
```

### Sinkhorn algorithm

To solve the regularized problem, we can use the Sinkhorn fixed point algorithm.
Let ``K \in \mathbb{R}^{n \times m}`` be the matrix defined by ``K_{ij} = \exp(-C_{ij} / \varepsilon)``.
Then the optimal coupling ``\hat{p}_{\varepsilon}(C)`` can be written as:
```math
\hat{p}_{\varepsilon}(C) = \mathrm{diag}(\hat{u}) ~ K ~ \mathrm{diag}(\hat{v}) \tag{1}
```
where ``\hat{u}`` and ``\hat{v}`` are the fixed points of the following Sinkhorn fixed point iteration:
```math
u^{t+1} = \frac{\mu}{Kv^t} \qquad \text{and} \qquad v^{t+1} = \frac{\nu}{K^\top u^t}
```

The implicit function theorem can be used to differentiate ``\hat{u}`` and ``\hat{v}`` with respect to ``C``, ``\mu`` and/or ``\nu``.
This can be combined with automatic differentiation of Equation (1) to find the Jacobian
```math
J = \frac{\partial ~ \mathrm{vec}(\hat{p}_{\varepsilon}(C))}{\partial ~ \mathrm{vec}(C)}
```
=#

# ## Implicit function wrapper

# For now, `ImplicitFunction` objects do not take multiple arguments, so we use non-constant global variables instead (even though we shouldn't)

d = 10
n = 3
m = 4

X = rand(d, n)
Y = rand(d, m)

μ = fill(1 / n, n)
ν = fill(1 / m, m)
C_vec = vec(pairwise(SqEuclidean(), X, Y, dims=2))

ε = 1.0;

#=
We now embed the Sinkhorn algorithm and optimality conditions inside an `ImplicitFunction` struct.
For technical reasons related to optimality checking, our forward procedure returns ``\hat{u}`` instead of ``\hat{p}_\varepsilon``.
=#

function forward(C_vec)
    C = reshape(C_vec, n, m)
    solver = OptimalTransport.build_solver(μ, ν, C, ε, SinkhornGibbs())
    OptimalTransport.solve!(solver)
    û = solver.cache.u
    return û
end

function conditions(C_vec, û)
    C = reshape(C_vec, n, m)
    K = exp.(.-C ./ ε)
    v̂ = ν ./ (K' * û)
    return û .- μ ./ (K * v̂)
end

implicit = ImplicitFunction(forward, conditions);

# ## Testing

# First, let us check that the forward pass works correctly

û, _ = forward(C_vec)

maximum(abs, conditions(C_vec, û))

# Using the implicit function defined above, we can build an autodiff-compatible implementation of `transportation_plan` which does not require backpropagating through the Sinkhorn iterations:

function transportation_plan(C_vec)
    C = reshape(C_vec, n, m)
    K = exp.(.-C ./ ε)
    û = implicit(C_vec)
    v̂ = ν ./ (K' * û)
    p̂_vec = vec(û .* K .* v̂')
    return p̂_vec
end;

# Let us compare with the result obtained using finite differences:

J_autodiff = Zygote.jacobian(transportation_plan, C_vec)[1]
J_finitediff = FiniteDiff.finite_difference_jacobian(transportation_plan, C_vec)
maximum(abs, J_autodiff - J_finitediff)

# The following tests are not included in the docs.  #src

@test maximum(abs, J_autodiff - J_finitediff) < 1e-7  #src


text:
#=

  Subset sum problem in Julia ConstraintSolver.jl

  From Katta G. Murty: "Optimization Models for Decision Making", page 340
  http://ioe.engin.umich.edu/people/fac/books/murty/opti_model/junior-7.pdf

  """
  Example 7.8.1

  A bank van had several bags of coins, each containing either
  16, 17, 23, 24, 39, or 40 coins. While the van was parked on the
  street, thieves stole some bags. A total of 100 coins were lost.
  It is required to find how many bags were stolen.
  """

  And we also solve some larger random problem instances.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using Distributions
const CS = ConstraintSolver
include("constraints_utils.jl")

function subset_sum_model(a,total,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,
                                            "all_solutions"=> all_solutions,
                                            # "all_optimal_solutions"=>all_solutions,
                                            "logging"=>[],

                                            # "traverse_strategy"=>:BFS,
                                            # "traverse_strategy"=>:DFS,
                                            "traverse_strategy"=>:DBFS,

                                            "branch_split"=>:Smallest,
                                            # "branch_split"=>:Biggest,
                                            # "branch_split"=>:InHalf,

                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                            "branch_strategy" => :IMPS, # default
                                            # "branch_strategy" => :ABS, # Activity Based Search
                                            # "activity.decay" => 0.999, # default 0.999
                                            # "activity.max_probes" => 10, # default, 10
                                            # "activity.max_confidence_deviation" => 20, # default 20

                                            # "simplify"=>false,
                                            # "simplify"=>true, # default

                                            "time_limit"=>6,

                                            # "backtrack" => false, # default true
                                            # "backtrack_sorting" => false, # default true

                                            # "lp_optimizer" => cbc_optimizer,
                                            "lp_optimizer" => glpk_optimizer,
                                            # "lp_optimizer" => ipopt_optimizer,
                                            ))

    if print_solutions
        println("a: $a")
        println("total: $total")
    end

    n = length(a)
    @variable(model, 0 <= x[1:n] <= n, Int)
    # @variable(model, 0 <= s <= n*maximum(a), Int)

    # scalar_product(model,x,a,total)

    @constraint(model,total == sum(a.*x))

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                # s_val = convert.(Integer,JuMP.value.(s; result=sol))
                println("x:$x_val")
                # println("s:$s_val")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

coins = [16, 17, 23, 24, 39, 40]
total = 100
@time subset_sum_model(coins,total,true,true)

#
# Some random instances
#
println("\nrandom instances")
for _ in 1:10
    # some random examples
    n = rand(DiscreteUniform(10,1_000))
    a = rand(DiscreteUniform(1,n),n)
    tot = rand(DiscreteUniform(1,3*sum(a)))
    println("\nn:$n total:$tot")
    # println("a:$a")
    @time subset_sum_model(a,tot,false,false)
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module ShaleDrillingLikelihood_InterpolationTest

using ShaleDrillingLikelihood
using Interpolations
using InteractiveUtils
using Test
using BenchmarkTools

using ShaleDrillingLikelihood:
    update_interpolation!, interpolation, scaled_interpolation,
    InPlaceInterp

using Interpolations: DimSpec # , scale

using Base: OneTo

const SHOW_CODE_WARNTYPE = false

# something to interpolate
A_x1 = 1 : 0.1 : 10
A_x2 = 1 : 0.5 : 20
A_x3 = OneTo(2)

f(x1, x2, x3) = log(x1+x2) + x3
A = [f(x1,x2,x3) for x1 in A_x1, x2 in A_x2, x3 in A_x3]

# for scaling
ranges = (A_x1, A_x2, A_x3)
splinetype(r::AbstractRange) = BSpline(Quadratic(Free(OnGrid())))
splinetype(r::AbstractUnitRange) = NoInterp()

# itpflag
it = splinetype.(ranges)
it isa DimSpec{BSpline} where {N}

# generate a test case
test_itp = interpolate(A, it)
test_sitp = scale(test_itp, ranges...)

# our object
inplace_itp = InPlaceInterp(A, it, ranges)
myitp = interpolation(inplace_itp)
mysitp = scaled_interpolation(inplace_itp)

# don't check out before updating coefs (prefilter)
@test !(myitp == test_itp)
@test !(mysitp == test_sitp)

# but aftewards, they do
update_interpolation!(inplace_itp)
@test myitp == test_itp
@test test_sitp == mysitp

if SHOW_CODE_WARNTYPE
    @code_warntype interpolation(inplace_itp)
    @code_warntype scaled_interpolation(inplace_itp)
    @code_warntype Interpolations.gradient(mysitp, 1.25, 1.25, 1)
end

vec = A_x2[2:end-1] .+ 0.1
jnk = zero(vec)

@inline function f1(grad, itp, vec)
    broadcast!(v -> itp(1,v,2), grad, vec)
end

@inline function f2(grad, itp, vec)
    grad .= itp(1,vec,2)
end

if SHOW_CODE_WARNTYPE
    @code_warntype f1(jnk, mysitp, vec)
    @code_warntype f2(jnk, mysitp, vec)
end

# timing to figure out most efficient way to udpate
@btime f1($jnk, $mysitp, $vec)
@btime f2($jnk, $mysitp, $vec)

# @code_typed optimize=true  f(jnk, mysitp, vec)

end


text:
module CyclicPackage

  using Boot

  include_folder(CyclicPackage, @__FILE__)

end


text:
function raster_plot!(ax, spikes::Array{Array{T,1},1}) where T
    for (trial,spike_train) in enumerate(spikes)
        ax.vlines(spike_train,trial,trial-1,
                   lw=0.3)
    end
end

function raster_hist(spikes)
    f,axs = subplots(2,1)
    raster_plot!(axs[1],spikes)
    axs[2].hist(reduce(vcat,spikes),bins=50)
end

function raster_hist(n::NeuronResponse,num_spikes,num_trials)
    spikes = [n.get_spikes(num_spikes) for _ in 1:num_trials]
    raster_hist(spikes)
end


text:
# using Ipaper
# using Test
# using PyPlot

@testset "Quantile" begin
  n = 20
  @test size(Quantile(rand(n, n, 30); dims=3)) == (n, n, 5)
  @test size(Quantile(rand(n, n, 30); dims=2)) == (n, 5, 30)
  @test size(Quantile(rand(n, n, 30); dims=1)) == (5, 20, 30)

  # Test for missing values
  @test Quantile([1, 2, missing, 4, 5, 10]) |> size == (5,)

  ## test for missing values
  # vector
  r1 = Quantile([missing, missing], [0.1, 0.2]; dims=1)
  @test all(isnan.(r1))

  # matrix
  x = rand(4, 4) |> to_missing
  x[1:2, :] .= missing
  r2 = Quantile(x, [0.1, 0.2]; dims=2)
  @test all(isnan.(r2[1:2, 1:2]))

  ## Quantile and nanquantile works
  x = rand(4, 4, 201)
  probs = [0.9, 0.99, 0.9999]

  r1 = Quantile(x, probs, dims=3)
  r2 = nanquantile(x, probs, dims=3)

  e_max = maximum(abs.(r1 - r2))
  @test e_max <= 1e-10
end


text:
mutable struct SVDDnegEps <: SVDDClassifier
    state::ModelState

    # model parameters
    C1::Float64
    C2::Float64
    kernel_fct::Kernel
    eps::Float64

    # training data
    data::Array{Float64,2}
    K::Array{Float64,2}
    adjust_K::Bool
    K_adjusted::Array{Float64,2}
    pools::Dict{Symbol, Array{Int64,1}}

    # fitted values
    alpha_values::Vector{Float64}
    const_term::Float64
    R::Float64

    function SVDDnegEps(data, pools; eps=0.0)
        m = new()
        m.C1 = 1.0
        m.C2 = 1.0
        m.eps = eps
        m.state = model_created
        m.data = data
        m.adjust_K = false
        m.pools = labelmap(pools)
        m.const_term = -Inf
        m.R = -Inf
        return m
    end
end

get_model_params(model::SVDDnegEps) = Dict(:C1 => model.C1, :C2 => model.C2, :eps => model.eps)

is_valid_param_value(model::SVDDnegEps, x::Type{Val{:C1}}, v) = 0 <= v <= 1
is_valid_param_value(model::SVDDnegEps, x::Type{Val{:C2}}, v) = is_valid_param_value(model, Val{:C1}, v)
is_valid_param_value(model::SVDDnegEps, x::Type{Val{:eps}}, v) = v >= 0

function set_C!(model::SVDDnegEps, C::Tuple{Number, Number})
    @assert 0 <= C[1] <= 1
    @assert 0 <= C[2] <= 1
    model.C1 = C[1]
    model.C2 = C[2]
    return nothing
end

set_C!(model::SVDDnegEps, C::Number) = set_C!(model, (C,C))

function set_eps!(model::SVDDnegEps, eps::Number)
    @assert eps >= 0
    model.eps = eps
    return nothing
end

function solve!(model::SVDDnegEps, solver::JuMP.OptimizerFactory)
    ULin = merge_pools(model.pools, :U, :Lin)
    length(ULin) > 0 || throw(ModelInvariantException("SVDDnegEps requires samples in pool :Lin or :U."))

    debug(LOGGER, "[SOLVE] Setting up QP for SVDDnegEps with $(is_K_adjusted(model) ? "adjusted" : "non-adjusted") kernel matrix.")
    QP = Model(solver)
    K = is_K_adjusted(model) ? model.K_adjusted : model.K

    @variable(QP, α[1:size(K,1)] >= 0)

    if haskey(model.pools, :Lout)
        @objective(QP, Max, sum(α[i] * model.eps for i in ULin) +
                             sum(α[i]*K[i,i] for i in ULin) -
                             sum(α[l]*K[l,l] for l in model.pools[:Lout]) -
                             sum(α[i]*α[j] * K[i,j] for i in ULin for j in ULin) +
                             2 * sum(α[l]*α[j] * K[l,j] for l in model.pools[:Lout] for j in ULin) -
                             sum(α[l]*α[m] * K[l,m] for l in model.pools[:Lout] for m in model.pools[:Lout]))

        @constraint(QP, sum(α[i] for i in ULin) - sum(α[l] for l in model.pools[:Lout]) == 1)
        @constraint(QP, α[ULin] .<= model.C1)
        @constraint(QP, α[model.pools[:Lout]] .<= model.C2)
    else # fall back to standard SVDD
        @objective(QP, Max, sum(α[i]*K[i,i] for i in ULin) -
                             sum(α[i]*α[j] * K[i,j] for i in ULin for j in ULin))
        @constraint(QP, sum(α) == 1)
        @constraint(QP, α[ULin] .<= model.C1)
    end
    debug(LOGGER, "[SOLVE] Solving QP with $(typeof(solver))...")
    JuMP.optimize!(QP)
    status = JuMP.termination_status(QP)
    debug(LOGGER, "[SOLVE] Finished with status: $(status).")
    model.alpha_values = JuMP.value.(α)
    return status
end

function get_alpha_prime(model::SVDDnegEps)
    if haskey(model.pools, :Lout)
        alpha = copy(model.alpha_values)
        alpha[model.pools[:Lout]] = - alpha[model.pools[:Lout]]
        return alpha
    else
        return model.alpha_values
    end
end

function get_support_vectors(model::SVDDnegEps)
    ULin = merge_pools(model.pools, :U, :Lin)
    length(ULin) > 0 || throw(ModelInvariantException("SVDDnegEps requires samples in pool :Lin or :U."))
    sv = filter!(x -> x in ULin, findall((model.alpha_values .> OPT_PRECISION) .& (model.alpha_values .< (model.C1 - OPT_PRECISION))))
    if haskey(model.pools, :Lout)
        sv = append!(sv, filter!(x -> x in model.pools[:Lout],
                                 findall((model.alpha_values .> OPT_PRECISION) .& (model.alpha_values .< (model.C2 - OPT_PRECISION)))))
    end
    return sv
end

function get_R_and_const_term(model::SVDDnegEps)
    @assert length(model.alpha_values) == size(model.K, 1)

    sv_idx = get_support_vectors(model)
    α = get_alpha_prime(model)

    const_term = sum(α[i] * α[j]*model.K[i,j] for i in eachindex(α) for j in eachindex(α))
    R_squared = 0.0
    tmp = 0.0
    # take the maximum R over all support vectors for numerical reasons
    for k in sv_idx
        tmp = model.K[k,k] - 2 * sum(α[i] * model.K[i,k] for i in eachindex(α)) + const_term + model.eps
        R_squared = tmp > R_squared ? tmp : R_squared
    end
    return (sqrt(R_squared), const_term)
end


text:
import Base: show
using Printf

# show(io::IO, x): One line show used for show, print, string etc.
# show(io::IO, ::MIME"text/plain", x): Multi line show used for display

# show with indent:
sshow(io::IO, x, indent) = print(io, repeat(' ', indent), x, '\n')

show(io::IO, ::MIME"text/plain", s::Block) = show(io, s)
show(io::IO, s::Block) = sshow(io, s, 0)
function sshow(io::IO, s::Block, indent::Int)
    sshow(io, (s.name === nothing ? "Block" : "$(s.name)"), indent)
    for l in s.layers
        sshow(io, l, indent+2)
    end
    #sshow(io, "end", indent)
end

show(io::IO, ::MIME"text/plain", s::Add) = show(io, s)
show(io::IO, s::Add) = sshow(io, s, 0)
function sshow(io::IO, r::Add, indent::Int)
    a = r.activation === nothing ? "" : "($(r.activation))"
    sshow(io, "Add$a", indent)
    sshow(io, r.blocks[1], indent+2)
    for i in 2:length(r.blocks)
        if r.blocks[i] == identity
            sshow(io, ".+ identity", indent+2)
        else
            sshow(io, ".+", indent+2)
            sshow(io, r.blocks[i], indent+2)
        end
    end
end

show(io::IO, ::MIME"text/plain", s::Mul) = show(io, s)
show(io::IO, s::Mul) = sshow(io, s, 0)
function sshow(io::IO, r::Mul, indent::Int)
    a = r.activation === nothing ? "" : "($(r.activation))"
    sshow(io, "Mul$a", indent)
    sshow(io, r.blocks[1], indent+2)
    for i in 2:length(r.blocks)
        if r.blocks[i] == identity
            sshow(io, ".* identity", indent+2)
        else
            sshow(io, ".*", indent+2)
            sshow(io, r.blocks[i], indent+2)
        end
    end
end

function sshow(io::IO, f::Function, indent::Int)
    name = string(f)
    if name[1] === '#'          # handle anonymous functions
        m = first(methods(f))
        lines = readlines(string(m.file))
        src = strip(lines[m.line], [' ', '\t', ',', ';'])
        sshow(io, src, indent)
    else
        sshow(io, name, indent)
    end
end

show(io::IO, ::MIME"text/plain", b::BatchNorm) = show(io, b)
show(io::IO, b::BatchNorm) = print(io, BatchNorm)

show(io::IO, ::MIME"text/plain", c::Conv) = show(io, c)
function show(io::IO, c::Conv)
    @printf(io, "Conv(%d×%d, %d⇒%d", c.wdims...)
    # Only print non-defaults
    if c.w !== nothing && eltype(c.w) !== Float32; print(io, ", ", eltype(c.w)); end
    if any(c.padding .!= 0); print(io, ", padding=$(c.padding)"); end
    if any(c.stride .!= 1); print(io, ", stride=$(c.stride)"); end
    if any(c.dilation .!= 1); print(io, ", dilation=$(c.dilation)"); end
    if c.groups != 1; print(io, ", groups=$(c.groups)"); end
    if !c.crosscorrelation; print(io, ", flipkernel"); end
    if c.channelmajor; print(io, ", channelmajor"); end
    if c.alpha != 1; print(io, ", alpha=$(c.alpha)"); end
    if c.beta != 0; print(io, ", beta=$(c.beta)"); end
    if c.bias !== nothing || c.binit !== nothing; print(io, ", bias"); end
    if c.normalization !== nothing; print(io, ", $(c.normalization)"); end
    if c.activation !== nothing; print(io, ", $(c.activation)"); end
    print(io, ")")
end

show(io::IO, ::MIME"text/plain", d::Linear) = show(io, d)
function show(io::IO, d::Linear)
    print(io, "Linear($(d.inputsize)=>$(d.outputsize)")
    if d.w !== nothing && eltype(d.w) !== Float32; print(io, ", ", eltype(d.w)); end
    if d.bias !== nothing || d.binit !== nothing; print(io, ", bias"); end
    if d.activation !== nothing; print(io, ", $(d.activation)"); end
    if d.dropout !== 0; print(io, ", dropout=$(d.dropout)"); end
    print(io, ")")
end

show(io::IO, ::MIME"text/plain", o::Op) = show(io, o)
function show(io::IO, o::Op)
    print(io, o.f)
    !isempty(o.args) && print(io, o.args)
    !isempty(o.kwargs) && print(io, (; o.kwargs...))
end

show(io::IO, ::MIME"text/plain", z::ZeroPad) = show(io, z)
show(io::IO, z::ZeroPad) = print(io, "ZeroPad$(z.padding)")


text:
#***********************************************************************
#  Copyright 2006 John A. Trangenstein
#
#  This software is made available for research and instructional use
#  only.
#  You may copy and use this software without charge for these
#  non-commercial purposes, provided that the copyright notice and
#  associated text is reproduced on all copies.
#  For all other uses (including distribution of modified versions),
#  please contact the author at
#    John A. Trangenstein
#    Department of Mathematics
#    Duke University
#    Durham, NC 27708-0320
#    USA
#  or
#    johnt@math.duke.edu
#
#  This software is made available "as is" without any assurance that it
#  is completely correct, or that it will work for your purposes.
#  Use the software at your own risk.
#***********************************************************************

function consdiff(fic::Int,lac::Int,fif::Int,laf::Int,fix::Int,lax::Int,
                  ifirst::Int,ilast::Int,
                  x::OffsetArray{Float64,1},
                  flux::OffsetArray{Float64,1},
                  conservd::OffsetArray{Float64,1})
#   ******************************************************************
#   update conservd to new time
#   ******************************************************************
    @inbounds for ic=ifirst:ilast
        conservd[ic] -= (flux[ic+1]-flux[ic]) / (x[ic+1]-x[ic])
    end
end

function stabledt(fc::Int,lc::Int,fm::Int,lm::Int,ifirst::Int,ilast::Int,
                  conserved::OffsetArray{Float64,1},
                  lambda_cell::OffsetArray{Float64,1},
                  x::OffsetArray{Float64,1})
#   ******************************************************************
#   compute stable timestep
#   ******************************************************************
    sdt=huge
    @inbounds for ic=ifirst:ilast
        abs_lambda=abs(lambda_cell[ic])
        dx = x[ic+1]-x[ic]
        if abs_lambda > roundoff*dx
            sdt=min(sdt,dx/abs_lambda)
        end
    end
    @assert sdt > 0.0
    sdt
end


text:
using FlexLinearAlgebra
using LinearAlgebra
using Test

v = FlexConvert([4,5,6])
@test v+v == 2v
@test v-v == FlexVector(1:3)
@test dot(v,v) == 77
v[2] = 4
v[3] = 4
@test v == 4*FlexOnes(1:3)
@test v[-1] == 0

x = collect(1:5)
v = FlexConvert(x)
@test Vector(v) == x
@test Set(keys(v)) == Set(1:5)
@test sum(values(v)) == sum(v)
@test keytype(v) == Int
@test valtype(0. * v) == Float64
@test length(v) == 5
@test haskey(v,1)
@test (-v) + v == 0v

v = FlexOnes(Int,1:6)
delete_entry!(v,5)
delete_entry!(v,5)
@test sum(v) == 5

M = FlexConvert(Matrix{Int}(I,3,3))
@test M[1,1]==1
S = Set(1:3)
@test Set(row_keys(M)) == S
@test Set(col_keys(M)) == S

A = [1 3; 4 5]
B = [-1.0 3; 5 -1.5]
v = [2; 3]
AA = FlexConvert(A)
BB = FlexConvert(B)
vv = FlexConvert(v)

@test valtype(BB) == Float64
S = Set(1:2)
@test Set(row_keys(BB)) == S
@test Set(col_keys(BB)) == S

@test Vector(AA*vv) == A*v
@test Matrix(AA*BB) == Matrix(A)*Matrix(B)
@test (AA*BB)*vv == AA*(BB*vv)
@test Set(row_keys(AA)) == Set(row_keys(BB))
AA[1,3] = 2
@test size(AA) == (2,3)

A = FlexMatrix{Int}(1:2,3:4)
A[1,4] = 5
@test sum(values(A)) == 5

A = FlexOnes(Int,1:5,1:5)
delete_row!(A,5)
delete_col!(A,5)
@test A == FlexOnes(Int,1:4,1:4)

A = FlexOnes(1:3,6:12) + im * FlexOnes(1:3,6:12)
@test Matrix(A') == Matrix(A)'

v = FlexOnes(Complex,1:4)
v[4] = 2im
@test dot(v,v) == (v'*v)[1]


text:
module MuscleSim
include("jules/Jules.jl")
js = Jules

allow_debug = false

function print_debug(message)
    if allow_debug
        print(message)
    end
end
type HillMuscleModel
    start_time::Float64
    end_time::Float64
    dt::Float64

    # Activation constants
    tau::Float64
    beta::Float64
    activ_lower_bound::Float64
    activ_upper_bound::Float64

    # Maximum muscle velocity and force
    V_max::Float64
    F_max::Float64
    L_optimal::Float64
    L_mt_initial::Float64

    # Tendon constants
    K_t::Float64
    L_st::Float64

    # function of time that describes the excitation
    excitation_func
end

# for attaching to compliant spring
type HillExternalModel
    L_load::Float64
    K_load::Float64
    L_total::Float64
end

type HillModelOutputs
    # Length and Velocity of the Muscle and Tendon
    L_mt::Vector{Float64} #
    V_mt::Vector{Float64} #

    # Force, normalized length, and Velocity of the muscle
    F_m::Vector{Float64} #
    F_mdot::Vector{Float64} #
    L_m::Vector{Float64} #
    V_m::Vector{Float64} #

    V_t::Vector{Float64} #
    L_t::Vector{Float64} #

    activation::Vector{Float64} #
    activation_dot::Vector{Float64}
    excitation::Vector{Float64} #

    time::Vector{Float64}
end

function CreateModel(;start_time = 0.0, end_time = 0.0, dt = 0.0, tau = 0.0, beta = 0.0, V_max = 0.0, F_max = 0.0,
                     L_optimal = 0.0, K_t = 0.0, K_sp = 0.0, L_st = 0.0, L_load = 0.0, L_mt_initial = 0.0,
                     L_total = 0.0, excitation_func = excite, activ_lower_bound = 0.1, activ_upper_bound = 1.0)

    model = HillMuscleModel(
        start_time,
        end_time,
        dt,

        tau,
        beta,
        activ_lower_bound,
        activ_upper_bound,

        V_max,
        F_max,
        L_optimal,
        L_mt_initial,

        K_t,
        L_st,
        excitation_func)

    external_model = HillExternalModel(
        L_load,
        K_sp,
        L_total)

    return (model, external_model)
end

#audot :: t -> (t -> b) -> b -> b
function audot(t, u, a; beta=BETA, tau=TAU)
    numerator = u(t) - (beta + (1 - beta) * u(t))*a
    return numerator/tau
end

function zero_excite(t::Float64)
    return 0.0
end

function excite(t::Float64)
    T0 = 1
    T1 = 6
    HIGH = 0.5
    LOW = 0.01
    ret = 0
    if t > T1
        return LOW
    end
    if t > T0
        return HIGH
    end
    return LOW
end

#norm_length_tension :: norm_length -> norm_force
function norm_length_tension(norm_length)
    A = -3.0508;
    B = 5.9758;
    C = -1.9597;
    if 0.42 <= norm_length <= 1.54
        return A.*(norm_length).^2 + B*norm_length + C;
    end

    return 0
end

#norm_inv_fv :: y -> norm_vel
function norm_inv_fv(y)
    if y > 0
        V_m = 0.995242 * exp(13.8817 * (y - 1.39214)) - 0.996815 * exp(-3.91442 * y)
        if V_m > 1.0
            V_m = 1.0
        end
    else
        V_m = -0.9968
    end
    return V_m
end

function calcL_t(model::HillMuscleModel, F_m)
    return model.L_st + F_m/model.K_t
end

function interp_activation(model::HillMuscleModel, t::Float64)
    return js.interp(model.time[1:length(model.activation)], model.activation, t);
end

function interp_length(model::HillMuscleModel, t::Float64)
    return js.interp(model.time, model.L_mt, t);
end

function interp_velocity(model::HillMuscleModel, t::Float64)
    return js.interp(model.time, model.V_mt, t);
end

function calcF_mdot(model::HillMuscleModel, t, F_m, V_mt)
    L_t = model.L_st + F_m/model.K_t
    L_m = (model.L_total - L_t - F_m/model.K_sp)
    L_mt = L_m + L_t

    norm_L_m = L_m/model.L_max

    y =
        F_m/(model.F_max * interp_activation(model, t) * norm_length_tension(norm_L_m))

    V_m = model.V_max*norm_inv_fv(y)

    V_t = V_mt - V_m

    #=
    print_debug("F_m $F_m\n")
    print_debug("numerator ")
    print_debug(interp_length(model, t) - F_m/model.K_t - model.L_st)
    print_debug("\n")
    print_debug("norm Lm $norm_L_m\n")
    print_debug("y $y\n")
    print_debug("norm_inf_fvy $(norm_inv_fv(y))\n")
    print_debug("V_m $V_m\n")
    print_debug("V_t $V_t\n")
    print_debug("\n")
    =#

    return F_mdot = model.K_t*V_t
end

# Model params -> activation -> muscle force -> mucle length -> muscle velocity
function calcMuscleVelocity(model::HillMuscleModel, activation, F_m, L_m)
    y = F_m/(model.F_max * activation * norm_length_tension(L_m/model.L_optimal))
    print_debug("calcMuscleVelocity ")
    print_debug("y $y ")
    print_debug("activation $activation ")
    print_debug("F_m $F_m ")
    print_debug("L_m $L_m\n")
    V_m = model.V_max*norm_inv_fv(y)
    return V_m
end

function calcFmdotnew(model::HillMuscleModel, V_m, V_mt)
    V_t = V_mt - V_m
    F_mdot = V_t*model.K_t

    print_debug("calcfmdot ")
    print_debug("V_m $V_m ")
    print_debug("V_mt $V_mt ")
    print_debug("V_t $V_t ")
    print_debug("F_mdot $F_mdot\n")
    return F_mdot
end

function calcV_mtExternal(external_model::HillExternalModel, model::HillMuscleModel, F_mdot)
    V_mt = -F_mdot/external_model.K_load # TODO: SHOULD THIS BE NEGATIVE?
    return V_mt
end

function calcL_mtExternal(external_model::HillExternalModel, model::HillMuscleModel, F_m)
    L_mt = external_model.L_total - (F_m/external_model.K_load + external_model.L_load)
    if (F_m/external_model.K_load + external_model.L_load) > external_model.L_total
        print("$F_m\n")
    end
    return L_mt
end

function gen_activation!(model::HillMuscleModel)
    model.excitation = map(model.excitation_func, model.time)
    a_dot = # Create anonymous function for activation.
        (t, a) ->
            audot(t, model.excitation_func, a, tau=model.tau, beta=model.beta)

    model.activation = # integrate to form activation
        foldl((acc, t) ->
            vcat(acc, js.RK4(acc[end], a_dot, t, model.dt)),
                model.excitation_func(model.time[1]), model.time) # integrate

    model.activation =
        collect(Iterators.take(model.activation, length(model.time)))
end

function simulate(model::HillMuscleModel, external_model::HillExternalModel)

    #F_m_dot = (t, F_m) -> calcF_mdot(model, t, F_m)

    #=
    model.F_m =
        foldl((acc, t) ->
            vcat(acc, js.RK4(acc[end], F_m_dot, t, model.dt)), model.time)
        =#

    initial_value_outputs = HillModelOutputs(
        zeros(length_time), #L_mt
        zeros(length_time), #V_mt

        zeros(length_time), #F_m
        zeros(length_time), #F_mdot
        zeros(length_time), #L_m
        zeros(length_time), #V_m

        zeros(length_time), #V_t
        zeros(length_time), #L_t

        zeros(length_time), #activation
        zeros(length_time), #activation_dot
        zeros(length_time), #excitation

        zeros(length_time) #time
        )

    initial_value_outputs.L_mt[1] = model.L_mt_init

    # take initial stab at calculating F_m
    temp_force = collect(0:0.01:2)
    temp_vel = map(norm_inv_fv, temp_force)
    FV = js.interp(temp_vel, temp_force, 0.0)
    initial_value_outputs.activation[1] = model.activ_lower_bound

    initial_value_outputsmodel.F_m[1] =
        model.F_max * activation * norm_length_tension((initial_value_outputsmodel.L_mt[1] - model.L_st)/model.L_max)*FV

    model.V_mt[1] = 0.0

    for i in 1:length(model.time)-1
    initial_value_outputs.F_m[i+1] =
        js.RK4(initial_value_outputs.F_m[i], (t, F_m) -> calcF_mdot(model, t, F_m, model.V_mt[i]), model.time[i], model.dt)

    model.V_mt[i+1] = calcF_mdot(model, model.time[i], model.F_m[i], model.V_mt[i])/model.K_sp

    print_debug(model.F_m[i])
    print_debug("\n")

    L_t = model.L_st + model.F_m[i]/model.K_t
    L_m = (model.L_total - L_t - model.F_m[i]/model.K_sp)
    model.L_mt[i+1] = L_m + L_t
    end
end

function simulateExternal(model::HillMuscleModel, external_model::HillExternalModel)

length_time = length(model.start_time:model.dt:model.end_time)

initial_value_outputs = HillModelOutputs(
    zeros(length_time), #L_mt
    zeros(length_time), #V_mt

    zeros(length_time), #F_m
    zeros(length_time), #F_mdot
    zeros(length_time), #L_m
    zeros(length_time), #V_m

    zeros(length_time), #V_t
    zeros(length_time), #L_t

    zeros(length_time), #activation
    zeros(length_time), #activation_dot
    zeros(length_time), #excitation

    zeros(length_time) #time
    )

    temp_excitation = model.excitation_func
    model.excitation_func = temp_excitation # t -> 0.1

    initial_value_outputs.L_mt[1] = model.L_mt_initial
    initial_value_outputs.L_t[1] = model.L_st
    initial_value_outputs.L_m[1] = model.L_mt_initial - model.L_st
    initial_value_outputs.activation[1] = js.constrain(model.excitation_func(model.start_time), model.activ_lower_bound, model.activ_upper_bound)

    # take initial stab at calculating F_m
    temp_force = collect(0:0.01:2)
    temp_vel = map(norm_inv_fv, temp_force)
    FV = js.interp(temp_vel, temp_force, 0.0)

    initial_value_outputs.F_m[1] =
        model.F_max * initial_value_outputs.activation[1] *
            norm_length_tension((initial_value_outputs.L_m[1])/model.L_optimal)*FV

    print_debug("F_m $(initial_value_outputs.F_m[1])\n")
    print_debug("L_t $(initial_value_outputs.L_t[1])\n")
    print_debug("L_m $(initial_value_outputs.L_m[1])\n")
    print_debug("L_mt $(initial_value_outputs.L_mt[1])\n")

# with low activation, simulate for a while to get initial F_m
    loopSimulationAB(model, external_model, initial_value_outputs)

#=
# the main simulation
    model.excitation_func = temp_excitation

outputs = HillModelOutputs(
    zeros(length_time), #L_mt
    zeros(length_time), #V_mt

    zeros(length_time), #F_m
    zeros(length_time), #F_mdot
    zeros(length_time), #L_m
    zeros(length_time), #V_m

    zeros(length_time), #V_t
    zeros(length_time), #L_t

    zeros(length_time), #activation
    zeros(length_time), #excitation

    zeros(length_time) #time
    )

    outputs.activation[1] = model.activ_lower_bound

# loop over simulation, then allow to settle
    loopSimulation(model, external_model, outputs)
    =#

return initial_value_outputs
end

function loopSimulation(model::HillMuscleModel, external_model::HillExternalModel, outputs::HillModelOutputs)
    time = model.start_time
    iteration = 1

    while time < model.end_time - model.dt

# calculate iteration + 1 values
        outputs.time[iteration] = time
        simulateStep(model, external_model, outputs, iteration, time)

        time += model.dt
        iteration += 1
    end
end

function loopSimulationAB(model::HillMuscleModel, external_model::HillExternalModel, outputs::HillModelOutputs)
    time = model.start_time
    iteration = 1

    for i in 1:4
        outputs.time[iteration] = time
        simulateStep(model, external_model, outputs, iteration, time)

        time += model.dt
        iteration += 1
    end

    while time < model.end_time - model.dt

# calculate iteration + 1 values
        outputs.time[iteration] = time
        simulateStepAB(model, external_model, outputs, iteration, time)

        time += model.dt
        iteration += 1
    end
end

function calcActivationAB(model::HillMuscleModel, previous_value, previous_adot, t)
#### based on excitation, calculate the next activation
# calculate current excitation

    a_dot = # Create anonymous function for activation.
        (t, a) ->
            audot(t, model.excitation_func, a, tau=model.tau, beta=model.beta)

    activation = # integrate to form activation
            js.adams_bashforth_moulton(previous_value, previous_adot, a_dot, t, model.dt)

    activation = js.constrain(activation, model.activ_lower_bound, model.activ_upper_bound)

    return (activation, a_dot(t, previous_value))
end

function simulateStepAB(model::HillMuscleModel, external_model::HillExternalModel,
    outputs::HillModelOutputs, iteration, time)

#### based on excitation, calculate the next activation
# calculate current excitation

    outputs.excitation[iteration] = model.excitation_func(time)

    if iteration == 1
        (outputs.activation[iteration + 1], outputs.activation_dot[iteration + 1]) =
        calcActivation(model,
            js.constrain(model.excitation_func(time), model.activ_lower_bound, model.activ_upper_bound), time)
    else
        (outputs.activation[iteration + 1], outputs.activation_dot[iteration + 1]) =
        calcActivation(model, outputs.activation[iteration], time)
    end

#=
 use the activation to calculate the next muscle model state
 this includes:
 F_m = F_t -
 F_mdot = F_tdot -
 V_m -
 V_t -
 L_t -
 L_m -
 =#

    print_debug("\nactual calc ")
    outputs.V_m[iteration + 1] =
        calcMuscleVelocity(model,
            outputs.activation[iteration], outputs.F_m[iteration], outputs.L_m[iteration])

    outputs.F_mdot[iteration + 1] = calcFmdotnew(model, outputs.V_m[iteration], outputs.V_mt[iteration])

    outputs.V_t[iteration + 1] = outputs.V_mt[iteration] - outputs.V_m[iteration]

    slope = (outputs.activation[iteration + 1] - outputs.activation[iteration])/model.dt
    interpolate_activation = t -> slope*t + outputs.activation[iteration] - slope*time

    print_debug("\nintegral calc ")
    outputs.F_m[iteration + 1] = js.adams_bashforth_moulton(
        outputs.F_m[iteration],
        outputs.F_mdot,
        (t, u) -> calcFmdotnew(
                    model,
                    calcMuscleVelocity(
                        model,
                        interpolate_activation(t),
                        u,
                        outputs.L_m[iteration]),
                    outputs.V_mt[iteration]),
        time,
        model.dt)

    outputs.F_m[iteration + 1] = js.constrain(outputs.F_m[iteration + 1], 0.0, model.F_max) #outputs.F_m[iteration + 1])

    F_m = outputs.F_m[iteration + 1]
    print_debug("\nFinal Fm $F_m \n")

    calcLengths(model, external_model, outputs, iteration, time)

#=
 interact with the outside world, calculating:
 L_mt
 V_mt
=#
    calcExternal(model, external_model, outputs, iteration, time)
end

function calcActivation(model::HillMuscleModel, previous_value, t)
#### based on excitation, calculate the next activation
# calculate current excitation

    a_dot = # Create anonymous function for activation.
        (t, a) ->
            audot(t, model.excitation_func, a, tau=model.tau, beta=model.beta)

    activation = # integrate to form activation
            js.RK4(previous_value, a_dot, t, model.dt)

    activation = js.constrain(activation, model.activ_lower_bound, model.activ_upper_bound)

    return (activation, a_dot(t, previous_value))
end

function simulateStep(model::HillMuscleModel, external_model::HillExternalModel,
    outputs::HillModelOutputs, iteration, time)

#### based on excitation, calculate the next activation
# calculate current excitation

    outputs.excitation[iteration] = model.excitation_func(time)

    if iteration == 1
        (outputs.activation[iteration + 1], outputs.activation_dot[iteration + 1]) =
            calcActivation(model,
                js.constrain(model.excitation_func(time), model.activ_lower_bound, model.activ_upper_bound), time)
    else
        (outputs.activation[iteration + 1], outputs.activation_dot[iteration + 1]) =
            calcActivation(model, outputs.activation[iteration], time)
    end

#=
 use the activation to calculate the next muscle model state
 this includes:
 F_m = F_t -
 F_mdot = F_tdot -
 V_m -
 V_t -
 L_t -
 L_m -
 =#

    print_debug("\nactual calc ")
    outputs.V_m[iteration + 1] =
        calcMuscleVelocity(model,
            outputs.activation[iteration], outputs.F_m[iteration], outputs.L_m[iteration])

    outputs.F_mdot[iteration + 1] = calcFmdotnew(model, outputs.V_m[iteration], outputs.V_mt[iteration])

    outputs.V_t[iteration + 1] = outputs.V_mt[iteration] - outputs.V_m[iteration]

    slope = (outputs.activation[iteration + 1] - outputs.activation[iteration])/model.dt
    interpolate_activation = t -> slope*t + outputs.activation[iteration] - slope*time

    print_debug("\nintegral calc ")
    outputs.F_m[iteration + 1] = js.RK4(
        outputs.F_m[iteration],
        (t, u) -> calcFmdotnew(
                    model,
                    calcMuscleVelocity(
                        model,
                        interpolate_activation(t),
                        u,
                        outputs.L_m[iteration]),
                    outputs.V_mt[iteration]),
        time,
        model.dt)

    outputs.F_m[iteration + 1] = js.constrain(outputs.F_m[iteration + 1], 0.0, outputs.F_m[iteration + 1])
    F_m = outputs.F_m[iteration + 1]
    print_debug("\nFinal Fm $F_m \n")

    calcLengths(model, external_model, outputs, iteration, time)
#=
 interact with the outside world, calculating:
 L_mt
 V_mt
=#
    calcExternal(model, external_model, outputs, iteration, time)
end

function calcLengths(model, external_model, outputs, iteration, time)

    outputs.L_t[iteration + 1] = calcL_t(model, outputs.F_m[iteration])
    Keq = external_model.K_load * model.K_t/(external_model.K_load + model.K_t)
    outputs.L_m[iteration + 1] =
        outputs.L_mt[iteration] - outputs.L_t[iteration]
        #external_model.L_total - (Keq * outputs.F_m[iteration] + model.L_st + external_model.L_load)

end

function calcExternal(model, external_model, outputs, iteration, time)

    outputs.V_mt[iteration + 1] =
        #0
        calcV_mtExternal(external_model, model, outputs.F_mdot[iteration])
        #calcvmt(model, time)

    outputs.L_mt[iteration + 1] =
        #model.L_mt_initial
        calcL_mtExternal(external_model, model, outputs.F_m[iteration])
        #calclmt(model, time)

end

delay = 2
end_ramp = 4
slope = 0.5*10.0^-2

function calclmt(model, time)
    if time < delay
        return model.L_mt_initial
    end

    if time < end_ramp
        return model.L_mt_initial +  (time - delay)*slope
    end

    return model.L_mt_initial + (end_ramp - delay) * slope
end

function calcvmt(model, time)
    if (time > delay) && (time < end_ramp)
        return slope
    end
    return 0.0
end

end


text:
using BioSequences, CSV, DataFrames

# Define custom function for nice imports
Base.parse(::Type{Vector{String}}, x::String) = Vector{String}(filter(x-> x != ", ", split(x, "\""))[2:end-1])
function Base.parse(::Type{Vector{Float64}}, x::String)
    number = split(split(x, "[")[end][1:end-1], ", ")
    number_list = Float64[]
    for num in number
        if num != ""
            push!(number_list, parse(Float64, num))
        else
            return push!(number_list, NaN)
        end
    end
    return number_list

end

function Base.parse(::Type{Tuple{Int64, Tuple{Int64, Int64}}}, x::T) where {T<:String}
    numbers = split(x, ", ")
    ind1 = parse(Int64, split(numbers[1], "(")[2])
    ind2 = parse(Int64, split(numbers[2], "(")[2])
    ind3 = parse(Int64, split(numbers[3], ")")[1])
    return (ind1, (ind2, ind3))
end

Base.parse(::Type{Vector{String}}, x::Missing) = String[]
Base.parse(::Type{Vector{Float64}}, x::Missing) = Float64[]

"""
    onehot_encoder(sequence::BioSequences.LongDNA)

Return One-Hot encoding of DNA sequence.
"""
function onehot_encoder(sequence::BioSequences.LongDNA)

    if occursin(ExactSearchQuery(dna"N"), sequence)
        throw(ArgumentError("Sequence cannot contain 'N'"))
    end

    encoder = Dict(DNA_A => 1, DNA_C => 2, DNA_G =>3, DNA_T => 4)
    encoded_seq = zeros(Int64, length(sequence), 4)
    for i in 1:length(sequence)
        encoded_seq[i, encoder[sequence[i]]] = 1
    end
    return encoded_seq
end

"""
    function eval_emat(
        sequence::BioSequences.LongDNA,
        emat;
        start_ind_seq::Int64=1,
        stop_ind_seq::Int64=length(sequence),
        start_ind_emat::Int64=1,
        stop_ind_emat::Int64=size(emat, 1)
    )

Evaluate energy matrix on given sequence.

# Parameters
- sequence: DNA sequence that is to be evaluated. Has to be either a string or BioSequences.LongDNA
- emat: Energy Matrix
- start_ind_seq: First index in the sequence to be evaluated. Default is 1.
- stop_ind_seq: Last index in the sequence to evaluated. Default is the last index of the sequence.
- start_index_seq: First row of energy matrix to be evaluated. Default is 1.
- stop_ind_seq: Last row of energy matrix to be evaluated. Default is the last row of the matrix.

"""
function eval_emat(
        sequence::Union{BioSequences.LongDNA, AbstractString},
        emat::Matrix;
        start_ind_seq::Int64=1,
        stop_ind_seq::Int64=length(sequence),
        start_ind_emat::Int64=1,
        stop_ind_emat::Int64=size(emat, 1)
    )

    if typeof(sequence) <: AbstractString
        sequence = BioSequences.LongDNA{4}(sequence)
    end

    if ismissing(stop_ind_emat)
        stop_ind_emat = size(emat)[1]
    end

    if (start_ind_seq < 0) || (start_ind_seq > stop_ind_seq) || (start_ind_seq >length(sequence))
        throw(ArgumentError("Starting index has to be within the sequence. Start Index: $start_ind_seq, Stop Index: $stop_ind_seq"))
    end

    if (start_ind_emat < 0) || (start_ind_emat > stop_ind_emat) || (start_ind_emat >size(emat)[1])
        throw(ArgumentError("Starting index has to be within the energy matrix. Start Index: $start_ind_emat, Stop Index: $stop_ind_emat"))
    end

    if (stop_ind_seq - start_ind_seq) != (stop_ind_emat-start_ind_emat)
        throw(ArgumentError("Length of sequence not equal to length of energy matrix"))
    end

    seq = sequence[start_ind_seq:stop_ind_seq]
    _emat = emat[start_ind_emat:stop_ind_emat, :]

    if size(_emat)[1] != length(seq)
        throw(ArgumentError("Energy matrix has to be the same length as sequence"))
    elseif size(_emat)[2] != 4
        throw(ArgumentError("Second dimension has to be of size 4 for energy matrix. Is $(size(emat)[2])"))
    end

    onehot_encoder(seq::BioSequences.LongDNA)
    encoded = onehot_encoder(seq)

    return sum(encoded .* _emat)
end

"""
    import_twist_order(filename)

Import csv file for twist order. Fixes types in dataframe.
"""
function import_twist_order(filename)
    df = CSV.read(
        filename,
        types=Dict(
            "fwd_primer" => String,
            "rev_primer1" => String,
            "rev_primer2" => String,
            "rev_primer3" => String,
        ),
        DataFrame
    )
    df.fwd_primer = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.fwd_primer)
    df.rev_primer1 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer1)
    df.rev_primer2 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer2)
    df.rev_primer3 = parse.(Tuple{Int64, Tuple{Int64, Int64}}, df.rev_primer3)

    df.genes = parse.(Vector{String}, df.genes)
    df.sequence = [LongDNA{4}(seq) for seq in df.sequence]
    return df
end


text:
using StochasticDiffEq, SparseArrays, Test

f(du,u,p,t) = (du.=1.01u)
function g(du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))

sol = solve(prob,EM(),dt=1/1000)

f(du,u,p,t) = (du.=1.01u)
function g(du,u,p,t)
  du[1,1] = 0.3
  du[1,2] = 0.6
  du[1,3] = 0.9
  du[1,4] = 0.12
  du[2,1] = 1.2
  du[2,2] = 0.2
  du[2,3] = 0.3
  du[2,4] = 1.8
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))

sol = solve(prob,SRA1())
sol = solve(prob,SRA2())
sol = solve(prob,SRA3())
sol = solve(prob,SOSRA())
sol = solve(prob,SOSRA2())
sol = solve(prob,SRA())

@test length(sol.W[1]) == 4

function g(du,u,p,t)
  @test typeof(du) <: SparseMatrixCSC
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=sprand(2,4,1.0))

sol = solve(prob,EM(),dt=1/1000)

@test length(sol.W[1]) == 4


text:
module QuandlAPI

using HTTP
using JSON
using Dates
using TimeSeries

include("requests.jl")
end


text:
using OptiMode
using LinearAlgebra, Statistics
using StaticArrays
using GeometryPrimitives
using ChainRules
using Zygote #, ForwardDiff, FiniteDifferences
using Optim, Interpolations
using Zygote: @ignore, dropgrad
using Rotations
using StaticArrays: Dynamic
using Rotations: RotY, MRP
using FFTW
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

function E_relpower_xyz(ms::ModeSolver{ND,T},ω²H) where {ND,T<:Real}
    E = 1im * ε⁻¹_dot( fft( kx_tc( reshape(ω²H[2],(2,size(ms.grid)...)),mn(ms),ms.M̂.mag), (2:1+ND) ), flat( ms.M̂.ε⁻¹ ))
    Es = reinterpret(reshape, SVector{3,Complex{T}},  E)
    Pₑ_xyz_rel = normalize([mapreduce((ee,epss)->(abs2(ee[a])*inv(epss)[a,a]),+,Es,ms.M̂.ε⁻¹) for a=1:3],1)
    return Pₑ_xyz_rel
end

TE_filter = (ms,ω²H)->E_relpower_xyz(ms,ω²H)[1]>0.7
TM_filter = (ms,ω²H)->E_relpower_xyz(ms,ω²H)[2]>0.7
oddX_filter = (ms,αX)->sum(abs2,𝓟x̄(ms.grid)*αX[2])>0.7
evenX_filter = (ms,αX)->sum(abs2,𝓟x(ms.grid)*αX[2])>0.7

LNx = rotate(MgO_LiNbO₃,Matrix(MRP(RotY(π/2))))
Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
# Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 256, 256, 1;
grid = Grid(Δx,Δy,Nx,Ny)
# p_pe = [
#        1.7,                #   top ridge width         `w_top`         [μm]
#        0.7,                #   ridge thickness         `t_core`        [μm]
#        0.5,                #   top layer etch fraction `etch_frac`     [1]
#        π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
#                ];
p_pe = [
    1.85,        # 700 nm top width of angle sidewall ridge
    0.7,        # 600 nm MgO:LiNbO₃ ridge thickness
    3.4 / 7.0,    # etch fraction (they say etch depth of 500 nm, full thickness 600 nm)
    0.5236,      # 30° sidewall angle in radians (they call it 60° , in our params 0° is vertical)
]
rwg(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNx,SiO₂,Δx,Δy)
rwg_pe(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNx,SiO₂,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).

geom_pe = rwg_pe(p_pe)
# ms = ModeSolver(0.9, rwg_pe(p_pe), grid)
p_pe_lower = [0.8, 0.3, 0., 0.01]
p_pe_upper = [2.5, 2., 1., π/4.]

# λs = reverse(1.45:0.02:1.65)
# ωs = 1 ./ λs

nω_jank = 20
ωs_jank = collect(range(1/2.25,1/1.95,length=nω_jank)) # collect(0.416:0.01:0.527)
λs_jank = 1 ./ ωs_jank

# ωs = ωs_jank #collect(range(0.6,0.7,length=10))
# λs = λs_jank #inv.(ωs)

# n1F,ng1F = solve_n(ms,ωs,rwg_pe(p_pe)); n1S,ng1S = solve_n(ms,2*ωs,rwg_pe(p_pe))
##

function sum_Δng²_FHSH(ωs,p)
    ms = Zygote.@ignore(ModeSolver(kguess(1/2.25,rwg_pe(p)), rwg(p), Grid(6.0,4.0,128,128); nev=4)) #ModeSolver(1.45, rwg_pe(p), gr)
	nω = length(ωs)
    # ngs_FHSH = solve_n(
	# 	ms,
	# 	vcat(ωs, 2*ωs ),
	# 	rwg_pe(p);
	# 	f_filter=(ms,ω²H)->E_relpower_xyz(ms,ω²H)[1]>0.6)[2]
	# ms = Zygote.@ignore ModeSolver(kguess(ωs[1],rwg_pe(p)), rwg_pe(p), grid; nev=4)

	nFS,ngFS,gvdFS,EFS = solve_n(
	    ms,
	    vcat(ωs,2*ωs),
	    rwg_pe(p);
	    f_filter=(ms,ω²H)->E_relpower_xyz(ms,ω²H)[1]>0.6)

    ngs_FH = ngFS[1:nω]
	ngs_SH = ngFS[nω+1:2*nω]
    Δng² = abs2.(ngs_SH .- ngs_FH)
    sum(Δng²)
end

# sum_Δng_FHSH(ωs,p_pe)

# warmup
println("warmup function runs")
p0 = copy(p_pe)
@show sum_Δng²_FHSH(ωs_jank,p_pe)
# @show vng0, vng0_pb = Zygote.pullback(x->var_ng(ωs,x),p0)
# @show grad_vng0 = vng0_pb(1)

# define function that computes value and gradient of function `f` to be optimized
# according to https://julianlsolvers.github.io/Optim.jl/stable/#user/tipsandtricks/
function fg!(F,G,x)
    value, value_pb = Zygote.pullback(x) do x
       # var_ng(ωs,x)
	   nω = 20
	   ωs = collect(range(1/2.25,1/1.95,length=nω))
	   sum_Δng²_FHSH(ωs,x)
    end
    if G != nothing
        G .= value_pb(1)[1]
    end
    if F != nothing
        # F = value
        return value
    end
end

# G0 = [0.,0.,0.,0.]
# @show fg!(0.,G0,p0)
# println("G0 = $G0")
##
rand_p0() = p_pe_lower .+ [rand()*(p_pe_upper[i]-p_pe_lower[i]) for i=1:4]

opts =  Optim.Options(
                        outer_iterations = 4,
                        iterations = 6,
                        # time_limit = 3*3600,
                        store_trace = true,
                        show_trace = true,
                        show_every = 1,
                        extended_trace = true,
                        x_tol = 1e-4, # Absolute tolerance in changes of the input vector x, in infinity norm. Defaults to 0.0.
                        f_tol = 1e-5, # Relative tolerance in changes of the objective value. Defaults to 0.0.
                        g_tol = 1e-5, # Absolute tolerance in the gradient, in infinity norm. Defaults to 1e-8. For gradient free methods, this will control the main convergence tolerance, which is solver specific.
                    )

println("########################### Opt 1 ##########################")

# p_opt1 = 0.62929, 0.71422,0.7658459,0.125366
res1 = optimize( Optim.only_fg!(fg!),
                p_pe_lower,
                p_pe_upper,
                rand_p0(),
                Fminbox(Optim.BFGS()),
                opts,
            )

println("########################### Opt 2 ##########################")

# res2 = optimize( Optim.only_fg!(fg!),
#                 rand_p0(),
#                 Optim.BFGS(),
#                 opts,
#             )

res2 = optimize( Optim.only_fg!(fg!),
                p_pe_lower,
                p_pe_upper,
                rand_p0(),
                Fminbox(Optim.BFGS();mu0=1e-7),
                opts,
            )

println("########################### Opt 3 ##########################")

# res3 = optimize( Optim.only_fg!(fg!),
#                 rand_p0(),
#                 Optim.BFGS(),
#                 opts,
#             )

res3 = optimize( Optim.only_fg!(fg!),
                p_pe_lower,
                p_pe_upper,
                rand_p0(),
                Fminbox(Optim.BFGS();mu0=1e-6),
                opts,
            )


text:
@doc raw"""
    Triangle3

An object `t` of the type [`Triangle3`](@ref) is a triangle in the
three-dimensional Euclidean space ``Ε^2``.

As the triangle is not a full-dimensional object, there is only a test whether
a point lies on the triangle or not.
"""
Triangle3

"""
    Triangle3(p::Point3, q::Point3, r::Point3)

Introduces a triangle `t` with vertices `p`, `q` and `r`.
"""
Triangle3(p::Point3, q::Point3, r::Point3)

"""
    ==(t₁::Triangle3, t₂::Triangle3)

Test for equality: two triangles `t₁` and `t₂` are equal, iff there exists a
cyclic permutation of the vertices of `t₂`, such that  they are equal to the
vertices of `t₁`.
"""
==(t₁::Triangle3, t₂::Triangle3)

"""
    vertex(t::Triangle3, i::Integer)

Returns the `i`ᵗʰ vertex modulo 3 of `t`.
"""
vertex(t::Triangle3, i::Integer)

"""
    supporting_plane(t::Triangle3)

Returns the supporting plane of `t`, with same orientation.
"""
supporting_plane(t::Triangle3)

"""
    is_degenerate(t::Triangle3)

Triangle `t` is degenerate, if the vertices are collinear.
"""
is_degenerate(t::Triangle3)

"""
    squared_area(t::Triangle3)

Returns a square of the area of `t`.
"""
squared_area(t::Triangle3)

"""
    bbox(t::Triangle3)

Returns a bounding box containing `t`.
"""
bbox(t::Triangle3)

"""
    transform(t::Triangle3, at::AffTransformation3)

Returns the triangle obtained by applying `at` on the three vertices of `t`.
"""
transform(t::Triangle3, at::AffTransformation3)


text:
using GeometricBase
using GeometricBase: test_interface
using GeometricIntegrators.Solutions
using GeometricProblems.HarmonicOscillator
using Test

nt = 10
Δt = 0.1
ni = 5

t0 = 0.0
x0 = rand(2)
y0 = rand(2)
q0 = rand(1)
p0 = rand(1)
z0 = rand(3)
λ0 = rand(1)
μ0 = rand(2)

t1 = 1.0
x1 = [rand(2) for j = 1:ni]
y1 = [rand(2) for j = 1:ni]
q1 = [rand(1) for j = 1:ni]
p1 = [rand(1) for j = 1:ni]
z1 = [rand(3) for j = 1:ni]
λ1 = [rand(1) for j = 1:ni]
μ1 = [rand(2) for j = 1:ni]

t2 = t1 + (t1 - t0)
x2 = [rand(2) for j = 1:ni]
y2 = [rand(2) for j = 1:ni]
q2 = [rand(1) for j = 1:ni]
p2 = [rand(1) for j = 1:ni]
z2 = [rand(3) for j = 1:ni]
λ2 = [rand(1) for j = 1:ni]
μ2 = [rand(2) for j = 1:ni]

tx = zero(x0)
ty = zero(y0)
tq = zero(q0)
tp = zero(p0)
tz = zero(z0)
tλ = zero(λ0)
tμ = zero(μ0)

xs = [rand(2) for i = 1:nt]
ys = [rand(2) for i = 1:nt]
qs = [rand(1) for i = 1:nt]
ps = [rand(1) for i = 1:nt]
zs = [rand(3) for i = 1:nt]
λs = [rand(1) for i = 1:nt]
μs = [rand(2) for i = 1:nt]

Xs = [rand(2) for i = 1:nt, j = 1:ni]
Ys = [rand(2) for i = 1:nt, j = 1:ni]
Qs = [rand(1) for i = 1:nt, j = 1:ni]
Ps = [rand(1) for i = 1:nt, j = 1:ni]
Zs = [rand(3) for i = 1:nt, j = 1:ni]
Λs = [rand(1) for i = 1:nt, j = 1:ni]
Ms = [rand(2) for i = 1:nt, j = 1:ni]

ode = harmonic_oscillator_ode()
dae = harmonic_oscillator_dae()
pode = harmonic_oscillator_pode()
pdae = harmonic_oscillator_pdae()

@testset "$(rpad("ODE Solution",80))" begin
    asol = AtomicSolution(ode)

    # test constructors and general functionality
    sol = Solution(ode, Δt, nt)
    @test typeof(sol) <: SolutionODE

    # test_interface(sol) # TODO reactivate

    sol0 = Solution(similar(ode, x0), Δt, nt)
    @test typeof(sol0) <: SolutionODE

    sol1 = Solution(similar(ode, x1), Δt, nt)
    @test typeof(sol1) <: SolutionODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0)
    get_initial_conditions!(sol, tx, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test asol.t == t0
    @test asol.q == x0

    δt, δx = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0

    set_initial_conditions!(sol1, similar(ode, t1, x2))
    for j = 1:ni
        get_initial_conditions!(sol1, tx, j)
        @test tx == x2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(ode, x0), Δt, nt)
    sol2 = Solution(similar(ode, x0), Δt, nt)
    for i = 1:nt
        tx .= xs[i]
        asol.q .= xs[i]
        set_solution!(sol1, tx, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == xs
    @test sol2.q[1:nt] == xs

    sol1 = Solution(similar(ode, x1), Δt, nt)
    sol2 = Solution(similar(ode, x1), Δt, nt)
    for i = 1:nt
        for j = 1:ni
            tx .= Xs[i, j]
            asol.q .= Xs[i, j]
            set_solution!(sol1, tx, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt, :] == Xs
    @test sol2.q[1:nt, :] == Xs

    # test nsave and nwrite parameters
    sol = Solution(ode, Δt, 20, nsave = 2)
    @test sol.nt == 10

    sol = Solution(ode, Δt, 20, nsave = 2, nwrite = 10)
    @test sol.nt == 5

    # test reset
    sol = Solution(ode, Δt, nt)
    reset!(sol)
    @test sol.t[0] == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt
end

@testset "$(rpad("PODE Solution",80))" begin
    asol = AtomicSolution(pode)

    # test constructors and general functionality
    sol = Solution(pode, Δt, nt)
    @test typeof(sol) <: SolutionPODE

    # test_interface(sol) # TODO reactivate

    sol0 = Solution(similar(pode, q0, p0), Δt, nt)
    @test typeof(sol0) <: SolutionPODE

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    @test typeof(sol1) <: SolutionPODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, q0, p0)
    get_initial_conditions!(sol, tq, tp, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tq == q0
    @test tp == p0
    @test asol.t == t0
    @test asol.q == q0
    @test asol.p == p0

    δt, δq, δp = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δq == q0
    @test δp == p0

    set_initial_conditions!(sol1, similar(pode, t1, q2, p2))
    for j = 1:ni
        get_initial_conditions!(sol1, tq, tp, j)
        @test tq == q2[j]
        @test tp == p2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(pode, q0, p0), Δt, nt)
    sol2 = Solution(similar(pode, q0, p0), Δt, nt)
    for i = 1:nt
        tq .= qs[i]
        tp .= ps[i]
        asol.q .= qs[i]
        asol.p .= ps[i]
        set_solution!(sol1, tq, tp, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == qs
    @test sol1.p[1:nt] == ps
    @test sol2.q[1:nt] == qs
    @test sol2.p[1:nt] == ps

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    sol2 = Solution(similar(pode, q1, p1), Δt, nt)
    for i = 1:nt
        for j = 1:ni
            tq .= Qs[i, j]
            tp .= Ps[i, j]
            asol.q .= Qs[i, j]
            asol.p .= Ps[i, j]
            set_solution!(sol1, tq, tp, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt, :] == Qs
    @test sol1.p[1:nt, :] == Ps
    @test sol2.q[1:nt, :] == Qs
    @test sol2.p[1:nt, :] == Ps

    # test nsave and nwrite parameters
    sol = Solution(pode, Δt, 20, nsave = 2)
    @test sol.nt == 10

    sol = Solution(pode, Δt, 20, nsave = 2, nwrite = 10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pode, Δt, nt)
    reset!(sol)
    @test sol.t[0] == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt
end

@testset "$(rpad("DAE Solution",80))" begin
    asol = AtomicSolution(dae)

    # test constructors and general functionality
    sol = Solution(dae, Δt, nt)
    @test typeof(sol) <: SolutionDAE

    # test_interface(sol) # TODO reactivate

    sol0 = Solution(similar(dae, z0, λ0), Δt, nt)
    @test typeof(sol0) <: SolutionDAE

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    @test typeof(sol1) <: SolutionDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, z0, λ0)
    get_initial_conditions!(sol, tz, tλ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tz == z0
    @test tλ == λ0
    @test asol.t == t0
    @test asol.q == z0
    @test asol.λ == λ0

    δt, δz, δλ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δz == z0
    @test δλ == λ0

    set_initial_conditions!(sol1, similar(dae, t1, z2, λ2))
    for j = 1:ni
        get_initial_conditions!(sol1, tz, tλ, j)
        @test tz == z2[j]
        @test tλ == λ2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(dae, z0, λ0), Δt, nt)
    sol2 = Solution(similar(dae, z0, λ0), Δt, nt)
    for i = 1:nt
        tz .= zs[i]
        tλ .= λs[i]
        asol.q .= zs[i]
        asol.λ .= λs[i]
        set_solution!(sol1, tz, tλ, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == zs
    @test sol1.λ[1:nt] == λs
    @test sol2.q[1:nt] == zs
    @test sol2.λ[1:nt] == λs

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    sol2 = Solution(similar(dae, z1, λ1), Δt, nt)
    for i = 1:nt
        for j = 1:ni
            tz .= Zs[i, j]
            tλ .= Λs[i, j]
            asol.q .= Zs[i, j]
            asol.λ .= Λs[i, j]
            set_solution!(sol1, tz, tλ, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt, :] == Zs
    @test sol1.λ[1:nt, :] == Λs
    @test sol2.q[1:nt, :] == Zs
    @test sol2.λ[1:nt, :] == Λs

    # test nsave and nwrite parameters
    sol = Solution(dae, Δt, 20, nsave = 2)
    @test sol.nt == 10

    sol = Solution(dae, Δt, 20, nsave = 2, nwrite = 10)
    @test sol.nt == 5

    # test reset
    sol = Solution(dae, Δt, nt)
    reset!(sol)
    @test sol.t[0] == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt
end

@testset "$(rpad("PDAE Solution",80))" begin
    asol = AtomicSolution(pdae)

    # test constructors and general functionality
    sol = Solution(pdae, Δt, nt)
    @test typeof(sol) <: SolutionPDAE

    # test_interface(sol) # TODO reactivate

    sol0 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    @test typeof(sol0) <: SolutionPDAE

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    @test typeof(sol1) <: SolutionPDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0, y0, μ0)
    get_initial_conditions!(sol, tx, ty, tμ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test ty == y0
    @test tμ == μ0
    @test asol.t == t0
    @test asol.q == x0
    @test asol.p == y0
    @test asol.λ == μ0

    tx .= 0
    ty .= 0
    get_initial_conditions!(sol, tx, ty, 1)
    @test tx == x0
    @test ty == y0

    δt, δx, δy, δμ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0
    @test δy == y0
    @test δμ == μ0

    set_initial_conditions!(sol1, similar(pdae, t1, x2, y2, μ2))
    for j = 1:ni
        get_initial_conditions!(sol1, tx, ty, tμ, j)
        @test tx == x2[j]
        @test ty == y2[j]
        @test tμ == μ2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol2 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol3 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    for i = 1:nt
        tx .= xs[i]
        ty .= ys[i]
        tμ .= μs[i]
        asol.q .= xs[i]
        asol.p .= ys[i]
        asol.λ .= μs[i]
        set_solution!(sol1, tx, ty, tμ, i)
        set_solution!(sol2, tx, ty, i)
        set_solution!(sol3, asol, i)
    end
    @test sol1.q[1:nt] == xs
    @test sol1.p[1:nt] == ys
    @test sol1.λ[1:nt] == μs
    @test sol2.q[1:nt] == xs
    @test sol2.p[1:nt] == ys
    @test sol3.q[1:nt] == xs
    @test sol3.p[1:nt] == ys
    @test sol3.λ[1:nt] == μs

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol2 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol3 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    for i = 1:nt
        for j = 1:ni
            tx .= Xs[i, j]
            ty .= Ys[i, j]
            tμ .= Ms[i, j]
            asol.q .= Xs[i, j]
            asol.p .= Ys[i, j]
            asol.λ .= Ms[i, j]
            set_solution!(sol1, tx, ty, tμ, i, j)
            set_solution!(sol2, tx, ty, i, j)
            set_solution!(sol3, asol, i, j)
        end
    end
    @test sol1.q[1:nt, :] == Xs
    @test sol1.p[1:nt, :] == Ys
    @test sol1.λ[1:nt, :] == Ms
    @test sol2.q[1:nt, :] == Xs
    @test sol2.p[1:nt, :] == Ys
    @test sol3.q[1:nt, :] == Xs
    @test sol3.p[1:nt, :] == Ys
    @test sol3.λ[1:nt, :] == Ms

    # test nsave and nwrite parameters
    sol = Solution(pdae, Δt, 20, nsave = 2)
    @test sol.nt == 10

    sol = Solution(pdae, Δt, 20, nsave = 2, nwrite = 10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pdae, Δt, nt)
    reset!(sol)
    @test sol.t[0] == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt
end


text:
@testset "Elements" begin
  @testset "Constructors" begin
    M = FlintZZ[1 2 3; 4 5 6]
    G = @inferred abelian_group(M)
    N = FlintZZ[1 2 3]
    a = @inferred GrpAbFinGenElem(G, N)
    @test parent(a) == G
    @test a.coeff == N

    G = @inferred abelian_group([3, 0])
    N = FlintZZ[1 1]
    a = @inferred GrpAbFinGenElem(G, N)
    @test @inferred parent(a) == G
    @test a.coeff == N
  end

  @testset "Generators" begin
    M = FlintZZ[1 2 3; 4 5 6]
    G = abelian_group(M)
    ge = @inferred gens(G)
    @test length(ge) == 3
    @test ge[1] == G[1]
    @test ge[2] == G[2]
    @test ge[3] == G[3]
  end

  @testset "Parent" begin
    G = @inferred abelian_group([3, 0])
    N = FlintZZ[1 1]
    a = @inferred GrpAbFinGenElem(G, N)
    @test @inferred parent(a) == G
  end

  @testset "String I/O" begin
    G = abelian_group([3, 0])
    N = FlintZZ[1 1]
    a = GrpAbFinGenElem(G, N)
    @test isa(string(a), String)
  end

  @testset "Hashing" begin
    G = abelian_group([3, 0])
    N = FlintZZ[1 1]
    a = GrpAbFinGenElem(G, N)
    @test isa(hash(a), UInt)
  end

  @testset "Indexing" begin
    G = abelian_group([3, 0])
    N = FlintZZ[1 2]
    a = GrpAbFinGenElem(G, N)
    @test @inferred a[1] == 1
    @test @inferred a[2] == 2
  end

  @testset "Comparison" begin
    G = abelian_group([3, 0])
    N = FlintZZ[1 2]
    a = GrpAbFinGenElem(G, N)
    b = GrpAbFinGenElem(G, deepcopy(N))
    @test @inferred a == b

    H = abelian_group([3, 0])
    c = GrpAbFinGenElem(H, N)
    @test_throws ErrorException a == c
  end

  @testset "Arithmetic" begin
    G = abelian_group([3, 3, 3])
    a = G[1]
    b = G[2]
    c = G([1, 1, 0])
    @test a + b == c
    @test -a == G([2, 0, 0])

    aa = @inferred(2 * a)
    @test aa == G([2, 0, 0])

    aa = @inferred(a * 2)
    @test aa == G([2, 0, 0])

    aa = @inferred(fmpz(2) * a)
    @test aa == G([2, 0, 0])
  end

  @testset "Neutral element" begin
    G = abelian_group([3, 3, 3])
    a = G[1]

    aa = @inferred(a * fmpz(2))
    @test aa == G([2, 0, 0])

    @test !iszero(a)

    c = G([0, 0, 0])
    @test iszero(c)
  end

  @testset "Parent object overloading" begin
    G = abelian_group([3, 3, 3])
    a = @inferred G(fmpz[1, 1, 1])
    @test parent(a) == G
    @test a.coeff == FlintZZ[1 1 1]

    a = @inferred G([1, 1, 1])
    @test parent(a) == G
    @test a.coeff == FlintZZ[1 1 1]

    M = FlintZZ[1 1 1]
    a = @inferred G(M)
    M[1, 1] = 3
    @test parent(a) == G
    @test a.coeff == FlintZZ[1 1 1]

    a = @inferred G[1]
    @test parent(a) == G
    @test a.coeff == FlintZZ[1 0 0]
  end

  @testset "Order" begin
    G = abelian_group([3, 3, 0])
    a = G[1]
    @test @inferred order(a) == 3

    G = abelian_group([3, 5, 0])
    a = G[1]
    @test @inferred order(a) == 3

    a = G[3]
    @test_throws ErrorException order(a)
  end

  @testset "Random elements" begin
    G = abelian_group([3, 5])
    a = @inferred rand(G)
    @test parent(a) == G

    G = abelian_group([3, 15])
    a = @inferred rand(G)
    @test parent(a) == G

    G = abelian_group([3, 0])
    @test_throws ErrorException rand(G)

    a = @inferred rand(G, 10)
    @test parent(a) == G
    @test -10 <= a[2] <= 10

    a = @inferred rand(G, fmpz(10))
    @test parent(a) == G
    @test -10 <= a[2] <= 10

    G = abelian_group([3, 0, 5, 0])
    @test_throws ErrorException rand(G)

    a = @inferred rand(G, 10)
    @test parent(a) == G
    @test -10 <= a[2] <= 10
    @test -10 <= a[4] <= 10

    a = @inferred rand(G, fmpz(10))
    @test parent(a) == G
    @test -10 <= a[2] <= 10
    @test -10 <= a[4] <= 10
  end

  @testset "Iterator" begin
    G = abelian_group([2, 0])
    @test_throws ErrorException begin for a in G end end
    G = abelian_group([fmpz(2)^100])
    @test_throws ErrorException begin for a in G end end

    G = abelian_group([3, 5, 10])
    @test length(G) == 3*5*10
    @test length(collect(G)) == 3*5*10

    G = abelian_group([3, 9, 27])
    @test length(G) == 3*9*27
    @test length(collect(G)) == 3*9*27
  end

  @testset "Helper" begin
    @testset "Reduce mod Hermite normal form" begin
      a = FlintZZ[21 32 43]
      H = FlintZZ[2 0 0 ; 0 3 0 ; 0 0 5]
      Hecke.reduce_mod_hnf_ur!(a, H)
      @test a == FlintZZ[1 2 3]

      a = FlintZZ[1 3 42]
      H = FlintZZ[1 1 14 ; 0 2 11 ; 0 0 17]
      Hecke.reduce_mod_hnf_ur!(a, H)
      @test a == FlintZZ[0 0 0]

      a = FlintZZ[0 0 1]
      H = FlintZZ[1 32 62 ; 0 45 90 ; 0 0 0]
      Hecke.reduce_mod_hnf_ur!(a, H)
      @test a == FlintZZ[0 0 1]
    end

    @testset "Smith normal form with transform" begin
      M = MatrixSpace(FlintZZ,1,1)([0])
      S = MatrixSpace(FlintZZ,1,1)([0])
      T,L,R = snf_with_transform(M, true, true)
      @test S == T
      @test L*M*R == T

      M = MatrixSpace(FlintZZ,1,1)([1])
      S = MatrixSpace(FlintZZ,1,1)([1])
      T,L,R = snf_with_transform(M, true, true)
      @test S == T
      @test L*M*R == T

      M = FlintZZ[834 599 214 915 ; 784 551 13 628 ; 986 5 649 100 ; 504 119 64 310 ]
      S = FlintZZ[1 0 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 0 36533330310]
      T,L,R = snf_with_transform(M, true, true)
      @test S == T
      @test L*M*R == T
      T,L,R = snf_with_transform(M, false, true)
      T,L,R = snf_with_transform(M, true, false)
      T,L,R = snf_with_transform(M, false, false)
    end
  end
end


text:
using IntUtils

function main()
    primes = sieve_of_eratosthenes(2_000_000)
    sum_primes = 0
    for i in 1:length(primes)
        if primes[i]
            sum_primes += i
        end
    end

    @printf("The sum of primes under 2000000 is %d\n", sum_primes)
end

main()


text:
@testset "785.is-graph-bipartite.jl" begin
    @test is_bipartite([[1,2,3], [0,2], [0,1,3], [0,2]]) == false
    @test is_bipartite([[1,3], [0,2], [1,3], [0,2]]) == true
    @test is_bipartite([[3, 6], [3, 6], [4, 6], [0, 1, 5], [2], [3], [0, 1, 2]]) == true
end


text:
export PTBLM
module PTBLM

    using DataDeps
    using ..MLDatasets: datafile, download_dep

    export

        traindata,
        testdata,

        download

    const DEPNAME = "PTBLM"
    const TRAINFILE = "ptb.train.txt"
    const TESTFILE  = "ptb.test.txt"

    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    traindata(; dir = nothing) = traindata(dir)
    testdata(; dir = nothing) = testdata(dir)

    function traindata(dir)
        path = datafile(DEPNAME, TRAINFILE, dir)
        xs = readdata(path)
        ys = makeys(xs)
        xs, ys
    end

    function testdata(dir)
        path = datafile(DEPNAME, TESTFILE, dir)
        xs = readdata(path)
        ys = makeys(xs)
        xs, ys
    end

    function readdata(path)
        lines = open(readlines, path)
        map(l -> Vector{String}(split(chomp(l))), lines)
    end

    function makeys(xs::Vector{Vector{String}})
        map(xs) do x
            y = copy(x)
            popfirst!(y)
            push!(y, "<eos>")
        end
    end

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: Penn Treebank sentences for language modeling
            Website: https://github.com/tomsercu/lstm

            -----------------------------------------------------
            WARNING: EXPERIMENTAL STATUS
            Please be aware that this dataset is from a secondary
            source. The provided interface by this package is not
            as developed as those for other datasets. We would
            welcome any contribution to provide this dataset in a
            more mature manner.
            ------------------------------------------------------

            The PTBLM dataset consists of Penn Treebank sentences
            for language modeling, available from tomsercu/lstm.
            The unknown words are replaced with <unk> so that the
            total vocaburary size becomes 10000.

            The files are available for download at the github
            repository linked above. Note that using the data
            responsibly and respecting copyright remains your
            responsibility.
            """,
            "https://raw.githubusercontent.com/tomsercu/lstm/master/data/" .* [TRAINFILE, TESTFILE],
            "218f4e6c7288bb5efeb03cc4cb8ae9c04ecd8462ebfba8e13e3549fab69dc25f",
        ))
    end
end


text:
using LoopVectorization
using StatsBase
using Tullio

const LIKELY_ERROR_CAUSES = """
1. Incorrect inputs -- check your program for bugs. If you provided an `r_eff` argument,
double check it is correct.
2. Your chains failed to converge. Check diagnostics.
3. You do not have enough posterior samples (ESS < ~25).
"""
const MIN_TAIL_LEN = 5  # Minimum size of a tail for PSIS to give sensible answers
const SAMPLE_SOURCES = ["mcmc", "vi", "other"]

export psis, PsisLoo, PsisLooMethod, Psis

###########################
###### RESULT STRUCT ######
###########################

"""
    Psis{R<:Real, AT<:AbstractArray{R, 3}, VT<:AbstractVector{R}}

A struct containing the results of Pareto-smoothed importance sampling.

# Fields

  - `weights`: A vector of smoothed, truncated, and normalized importance sampling weights.
  - `pareto_k`: Estimates of the shape parameter `k` of the generalized Pareto distribution.
  - `ess`: Estimated effective sample size for each LOO evaluation, based on the variance of
    the weights.
  - `sup_ess`: Estimated effective sample size for each LOO evaluation, based on the
    supremum norm, i.e. the size of the largest weight. More likely than `ess` to warn when
    importance sampling has failed. However, it can have a high variance.
  - `r_eff`: The relative efficiency of the MCMC chain, i.e. ESS / posterior sample size.
  - `tail_len`: Vector indicating how large the "tail" is for each observation.
  - `posterior_sample_size`: How many draws from an MCMC chain were used for PSIS.
  - `data_size`: How many data points were used for PSIS.
"""
struct Psis{
    RealType <: Real,
    AT <: AbstractArray{RealType, 3},
    VT <: AbstractVector{RealType},
}
    weights::AT
    pareto_k::VT
    ess::VT
    sup_ess::VT
    r_eff::VT
    tail_len::Vector{Int}
    posterior_sample_size::Int
    data_size::Int
end

function Base.show(io::IO, ::MIME"text/plain", psis_object::Psis)
    table = hcat(psis_object.pareto_k, psis_object.ess, psis_object.sup_ess)
    post_samples = psis_object.posterior_sample_size
    data_size = psis_object.data_size
    println("Results of PSIS with $post_samples posterior samples and $data_size cases.")
    _throw_pareto_k_warning(psis_object.pareto_k)
    return pretty_table(
        table;
        compact_printing=false,
        header=[:pareto_k, :ess, :sup_ess],
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end

###########################
####### PSIS FNCTNS #######
###########################

"""
    psis(
        log_ratios::AbstractArray{T<:Real},
        r_eff::AbstractVector;
        source::String="mcmc"
    ) -> Psis

Implements Pareto-smoothed importance sampling (PSIS).

# Arguments
## Positional Arguments
  - `log_ratios::AbstractArray`: A 2d or 3d array of (unnormalized) importance ratios on the
    log scale. Indices must be ordered as `[data, step, chain]`. The chain index can be left
    off if there is only one chain, or if keyword argument `chain_index` is provided.
  - $R_EFF_DOC

## Keyword Arguments

  - $CHAIN_INDEX_DOC
  - `source::String="mcmc"`: A string or symbol describing the source of the sample being
    used. If `"mcmc"`, adjusts ESS for autocorrelation. Otherwise, samples are assumed to be
    independent. Currently permitted values are $SAMPLE_SOURCES.

See also: [`relative_eff`]@ref, [`psis_loo`]@ref, [`psis_ess`]@ref.
"""
function psis(
    log_ratios::AbstractArray{<:Real, 3};
    r_eff::AbstractVector{<:Real}=similar(log_ratios, 0),
    source::Union{AbstractString, Symbol}="mcmc"
)

    source = lowercase(String(source))
    dims = size(log_ratios)

    data_size = dims[1]
    post_sample_size = dims[2] * dims[3]

    # Reshape to matrix (easier to deal with)
    log_ratios = reshape(log_ratios, data_size, post_sample_size)
    r_eff = _generate_r_eff(log_ratios, dims, r_eff, source)
    weights = similar(log_ratios)
    # Shift ratios by maximum to prevent overflow
    @tturbo @. weights = exp(log_ratios - $maximum(log_ratios; dims=2))

    _check_input_validity_psis(reshape(log_ratios, dims), r_eff)

    tail_length = Vector{Int}(undef, data_size)
    ξ = similar(r_eff)
    @inbounds Threads.@threads for i in eachindex(tail_length)
        tail_length[i] = _def_tail_length(post_sample_size, r_eff[i])
        ξ[i] = @views ParetoSmooth._do_psis_i!(weights[i, :], tail_length[i])
    end

    @tullio norm_const[i] := weights[i, j]
    @tturbo weights .= weights ./ norm_const
    ess = psis_ess(weights, r_eff)
    inf_ess = sup_ess(weights, r_eff)

    weights = reshape(weights, dims)

    return Psis(
        weights,
        ξ,
        ess,
        inf_ess,
        r_eff,
        tail_length,
        post_sample_size,
        data_size
    )
end

function psis(
    log_ratios::AbstractMatrix{<:Real};
    chain_index::AbstractVector{<:Integer}=_assume_one_chain(log_ratios),
    kwargs...,
)
    new_log_ratios = _convert_to_array(log_ratios, chain_index)
    return psis(new_log_ratios; kwargs...)
end

"""
    _do_psis_i!(is_ratios::AbstractVector{Real}, tail_length::Integer) -> T

Do PSIS on a single vector, smoothing its tail values.

# Arguments

  - `is_ratios::AbstractVector{Real}`: A vector of importance sampling ratios,
    scaled to have a maximum of 1.

# Returns

  - `T<:Real`: ξ, the shape parameter for the GPD; big numbers indicate thick tails.
"""
function _do_psis_i!(is_ratios::AbstractVector{T}, tail_length::Integer) where {T <: Real}

    len = length(is_ratios)
    tail_start = len - tail_length + 1  # index of smallest tail value

    # sort is_ratios and also get results of sortperm() at the same time
    ratio_index = collect(zip(is_ratios, Base.OneTo(len)))
    partialsort!(ratio_index, (tail_start-1):len; by=first)
    is_ratios .= first.(ratio_index)
    @views tail = is_ratios[tail_start:len]
    _check_tail(tail)

    # Get value just before the tail starts:
    cutoff = is_ratios[tail_start - 1]
    ξ = _psis_smooth_tail!(tail, cutoff)

    # truncate at max of raw weights (1 after scaling)
    clamp!(is_ratios, 0, 1)
    # unsort the ratios to their original position:
    invpermute!(is_ratios, last.(ratio_index))

    return ξ::T
end

"""
    _def_tail_length(log_ratios::AbstractVector, r_eff::Real) -> tail_len::Integer

Define the tail length as in Vehtari et al. (2019).
"""
function _def_tail_length(length::Int, r_eff::Real)
    return Int(ceil(min(length / 5, 3 * sqrt(length / r_eff))))
end

"""
    _psis_smooth_tail!(tail::AbstractVector{T}, cutoff::T) where {T<:Real} -> ξ::T

Takes an *already sorted* vector of observations from the tail and smooths it *in place*
with PSIS before returning shape parameter `ξ`.
"""
function _psis_smooth_tail!(tail::AbstractVector{T}, cutoff::T) where {T <: Real}
    len = length(tail)
    if any(isinf.(tail))
        return ξ = Inf
    else
        @turbo @. tail = tail - cutoff

        # save time not sorting since tail is already sorted
        ξ, σ = gpdfit(tail)
        @turbo @. tail = gpd_quantile(($(1:len) - 0.5) / len, ξ, σ) + cutoff
    end
    return ξ
end

##########################
#### HELPER FUNCTIONS ####
##########################

"""
Generate the relative effective sample size if not provided by the user.
"""
function _generate_r_eff(
    weights::AbstractArray{R},
    dims::Base.AbstractVecOrTuple,
    r_eff::T,
    source::String,
)::T where {R<:Real, T<:AbstractVector{R}}
    output::T = similar(r_eff, dims[1])
    if isempty(r_eff)
        if source == "mcmc"
            @info "Adjusting for autocorrelation. If the posterior samples are not " *
                  "autocorrelated, specify the source of the posterior sample using the " *
                  "keyword argument `source`. MCMC samples are always autocorrelated; VI " *
                  "samples are not."
            output .= relative_eff(reshape(weights, dims))
        elseif source ∈ SAMPLE_SOURCES
            @info "Samples have not been adjusted for autocorrelation. If the posterior " *
                  "samples are autocorrelated, as in MCMC methods, ESS estimates will be " *
                  "upward-biased, and standard error estimates will be downward-biased. " *
                  "MCMC samples are always autocorrelated; VI samples are not."
            return output .= ones(R, dims[1])
        else
            throw(
                ArgumentError(
                    "$source is not a valid source. Valid sources are $SAMPLE_SOURCES."
                ),
            )
            return output .= ones(R, dims[1])
        end
    else
        return r_eff
    end
    if any(_invalid_number, r_eff)
        throw(
            ArgumentError(
                "PSIS-LOO has encountered an error calculating ESS values for your " *
                "Markov chains. Please check your inputs. $LIKELY_ERROR_CAUSES"
            )
        )
    end
    return output::T
end

"""
Make sure all inputs to `psis` are valid.
"""
function _check_input_validity_psis(
    log_ratios::AbstractArray{T, 3}, r_eff::AbstractVector{T}
) where {T <: Real}
    if any(_invalid_number, log_ratios)
        throw(DomainError("Invalid input for `log_ratios` (contains NaN  or inf values)."))
    elseif isempty(log_ratios)
        throw(ArgumentError("Invalid input for `log_ratios` (array is empty)."))
    elseif length(r_eff) ≠ size(log_ratios, 1)
        throw(ArgumentError("Size of `r_eff` does not equal the number of data points."))
    end
    return nothing
end

"""
Check the tail to make sure a GPD fit is possible.
"""
function _check_tail(tail::AbstractVector{T}) where {T <: Real}
    if maximum(tail) ≈ minimum(tail)
        throw(
            ArgumentError(
                "Unable to fit generalized Pareto distribution: all tail values are the " *
                "same. Likely causes are: \n$LIKELY_ERROR_CAUSES",
            ),
        )
    elseif length(tail) < MIN_TAIL_LEN
        throw(
            ArgumentError(
                "Unable to fit generalized Pareto distribution: tail length was too " *
                "short. Likely causese are: \n$LIKELY_ERROR_CAUSES",
            ),
        )
    end
    return nothing
end

"""
Check if an input is invalid.
"""
function _invalid_number(x::Real)
    return isinf(x) || isnan(x)
end


text:
# TODO, add `<:` for concrete dispatch when using StringThemSamplableBeliefs
StringThemSamplableBeliefs = Union{<:Uniform, <:Normal, <:MvNormal, <:ZeroMeanDiagNormal, <:Categorical, <:DiscreteNonParametric, <:Rayleigh, <:BallTreeDensity, <:ManifoldKernelDensity, <:AliasingScalarSampler, <:HeatmapGridDensity, <:LevelSetGridNormal}

## TODO, TBD
# Base.@kwdef struct PackedDiscreteNonParametric <: PackedSamplableBelief
#   _type::String        = "IncrementalInference.PackedDiscreteNonParametric"
# end

Base.@kwdef struct PackedCategorical <: PackedSamplableBelief
  _type::String        = "IncrementalInference.PackedCategorical"
  p::Vector{Float64}   = [1.0;]
end

Base.@kwdef mutable struct PackedUniform <: PackedSamplableBelief
  _type::String  = "IncrementalInference.PackedUniform"
  a::Float64     = 0.0
  b::Float64     = 1.0
  PackedSamplableTypeJSON::String = "IncrementalInference.PackedUniform"
end

Base.@kwdef struct PackedNormal <: PackedSamplableBelief
  _type::String  = "IncrementalInference.PackedNormal"
  mu::Float64    = 0.0
  sigma::Float64 = 1.0
end

Base.@kwdef struct PackedZeroMeanDiagNormal <: PackedSamplableBelief
  _type::String        = "IncrementalInference.PackedZeroMeanDiagNormal"
  diag::Vector{Float64}= ones(1)
end

Base.@kwdef struct PackedZeroMeanFullNormal <: PackedSamplableBelief
  _type::String        = "IncrementalInference.PackedZeroMeanFullNormal"
  cov::Vector{Float64} = ones(1)
end

Base.@kwdef mutable struct PackedDiagNormal <: PackedSamplableBelief
  _type::String        = "IncrementalInference.PackedDiagNormal"
  mu::Vector{Float64}  = zeros(1)
  diag::Vector{Float64}= ones(1)
end

Base.@kwdef struct PackedFullNormal <: PackedSamplableBelief
  _type::String        = "IncrementalInference.PackedFullNormal"
  mu::Vector{Float64}  = zeros(1)
  cov::Vector{Float64} = ones(1)
end

Base.@kwdef struct PackedRayleigh <: PackedSamplableBelief
  _type::String = "IncrementalInference.PackedRayleigh"
  sigma::Float64 = 1.0
end

#


text:
export lambda
"""
    @lambda e

Convert `e` to an anonymous function with `_` as the input variable.

# Examples
```julia
lambda_function = @lambda vcat(_, 2)
@test lambda_function(1) == vcat(1, 2)
```
"""
lambda(e) = Expr(:->, :_, e)

"""
    @lambda(f, e, v = :\_)

[`lambda`](@ref) `e` then call `f` on `e` and `v`.

# Examples
```julia
f = :map
e = :(vcat(_, 1))
lambda(f, e)

_ = [1, 2]
@test map(_ -> vcat(_, 1), _) == @lambda map vcat(_, 1)
```
"""
lambda(f, e, v = :_) = @chain begin
    e
    lambda
    Expr(:call, f, _, v)
end

"""
    @lambda(f::Expr, e, v = \_)

[`lambda`](@ref) `e` then insert `e` as as the first and `v` as the last
argument to `f`.

`f` must be a call.

# Examples
```julia
e = :(_ + 1)
f = :(NullableArrays.map(lift = true))

lambda(f, e)

_ = NullableArrays.NullableArray([1, 2], [false, true])

result = @lambda NullableArrays.broadcast(lift = true) _ + 1

@test result.values[1] == 2
@test result.isnull == [false, true]

# `f` must be a call
@test_throws ErrorException lambda(:(import ChainMap), :(_ + 1) )
```
"""
lambda(f::Expr, e, v = :_) = MacroTools.@match f begin
    function_call_(arguments__) => @chain begin
        e
        lambda
        Expr(:call, function_call, _, arguments..., v)
    end
    f_=> error("`f` must be a call")
end

@nonstandard lambda
export @lambda

"""
    @map(e, v = \_)

A convenience macro for [`lambda`](@ref) where `f` = `map`

# Examples
```julia
_ = [1, 2]
@test map(_ -> vcat(_, 1), _) == @map vcat(_, 1)
```
"""
macro map(e, v = :_)
    esc(lambda(:map, e, v))
end
export @map


text:
export PeriodicSegment

"""
	PeriodicSegment(a,b)

represents a periodic interval from `a` to `b`, that is, the point
`b` is identified with `a`.
"""
struct PeriodicSegment{T} <: PeriodicDomain{T}
    a::T
    b::T
    PeriodicSegment{T}() where {T} = new{T}(0,2convert(T,π))
    PeriodicSegment{T}(a,b) where {T} = new{T}(a,b)
end

PeriodicSegment()=PeriodicSegment{Float64}()
PeriodicSegment(a::Int,b::Int) = PeriodicSegment(Float64(a),Float64(b)) #convenience method
PeriodicSegment(a,b) = PeriodicSegment{promote_type(typeof(a),typeof(b))}(a,b)
PeriodicSegment(a::Tuple,b::Tuple) = Interval(Vec(a...),Vec(b...))

function convert(::Type{PeriodicSegment}, d::ClosedInterval)
	a,b = d.left,d.right
    @assert isfinite(a) && isfinite(b)
    PeriodicSegment(a,b)
end

Segment(d::PeriodicSegment) = Segment(leftendpoint(d),rightendpoint(d))
Interval(d::PeriodicSegment) = Interval(leftendpoint(d),rightendpoint(d))
PeriodicSegment(d::Segment) = PeriodicSegment(leftendpoint(d),rightendpoint(d))

convert(::Type{PeriodicSegment{T}}, d::PeriodicSegment) where {T<:Number} = PeriodicSegment{T}(leftendpoint(d),rightendpoint(d))

isambiguous(d::PeriodicSegment) = all(isnan(leftendpoint(d))) && all(isnan(rightendpoint(d)))
convert(::Type{PeriodicSegment{T}},::AnyDomain) where {T<:Number} = PeriodicSegment{T}(NaN,NaN)
convert(::Type{PeriodicSegment{Vec{d,T}}},::AnyDomain) where {d,T} = PeriodicSegment(Vec(fill(NaN,d)...),Vec(fill(NaN,d)...))
convert(::Type{PeriodicSegment{T}},::AnyDomain) where {T} = PeriodicSegment(nan(T),nan(T))
PeriodicSegment{T}(d) where T = convert(PeriodicSegment{T}, d)
PeriodicSegment(d) = convert(PeriodicSegment, d)

## Information
leftendpoint(d::PeriodicSegment) = d.a
rightendpoint(d::PeriodicSegment) = d.b

first(d::PeriodicSegment) = leftendpoint(d)

issubset(a::PeriodicSegment,b::PeriodicSegment) = first(a)∈b && (a.b∈b || a.b == b.b)

# we disable last since the domain is "periodic"

## Map periodic interval

tocanonical(d::PeriodicSegment{T},x) where {T} = π*(tocanonical(Segment(d),x)+1)
tocanonicalD(d::PeriodicSegment{T},x) where {T} = π*tocanonicalD(Segment(d),x)
fromcanonical(d::PeriodicSegment,θ) = fromcanonical(Segment(d),θ/π-1)
fromcanonicalD(d::PeriodicSegment,θ) = fromcanonicalD(Segment(d),θ/π-1)/π

arclength(d::PeriodicSegment) = norm(complexlength(d))
angle(d::PeriodicSegment) = angle(complexlength(d))
complexlength(d::PeriodicSegment) = rightendpoint(d) - leftendpoint(d)
reverseorientation(d::PeriodicSegment) = PeriodicSegment(rightendpoint(d), leftendpoint(d))

==(d::PeriodicSegment,m::PeriodicSegment) = leftendpoint(d) == m.a && rightendpoint(d) == m.b

## algebra

for op in (:*,:+,:-)
    @eval begin
        $op(c::Number,d::PeriodicSegment) = PeriodicSegment($op(c,leftendpoint(d)),$op(c,rightendpoint(d)))
        $op(d::PeriodicSegment,c::Number) = PeriodicSegment($op(leftendpoint(d),c),$op(rightendpoint(d),c))
    end
end

@eval /(d::PeriodicSegment,c::Number) = PeriodicSegment(/(leftendpoint(d),c),/(rightendpoint(d),c))

+(d1::PeriodicSegment,d2::PeriodicSegment) = PeriodicSegment(d1.a+d2.a,d1.b+d2.b)

# These are needed for spaces to auto-convert [a,b] to Interval
function convert(::Type{PeriodicDomain},d::ClosedInterval)
    a,b=d.left,d.right
    a,b = float(a),float(b)
    if isinf(norm(a)) && isinf(norm(b))
        PeriodicLine(d)
    elseif isinf(norm(a)) || isinf(norm(b))
        error("PeriodicRay not implemented")
    else
        PeriodicSegment(d)
    end
end

issubset(a::PeriodicSegment, b::IntervalOrSegment) = Segment(endpoints(a)...)⊆b
issubset(a::IntervalOrSegment, b::PeriodicSegment) = PeriodicSegment(endpoints(a)...)⊆b

boundary(d::ProductDomain{Tuple{A,B}}) where {A<:IntervalOrSegment,B<:PeriodicSegment} =
    UnionDomain((PeriodicSegment(Vec(rightendpoint(factor(d,1)),leftendpoint(factor(d,2))),Vec(rightendpoint(factor(d,1)),rightendpoint(factor(d,2)))),
        PeriodicSegment(Vec(leftendpoint(factor(d,1)),rightendpoint(factor(d,2))),Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2))))))
boundary(d::ProductDomain{Tuple{A,B}}) where {A<:PeriodicSegment,B<:IntervalOrSegment} =
    UnionDomain((PeriodicSegment(Vec(leftendpoint(factor(d,1)),leftendpoint(factor(d,2))),Vec(rightendpoint(factor(d,1)),leftendpoint(factor(d,2)))),
        PeriodicSegment(Vec(rightendpoint(factor(d,1)),rightendpoint(factor(d,2))),Vec(leftendpoint(factor(d,1)),rightendpoint(factor(d,2))))))
boundary(d::ProductDomain{Tuple{A,B}}) where {A<:PeriodicSegment,B<:PeriodicSegment} = EmptyDomain()

union_rule(A::SumSpace{<:Any,<:PeriodicSegment}, B::Space{<:IntervalOrSegment}) =
    union(Space(Interval(domain(A))), B)

Base.setdiff(d::PeriodicSegment, ptsin::UnionDomain{AS}) where {AS <: AbstractVector{P}} where {P <: Point} =
    affine_setdiff(d, ptsin)


text:
function onedMesh(xb, order::Int, numElem::Int)
#= """
  onedMesh    - Generate a one dimensional mesh with numElem elements.
                This routine returns elements of the same type as
                the inputs xb, eConnb

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    x, eConn, indexU, indexC = onedMesh(xb,order,numElem)
  ```

  Arguments:
  - `xb`: Endpoints of the interval
  - `order`: The order of the finite element mesh, 1-linear, 2-quadratic, etc.
  - `numElem`: Number of elements

  Output arguments:
  - `x`: Nodal coordinates of the mesh
  - `eConn`: Element connectivity of adapted mesh
  - `indexU`: Node numbers of unknowns
  - `indexC`: Node numbers of boundaries
""" =#

  #  Create the element connectivity
  eConn = zeros(Int64,numElem,order+1)
  for k=1:numElem
    if (order == 1)
      eConn[k,:] = [k, k+1]
    elseif (order == 2)
      eConn[k,:] = [2*k-1, 2*k, 2*k+1]
    elseif (order == 3)
      eConn[k,:] = [3*k-2, 3*k-1, 3*k, 3*k+1]
    else
      eConn[k,:] = order*k-(order-1):order*k+1
    end
  end

  nNodes = eConn[end,end]
  x = collect(range(xb[1],xb[2],length=nNodes))

  indexU = [2:nNodes-1;]
  indexC = [1;nNodes]

  return x, eConn, indexU, indexC
end