text:
using Base.Threads
using BenchmarkTools
using LinearAlgebra
using LoopVectorization
using MPI
using Strided
using UltraDark: psi_whole_step!, Grids, PencilGrids

if ~MPI.Initialized()
    MPI.Init()
end

n = 128

Δt = 1.0

function step_prim_for!(Δt, psi, phi, nothing)
    for i in eachindex(psi)
        psi[i] *= exp(- im * Δt / 1 * phi[i])
    end
end

function step_prim_for_threads!(Δt, psi, phi, nothing)
    @threads for i in eachindex(psi)
        psi[i] *= exp(- im * Δt / 1 * phi[i])
    end
end

function step_for!(Δt, grids, nothing)
    for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_cartesian!(Δt, grids, nothing)
    for i in CartesianIndices(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_threads!(Δt, grids, nothing)
    @threads for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_inbounds_threads!(Δt, grids, nothing)
    @inbounds @threads for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_fastmath_inbounds_threads!(Δt, grids, nothing)
    @fastmath @inbounds @threads for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_avx!(Δt, grids, nothing)
    @avx for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_avxt!(Δt, grids, nothing)
    @avxt for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_turbo!(Δt, grids, nothing)
    @turbo for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_for_tturbo!(Δt, grids, nothing)
    @tturbo for i in eachindex(grids.ψx)
        grids.ψx[i] *= exp(- im * Δt / 1 * grids.Φx[i])
    end
end

function step_broadcast!(Δt, grids, nothing)
    @. grids.ψx *= exp(- im * Δt / 1 * grids.Φx)
end

function step_broadcast_strided!(Δt, grids, nothing)
    @strided @. grids.ψx *= exp(- im * Δt / 1 * grids.Φx)
end

# function step_broadcast_turbo!(Δt, grids, nothing)
#     @turbo @. grids.ψx *= exp(- im * Δt / 1 * grids.Φx)
# end

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
suite = BenchmarkGroup()

for (label, grids) in [("Grids", Grids(1.0, n)), ("PencilGrids", PencilGrids(1.0, n))]

    grids.ψx .= randn(size(grids.ψx))
    grids.ψk .= randn(size(grids.ψk))

    suite[label] = BenchmarkGroup()

    for f in [
              step_prim_for!,
              step_prim_for_threads!,
             ]
        suite[label]["$f"] = @benchmarkable $f($Δt, $(grids.ψx), $(grids.Φx), $nothing)
    end

    for f in [
              psi_whole_step!,
              step_for!,
              step_for_cartesian!,
              step_for_threads!,
              step_for_inbounds_threads!,
              step_for_fastmath_inbounds_threads!,
              step_for_avx!,
              step_for_avxt!,
              step_for_turbo!,
              step_broadcast!,
              # step_broadcast_strided!,
              # step_broadcast_turbo!,
             ]
        suite[label]["$f"] = @benchmarkable $f($Δt, $grids, $nothing)
    end

end

@info "Tuning"
@show tune!(suite)

@info "Running benchmarks"
results = run(suite)

@info "Results"
@show res_min = minimum(results)
# @show res_med = median(results)

@show judge(res_min["Grids"], res_min["PencilGrids"])


text:
using PQEcolaPoint
using Dates

# what is the contract?
short_put_contract = PutContractModel()
short_put_contract.ticker = "XYZ"
short_put_contract.expiration_date = Date(2022, 04, 14)
short_put_contract.strike_price = 105.0
short_put_contract.premium = 0.0
short_put_contract.number_of_contracts = 1
short_put_contract.direction = 1

# other version -
val_vega = vega(short_put_contract; Sₒ=113.46, number_of_levels=80, σ=0.5212, T=(27.0 / 365), μ=0.0047)

ext: jl


text:
# addition of positive definite matrices

using PDMats
using Base.Test

for T in [Float64,Float32]

  print_with_color(:blue, "Testing addition with eltype = $T\n")
  M = convert(Array{T,2},[4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
  V = convert(Array{T,1},[1.5, 2.5, 2.0])
  X = convert(T,2.0)

  pm1 = PDMat(M)
  pm2 = PDiagMat(V)
  pm3 = ScalMat(3,X)
  pm4 = X*I
  #do not test PDSparseMat for Float32 on Julia 0.3 - impossible because of issue #14076
  pm5 = (T==Float64 || VERSION > v"0.4.2")?PDSparseMat(sparse(M)):PDMat(M)

  pmats = Any[pm1, pm2, pm3, pm5]

  for p1 in pmats, p2 in pmats
	  pr = p1 + p2
	  @test size(pr) == size(p1)
	  @test_approx_eq full(pr) full(p1) + full(p2)

	  pr = pdadd(p1, p2, convert(T,1.5))
	  @test size(pr) == size(p1)
	  @test_approx_eq full(pr) full(p1) + full(p2) * convert(T,1.5)
  end

  for p1 in pmats
        pr = p1 + pm4
        @test size(pr) == size(p1)
        @test_approx_eq full(pr) full(p1) + pm4
  end
end


text:
using RandomLinearAlgebraSolvers
using LinearAlgebra, LLSModels, SparseArrays, Stopping
using Random, Test

using MatrixDepot

Random.seed!(1234)

# algo_list = [:RandomVectorSketch, :randomInv]
algo_list = [:RandomizedKaczmarz, :RandomizedBlockKaczmarz, :RandomizedCD]
algo_pd_list = [:RandomizedCD2, :RandomizedNewton]

@testset "RandomLinearAlgebraSolvers.jl" begin
  for T in [Float16, Float32, Float64], algo in union(algo_list, algo_pd_list)
    A = T[1 0; 0 1]
    b = T[1; 1]
    stp = eval(algo)(RLAStopping(A, b))
    sol = stp.current_state.x
    @test eltype(sol) == T
  end
end

function overdetermined_with_random_matrices(T = Float64, m = 10, n = 5)
  A = rand(T, m, n)
  sol = rand(T, n)
  b = A * sol
  return A, b, sol
end

function overdetermined_with_sparse_random_matrices(T = Float64, m = 10, n = 5)
  A = sprand(T, m, n, 1 / log(n * m))
  sol = rand(T, n)
  b = A * sol
  return A, b, sol
end

function overdetermined_with_random_pd_matrices(T = Float64, n = 10)
  A = rand(T, n, n)
  A = A' * A
  A = (A + A') / 2
  sol = rand(T, n)
  b = A * sol
  return A, b, sol
end

function hilbert(T = Float64, n = 4)
  A = matrixdepot("hilb", n)
  sol = rand(T, n)
  b = A * sol
  return A, b, sol
end

pb_generators =
  [:hilbert, :overdetermined_with_random_matrices, :overdetermined_with_sparse_random_matrices]

pb_pd_generators = [:overdetermined_with_random_pd_matrices]

for meth in algo_list, fpb in pb_generators
  @testset "$meth: Test $(fpb)" begin
    A, b, sol = eval(fpb)()
    stp = RLAStopping(A, b, max_iter = 1000000, max_eval = 100000, rtol = 1e-2, atol = 0.0)
    eval(meth)(stp)
    @test status(stp) == :Optimal
    @show @allocated eval(meth)(stp)
  end
end

for meth in union(algo_list, algo_pd_list), fpb in pb_pd_generators
  @testset "$meth: Test $(fpb)" begin
    A, b, sol = eval(fpb)()
    stp = RLAStopping(A, b, max_iter = 1000000, max_eval = 100000, rtol = 1e-2, atol = 0.0)
    eval(meth)(stp)
    @test status(stp) == :Optimal
    @show @allocated eval(meth)(stp)
  end
end


text:
"""
    mutable struct Reaction
        id::String
        metabolites::Dict{String,Float64}
        lb::Float64
        ub::Float64
        grr::Maybe{GeneAssociation}
        subsystem::Maybe{String}
        notes::Notes
        annotations::Annotations
        objective_coefficient::Float64
    end

A structure for representing a single reaction in a [`StandardModel`](@ref).
"""
mutable struct Reaction
    id::String
    metabolites::Dict{String,Float64}
    lb::Float64
    ub::Float64
    grr::Maybe{GeneAssociation}
    subsystem::Maybe{String}
    notes::Notes
    annotations::Annotations
    objective_coefficient::Float64
end

"""
    Reaction(
        id = "";
        metabolites = Dict{String,Float64}(),
        lb = -_constants.default_reaction_bound,
        ub = _constants.default_reaction_bound,
        grr = nothing,
        subsystem = nothing,
        notes = Notes(),
        annotations = Annotations(),
        objective_coefficient = 0.0,
    )

A constructor for Reaction that only takes a reaction `id` and
assigns default/uninformative values to all the fields that are not
explicitely assigned.
"""
function Reaction(
    id = "";
    metabolites = Dict{String,Float64}(),
    lb = -_constants.default_reaction_bound,
    ub = _constants.default_reaction_bound,
    grr = nothing,
    subsystem = nothing,
    notes = Notes(),
    annotations = Annotations(),
    objective_coefficient = 0.0,
)
    mets = Dict(k => float(v) for (k, v) in metabolites)
    return Reaction(
        id,
        mets,
        lb,
        ub,
        grr,
        subsystem,
        notes,
        annotations,
        objective_coefficient,
    )
end

"""
    Reaction(
        id::String,
        metabolites::Dict{String,Union{Int, Float64}},
        dir = :bidirectional;
        default_bound = _constants.default_reaction_bound,
    )

Convenience constructor for `Reaction`. The reaction equation is specified using
`metabolites`, which is a dictionary mapping metabolite ids to stoichiometric
coefficients. The direcion of the reaction is set through `dir` which can take
`:bidirectional`, `:forward`, and `:reverse` as values. Finally, the
`default_bound` is the value taken to mean infinity in the context of constraint
based models, often this is set to a very high flux value like 1000.
"""
function Reaction(
    id::String,
    metabolites,
    dir = :bidirectional;
    default_bound = _constants.default_reaction_bound,
)
    if dir == :forward
        lb = 0.0
        ub = default_bound
    elseif dir == :reverse
        lb = -default_bound
        ub = 0.0
    else
        lb = -default_bound
        ub = default_bound
    end
    Reaction(id; metabolites = metabolites, lb = lb, ub = ub)
end


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(The BTG Program)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#################################################################################
#copy this box into a file to use the BTG program
using LinearAlgebra
using Test
using Distributions
using DataFrames
using CSV
using Polynomials
using PyPlot
include("../validation/det_loocv.jl")
include("../validation/loocv.jl")
include("../quadrature/quadrature.jl")
include("../transforms/transforms.jl")
include("../priors/priors.jl")
include("../various_mle_things/incremental.jl")
include("../kernels/kernel.jl")
include("../datastructs.jl")
include("../computation/buffers0.jl") #datastruct, kernel, incremental, quadrature
include("../model0.jl") #buffers, datastructs, several auxiliary
include("../computation/finitedifference.jl")
include("../computation/tdist.jl") #model0 and buffer0
#################################################################################

df = DataFrame(CSV.File("../datasets/abalone.csv"))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age
normalizing_constant = maximum(target)
target = target/normalizing_constant #normalization

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~(btg1)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Attributes:
# - single length scale
# - 2-dimensional covariates
# - 3-multidimensional location vectors

#ind = 120:125
#ind = [1:9;11:25]
#ind = 1:30
ind = 100:150
#ind = 45:65
posx = 1:4 #
posc = 1:1
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean

Fx = data[ind, posc]
y = float(target[ind])
#pind = 10:10 #prediction index
pind = 27:27
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
rangeθ = [100.0 200]
#rangeθ = [50.0 1000]
rangeλ = [0.1 1.5] #we will always used 1 range scale for lambda
btg1 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
#θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node
#λ1 = btg1.nodesWeightsλ.nodes[3]
##################################################################################################

println("Test point attributes")
println("Location:")
println(data[pind, posx])
println("Covariates")
println(data[pind, posc])
println("Label")
println(target[pind])

#(_, _, _) = solve(btg1)
(pdf, cdf, dpdf) = solve(btg1);
    a = y0 -> dpdf(x0, Fx0, y0);
    b = y0 -> pdf(x0, Fx0, y0);
    c = y0 -> cdf(x0, Fx0, y0);

#println("Checking derivative of btg_pdf...")
(_, _, plt1, pol1) = checkDerivative(b, a, 0.5, nothing, 4, 12, 10); #first arg is function, second arg is derivative
(_, _, plt2, pol2) = checkDerivative(c, b, 0.5, nothing, 4, 12, 10); #first arg is function, second arg is derivative

PyPlot.figure(1)
(x, y) = plt_data(b, .0001, 2, 100)
(x1, y1) = plt_data(c, .0001, 2, 100)

PyPlot.plot(x, y)
#PyPlot.plot(x1, y1)
#plt(b, .01, 1; label = "pdf")
#plt!(c, .01, 1; label = "cdf")
#plt(target[pind])

for i = 151:156
    pind = i:i
    #prediction data
    Fx0 = reshape(data[pind, posc], 1, length(posc))
    x0 = reshape(data[pind, posx], 1, length(posx))
    b = y0 -> pdf(x0, Fx0, y0);
    c = y0 -> cdf(x0, Fx0, y0);
    figure(i)
    (x, y) = plt_data(b, .01, 2, 400)
    (x1, y1) = plt_data(c, .01, 2, 400)
    PyPlot.plot(x, y)
    PyPlot.plot(x1, y1)
    PyPlot.axvline(x=target[i])
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@printf("\nTesting BlockArrays\n")

using AbstractOperators.BlockArrays
using LinearAlgebra
using Printf

T = Float64
x = [2.0; 3.0]
xb = ([2.0; 3.0], Complex{T}[4.0; 5.0; 6.0], [1.0 2.0 3.0; 4.0 5.0 6.0])

x2 = randn(2)
x2b = (randn(2), randn(3)+im.*randn(3), randn(2,3))

@test blocksize(x) == (2,)
@test blocksize(xb) == ((2,),(3,),(2,3))

@test blockeltype(x) == T
@test blockeltype(xb) == (T, Complex{T},T)

@test blocklength(x) == 2
@test blocklength(xb) == 2+3+6

@test blockvecnorm(x) == norm(x)
@test blockvecnorm(xb) == sqrt(norm(xb[1])^2+norm(xb[2])^2+norm(xb[3])^2)

@test blockmaxabs(x) == 3.0
@test blockmaxabs(xb) == 6.0

@test typeof(blocksimilar(x)) == typeof(x)
@test typeof(blocksimilar(xb)) == typeof(xb)

@test blockcopy(x) == x
@test blockcopy(xb) == xb

y = blocksimilar(x)
yb = blocksimilar(xb)
blockcopy!(y,x)
blockcopy!(yb,xb)

@test y == x
@test yb == xb

y = blocksimilar(x)
yb = blocksimilar(xb)
blockset!(y,x)
blockset!(yb,xb)

@test y == x
@test yb == xb

z = blockvecdot(x,x2)
zb = blockvecdot(xb,x2b)
@test z == dot(x,x2)
@test zb == dot(xb[1],x2b[1])+dot(xb[2],x2b[2])+dot(xb[3],x2b[3])

y = blockzeros(x)
yb = blockzeros(xb)
@test y == zeros(2)
@test yb == (zeros(2),zeros(3)+im*zeros(3),zeros(2,3))

y = blockzeros(blocksize(x))
yb = blockzeros(blocksize(xb))
@test y == zeros(2)
@test yb == (zeros(2),zeros(3)+im*zeros(3),zeros(2,3))

y = blockzeros(blockeltype(x),blocksize(x))
yb = blockzeros(blockeltype(xb),blocksize(xb))
@test y == zeros(2)
@test yb == (zeros(2),zeros(3)+im*zeros(3),zeros(2,3))

y = blockones(x)
yb = blockones(xb)
@test y == ones(2)
@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))

y = blockones(blocksize(x))
yb = blockones(blocksize(xb))
@test y == ones(2)
@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))

y = blockones(blockeltype(x),blocksize(x))
yb = blockones(blockeltype(xb),blocksize(xb))
@test y == ones(2)
@test yb == (ones(2),ones(3)+im*zeros(3),ones(2,3))

blockscale!(y,2,x2)
blockscale!(yb,2,x2b)

@test y ==  2 .*x2
@test yb == (2 .*x2b[1], 2 .*x2b[2], 2 .*x2b[3])

blockcopy!(y,x)
blockcopy!(yb,xb)

blockcumscale!(y,2,x2)
blockcumscale!(yb,2,x2b)

@test y == x .+ 2 .*x2
@test yb == (xb[1] .+ 2 .*x2b[1], xb[2] .+  2 .*x2b[2], xb[3] .+  2 .*x2b[3])

blockaxpy!(y,x,2,x2)
blockaxpy!(yb,xb,2,x2b)

@test y == x .+ 2 .*x2
@test yb == (xb[1] .+ 2 .*x2b[1], xb[2] .+ 2 .*x2b[2], xb[3] .+ 2 .*x2b[3])

x = (ones(Float64, 5), zeros(Float64, 2, 3))
@test blockiszero(x) == false

y = (zeros(Complex{Float64}, 3, 2), ones(Complex{Float64}, 5))
@test blockiszero(x) == false

x = (zeros(Float64, 5), zeros(Complex{Float64}, 2, 3))
@test blockiszero(x) == true

y = (zeros(Float64, 2, 3), zeros(Complex{Float64}, 5))
@test blockiszero(x) == true


text:
module Heat1D

using LinearAlgebra
using SparseArrays

include("structures.jl")
export FDHeatProblem, PeriodicBC, BEuler
include("constructors.jl")
export FDHeatProblem
include("assembly.jl")
export assemble_system!
include("solve.jl")
export solve_heat

end # module


text:
module StructuralBase

using Optim, MDPTools, Reexport
@reexport using StatsBase
const AbstractUtility = MDPTools.AbstractUtility
const AbstractState = MDPTools.AbstractState
export EstimationMethod, EstimationResults, AbstractUtility,
        AbstractTrace, ConvergenceInfo, LinearUtility,
        AbstractState, hessian, gradient, loglikelihood,
        tstat

import StatsBase: loglikelihood
import Base: gradient

abstract EstimationMethod
abstract ConvergenceInfo

# --- AbstractTrace ---

abstract AbstractTrace
abstract Trace

type EstimationResults{T<:EstimationMethod, Tf<:Real}
    E::T
    loglikelihood::Tf
    ∇loglikelihood::Vector{Tf}
    ∇²loglikelihood::Matrix{Tf}
    tdf::Optim.TwiceDifferentiable #TwiceDifferentiableFunction
    coef::Vector{Tf}
    conv::ConvergenceInfo
    trace::Any#::Optim.MultivariateOptimizationResults FIXME should be MultivariateOptimizationResults
    nobs
    meta
end

# computational details
convinfo(res::EstimationResults) = res.conv
outer_iterations(res::EstimationResults) = convinfo(res).outer
inner_iterations(res::EstimationResults) = inner_iterations(convinfo(res))

newton_iterations(res::EstimationResults) = newton_iterations(convinfo(res))
contraction_iterations(res::EstimationResults) = contraction_iterations(convinfo(res))

# loglikelihood
nobs(res) = res.nobs
loglikelihood(res) = res.loglikelihood
loglikelihood(res, x) = -res.tdf.f(x)*nobs(res)
gradient(res) = res.∇loglikelihood
function gradient(res, x)
    g = similar(x)
    res.tdf.g!(x, g)
    -g*nobs(res)
end
hessian(res::EstimationResults) = res.∇²loglikelihood
function hessian(res::EstimationResults, x)
    n = length(x)
    h = zeros(n, n)
    res.tdf.h!(x, h)
    -h*nobs(res)
end

# coefficients and inference
coef(res::EstimationResults) = res.coef
coef(res, i) = coef(res)[i]
vcov(res::EstimationResults) = inv(-hessian(res))
stderr(res::EstimationResults) = sqrt.(diag(vcov(res)))
stderr(res::EstimationResults, i) = sqrt.(diag(vcov(res)))[i]
tstat(res, i, x) = (coef(res, i)-x)/(stderr(res, i))
tstat(res, i) = tstat(res, i, 0)
tstat(res, x::Vector) = [tstat(res, i, x[i]) for i = 1:length(x)]
tstats(res) = tstat(res, zeros(coef(res)))
end # module


text:
FECellBasisStyle(::Type{T}) where T = Val{false}()

FECellBasisStyle(cell_basis) = FECellBasisStyle(typeof(cell_basis))

FECellBasisStyle(::Type{T}) where T<:CellField = Val(is_basis(T))

"""
"""
function is_a_fe_cell_basis(cell_basis)
  v = FECellBasisStyle(cell_basis)
  get_val_parameter(v)
end

function is_a_fe_cell_basis(::Type)
  @unreachable "is_a_fe_cell_basis cannot be called on types"
end


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.
#***********************************************************************

using OffsetArrays

function initsl(ncells::Int, fc::Int, lc::Int, fx::Int, lx::Int, ifirst::Int, ilast::Int,
                conservd::OffsetArray{Float64,1},
                x::OffsetArray{Float64,1})
#   ******************************************************************
#   interior values only; others defined by boundary conditions
#   ******************************************************************
    dx=(x_right-x_left)/ncells
    @inbounds for ie=ifirst:ilast+1
        x[ie]=x_left+ie*dx
    end
    ijump=max(ifirst-1,min(convert(Int,round(ncells*(jump-x_left)/(x_right-x_left))),ilast+1))
    @inbounds for ic=ifirst:ijump-1
        conservd[ic]=statelft
    end
    frac=(jump-x_left-ijump*dx)/(x_right-x_left)
    conservd[ijump]=statelft*frac+statergt*(1.0-frac)
    @inbounds for ic=ijump+1:ilast
        conservd[ic]=statergt
    end
end

function bcmesh(fim::Int,lam::Int,ncells::Int, x::OffsetArray{Float64,1})
#   ******************************************************************
#   right side
#   ******************************************************************
    if lam >= ncells
        it=ncells
        dx=x(it)-x(it-1)
        @inbounds for ie=it+1:lam+1
            x[ie]=x[it]+dx*(ie-it)
        end
    end
#   ******************************************************************
#   left side
#   ******************************************************************
    if fim < 0
        dx=x(1)-x(0)
        @inbounds for ie=fim:-1
            x[ie]=x[0]+dx*(ie)
        end
    end
end

function bccells(fic::Int,lac::Int,ncells::Int, conservd::OffsetArray{Float64,1})
#   ******************************************************************
#   right side
#   ******************************************************************
    if lac >= ncells
        # println("treating right side")
        @inbounds for ic=ncells:lac
            # outgoing wave:
            conservd[ic]=conservd[ncells-1]
        end
    end
#     ******************************************************************
#     left side
#     ******************************************************************
    if fic < 0
        # println("treating left side")
        @inbounds for ic=fic:-1
            # outgoing wave:
            # conservd[ic]=conservd[0]

            # specified value
            conservd[ic]=statelft
        end
    end
end


text:
using Test
using MLJText

@testset "abstract text transformer" begin
    include("abstract_text_transformer.jl")
end

@testset "scitypes" begin
    include("scitypes.jl")
end


text:
# ----------------------------------------
#     Non-linear activation functions
# ----------------------------------------
import Base.abs
import Base.reshape
import Base.vcat
import Base.sind
import Base.sinpi

import Base.exp2
import Base.exp10

import Base.log2
import Base.log10

export abs
export reshape
export sind
export sinpi

export exp2
export exp10

export log2
export log10

export min2max, min2max!
## -------------------------------------------------------- min2max
"""
    min2max!(x::AbstractArray; lower=0.0, upper=1.0) -> x

limit the scope of the data, i.e. ⤦\n
    @. x = min(max(x, lower), upper)
"""
function min2max!(x::AbstractArray; lower=0.0, upper=1.0)
    T = eltype(x)
    L = T(lower)
    U = T(upper)
    @. x = min(max(x, L), U)
end

"""
    min2max(x::AbstractArray; lower=0.0, upper=1.0) -> y

limit the scope of the data, i.e. ⤦\n
    y = min.(max.(x, lower), upper)
"""
function min2max(x::AbstractArray; lower=0.0, upper=1.0)
    T = eltype(x)
    L = T(lower)
    U = T(upper)
    return min.(max.(x, L), U)
end

"""
    min2max!(x::Variable{S}; lower=0.0, upper=1.0) -> y::Variable{S}

limit the scope of the data, i.e. ⤦\n
    y = Variable{S}(min2max!(ᵛ(x), lower=lower, upper=upper), x.backprop)
"""
function min2max!(x::Variable{S}; lower=0.0, upper=1.0) where S
    y = Variable{S}(min2max!(ᵛ(x), lower=lower, upper=upper), x.backprop)
    if x.backprop
        function min2maxBackward()
            if need2computeδ!(x)
                T = eltype(S)
                L = T(lower)
                U = T(upper)
                ∇ = L .< ᵛ(x) .< U
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, min2maxBackward)
    end
    return y
end

"""
    min2max(x::Variable{S}; lower=0.0, upper=1.0) where S -> y::Variable{S}

limit the scope of the data, i.e. ⤦\n
    y = Variable{S}(min2max(ᵛ(x), lower=lower, upper=upper), x.backprop)
"""
function min2max(x::Variable{S}; lower=0.0, upper=1.0) where S
    y = Variable{S}(min2max(ᵛ(x), lower=lower, upper=upper), x.backprop)
    if x.backprop
        function min2maxBackward()
            if need2computeδ!(x)
                T = eltype(S)
                L = T(lower)
                U = T(upper)
                ∇ = L .< ᵛ(x) .< U
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, min2maxBackward)
    end
    return y
end

export relu, relu!
## -------------------------------------------------------- relu
function relu!(x::AbstractArray)
    @. x = max(0.0, x)
end

function relu(x::AbstractArray)
    𝟘 = eltype(x)(0.0)
    return max.(𝟘, x)
end

function relu!(x::Variable{T}) where T
    y = Variable{T}(relu!(ᵛ(x)), x.backprop)
    if x.backprop
        function reluBackward()
            if need2computeδ!(x)
                ∇ = ᵛ(x) .> 0.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, reluBackward)
    end
    return y
end

function relu(x::Variable{T}) where T
    y = Variable{T}(relu(ᵛ(x)), x.backprop)
    if x.backprop
        function reluBackward()
            if need2computeδ!(x)
                ∇ = ᵛ(x) .> 0.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, reluBackward)
    end
    return y
end

export relu1, relu1!
## -------------------------------------------------------- relu1
function relu1!(x::AbstractArray)
    @. x = min(max(x, 0.0), 1.0)
end

function relu1(x::AbstractArray)
    T = eltype(x)
    𝟘 = T(0.0)
    𝟙 = T(1.0)
    return min.(max.(x, 𝟘), 𝟙)
end

function relu1!(x::Variable{T}) where T
    y = Variable{T}(relu1!(ᵛ(x)), x.backprop)
    if x.backprop
        function relu1Backward()
            if need2computeδ!(x)
                ∇ = 0.0 .< ᵛ(x) .< 1.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, relu1Backward)
    end
    return y
end

function relu1(x::Variable{T}) where T
    y = Variable{T}(relu1(ᵛ(x)), x.backprop)
    if x.backprop
        function relu1Backward()
            if need2computeδ!(x)
                ∇ = 0.0 .< ᵛ(x) .< 1.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, relu1Backward)
    end
    return y
end

export relu6, relu6!
## -------------------------------------------------------- relu6
function relu6!(x::AbstractArray)
    @. x = min(max(x, 0.0), 6.0)
end

function relu6(x::AbstractArray)
    T = eltype(x)
    𝟘 = T(0.0)
    𝟞 = T(6.0)
    return min.(max.(x, 𝟘), 𝟞)
end

function relu6!(x::Variable{T}) where T
    y = Variable{T}(relu6!(ᵛ(x)), x.backprop)
    if x.backprop
        function relu6Backward()
            if need2computeδ!(x)
                ∇ = 0.0 .< ᵛ(x) .< 6.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, relu6Backward)
    end
    return y
end

function relu6(x::Variable{T}) where T
    y = Variable{T}(relu6(ᵛ(x)), x.backprop)
    if x.backprop
        function relu6Backward()
            if need2computeδ!(x)
                ∇ = 0.0 .< ᵛ(x) .< 6.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, relu6Backward)
    end
    return y
end

export hardtanh, hardtanh!
## -------------------------------------------------------- hardtanh
function hardtanh!(x::AbstractArray)
    T  = eltype(x)
    𝟙₋ = T(-1.0)
    𝟙₊ = T( 1.0)
    @. x = min(max(x, 𝟙₋), 𝟙₊)
end

function hardtanh(x::AbstractArray)
    T = eltype(x)
    𝟙₋ = T(-1.0)
    𝟙₊ = T( 1.0)
    return min.(max.(x, 𝟙₋), 𝟙₊)
end

function hardtanh!(x::Variable{T}) where T
    y = Variable{T}(hardtanh!(ᵛ(x)), x.backprop)
    if x.backprop
        function hardtanhBackward()
            if need2computeδ!(x)
                ∇ = abs(ᵛ(x)) .< 1.0f0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, hardtanhBackward)
    end
    return y
end

function hardtanh(x::Variable{T}) where T
    y = Variable{T}(hardtanh(ᵛ(x)), x.backprop)
    if x.backprop
        function hardtanhBackward()
            if need2computeδ!(x)
                ∇ = abs(ᵛ(x)) .< 1.0
                δ(x) .+= δ(y) .* ∇
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, hardtanhBackward)
    end
    return y
end

export leakyrelu, leakyrelu!
## -------------------------------------------------------- leakyrelu
function leakyrelu!(x::AbstractArray)
    ZPONE = eltype(x)(0.1)
    @. x = max(ZPONE * x, x)
end

function leakyrelu(x::AbstractArray)
    ZPONE = eltype(x)(0.1)
    return max.(ZPONE .* x, x)
end

function leakyrelu!(x::Variable{T}) where T
    ZPONE = eltype(x)(0.1)
    tempv = ᵛ(x) .* ZPONE
    ᵛ(x) .= max.(ᵛ(x), tempv)
    y = Variable{T}(ᵛ(x), x.backprop)
    if x.backprop
        mask1 = ᵛ(x) .> tempv
        mask2 = .!mask1
        function leakyreluBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (mask1 .+ ZPONE .* mask2)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, leakyreluBackward)
    end
    return y
end

function leakyrelu(x::Variable{T}) where T
    ZPONE = eltype(x)(0.1)
    tempv = ᵛ(x) .* ZPONE
    mask1 = ᵛ(x) .> tempv
    mask2 = .!mask1
    y = Variable{T}(max.(tempv, ᵛ(x)), x.backprop)
    if x.backprop
        function leakyreluBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (mask1 + ZPONE .* mask2)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, leakyreluBackward)
    end
    return y
end

export sigmoid, sigmoid!
## -------------------------------------------------------- sigmoid
function sigmoid!(x::AbstractArray)
    𝟙 = eltype(x)(1.0)
    @. x = 𝟙 / (𝟙 + exp(-x))
end

function sigmoid(x::AbstractArray)
    𝟙 = eltype(x)(1.0)
    return 𝟙 ./ (𝟙 .+ exp.(-x))
end

function sigmoid!(x::Variable{T}) where T
    y = Variable{T}(sigmoid!(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙 = eltype(x)(1)
        function sigmoidBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y) .* (𝟙 .- ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sigmoidBackward)
    end
    return y
end

function sigmoid(x::Variable{T}) where T
    y = Variable{T}(sigmoid(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙 = eltype(x)(1.0)
        function sigmoidBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y) .* (𝟙 .- ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sigmoidBackward)
    end
    return y
end

export swish, swish!
## -------------------------------------------------------- swish
function swish!(x::AbstractArray)
    𝟙 = eltype(x)(1.0)
    @. x = x / (𝟙 + exp(-x))
end

function swish(x::AbstractArray)
    𝟙 = eltype(x)(1.0)
    return  x ./ (𝟙 .+ exp.(-x))
end

function swish!(x::Variable{T}) where T
    return dotMul(sigmoid(x), x)
end

function swish(x::Variable{T}) where T
    return dotMul(sigmoid(x), x)
end

export softmax
## -------------------------------------------------------- softmax
function softmax(x::AbstractArray; dims::Union{Int,NTuple{N,Int}}) where N
    y = exp.(x .- maximum(x, dims=dims))
    Σ = eltype(x)(1.0) ./ sum(y, dims=dims)
    return y .* Σ
end

function softmax(x::Variable{T}; dims::Union{Int,NTuple{N,Int}}) where {T,N}
    y = Variable{T}(softmax(ᵛ(x); dims=dims), x.backprop)
    if x.backprop
        function softmaxBackward()
            if need2computeδ!(x)
                ẏy = δ(y) .* ᵛ(y);
                δ(x) .+= ẏy .- ᵛ(y) .* sum(ẏy, dims=dims);
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, softmaxBackward)
    end
    return y
end

# -----------------
# 不常用激活函数....
# -----------------
export softplus, softplus!
function softplus!(x::AbstractArray)
    @. x = log(1.0 + exp(x))
end

function softplus(x::AbstractArray)
    𝟙 = eltype(x)(1.0)
    return log.( 𝟙 .+ exp.(x) )
end

function softplus!(x::Variable{T}) where T
    y = Variable{T}(softplus(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙 = eltype(x)(1.0)
        function softplusBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (𝟙 .+ exp.( - ᵛ(x) ))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, softplusBackward)
    end
    return y
end

function softplus(x::Variable{T}) where T
    y = Variable{T}(softplus(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙 = eltype(x)(1.0)
        function softplusBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (𝟙 .+ exp.( - ᵛ(x) ))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, softplusBackward)
    end
    return y
end

export exp!
function exp!(x::Variable{T}) where T
    y = Variable{T}(exp!(ᵛ(x)), x.backprop)
    if x.backprop
        function expBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, expBackward)
    end
    return y
end

function Base.:exp(x::Variable{T}) where T
    y = Variable{T}(exp(ᵛ(x)), x.backprop)
    if x.backprop
        function expBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, expBackward)
    end
    return y
end

export log!
function log!(x::Variable{T}) where T
    y = Variable{T}(log(ᵛ(x)), x.backprop)
    if x.backprop
        function logBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ ᵛ(x)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, logBackward)
    end
    return y
end

function Base.:log(x::Variable{T}) where T
    y = Variable{T}(log(ᵛ(x)), x.backprop)
    if x.backprop
        function logBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ ᵛ(x)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, logBackward)
    end
    return y
end

export abs!
function abs!(x::AbstractArray)
    @. x = abs(x)
end

function Base.:abs(x::AbstractArray)
    return abs.(x)
end

function abs!(x::Variable{T}) where T
    y = Variable{T}(abs(ᵛ(x)), x.backprop)
    if x.backprop
        function absBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* sign.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, absBackward)
    end
    return y
end

function Base.:abs(x::Variable{T}) where T
    y = Variable{T}(abs(ᵛ(x)), x.backprop)
    if x.backprop
        function absBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* sign.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, absBackward)
    end
    return y
end

function Base.:reshape(x::Variable{T}, newsize) where T
    y = Variable{T}( reshape(ᵛ(x), newsize), x.backprop )
    if x.backprop
        function reshapeBackward()
            if need2computeδ!(x)
                δ(x) .+= reshape(δ(y), x.shape)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, reshapeBackward)
    end
    return y
end

function exp2!(x::AbstractArray)
    @. x = exp2(x)
end

function Base.:exp2(x::AbstractArray)
    return exp2.(x)
end

function exp2!(x::Variable{T}) where T
    # exp2 represents y = 2^x
    y = Variable{T}(exp2!(ᵛ(x)), x.backprop)
    if x.backprop
        𝟚 = eltype(x)(2)
        function exp2Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* log(𝟚) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, exp2Backward)
    end
    return y
end

function Base.:exp2(x::Variable{T}) where T
    # EXP2 represents y = 2^x
    y = Variable{T}(exp2(ᵛ(x)), x.backprop)
    if x.backprop
        𝟚 = eltype(x)(2)
        function exp2Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* log(𝟚) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, exp2Backward)
    end
    return y
end

function exp10!(x::AbstractArray)
    @. x = exp10(x)
end

function Base.:exp10(x::AbstractArray)
    return exp10.(x)
end

function exp10!(x::Variable{T}) where T
    # EXP10 represents y = 10^x
    y = Variable{T}(exp10!(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙𝟘 = eltype(x)(10)
        function exp10Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* log(𝟙𝟘) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, exp10Backward)
    end
    return y
end

function Base.:exp10(x::Variable{T}) where T
    # EXP10 represents y = 10^x
    y = Variable{T}(exp10(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙𝟘 = eltype(x)(10)
        function exp10Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* log(𝟙𝟘) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, exp10Backward)
    end
    return y
end

function log2!(x::AbstractArray)
    @. x = log2(x)
end

function Base.:log2(x::AbstractArray)
    return log2.(x)
end

function log2!(x::Variable{T}) where T
    y = Variable{T}(log2(ᵛ(x)), x.backprop)
    if x.backprop
        𝟚 = eltype(x)(2)
        function log2Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (log(𝟚) .* ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, log2Backward)
    end
    return y
end

function Base.:log2(x::Variable{T}) where T
    y = Variable{T}(log2(ᵛ(x)), x.backprop)
    if x.backprop
        𝟚 = eltype(x)(2)
        function log2Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (log(𝟚) .* ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, log2Backward)
    end
    return y
end

function log10!(x::AbstractArray)
    @. x = log10(x)
end

function Base.:log10(x::AbstractArray)
    return log10.(x)
end

function log10!(x::Variable{T}) where T
    y = Variable{T}(log10(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙𝟘 = eltype(x)(10)
        function log10Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (log(𝟙𝟘) .* ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, log10Backward)
    end
    return y
end

function Base.:log10(x::Variable{T}) where T
    y = Variable{T}(log10(ᵛ(x)), x.backprop)
    if x.backprop
        𝟙𝟘 = eltype(x)(10)
        function log10Backward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (log(𝟙𝟘) .* ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, log10Backward)
    end
    return y
end

function sec!(x::Variable{T}) where T
    # SEC represents y = sec(x)
    y = Variable{T}(sec(ᵛ(x)), x.backprop)
    if x.backprop
        function secBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y) .* tan.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, secBackward)
    end
    return y
end

function Base.:sec(x::Variable{T}) where T
    # SEC represents y = sec(x)
    y = Variable{T}(sec(ᵛ(x)), x.backprop)
    if x.backprop
        function secBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* ᵛ(y) .* tan.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, secBackward)
    end
    return y
end

function sqrt!(x::Variable{T}) where T
    y = Variable{T}(sqrt!(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟚 = TOO(2.000)
        ϵ = TOO(1e-38)
        function sqrtBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (𝟚 .* (ᵛ(y) .+ ϵ))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sqrtBackward)
    end
    return y
end

function Base.:sqrt(x::Variable{T}) where T
    y = Variable{T}(sqrt(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟚 = TOO(2.000)
        ϵ = TOO(1e-38)
        function sqrtBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) ./ (𝟚 .* (ᵛ(y) .+ ϵ))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sqrtBackward)
    end
    return y
end

# -- tan serials --
export tan!
function tan!(x::Variable{T}) where T
    y = Variable{T}(tan!(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tanBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (𝟙 .+ ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tanBackward)
    end
    return y
end

function Base.:tan(x::Variable{T}) where T
    y = Variable{T}(tan(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tanBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (𝟙 .+ ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tanBackward)
    end
    return y
end

export tand!
function tand!(x::AbstractArray)
    @. x = tand(x)
end

function Base.:tand(x::AbstractArray)
    return tand.(x)
end

function tand!(x::Variable{T}) where T
    y = Variable{T}(tand!(ᵛ(x)), x.backprop)
    if x.backprop
        DPI = TOO(pi/180)
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tandBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* DPI .* (𝟙 .+ ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tandBackward)
    end
    return y
end

function Base.:tand(x::Variable{T}) where T
    y = Variable{T}(tand(ᵛ(x)), x.backprop)
    if x.backprop
        DPI = TOO(pi/180)
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tandBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* DPI .* (𝟙 .+ ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tandBackward)
    end
    return y
end

export tanh!
function tanh!(x::Variable{T}) where T
    y = Variable{T}(tanh!(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tanhBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (𝟙 .- ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tanhBackward)
    end
    return y
end

function Base.:tanh(x::Variable{T}) where T
    y = Variable{T}(tanh(ᵛ(x)), x.backprop)
    if x.backprop
        TOO = eltype(x)
        𝟙 = TOO(1.0)
        𝟚 = TOO(2.0)
        function tanhBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* (𝟙 .- ᵛ(y).^𝟚)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, tanhBackward)
    end
    return y
end

export tanhshrink, tanhshrink!
function tanhshrink!(x::AbstractArray)
    @. x = x - tanh(x)
end

function tanhshrink(x::AbstractArray)
    return  x - tanh(x)
end

function tanhshrink!(x::Variable{T}) where T
    return x - tanh(x)
end

function tanhshrink(x::Variable{T}) where T
    return x - tanh(x)
end

# # -- sin serials --
export sin!
function sin!(x::Variable{T}) where T
    y = Variable{T}(sin(ᵛ(x)), x.backprop)
    if x.backprop
        function sinBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* cos.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sinBackward)
    end
    return y
end

function Base.:sin(x::Variable{T}) where T
    y = Variable{T}(sin(ᵛ(x)), x.backprop)
    if x.backprop
        function sinBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* cos.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sinBackward)
    end
    return y
end

export sinc!
function sinc!(x::AbstractArray)
    @. x = sinc(x)
end

function Base.:sinc(x::AbstractArray)
    return sinc.(x)
end

function sinc!(x::Variable{T}) where T
    # sinc represents y = sin(pi*x)/(pi*x)
    y = Variable{T}(sinc(ᵛ(x)), x.backprop)
    if x.backprop
        function sincBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* cosc.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sincBackward)
    end
    return y
end

function Base.:sinc(x::Variable{T}) where T
    # sinc represents y = sin(pi*x)/(pi*x)
    y = Variable{T}(sinc(ᵛ(x)), x.backprop)
    if x.backprop
        function sincBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* cosc.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sincBackward)
    end
    return y
end

export sind!
function sind!(x::AbstractArray)
    @. x = sind(x)
end

function Base.:sind(x::AbstractArray)
    return sind.(x)
end

function sind!(x::Variable{T}) where T
    y = Variable{T}(sind(ᵛ(x)), x.backprop)
    if x.backprop
        DPI = eltype(x)(pi/180)
        function sindBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* DPI .* cosd.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sindBackward)
    end
    return y
end

function Base.:sind(x::Variable{T}) where T
    y = Variable{T}(sind(ᵛ(x)), x.backprop)
    if x.backprop
        DPI = eltype(x)(pi/180) # 1 rad⁻¹
        function sindBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* DPI .* cosd.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sindBackward)
    end
    return y
end

export sinpi!
function sinpi!(x::AbstractArray)
    @. x = sinpi(x)
end

function Base.:sinpi(x::AbstractArray)
    return sinpi.(x)
end

function sinpi!(x::Variable{T}) where T
    y = Variable{T}(sinpi(ᵛ(x)), x.backprop)
    if x.backprop
        𝝅 = eltype(x)(pi)
        function sinpiBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* 𝝅 .* cospi.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sinpiBackward)
    end
    return y
end

function Base.:sinpi(x::Variable{T}) where T
    y = Variable{T}(sinpi(ᵛ(x)), x.backprop)
    if x.backprop
        𝝅 = eltype(x)(pi)
        function sinpiBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y) .* 𝝅 .* cospi.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, sinpiBackward)
    end
    return y
end

export linearsin,linearsin!
function linearsin!(x::AbstractArray)
    @. x = sin(x) + x
end

function linearsin(x::AbstractArray)
    return sin(x) + x
end

function linearsin!(x::Variable{T}) where T
    return sin(x) + x
end

function linearsin(x::Variable{T}) where T
    return sin(x) + x
end

export cos!
function cos!(x::Variable{T}) where T
    y = Variable{T}(cos(ᵛ(x)), x.backprop)
    if x.backprop
        function cosBackward()
            if need2computeδ!(x)
                δ(x) .+= - δ(y) .* sin.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, cosBackward)
    end
    return y
end

function Base.:cos(x::Variable{T}) where T
    y = Variable{T}(cos(ᵛ(x)), x.backprop)
    if x.backprop
        function cosBackward()
            if need2computeδ!(x)
                δ(x) .+= - δ(y) .* sin.(ᵛ(x))
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, cosBackward)
    end
    return y
end

export inv!
function inv!(x::Variable{T}) where T
    y = Variable{T}(inv!(ᵛ(x)), x.backprop)
    if x.backprop
        function invBackward()
            if need2computeδ!(x)
                δ(x) .+= - δ(y) .* ᵛ(y) .* ᵛ(y);
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, invBackward)
    end
    return y
end

function inv(x::Variable{T}) where T
    y = Variable{T}(inv(ᵛ(x)), x.backprop)
    if x.backprop
        function invBackward()
            if need2computeδ!(x)
                δ(x) .+= - δ(y) .* ᵛ(y) .* ᵛ(y)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, invBackward)
    end
    return y
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using MultilinearOpt
using JuMP
using Gurobi
using Test
using Random
using LinearAlgebra
using MathOptInterface

const MOI = MathOptInterface

function test_constraint_violation(f, xbounds::NTuple{2}, ybounds::NTuple{2};
        method::Symbol, disc_level::Integer, rng::AbstractRNG, atol::Real, num_tests::Integer,
        sense::MOI.OptimizationSense)
    m = Model(with_optimizer(Gurobi.Optimizer, OutputFlag=0))
    @variable m x lower_bound=xbounds[1] upper_bound=xbounds[2]
    @variable m y lower_bound=ybounds[1] upper_bound=ybounds[2]
    @variable m z
    @constraint m z == f(x, y)
    relaxbilinear!(m, method=method, disc_level=disc_level)
    for i = 1 : num_tests
        if sense != MOI.FEASIBILITY_SENSE
            x̄ = xbounds[1] + rand(rng) * (xbounds[2] - xbounds[1])
            ȳ = ybounds[1] + rand(rng) * (ybounds[2] - ybounds[1])
            z̄ = f(x̄, ȳ) + rand(rng) - 0.5
            min_objective = sum(x -> x^2, [x, y, z] - [x̄, ȳ, z̄])
            if sense == MOI.MIN_SENSE
                @objective m Min min_objective
            else
                @objective m Max -min_objective
            end
        end
        optimize!(m)
        bounds_tol = 1e-8
        @test xbounds[1] - bounds_tol <= value(x) <= xbounds[2] + bounds_tol
        @test ybounds[1] - bounds_tol <= value(y) <= ybounds[2] + bounds_tol
        violation = f(value(x), value(y)) - value(z)
        # @show violation
        # @show value(x) value(y) value(z)
        # println()
        @test violation ≈ 0 atol=atol
    end
end

@testset "constraint violation" begin
    rng = MersenneTwister(1)
    atol = 1.1e-2
    num_tests = 10
    xbounds = (-0.5, 0.5)
    ybounds = (-0.4, 0.6)
    for (method, disc_level) in  [
            (:Logarithmic1D, 31),
            (:Logarithmic2D, 6),
            (:ZigZag1D, 31),
            (:ZigZag2D, 6),
            # (:Unary, 31) # broken
            (:MisenerLinear, 31),
            (:MisenerLog1, 31),
            (:MisenerLog2, 31)]

        for sense in [MOI.MIN_SENSE, MOI.MAX_SENSE, MOI.FEASIBILITY_SENSE]
            test_constraint_violation(*, xbounds, ybounds;
                method=method, disc_level=disc_level,
                atol=atol, rng=rng, num_tests=num_tests, sense=sense)

            test_constraint_violation((x, y) -> x * y + 3, xbounds, ybounds;
                method=method, disc_level=disc_level,
                atol=atol, rng=rng, num_tests=num_tests, sense=sense)

            test_constraint_violation((x, y) -> x^2, xbounds, ybounds;
                method=method, disc_level=disc_level,
                atol=atol, rng=rng, num_tests=num_tests, sense=sense)

            test_constraint_violation((x, y) -> x^2 + x * y - y, xbounds, ybounds;
                method=method, disc_level=disc_level,
                atol=2 * atol, rng=rng, num_tests=num_tests, sense=sense)
        end
    end
end


text:
###
### Particle Filtering and Particle MCMC Samplers.
###

#######################
# Particle Transition #
#######################

"""
    ParticleTransition{T, F<:AbstractFloat}

Fields:
- `θ`: The parameters for any given sample.
- `lp`: The log pdf for the sample's parameters.
- `le`: The log evidence retrieved from the particle.
- `weight`: The weight of the particle the sample was retrieved from.
"""
struct ParticleTransition{T, F<:AbstractFloat}
    θ::T
    lp::F
    le::F
    weight::F
end

function additional_parameters(::Type{<:ParticleTransition})
    return [:lp,:le, :weight]
end

DynamicPPL.getlogp(t::ParticleTransition) = t.lp

####
#### Generic Sequential Monte Carlo sampler.
####

"""
$(TYPEDEF)

Sequential Monte Carlo sampler.

# Fields

$(TYPEDFIELDS)
"""
struct SMC{space, R} <: ParticleInference
    resampler::R
end

"""
    SMC(space...)
    SMC([resampler = ResampleWithESSThreshold(), space = ()])
    SMC([resampler = resample_systematic, ]threshold[, space = ()])

Create a sequential Monte Carlo sampler of type [`SMC`](@ref) for the variables in `space`.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
function SMC(resampler = Turing.Core.ResampleWithESSThreshold(), space::Tuple = ())
    return SMC{space, typeof(resampler)}(resampler)
end

# Convenient constructors with ESS threshold
function SMC(resampler, threshold::Real, space::Tuple = ())
    return SMC(Turing.Core.ResampleWithESSThreshold(resampler, threshold), space)
end
SMC(threshold::Real, space::Tuple = ()) = SMC(resample_systematic, threshold, space)

# If only the space is defined
SMC(space::Symbol...) = SMC(space)
SMC(space::Tuple) = SMC(Turing.Core.ResampleWithESSThreshold(), space)

mutable struct SMCState{V<:VarInfo, F<:AbstractFloat} <: AbstractSamplerState
    vi                   ::   V
    # The logevidence after aggregating all samples together.
    average_logevidence  ::   F
    particles            ::   ParticleContainer
end

function SMCState(model::Model)
    vi = VarInfo(model)
    particles = ParticleContainer(Trace[])

    return SMCState(vi, 0.0, particles)
end

function Sampler(alg::SMC, model::Model, s::Selector)
    dict = Dict{Symbol, Any}()
    state = SMCState(model)
    return Sampler(alg, dict, s, state)
end

function AbstractMCMC.sample_init!(
    ::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:SMC},
    N::Integer;
    kwargs...
)
    # set the parameters to a starting value
    initialize_parameters!(spl; kwargs...)

    # reset the VarInfo
    vi = spl.state.vi
    reset_num_produce!(vi)
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)
    empty!(vi)

    # create a new set of particles
    T = Trace{typeof(spl),typeof(vi),typeof(model)}
    particles = T[Trace(model, spl, vi) for _ in 1:N]

    # create a new particle container
    spl.state.particles = pc = ParticleContainer(particles)

    # Perform particle sweep.
    logevidence = sweep!(pc, spl.alg.resampler)
    spl.state.average_logevidence = logevidence

    return
end

function AbstractMCMC.step!(
    ::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:SMC},
    ::Integer,
    transition;
    iteration=-1,
    kwargs...
)
    # check that we received a real iteration number
    @assert iteration >= 1 "step! needs to be called with an 'iteration' keyword argument."

    # grab the weight
    pc = spl.state.particles
    weight = getweight(pc, iteration)

    # update the master vi
    particle = pc.vals[iteration]
    params = tonamedtuple(particle.vi)

    # This is pretty useless since we reset the log probability continuously in the
    # particle sweep.
    lp = getlogp(particle.vi)

    return ParticleTransition(params, lp, spl.state.average_logevidence, weight)
end

####
#### Particle Gibbs sampler.
####

"""
$(TYPEDEF)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

# Fields

$(TYPEDFIELDS)
"""
struct PG{space,R} <: ParticleInference
    """Number of particles."""
    nparticles::Int
    """Resampling algorithm."""
    resampler::R
end

isgibbscomponent(::PG) = true

"""
    PG(n, space...)
    PG(n, [resampler = ResampleWithESSThreshold(), space = ()])
    PG(n, [resampler = resample_systematic, ]threshold[, space = ()])

Create a Particle Gibbs sampler of type [`PG`](@ref) with `n` particles for the variables
in `space`.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
function PG(
    nparticles::Int,
    resampler = Turing.Core.ResampleWithESSThreshold(),
    space::Tuple = (),
)
    return PG{space, typeof(resampler)}(nparticles, resampler)
end

# Convenient constructors with ESS threshold
function PG(nparticles::Int, resampler, threshold::Real, space::Tuple = ())
    return PG(nparticles, Turing.Core.ResampleWithESSThreshold(resampler, threshold), space)
end
function PG(nparticles::Int, threshold::Real, space::Tuple = ())
    return PG(nparticles, resample_systematic, threshold, space)
end

# If only the number of particles and the space is defined
PG(nparticles::Int, space::Symbol...) = PG(nparticles, space)
function PG(nparticles::Int, space::Tuple)
    return PG(nparticles, Turing.Core.ResampleWithESSThreshold(), space)
end

mutable struct PGState{V<:VarInfo, F<:AbstractFloat} <: AbstractSamplerState
    vi                   ::   V
    # The logevidence after aggregating all samples together.
    average_logevidence  ::   F
end

function PGState(model::Model)
    vi = VarInfo(model)
    return PGState(vi, 0.0)
end

const CSMC = PG # type alias of PG as Conditional SMC

"""
    Sampler(alg::PG, model::Model, s::Selector)

Return a `Sampler` object for the PG algorithm.
"""
function Sampler(alg::PG, model::Model, s::Selector)
    info = Dict{Symbol, Any}()
    state = PGState(model)
    return Sampler(alg, info, s, state)
end

function AbstractMCMC.step!(
    ::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:PG},
    ::Integer,
    transition;
    kwargs...
)
    # obtain or create reference particle
    vi = spl.state.vi
    ref_particle = isempty(vi) ? nothing : forkr(Trace(model, spl, vi))

    # reset the VarInfo before new sweep
    reset_num_produce!(vi)
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)

    # create a new set of particles
    num_particles = spl.alg.nparticles
    T = Trace{typeof(spl),typeof(vi),typeof(model)}
    if ref_particle === nothing
        particles = T[Trace(model, spl, vi) for _ in 1:num_particles]
    else
        particles = Vector{T}(undef, num_particles)
        @inbounds for i in 1:(num_particles - 1)
            particles[i] = Trace(model, spl, vi)
        end
        @inbounds particles[num_particles] = ref_particle
    end

    # create a new particle container
    pc = ParticleContainer(particles)

    # Perform a particle sweep.
    logevidence = sweep!(pc, spl.alg.resampler)

    # pick a particle to be retained.
    Ws = getweights(pc)
    indx = randcat(Ws)

    # extract the VarInfo from the retained particle.
    params = tonamedtuple(vi)
    spl.state.vi = pc.vals[indx].vi

    # This is pretty useless since we reset the log probability continuously in the
    # particle sweep.
    lp = getlogp(spl.state.vi)

    # update the master vi.
    return ParticleTransition(params, lp, logevidence, 1.0)
end

function AbstractMCMC.sample_end!(
    ::AbstractRNG,
    ::AbstractModel,
    spl::Sampler{<:ParticleInference},
    N::Integer,
    ts::Vector{<:ParticleTransition};
    resume_from = nothing,
    kwargs...
)
    # Exponentiate the average log evidence.
    # loge = exp(mean([t.le for t in ts]))
    loge = mean(t.le for t in ts)

    # If we already had a chain, grab the logevidence.
    if resume_from isa MCMCChains.Chains
        # pushfirst!(samples, resume_from.info[:samples]...)
        pre_loge = resume_from.logevidence
        # Calculate new log-evidence
        pre_n = length(resume_from)
        loge = (pre_loge * pre_n + loge * N) / (pre_n + N)
    elseif resume_from !== nothing
        error("keyword argument `resume_from` has to be `nothing` or a `MCMCChains.Chains` object")
    end

    # Store the logevidence.
    spl.state.average_logevidence = loge
end

function DynamicPPL.assume(
    rng,
    spl::Sampler{<:Union{PG,SMC}},
    dist::Distribution,
    vn::VarName,
    ::Any
)
    vi = current_trace().vi
    if inspace(vn, spl)
        if ~haskey(vi, vn)
            r = rand(rng, dist)
            push!(vi, vn, r, dist, spl)
        elseif is_flagged(vi, vn, "del")
            unset_flag!(vi, vn, "del")
            r = rand(rng, dist)
            vi[vn] = vectorize(dist, r)
            setgid!(vi, spl.selector, vn)
            setorder!(vi, vn, get_num_produce(vi))
        else
            updategid!(vi, vn, spl)
            r = vi[vn]
        end
    else # vn belongs to other sampler <=> conditionning on vn
        if haskey(vi, vn)
            r = vi[vn]
        else
            r = rand(rng, dist)
            push!(vi, vn, r, dist, Selector(:invalid))
        end
        lp = logpdf_with_trans(dist, r, istrans(vi, vn))
        acclogp!(vi, lp)
    end
    return r, 0
end

function DynamicPPL.observe(spl::Sampler{<:Union{PG,SMC}}, dist::Distribution, value, vi)
    produce(logpdf(dist, value))
    return 0
end

####
#### Resampling schemes for particle filters
####

# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# Default resampling scheme
function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
    return resample_systematic(w, num_particles)
end

# More stable, faster version of rand(Categorical)
function randcat(p::AbstractVector{<:Real})
    T = eltype(p)
    r = rand(T)
    s = 1
    for j in eachindex(p)
        r -= p[j]
        if r <= zero(T)
            s = j
            break
        end
    end
    return s
end

function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
    return rand(Distributions.sampler(Categorical(w)), num_particles)
end

function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)
    # Pre-allocate array for resampled particles
    indices = Vector{Int}(undef, num_particles)

    # deterministic assignment
    residuals = similar(w)
    i = 1
    @inbounds for j in 1:length(w)
        x = num_particles * w[j]
        floor_x = floor(Int, x)
        for k in 1:floor_x
            indices[i] = j
            i += 1
        end
        residuals[j] = x - floor_x
    end

    # sampling from residuals
    if i <= num_particles
        residuals ./= sum(residuals)
        rand!(Categorical(residuals), view(indices, i:num_particles))
    end

    return indices
end

"""
    resample_stratified(weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by stratified resampling.

In stratified resampling `n` ordered random numbers `u₁`, ..., `uₙ` are generated, where
``uₖ \\sim U[(k - 1) / n, k / n)``. Based on these numbers the samples `x₁`, ..., `xₙ`
are selected according to the multinomial distribution defined by the normalized `weights`,
i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_stratified(weights::AbstractVector{<:Real}, n::Integer)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]

    # generate all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # sample next `u` (scaled by `n`)
        u = oftype(v, i - 1 + rand())

        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample
    end

    return samples
end

"""
    resample_systematic(weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by systematic resampling.

In systematic resampling a random number ``u \\sim U[0, 1)`` is used to generate `n` ordered
numbers `u₁`, ..., `uₙ` where ``uₖ = (u + k − 1) / n``. Based on these numbers the samples
`x₁`, ..., `xₙ` are selected according to the multinomial distribution defined by the
normalized `weights`, i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_systematic(weights::AbstractVector{<:Real}, n::Integer)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]
    u = oftype(v, rand())

    # find all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample

        # update `u`
        u += one(u)
    end

    return samples
end


text:
#=

This file is the entry-point for the Absorption module.

It includes this module's source files and exports the relevant keywords.

=#

module Absorption

using SpecialFunctions          # For complex error functions
using Parameters                # For constructing HitranTable with keywords
using DocStringExtensions       # For simplifying docstring
using Interpolations            # For interpolating in lookup tables and interpolating CS
using JLD2                      # For saving and loading the interpolator
using ProgressMeter             # For showing progress, especially in creating interpolator
using KernelAbstractions        # For heterogeneous (GPU+CPU) programming
using CUDAKernels               # Access to CUDADevice
using CUDA                      # For GPU programming
using ForwardDiff, DiffResults  # For auto-differentiation
using NetCDF                    # For loading NetCDF files with constants
using ..Architectures           # For GPU/CPU convenience
using ..Architectures: CPU, GPU # Again for GPU/CPU convenience
import DataInterpolations: CubicSpline as DI_CS # For use in qoft

include("constants/constants.jl")                   # Scientific and mathematical constants
include("constants/mol_weights.jl")                 # Molecular weights
include("constants/TIPS_2017.jl")                   # Partition sums data
include("types.jl")                                 # All types used in this module
include("read_hitran.jl")                                # HITRAN file-related functions
include("complex_error_functions.jl")               # CEFs used in line broadening
include("make_model_helpers.jl")                  # CS interpolator functions
include("compute_absorption_cross_section.jl")      # Cross-section from HITRAN
include("autodiff_helper.jl")                       # Auto-differentiation

# Export the Cross Section models
export AbstractCrossSectionModel, HitranModel, InterpolationModel

# Export the absorption_cross_section functions
export compute_absorption_cross_section, absorption_cross_section

# Export the hitran functions from hitran.jl
export read_hitran, make_hitran_model

# Export the broadening function types
export Doppler, Lorentz, Voigt, AbstractBroadeningFunction

# Export the complex error functions
export AbstractComplexErrorFunction, HumlicekErrorFunction,
       HumlicekWeidemann32VoigtErrorFunction, HumlicekWeidemann32SDErrorFunction,
       CPF12ErrorFunction, ErfcHumliErrorFunctionVoigt, ErfcHumliErrorFunctionSD,
       ErfcErrorFunction

# Export the hitran table struct type
export HitranTable

# Export the interpolator functions
export make_interpolation_model, save_interpolation_model, load_interpolation_model

end


text:
using ProximalOperators
using ProximalOperators: Zero

ProximalOperators.Conjugate(_::Zero) = IndZero()
ProximalOperators.Conjugate(_::IndZero) = Zero()
ProximalOperators.Conjugate(f::SqrNormL2) = SqrNormL2(1.0 / f.lambda)

# TODO: Add other useful functions and calculus rules such as translation


text:
"""
    stack(x::NamedTuple) -> NamedTuple
    stack(x::AbstractArray{NamedTuple}) -> NamedTuple
    stack(x::AbstractArray{AbstractArray{<:NamedTuple}}) -> NamedTuple

Given a container of `NamedTuple`s, concatenate them, using the container dimensions as the
dimensions of the resulting arrays.

# Examples

```@example
using ArviZ
nchains, ndraws = 4, 100
data = [(x = rand(), y = randn(2), z = randn(2, 3)) for _ in 1:nchains, _ in 1:ndraws];
stacked_data = ArviZ.stack(data);
```
"""
stack(x) = x
stack(x::AbstractArray{T}) where {T<:Number} = Array(x)
stack(x::AbstractArray) = stack(stack.(x))
stack(x::NamedTuple) = (; (k => stack(v) for (k, v) in pairs(x))...)
function stack(x::AbstractArray{S}) where {T<:Number,N,S<:AbstractArray{T,N}}
    ret = Array{T}(undef, (size(x)..., size(x[1])...))
    @simd for k in keys(x)
        @inbounds setindex!(ret, x[k], k, (Colon() for _ in 1:N)...)
    end
    return ret
end
function stack(x::AbstractArray{<:NamedTuple{K}}) where {K}
    length(x) == 0 && return
    @inbounds x1 = x[1]
    ret = NamedTuple()
    for k in K
        v = replacemissing.(stack.(getproperty.(x, k)))
        ret = merge(ret, (k => stack(v),))
    end
    return ret
end

@doc doc"""
    from_namedtuple(posterior::NamedTuple; kwargs...) -> InferenceData
    from_namedtuple(posterior::Vector{<:NamedTuple}; kwargs...) -> InferenceData
    from_namedtuple(posterior::Matrix{<:NamedTuple}; kwargs...) -> InferenceData
    from_namedtuple(posterior::Vector{Vector{<:NamedTuple}}; kwargs...) -> InferenceData
    from_mcmcchains(
        posterior::NamedTuple,
        sample_stats::Any,
        posterior_predictive::Any,
        predictions::Any,
        log_likelihood::Any;
        kwargs...
    ) -> InferenceData

Convert a `NamedTuple` or container of `NamedTuple`s to an `InferenceData`.

If containers are passed, they are flattened into a single `NamedTuple` with array elements
whose first dimensions correspond to the dimensions of the containers.

# Arguments

- `posterior`: The data to be converted. It may be of the following types:
    + `::NamedTuple`: The keys are the variable names and the values are arrays with
        dimensions `(nchains, ndraws, sizes...)`.
    + `::Vector{<:NamedTuple}`: Each element is a `NamedTuple` from a chain with
        `Array`/`MonteCarloMeasurements.Particle` values with dimensions
        `(ndraws, sizes...)`.
    + `::Matrix{<:NamedTuple}`: Each element is a single draw from a single chain, with
        array/scalar values with dimensions `sizes`. The dimensions of the matrix container
        are `(nchains, ndraws)`
    + `::Vector{Vector{<:NamedTuple}}`: The same as the above case.

# Keywords

- `posterior_predictive::Any=nothing`: Draws from the posterior predictive distribution
- `sample_stats::Any=nothing`: Statistics of the posterior sampling process
- `predictions::Any=nothing`: Out-of-sample predictions for the posterior.
- `prior::Any=nothing`: Draws from the prior
- `prior_predictive::Any=nothing`: Draws from the prior predictive distribution
- `sample_stats_prior::Any=nothing`: Statistics of the prior sampling process
- `observed_data::Dict{String,Array}=nothing`: Observed data on which the `posterior` is
     conditional. It should only contain data which is modeled as a random variable. Keys
     are parameter names and values.
- `constant_data::Dict{String,Array}=nothing`: Model constants, data included in the model
     which is not modeled as a random variable. Keys are parameter names and values.
- `predictions_constant_data::Dict{String,Array}=nothing`: Constants relevant to the model
     predictions (i.e. new `x` values in a linear regression).
- `log_likelihood::Any=nothing`: Pointwise log-likelihood for the data. It is recommended
     to use this argument as a dictionary whose keys are observed variable names and whose
     values are log likelihood arrays.
- `library=nothing`: Name of library that generated the draws
- `coords::Dict{String,Vector}=nothing`: Map from named dimension to named indices
- `dims::Dict{String,Vector{String}}=nothing`: Map from variable name to names of its
     dimensions

# Returns

- `InferenceData`: The data with groups corresponding to the provided data

# Examples

```@example
using ArviZ
nchains, ndraws = 2, 10

data1 = (
    x = rand(nchains, ndraws),
    y = randn(nchains, ndraws, 2),
    z = randn(nchains, ndraws, 3, 2),
)
idata1 = from_namedtuple(data1)

data2 = [(x = rand(ndraws), y = randn(ndraws, 2), z = randn(ndraws, 3, 2)) for _ = 1:nchains];
idata2 = from_namedtuple(data2)

data3 = [(x = rand(), y = randn(2), z = randn(3, 2)) for _ = 1:nchains, _ = 1:ndraws];
idata3 = from_namedtuple(data3)

data4 = [[(x = rand(), y = randn(2), z = randn(3, 2)) for _ = 1:ndraws] for _ = 1:nchains];
idata4 = from_namedtuple(data4)
```
"""
from_namedtuple

function from_namedtuple(
    posterior,
    posterior_predictive,
    sample_stats,
    predictions,
    log_likelihood;
    library = nothing,
    kwargs...,
)
    all_idata = InferenceData()
    post_dict = posterior === nothing ? nothing : convert(Dict, posterior)
    for (group, group_data) in [
        :posterior_predictive => posterior_predictive,
        :sample_stats => sample_stats,
        :predictions => predictions,
        :log_likelihood => log_likelihood,
    ]
        group_data === nothing && continue
        if post_dict !== nothing
            if group_data isa Union{Symbol,String}
                group_data = [Symbol(group_data)]
            end
            if group_data isa Union{AbstractVector{Symbol},NTuple{N,Symbol} where {N}}
                group_data = popsubdict!(post_dict, group_data)
            end
            isempty(group_data) && continue
        end
        group_dataset = convert_to_dataset(group_data; kwargs...)
        if library !== nothing
            setattribute!(group_dataset, "inference_library", string(library))
        end
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    (post_dict === nothing || isempty(post_dict)) && return all_idata

    group_dataset = convert_to_dataset(post_dict; kwargs...)
    if library !== nothing
        setattribute!(group_dataset, "inference_library", string(library))
    end
    concat!(all_idata, InferenceData(posterior = group_dataset))

    return all_idata
end
function from_namedtuple(
    posterior::Union{NamedTuple,Nothing} = nothing;
    posterior_predictive = nothing,
    sample_stats = nothing,
    predictions = nothing,
    prior = nothing,
    prior_predictive = nothing,
    sample_stats_prior = nothing,
    observed_data = nothing,
    constant_data = nothing,
    predictions_constant_data = nothing,
    log_likelihood = nothing,
    library = nothing,
    kwargs...,
)
    all_idata = from_namedtuple(
        posterior,
        posterior_predictive,
        sample_stats,
        predictions,
        log_likelihood;
        library = library,
        kwargs...,
    )

    if any(x -> x !== nothing, [prior, prior_predictive, sample_stats_prior])
        pre_prior_idata = convert_to_inference_data(
            prior;
            posterior_predictive = prior_predictive,
            sample_stats = sample_stats_prior,
            library = library,
            kwargs...,
        )
        prior_idata = rekey(
            pre_prior_idata,
            Dict(
                :posterior => :prior,
                :posterior_predictive => :prior_predictive,
                :sample_stats => :sample_stats_prior,
            ),
        )
        concat!(all_idata, prior_idata)
    end

    for (group, group_data) in [
        :observed_data => observed_data,
        :constant_data => constant_data,
        :predictions_constant_data => predictions_constant_data,
    ]
        group_data === nothing && continue
        group_dict = convert(Dict, group_data)
        group_dataset =
            convert_to_constant_dataset(group_dict; library = library, kwargs...)
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    return all_idata
end
function from_namedtuple(data::AbstractVector{<:NamedTuple}; kwargs...)
    return from_namedtuple(stack(data); kwargs...)
end
function from_namedtuple(data::AbstractMatrix{<:NamedTuple}; kwargs...)
    return from_namedtuple(stack(data); kwargs...)
end
function from_namedtuple(data::AbstractVector{<:AbstractVector{<:NamedTuple}}; kwargs...)
    return from_namedtuple(stack(data); kwargs...)
end

"""
    convert_to_inference_data(obj::NamedTuple; kwargs...) -> InferenceData
    convert_to_inference_data(obj::Vector{<:NamedTuple}; kwargs...) -> InferenceData
    convert_to_inference_data(obj::Matrix{<:NamedTuple}; kwargs...) -> InferenceData
    convert_to_inference_data(obj::Vector{Vector{<:NamedTuple}}; kwargs...) -> InferenceData

Convert `obj` to an [`InferenceData`](@ref). See [`from_namedtuple`](@ref) for a description
of `obj` possibilities and `kwargs`.
"""
convert_to_inference_data(data::NamedTuple; kwargs...) = from_namedtuple(data; kwargs...)
function convert_to_inference_data(data::AbstractVector{<:NamedTuple}; kwargs...)
    return from_namedtuple(data; kwargs...)
end
function convert_to_inference_data(data::AbstractMatrix{<:NamedTuple}; kwargs...)
    return from_namedtuple(data; kwargs...)
end
function convert_to_inference_data(
    data::AbstractVector{<:AbstractVector{<:NamedTuple}};
    kwargs...,
)
    return from_namedtuple(data; kwargs...)
end


text:
module SimpleCarModels

using StaticArrays
using DifferentialDynamicsModels
using SpecialFunctions
import DifferentialDynamicsModels: mod2piF, adiff

include("math.jl")
include("models.jl")
include("dubins.jl")
include("reedsshepp.jl")
include("dubinsCC.jl")
include("elementary.jl")

end # module


text:
"""
Based on
https://hyperopt.github.io/hyperopt/
"""
module TestQuantLogUniform

using Test
using TreeParzen

@testset "quantised log uniform" begin

    qlu = HP.QuantLogUniform(:qlu, log(1.01), log(20.0), 2.0)
    N = 10_000

    qlu_samples = [TreeParzen.Resolve.node(qlu, TreeParzen.Trials.ValsDict()) for i in 1:N]

    @test all(qlu_samples .% 2 .== 0)
    @test minimum(qlu_samples) == 2
    @test maximum(qlu_samples) == 20

    # get the distinct samples, we already asserted they are all mod 2
    # and their max and min so we don't really need to check against a set
    sample_vals = sort(unique(qlu_samples))'
    counts = sum(sample_vals .== qlu_samples; dims=1)
    # reverse sort because we expect biggest first and smallest last
    ordered_indices = sortperm(dropdims(counts; dims=1); rev=true)

    # check that its always ordered from bigger to smaller;
    # it should be fairly rare for one to get bigger than one preceeding it
    # but we allow it by adding [[1, -1, 2]] to the diffcheck, which allows
    # elements to be swapped in 1 position only
    @test all(in.(diff(ordered_indices), [[1, -1, 2]]))

end

end # module TestQuantLogUniform
true


text:
module SingularSpectrumAnalysis

using LinearAlgebra, Statistics, RecipesBase, Requires

export hankel, hankelize, elementary, reconstruct, hsvd, autogroup, analyze

export sigmaplot, pairplot

export PredictionData, trendorder, fit_trend, pred, trend, seasons

@recipe function sigmaplot(USV::SVD; cumulative=false)
    seriestype := :scatter
    title --> "Normalized Singular Value Plot"
    ylabel --> "\$\\sigma / \\sum\\sigma_i\$"
    S = USV.S./sum(USV.S)
    cumulative ? cumsum(S) : S
end

@userplot Pairplot
"""
pairplot(USV, groupings)
Usage:
```julia
USV = hsvd(data,L)
seasonal_groupings = [1:2, 4:5]
pairplot(USV,seasonal_groupings)
```
"""
pairplot

@recipe function pairplot(h::Pairplot)
    USV, groupings = h.args[1:2]
    M = length(groupings)
    layout := M
    legend --> false
    for m = 1:M
        i = groupings[m]
        @assert length(i) == 2 "pairplot: All groupings have to be pairs"
        elements = USV.U[:,i].*sqrt.(USV.S[i])'
        subplot := m
        @series elements[:,1], elements[:,2]
    end
    nothing
end

"""
X = hankel(x,window_size)
Form the trajectory matrix `X` of size KxL, K = N-L+1
x can be a vector or a matrix For multivariate SSA
"""
function hankel(x,L)
    N = size(x,1)
    D = isa(x,AbstractVector) ? 1 : size(x,2)
    @assert L <= N/2 "L has to be less than N/2 = $(N/2)"
    K = N-L+1
    X = zeros(K,L*D)
    for d = 1:D
        for j = 1:L, i = 1:K
            k = i+j-1
            X[i,j+(d-1)*L] = x[k,d]
        end
    end
    X
end

"""
Ui = elementary(USV,I)
Computes the sum Uᵢ*Sᵢ*Vᵢ' for all i in I
If I is 1:L, the returned matrix is identical to U*S*V'
"""
function elementary(USV,I)
    sum(USV.U[:,i]*USV.S[i]*USV.Vt[i,:]' for i in I)
end

"""
means = hankelize(X)
Computes a timeseries from an approximate Hankel matrix by diagoal averaging (Hankelization). Note that the returned value is a vector and not a hankel matrix. A hankel matrix is easily obtainable by `hankel(hankelize(X), size(X,2))`
"""
function hankelize(X)
    K,L = size(X)
    @assert K >= L "The input matrix must be a tall matrix"
    N = K+L-1
    means = zeros(N)
    for n = 1:N
        if n <= K
            rangei = n:-1:max(1,n-L+1)
            rangej = 1:min(n,L)
        else
            rangei = K:-1:max(n-L+1,n-K+2)
            rangej = n-K+1:L
        end
        means[n] = mean(X[rangei[k],rangej[k]] for k in eachindex(rangei))
    end
    means
end
"""
USV = hsvd(y,L)
Form a trajectory hankel matrix from data `y` and compute svd on this Hankelmatrix
"""
function hsvd(y,L)
    X = hankel(y,L) # Form trajectory matrix
    USV = svd(X)
end

"""
yrt,yrs = reconstruct(USV::SVD, trends, seasonal::AbstractArray)
yr      = reconstruct(USV::SVD, groupings::AbstractArray)

Compute a reconstruction of the time-series based on an SVD object obtained from `hsvd` and user selected groupings. See also `?SingularSpectrumAnalysis`
"""
function reconstruct(USV, trends, seasonal::AbstractArray)
    if isa(trends,Number)
        trends = [trends]
    end
    yrt = reconstruct(USV, trends)
    yrs = reconstruct(USV, seasonal)
    size(yrt,2) == 1 && (yrt = vec(yrt))
    yrt, yrs
end

function reconstruct(USV, groupings::AbstractArray)
    M = length(groupings)
    K,L = size(USV.U)
    N = K+L-1
    yr = zeros(N,M)
    for m = 1:M
        X = elementary(USV,groupings[m])
        yr[:,m] = hankelize(X)
    end
    yr
end

"""
trend, seasonal_groupings = autogroup(USV, th = 0.95)
Try to automatically group singular values. `th ∈ (0,1)` determins the percentage of variance to explain.
"""
function autogroup(USV, th = 0.95)
    nS = USV.S .- minimum(USV.S)
    nS ./= sum(nS)
    cs = cumsum(nS)
    ind = findfirst(x->x > th, cs)
    iseven(ind) && (ind -= 1)
    t = findtrend(nS[1:ind])
    seasonal_groupings = UnitRange[]
    for i = [1:2:t-1; t+1:2:ind]
        push!(seasonal_groupings, i:i+1)
    end
    t, seasonal_groupings
end

function findtrend(S)
    function eval_trend(ind)
        iseven(ind) && return Inf
        cost = 0.0
        for i = [1:2:ind-1; ind+1:2:length(S)]
            cost += S[i] - S[i+1]
        end
        cost
    end
    costs = eval_trend.(1:length(S))
    _,ind = findmin(costs)
    ind
end

function analyze(y,L)
    USV = hsvd(y,L)
    trend, seasonal_groupings = autogroup(USV)
    reconstruct(USV, trend, seasonal_groupings)
end

mutable struct PredictionData
    trend_parameters
    trend_regressor
    seasonal_models
    seasonal_predictions
end
trendorder(pd) = length(pd.trend_parameters)-1

trend(pd::PredictionData) = pd.trend_regressor*pd.trend_parameters
trend(pd::PredictionData, i) = pd.trend_regressor[i,:]'*pd.trend_parameters
seasons(pd::PredictionData) = pd.seasonal_predictions
seasons(pd::PredictionData,i) = getindex.(pd.seasonal_predictions,i)
Base.getindex(pd::PredictionData) = trend(pd,i), seasons(pd,i)

@recipe function plotpd(pd::PredictionData)
    layout := 1+length(pd.seasonal_models)
    legend --> false
    @series begin
        subplot := 1
        title --> "Estimated trend"
        pd.trend_regressor*pd.trend_parameters
    end
    for (i,s) in enumerate(pd.seasonal_predictions)
        @series begin
            subplot := 1+i
            title --> "Seasonal predictions $i"
            s
        end
    end
end

"""
A,x = fit_trend(yt::AbstractVector, order=1)
Fit an n:th order polynomial to the trend signal `yt`
Returns the regressor matrix `A` and the coefficients `x` with the first coefficient being the 0:th order component (mean).
"""
function fit_trend(yt::AbstractVector, order=1)
    N = length(yt)
    A = (1:N).^(0:order)'
    x = A\yt
    A,x
end

function __init__()
    @require ControlSystemIdentification="3abffc1c-5106-53b7-b354-a47bfc086282" begin
    CSI = ControlSystemIdentification

    function PredictionData(yt,ys; trend_order=1, ar_order=2)
        A,x = fit_trend(yt, trend_order)  # Fits a polynomial of order 1 (line)
        ns = size(ys,2) # number of seasonal components
        models = [CSI.ar(1, ys[:,i], ar_order)[1] for i = 1:ns] # Fit one model per component
        # We can use these models to for one-step predictions
        seasonal_predictions = [CSI.predict(models[i], ys[:,i])  for i = 1:ns]
        # for s in seasonal_predictions
        #     for i = 1:ar_order insert!(s,1,0.) end
        # end
        PredictionData(x,A,models,seasonal_predictions)
    end

    function _pred(pd::PredictionData)
        A,x,ys,models = pd.trend_regressor, pd.trend_parameters, pd.seasonal_predictions, pd.seasonal_models
        @assert size(A,1) >= 3
        @assert (ys isa AbstractVecOrMat{<:Number} && size(ys,1) >= 3) || size(ys[1],1) >= 3
        ns = length(models)
        ysh = [CSI.predict(models[i], ys[i][end-2:end])  for i = 1:ns]
        ysh = reduce(hcat, ysh)
        yth = A[end,:]'*x
        yth, ysh
    end

    """
    pd = pred(pd::PredictionData, n=1)
    Form `n` step prediction. This requires you to have constructed a `PredictionData` object according to
    ```julia
    yt, ys = analyze(yn, L) # trend and seasons
    pd = fit_trend(yt, nt)  # Fits a polynomial of order nt (nt = 1 -> line)
    na = 2 # Autoregressive order
    # Next line fits one AR model for each seasonal component
    ns = size(ys,2) # number of seasonal components
    using ControlSystemIdentification
    models = [ar(1, ys[:,i], na)[1] for i = 1:ns] # Fit one model per component
    pd.seasonal_models = models
    seasonal_predictions = [predict(models[i], ys[:,i])  for i = 1:ns]
    pd.seasonal_predictions = seasonal_predictions
    ```
    """
    function pred(pd::PredictionData, n=1)
        n == 0 && (return pd)
        A,x,ys,models = pd.trend_regressor, pd.trend_parameters, pd.seasonal_predictions, pd.seasonal_models
        yth, ysh = _pred(pd)
        pd.trend_regressor = [A; (size(A,1)+1).^(0:trendorder(pd))']
        for i = eachindex(ys)
            pd.seasonal_predictions[i] = [ys[i]; ysh[i]]
        end
        pred(pd,n-1)
    end
end
end
end


text:
"""
    partition(s, tokenizer, blockifier, tokens; disable, postproc)

Go through a piece of text, either with an existing tokenization or an empty
one, tokenize if needed with the given tokenizer, blockify with the given
blockifier, and return a partition of the text into a vector of Blocks.

## Args

## KwArgs

    * disable:  list of token names to ignore (e.g. if want to allow math)
    * postproc: postprocessing to
"""
function partition(
            s::SS,
            tokenizer::Function,
            blockifier::Function;
            tokens::SubVector{Token}=EMPTY_TOKEN_SVEC,
            disable::Vector{Symbol}=Symbol[],
            postproc::Function=identity
            )::Vector{Block}

    parts = Block[]
    isempty(s) && return parts
    if isempty(tokens)
        tokens = tokenizer(s)
    end
    # NOTE: we need to be explicit here as in the recursive case, when partitioning
    # a block, there will not be a LR and EOS token, we're just getting the blocks'
    # inner token (and so checking the length of the tokens is insufficient)
    if length(tokens) == 2 && tokens[1].name == :LINE_RETURN && tokens[2].name == :EOS
        return [TextBlock(s)]
    end

    # disable tokens if desired
    isempty(disable) || filter!(t -> t.name ∉ disable, tokens)

    # form Blocks
    blocks = blockifier(tokens)
    # discard first block if it's a 0-length P_BREAK
    if !isempty(blocks) && iszero(to(blocks[1]))
        deleteat!(blocks, 1)
    end
    isempty(blocks) && return [TextBlock(s, tokens)]

    # disable additional blocks if desired
    isempty(disable) || filter!(t -> t.name ∉ disable, blocks)

    # Form a full partition with text blocks and blocks.
    parent = parent_string(s)
    first_block = blocks[1]
    last_block  = blocks[end]

    # add Text at beginning if first block is not there
    if from(s) < from(first_block)
        inter = subs(parent, from(s), prev_index(first_block))
        tb    = TextBlock(inter, tokens)
        push!(parts, tb)
    end

    # Go through blocks and add text with what's between them
    for i in 1:length(blocks)-1
        bi   = blocks[i]
        bip1 = blocks[i+1]
        push!(parts, blocks[i])
        inter = subs(parent, next_index(bi), prev_index(bip1))
        isempty(inter) || push!(parts, TextBlock(inter, tokens))
    end
    push!(parts, last_block)

    # add Text at the end if last block is not there
    if to(s) > to(last_block)
        inter = subs(parent, next_index(last_block), to(s))
        push!(parts, TextBlock(inter, tokens))
    end

    # Postprocessing (e.g. forming blockquotes, lists etc)
    return postproc(parts)
end
@inline partition(s::String, a...; kw...) = partition(subs(s), a...; kw...)
@inline partition(b::Block, a...; kw...)  = partition(content(b), a...; tokens=b.inner_tokens)

"""
    tokenizer_factory(; templates, postproc)

Arguments:
----------
    templates: a dictionary or matchers to find tokens.
    postproc: a function to apply on tokens after they've been found e.g. to merge
        them or filter them etc.

Returns:
--------
    A function that takes a string and returns a vector of tokens.
"""
@inline function tokenizer_factory(;
            templates::LittleDict=MD_TOKENS
            )::Function
    return s -> find_tokens(s, templates)
end

default_md_tokenizer   = tokenizer_factory()
default_math_tokenizer = tokenizer_factory(templates=MD_MATH_TOKENS)
default_html_tokenizer = tokenizer_factory(templates=HTML_TOKENS)

default_md_blockifier   = t -> find_blocks(subv(t), is_md=true)
default_html_blockifier = t -> find_blocks(subv(t), is_md=false)

@inline md_partition(e; kw...) =
    partition(e, default_md_tokenizer, default_md_blockifier;
              postproc=default_md_postproc!, kw...)

@inline math_partition(e; kw...) =
    partition(e, default_math_tokenizer, default_md_blockifier; kw...)

@inline html_partition(e; kw...) =
    partition(e, default_html_tokenizer, default_html_blockifier; kw...)

function default_md_postproc!(blocks::Vector{Block})
    form_blockquotes!(blocks)
    form_lists!(blocks)
    form_tables!(blocks)
    form_refs!(blocks)
    remove_inner!(blocks)
    return blocks
end

"""
    md_grouper(blocks)

Form begin-end spans keeping track of tokens and group text and inline blocks
after partition, this helps in forming paragraphs.
"""
function md_grouper(blocks::Vector{Block})::Vector{Group}

    groups   = Group[]
    cur_role = :NONE
    cur_head = 0
    i        = 1
    n_blocks = length(blocks)

    @inbounds while i <= n_blocks
        bi = blocks[i]
        br = bi.name in INLINE_BLOCKS ? :PARAGRAPH : bi.name

        if br != :PARAGRAPH
            _close_open_paragraph!(groups, blocks, cur_head, i)
            if br == :LX_BEGIN
                blocks[i+1].name == :CU_BRACKETS || enverr()
                env_name = content(blocks[i+1])
                # look ahead trying to find the proper closing \end{...}
                open_depth    = 1
                closing_index = -1
                probe = blocks[i+1]
                for j in i+1:n_blocks-1
                    cand  = probe.name
                    probe = blocks[j+1]
                    if (
                        cand == :LX_END &&
                        probe.name == :CU_BRACKETS &&
                        content(probe) == env_name
                        )
                        open_depth -= 1
                    elseif (
                            cand == :LX_BEGIN &&
                            probe.name == :CU_BRACKETS &&
                            content(probe) == env_name
                            )
                        open_depth += 1
                    end
                    if open_depth == 0
                        closing_index = j+1  # take the brace
                        break
                    end
                end
                closing_index == -1 && enverr("{$env_name}")
                # here we have a range of blocks from i -> closing which form a
                # complete begin{...} --> end{...}
                # it can't be empty as first and last are necessarily begin/end
                br = Symbol("ENV_$(env_name)")
                push!(groups, Group(blocks[i:closing_index]; role=br))
                # move the head to the end
                i = closing_index
            else
                push!(groups, Group(bi; role=br))
            end
            cur_head = 0
            cur_role = br

        elseif i == length(blocks)
            cur_head = ifelse(cur_head == 0, i, cur_head)
            _close_open_paragraph!(groups, blocks, cur_head, i+1)

        else
            cur_head = ifelse(cur_head == 0, i, cur_head)
        end
        i += 1
    end

    # finalise by removing P_BREAK
    filter!(g -> g.role != :P_BREAK, groups)
    return groups
end

function _close_open_paragraph!(groups, blocks, cur_head, i)
    if cur_head != 0
        push!(groups, Group(blocks[cur_head:i-1]; role=:PARAGRAPH))
    end
    return
end

enverr(n="") = parser_exception(:BlockNotClosed, """
    An opening \\begin$n was found but either not followed by a brace or
    left open.
    """)

"""
    split_args(s)

Take a string like 'foo "bar baz" 1' and return a string that is split along
whitespaces preserving quoted strings. So ["foo", "\"bar baz\"", "1"].
"""
function split_args(s::SS)
    # 1. find single-quoted / triply-quoted strings
    # 2. split the string outside of the quoted strings
    # 3. return the list of strings
    #
    # Ex:   foo "bar baz" 1
    #
    # expected output is ["foo", "\"bar baz\"", "1"]
    #
    # (specific parsing/processing is then left to the user, the
    # contract is that the user can then join(output, " ") and
    # get something equivalent in terms of how Franklin parses it)
    #
    parts = partition(
        s,
        _s -> find_tokens(_s, ARGS_TOKENS),
        _t -> (b = Block[]; _find_blocks!(b, subv(_t), ARGS_BLOCKS); b),
    )

    # form a dummy string with |__STR__| --> "foo |__STR__| 1"
    dummy   = IOBuffer()
    insert  = "|__STR__|"
    strings = String[]
    for p in parts
        if p.name == :TEXT
            write(dummy, content(p))
        else
            write(dummy, insert)
            push!(strings, content(p))
        end
    end
    # split the dummy string along white spaces
    splits = split(String(take!(dummy)))

    # reform the arguments in which |__STR__| appears.
    i = 1
    args = String[]
    for sp in splits
        if occursin(insert, sp)
            push!(args, replace(sp, insert => "\"" * strings[i] * "\""))
            i += 1
        else
            push!(args, sp)
        end
    end
    return args
end
split_args(s::String) = split_args(subs(s))