text:
raw"""
All the inputs are consumed and the difference of the output sums goes to the
miner

Input hashes and output index point to the original transaction
The unlocking script must match the locking script

Transaction:
  - Version
  - Input Counter
  - Inputs
    - hash
        The transaction origin id
    - output index
        The transacion output idx
    - unlocking script size
    - unlocking script
    - sequence number
  - Outputs
    - amount
    - locking script size
    - locking script

Hexdump of the genesis block area blk00000.dat:

f9be b4d9 1d01 0000 0100 0000 0000 0000  ................
^^^^^^^^^ magic bytes
          ^^^^^^^^^ block size
                    ^^^^^^^^^ version
                              ^^^^^^^^^
0000 0000 0000 0000 0000 0000 0000 0000  ................
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0000 0000 0000 0000 0000 0000 3ba3 edfd  ............;...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ previous hash
                              ^^^^^^^^^
7a7b 12b2 7ac7 2c3e 6776 8f61 7fc8 1bc3  z{..z.,>gv.a....
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
888a 5132 3a9f b8aa 4b1e 5e4a 29ab 5f49  ..Q2:...K.^J)._I
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ merkle root
                              ^^^^^^^^^ time stamp
ffff 001d 1dac 2b7c 0101 0000 0001 0000  ......+|........
^^^^^^^^^ difficulty
          ^^^^^^^^^ nonce
                    ^^ number of transactions (varint)
                      ^^^^^^^^^^ version
                                ^^ number of inputs
                                   ^^^^
0000 0000 0000 0000 0000 0000 0000 0000  ................
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0000 0000 0000 0000 0000 0000 0000 ffff  ................
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ originating transaction hash
                                   ^^^^
ffff 4d04 ffff 001d 0104 4554 6865 2054  ..M.......EThe T
^^^^ output index of the originating transaction
     ^^ script length (varint)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
696d 6573 2030 332f 4a61 6e2f 3230 3039  imes 03/Jan/2009
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2043 6861 6e63 656c 6c6f 7220 6f6e 2062   Chancellor on b
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
7269 6e6b 206f 6620 7365 636f 6e64 2062  rink of second b
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
6169 6c6f 7574 2066 6f72 2062 616e 6b73  ailout for banks
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ unlocking script
ffff ffff 0100 f205 2a01 0000 0043 4104  ........*....CA.
^^^^^^^^^ sequence number
          ^^ number of outputs
            ^^^^^^^^^^^^^^^^^^^^ amount (Big endian Float64)
                                ^^ script length
                                   ^^^^
678a fdb0 fe55 4827 1967 f1a6 7130 b710  g....UH'.g..q0..
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
5cd6 a828 e039 09a6 7962 e0ea 1f61 deb6  \..(.9..yb...a..
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
49f6 bc3f 4cef 38c4 f355 04e5 1ec1 12de  I..?L.8..U......
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
5c38 4df7 ba0b 8d57 8a4c 702b 6bf1 1d5f  \8M....W.Lp+k.._
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ac00 0000 00f9 beb4 d9d7 0000 0001 0000  ................
^^ locking script
  ^^^^^^^^^^ lock time
            ^^^^^^^^^^ magic bytes
006f e28c 0ab6 f1b3 72c1 a6a2 46ae 63f7  .o......r...F.c.
4f93 1e83 65e1 5a08 9c68 d619 0000 0000  O...e.Z..h......
0098 2051 fd1e 4ba7 44bb be68 0e1f ee14  .. Q..K.D..h....
677b a1a3 c354 0bf7 b1cd b606 e857 233e  g{...T.......W#>
0e61 bc66 49ff ff00 1d01 e362 9901 0100  .a.fI......b....
0000 0100 0000 0000 0000 0000 0000 0000  ................
0000 0000 0000 0000 0000 0000 0000 0000  ................
0000 00ff ffff ff07 04ff ff00 1d01 04ff  ................
ffff ff01 00f2 052a 0100 0000 4341 0496  .......*....CA..
"""

"""
    TransactionInput

Data Structure Storing Transaction Inputs
"""
struct TransactionInput
    # TODO: remove unlocking_script_size
    # the hash of the transaction where the output comes from
    hash                  :: UInt256
    output_index          :: UInt32
    unlocking_script_size :: UInt64
    unlocking_script      :: Vector{UInt8}
    sequence_number       :: UInt32
end

"""
    TransactionOutput

Data Structure Storing Transaction Outputs
"""
struct TransactionOutput
    # TODO: remove locking_script_size
    amount              :: UInt64 # in satoshis = 1e-8 bitcoins
    locking_script_size :: UInt64
    locking_script      :: Vector{UInt8}
end

"""
    Witness

Data Structure Storing The Witness for SegWit transactions
"""
struct Witness
    data :: Vector{Vector{UInt8}}
end

"""
    Transaction

Data Structure Storing Transactions

To get the hash of a transaction

```julia
double_sha256(tx)
```
"""
struct Transaction
    # TODO: remove *_counter
    version         :: UInt32
    marker          :: UInt8
    flag            :: UInt8
    input_counter   :: UInt64
    inputs          :: Vector{TransactionInput}
    output_counter  :: UInt64
    outputs         :: Vector{TransactionOutput}
    witness_counter :: UInt64
    witnesses       :: Vector{Witness}
    lock_time       :: UInt32
end

function TransactionInput(io)

    # TODO: Cannot read UInt256 directly from io here:
    in_hash = read!(io, Array{UInt256}(undef, 1))[1]
    # in_hash = read(io, UInt256)
    output_index = read(io, UInt32)
    unlocking_script_size = signed(read_varint(io))

    unlocking_script = read!(io, Array{UInt8}(undef, unlocking_script_size))
    sequence_number = read(io, UInt32)

    TransactionInput(
        in_hash,
        output_index,
        unlocking_script_size,
        unlocking_script,
        sequence_number
    )
end

function TransactionOutput(io)

    out_amount = read(io, UInt64)
    out_script_size = signed(read_varint(io))
    locking_script = read!(io, Array{UInt8}(undef, out_script_size))

    TransactionOutput(out_amount,
                      out_script_size,
                      locking_script)
end

function Witness(io::IO)
    n_items = read_varint(io)
    data = Vector{Vector{UInt8}}(undef, n_items)

    for i in 1:n_items
        l = read_varint(io)
        data[i] = read!(io, Array{UInt8, 1}(undef, l))
    end

    Witness(data)
end

function Transaction(io)

    version = read(io, UInt32)

    marker_or_n_in = read_varint(io)

    is_segwit = marker_or_n_in == zero(marker_or_n_in)
    if is_segwit
        marker = UInt8(marker_or_n_in)
        @assert marker == 0x00
        flag = read(io, UInt8)
        @assert flag == 0x01
        n_in = read_varint(io)
    else
        # these two don't matter
        marker = 0xff
        flag = 0xff
        n_in = marker_or_n_in
    end

    @assert n_in > 0
    inputs = TransactionInput[
        TransactionInput(io)
        for i in 1:n_in
    ]

    n_out = read_varint(io)
    @assert n_out > 0
    outputs = TransactionOutput[
        TransactionOutput(io)
        for i in 1:n_out
    ]

    if is_segwit
        n_witness = n_in
        @assert n_witness > 0
        witnesses = Witness[
            Witness(io)
            for i in 1:n_witness
        ]
    else
        n_witness = zero(UInt64)
        witnesses = Witness[]
    end

    locktime = read(io, UInt32)

    Transaction(version,
                marker,
                flag,
                n_in,
                inputs,
                n_out,
                outputs,
                n_witness,
                witnesses,
                locktime)
end

# TODO: is this the same as bytes2hex?
hexarray(x::Array{UInt8}) = mapreduce(x -> string(x, base = 16), *, x)

function showcompact(io::IO, input::TransactionInput)
    # TODO: do something useful here
    println(io, "")
end

function Base.show(io::IO, input::TransactionInput)
    if !get(io, :compact, false)
        println(io, "Transaction input:")
        println(io, "  Hash:                  " * string(input.hash,            base = 16))
        println(io, "  Output index:          " * string(input.output_index,    base = 10))
        println(io, "  Input Sequence:        " * string(input.sequence_number, base = 10))
    end
end

# function Base.showall(io::IO, input::TransactionInput)
#     println(io, "Transaction input:")
#     println(io, "  Hash:                  " * input.hash)
#     println(io, "  Output index:          " * input.output_index)
#     println(io, "  Unlocking script size: " * input.unlocking_script_size)
#     println(io, "  Unlocking script:      " * hexarray(input.unlocking_script))
#     println(io, "  Input Sequence:        " * input.sequence_number)
# end

function Base.show(io::IO, output::TransactionOutput)
    println(io, "Transaction output: " * string(output.amount, base = 10))
end
# function Base.showall(io::IO, output::TransactionOutput)
#     println(io, "Transaction output: "    * string(output.amount,              base = 10))
#     println(io, "  Locking script size: " * string(output.locking_script_size, base = 10))
#     println(io, "  Locking script:      " * hexarray(output.locking_script))
# end

function showcompact(io::IO, transaction::Transaction)
    # TODO: add id here
    println(io, "Transaction: ")
end

function Base.show(io::IO, transaction::Transaction)
    # TODO: add id here
    if !get(io, :compact, false)
        println(io, "Transaction: ")
        println(io, "  Version:        " * string(transaction.version, base = 10))
        println(io, "  Input counter:  " * string(transaction.input_counter, base = 10))
        println(io, "  Output counter: " * string(transaction.output_counter, base = 10))
        println(io, "  Lock time:      " * string(transaction.lock_time, base = 10))
    end
end

# function Base.showall(io::IO, transaction::Transaction)
#     # TODO: add id here
#     println(io, "Transaction: ")
#     println(io, "  Version:        " * string(transaction.version, base = 10))
#     println(io, "  Input counter:  " * string(transaction.input_counter, base = 10))
#     for i in 1:transaction.input_counter
#         show(transaction.inputs[i])
#     end
#     println(io, "  Output counter: " * string(transaction.output_counter, base = 10))
#     for i in 1:transaction.output_counter
#         show(transaction.outputs[i])
#     end
#     println(io, "  Lock time:      " * string(transaction.lock_time, base = 10))
# end

# TODO: little endian only:
function sha256(tx::Transaction)

    ctx = SHA.SHA256_CTX()

    SHA.update!(ctx, to_byte_tuple(tx.version))
    SHA.update!(ctx, to_varint(tx.input_counter))
    for i in 1:tx.input_counter
        SHA.update!(ctx, to_byte_tuple(tx.inputs[i].hash))
        SHA.update!(ctx, to_byte_tuple(tx.inputs[i].output_index))
        SHA.update!(ctx, to_varint(tx.inputs[i].unlocking_script_size))
        SHA.update!(ctx, tx.inputs[i].unlocking_script)
        SHA.update!(ctx, to_byte_tuple(tx.inputs[i].sequence_number))
    end
    SHA.update!(ctx, to_varint(tx.output_counter))
    for i in 1:tx.output_counter
        SHA.update!(ctx, to_byte_tuple(tx.outputs[i].amount))
        SHA.update!(ctx, to_varint(tx.outputs[i].locking_script_size))
        SHA.update!(ctx, tx.outputs[i].locking_script)
    end
    SHA.update!(ctx, to_byte_tuple(tx.lock_time))

    return SHA.digest!(ctx)
end

function double_sha256(tx::Transaction)::UInt256
    tx |> sha256 |> sha256 |> x -> reinterpret(UInt256, x)[1]
end

# TODO: copyto! is save, right?
# TODO: find an allocation free version for these: copyto!(x, i, to_byte_tuple(y))
"""
mutates data in place, advances idx and returns it.
"""
function dump_txin_data!(data, idx, txin::TransactionInput)
    copyto!(data, idx, to_byte_tuple(txin.hash))
    idx += sizeof(txin.hash)

    copyto!(data, idx, to_byte_tuple(txin.output_index))
    idx += sizeof(txin.output_index)

    unlocking_script_size = to_varint(txin.unlocking_script_size)
    copyto!(data, idx, unlocking_script_size)
    idx += sizeof(unlocking_script_size)

    copyto!(data, idx, txin.unlocking_script)
    idx += sizeof(txin.unlocking_script)

    copyto!(data, idx, to_byte_tuple(txin.sequence_number))
    idx += sizeof(txin.sequence_number)

    return idx
end

# TODO: copyto! is save, right?
# TODO: find an allocation free version for these: copyto!(x, i, to_byte_tuple(y))
"""
mutates data in place, advances idx and returns it.
"""
function dump_txout_data!(data, idx, txout::TransactionOutput)
    copyto!(data, idx, to_byte_tuple(txout.amount))
    idx += sizeof(txout.amount)

    locking_script_size = to_varint(txout.locking_script_size)
    copyto!(data, idx, locking_script_size)
    idx += sizeof(locking_script_size)

    copyto!(data, idx, txout.locking_script)
    idx += sizeof(txout.locking_script)

    return idx
end

# TODO: copyto! is save, right?
# TODO: find an allocation free version for these: copyto!(x, i, to_byte_tuple(y))
"""
mutates data in place, advances idx and returns it.
"""
function dump_tx_data!(data, idx, tx::Transaction)
    copyto!(data, idx, to_byte_tuple(tx.version))
    idx += sizeof(tx.version)

    is_segwit = tx.marker == 0x00
    if is_segwit
        data[idx] = tx.marker
        idx += 1
        data[idx] = tx.flag
        idx += 1
    end

    input_counter = to_varint(tx.input_counter)
    copyto!(data, idx, input_counter)
    idx += sizeof(input_counter)
    for i in 1:tx.input_counter
        idx = dump_txin_data!(data, idx, tx.inputs[i])
    end

    output_counter = to_varint(tx.output_counter)
    copyto!(data, idx, output_counter)
    idx += sizeof(output_counter)
    for i in 1:tx.output_counter
        idx = dump_txout_data!(data, idx, tx.outputs[i])
    end

    if is_segwit
        for i in 1:tx.input_counter
            n_items = length(tx.witnesses[i].data)
            n_items_bytes = to_varint(n_items)
            copyto!(data, idx, n_items_bytes)
            idx += length(n_items_bytes)

            for j in 1:n_items
                l = length(tx.witnesses[i].data[j])
                l_bytes = to_varint(l)
                copyto!(data, idx, l_bytes)
                idx += length(l_bytes)

                copyto!(data, idx, tx.witnesses[i].data[j])
                idx += l
            end
        end
    end

    copyto!(data, idx, to_byte_tuple(tx.lock_time))
    idx += sizeof(tx.lock_time)

    return idx
end

# TODO: common abstraction with "Link"?
struct FilePointer
    file_number   :: UInt64
    file_position :: UInt64
end

function Transaction(fp::FilePointer)
    fp.file_number |>
        get_block_chain_file_path |>
        x -> open(x) do fh
            seek(fh, fp.file_position)
            Transaction(fh)
        end
end

total_output(tx::Transaction) = sum(x -> x.amount, tx.outputs)

is_coinbase(txin::TransactionInput) = txin.hash == zero(UInt256)


text:
using DiffEqBase: remake, AbstractODEProblem
using ...Bifurcations: FixedPointBifurcationProblem, DiffEqWrapper
using ..ArrayUtils: canonicalize
using ..Codim1
using ..Codim1: Codim1Solver, AbstractSpecialPoint, resolve_point

abstract type Codim2Problem{skind, tkind} <: BifurcationProblem{skind, tkind} end

"""
    BifurcationProblem(point::AbstractSpecialPoint,
                       solver::Codim1Solver,
                       param_axis2::Lens,
                       t2_domain::Tuple)

Construct codimension-2 bifurcation problem given a bifurcation
`point`.
"""
function BifurcationProblem(point::AbstractSpecialPoint,
                            solver::Codim1Solver,
                            param_axis2::Lens,
                            t2_domain::Tuple;
                            kwargs...)

    cd1_prob = solver.prob :: FixedPointBifurcationProblem
    de_wrapper = cd1_prob.p :: DiffEqWrapper
    de_prob = de_wrapper.de_prob :: AbstractODEProblem

    param_axis1 = de_wrapper.param_axis
    t_domain = (
        (cd1_prob.t_domain[1], t2_domain[1]),
        (cd1_prob.t_domain[2], t2_domain[2]),
    )

    # Manually cast.  The fact that `resolved.u` is not of type xtype
    # indicates type instability somewhere.  Until it is fixed, let's
    # silently cast always.  After type stability is established, it
    # has to be changed to emit warning before cast.
    xtype = typeof(de_prob.u0)

    resolved = resolve_point(point, solver)
    x0 = xtype(resolved.u[1:end - 1])
    t0 = SVector(resolved.u[end], get(de_prob.p, param_axis2))

    vals, vecs = eigen(Codim1.ds_jacobian(resolved.J))
    # TOOD: use eigs (depending on size(J)?)
    if timekind(point) isa Discrete
        val, idx = findmax(@. abs(abs(vals) - 1))
        # @assert val ≈ 1
    else
        val, idx = findmin(@. abs(real(vals)))
        # @assert val ≈ 0
    end
    v0 = cast_container(xtype, vecs[:, idx])
    v0 = canonicalize(v0)
    w0 = imag(vals[idx])
    if eltype(v0) <: Complex && w0 < 0
        w0 = -w0
        v0 = conj(v0)
    end

    # Sanity checks:
    if point.point_type == Codim1.PointTypes.saddle_node
        @assert eltype(v0) <: Real
    elseif point.point_type == Codim1.PointTypes.hopf
        @assert eltype(v0) <: Complex
    end
    if xtype <: SVector
        v0 :: SVector
    end

    return DiffEqCodim2Problem(
        de_prob,
        param_axis1,
        param_axis2,
        t_domain;
        x0 = x0,
        v0 = v0,
        w0 = w0,
        t0 = t0,
        augmented_system = preferred_augsys(point),
        kwargs...)
end

ext: jl


text:
"""
 WeightSpace represents a space that weights another space.
 Overload weight(S,x).
"""
abstract type WeightSpace{S,DD,RR} <: Space{DD,RR} end

domain(S::WeightSpace) = domain(S.space)

points(sp::WeightSpace,n) = points(sp.space,n)

struct WeightSpacePlan{S,P,T,V}
    space::S
    plan::P
    points::Vector{T}
    weights::Vector{V}
end

struct IWeightSpacePlan{S,P,T,V}
    space::S
    plan::P
    points::Vector{T}
    weights::Vector{V}
end

function plan_transform(S::WeightSpace,vals::AbstractVector)
    pts=points(S,length(vals))
    WeightSpacePlan(S,plan_transform(S.space,vals),pts,weight.(Ref(S),pts))
end

function plan_itransform(S::WeightSpace,vals::AbstractVector)
    pts=points(S,length(vals))
    IWeightSpacePlan(S,plan_itransform(S.space,vals),pts,weight.(Ref(S),pts))
end

*(P::WeightSpacePlan,vals::AbstractVector) = P.plan*(vals./P.weights)
*(P::IWeightSpacePlan,cfs::AbstractVector) = P.weights.*(P.plan*cfs)

# used for ProductFun
transform(sp::WeightSpace,vals::AbstractVector,plan::WeightSpacePlan) =
    transform(sp.space,vals./(sp==plan.space ? plan.weights : weight.(sp,plan.points)),plan.plan)
itransform(sp::WeightSpace,cfs::AbstractVector,plan::WeightSpacePlan) =
    itransform(sp.space,cfs,plan.plan).*(sp==plan.space ? plan.weights : weight.(sp,plan.points))

evaluate(f::AbstractVector,S::WeightSpace,x) = weight(S,x)*evaluate(f,S.space,x)

function evaluate(f::AbstractVector,S::WeightSpace,x...)
    fv=evaluate(f,S.space,x...)
    weight(S,x...).*fv
end

# recurrence is inhereted
for FUNC in (:recα,:recβ,:recγ)
    @eval $FUNC(T,ws::WeightSpace,k) = $FUNC(T,ws.space,k)
end

include("divide_singularity.jl")
include("JacobiWeight.jl")
include("JacobiWeightOperators.jl")
include("JacobiWeightChebyshev.jl")
include("HeavisideSpace.jl")
include("DiracSpace.jl")
include("LogWeight.jl")
include("ExpWeight.jl")


text:
# The order of the numerator polynomial in z = \omega^2 that we need to
# see if it has any zeros:
using Polynomials

function sturms_theorem(x)
n_lor = round(Int, length(x)/4)
pord = 2*(n_lor-1) + 1
# Compute coefficients in the numerator & denominator of PSD:
# The PSD is given by:
#
# \sum_i (a_i z +b_i)/(z^2 + c_i*z + d_i)
#
# where z = \omega^2
#
# We are computing the number of positive real zeros of this equation.
#

# Now, loop over coefficients:

a = zeros(n_lor)
b = zeros(n_lor)
c = zeros(n_lor)
d = zeros(n_lor)

for i=1:n_lor
  ar_j = x[(i-1)*4+1]
  ai_j = x[(i-1)*4+2]
  br_j = x[(i-1)*4+3]
  bi_j = x[(i-1)*4+4]
  a[i] = ar_j*br_j-ai_j*bi_j
  b[i] = (bi_j^2+br_j^2)*(ai_j*bi_j+ar_j*br_j)
  c[i] = 2*(br_j^2-bi_j^2)
  d[i] = (bi_j^2+br_j^2)^2
end

# Initialize a polynomial:
p0 = Poly(zeros(pord+1))
for i=1:n_lor
# The polynomial for the current Lorentzian term in the common-denominator expression:
  pcur = Poly([b[i],a[i]])
  for j=1:n_lor
# Only multiply by the denominators from the other Lorentzians:
    if j != i
      pcur *= Poly([d[j],c[j],1])
    end
  end
  p0 += pcur
end

# Compute the roots of this polynomial:
#poly_root = roots(p0)

# Now that we've computed coefficients of the polynomial, we just need
# to apply Sturm's theorem!

# Set up an array to hold signs of polynomial at zero & infinity:
f_of_0 = zeros(pord+1)
f_of_inf = zeros(pord+1)

# Take the derivative of the polynomial:
p1 = polyder(p0)
#println(p0)
#println(p1)
# Insert the coefficients of the z^0 term:
f_of_0[1] = p0(0)
f_of_0[2] = p1(0)
# Insert the coefficient of the z^(p_ord-i) term:
f_of_inf[1] = p0[pord]
f_of_inf[2] = p1[pord-1]

# Now, loop over the Sturm polynomial series:
for i=3:pord+1
  p2 = -rem(p0,p1)
# Check that round-off error hasn't left us with a polynomial
# that is the same order as p0 or p1, but with a small coefficient:
  if length(p2) >= length(p1)
    coeff = zeros(pord-i+2)
    for j=0:pord-i+1
      coeff[j+1]=p2[j]
    end
    p2 = Poly(coeff)
  end
# Insert the z^0 term:
  f_of_0[i]=p2(0)
# Insert the z -> \infty term:
  f_of_inf[i]=p2[pord-i+1]
# Now move promote these polynomials, readying them for recursion in the next step:
  p0=copy(p1)
  p1=copy(p2)
end

# Now we'll compute the number of sign changes at z=0:
sig_0 = 0
for i=1:pord
#  if f_of_0[i+1]*f_of_0[i] < 0 && abs(f_of_0[i+1]) > eps() && abs(f_of_0[i]) > eps()
  if sign(f_of_0[i+1]) != sign(f_of_0[i])
    sig_0 +=1
  end
end
# Next, compute the number of sign changes at z=\infty:
sig_inf = 0
for i=1:pord
#  if f_of_inf[i+1]*f_of_inf[i] < 0 && abs(f_of_inf[i+1]) > eps() && abs(f_of_inf[i]) > eps()
  if sign(f_of_inf[i+1]) != sign(f_of_inf[i])
    sig_inf +=1
  end
end

# The difference between z=0 & z=\infty gives the number of positive, real roots:
#println("Number of positive, real roots: ",sig_0-sig_inf)

# Print out the actual roots for inspection:
#println("Roots: ",poly_root)

# Now, determine which roots are positive & real from the numerial solver:
n_pos_real = sig_0-sig_inf
return n_pos_real
end


text:
using EffectiveWaves, Test

# this is for high volume and frequency
@testset "compare wienger-hopf and match method reflection" begin

    medium = Medium(1.0,1.0+0.0im)

    ω = 1.3

    specie = Specie(ρ=0.0, r=0.5, c=0.0, volfrac = 0.2)
    radius_multiplier = 1.01

    tol = 1e-9
    hankel_order=0

    ws = effective_waves(ω, medium, [specie];
        radius_multiplier = radius_multiplier,
        tol=tol, hankel_order = hankel_order,
        apply_meshing = false,
        num_wavenumbers = 20,
        mesh_points = 30, mesh_size = 2.);

    function Rerror(θ)

        match_ws = MatchWave(ω, medium, specie;
            radius_multiplier = radius_multiplier,
            hankel_order = hankel_order,
            tol = tol, wave_effs = ws[1:15],
            θin=θ,
            max_size = 600,
        );

        Rm = reflection_coefficient(ω, match_ws, medium, specie; θin = θ)
        Rw = wienerhopf_reflection_coefficient(ω, medium, [specie];
                θin = θ,
                tol=tol,
                radius_multiplier = radius_multiplier,
                hankel_order = hankel_order
        )
        return abs(Rm-Rw)/abs(Rw)
    end

    Rerror(0.0) # 3.771528229242978e-5
    Rerror(0.5) # 5.2531420000213484e-5
    Rerror(1.4) # 0.0005018041373008065

    @test Rerror(0.0) < 1e-4
    @test Rerror(0.5) < 1e-4
    @test Rerror(1.4) < 1e-3

end

# Note the wiener-hopf method is inaccurate or takes along time to calculate a very low-frequency solution
@testset "compare wienger-hopf for low-ish frequency" begin

    medium = Medium(1.0,1.0+0.0im)
    ω = 0.01
    hankel_order=0

    # Always Dirichlet boundary conditions for hankel_order=0!
    specie = Specie(ρ = 0.0, r=0.01, c=6.0, volfrac = 0.3)

    tol = 1e-7

    eff_medium = effective_medium(medium, [specie])

    function Rerror(θ)
        ws = effective_waves(ω, medium, [specie];
            tol=tol, θin = θ,
            hankel_order = hankel_order,
            num_wavenumbers = 2,
            extinction_rescale = true);
        R1 = reflection_coefficient(ω, ws[1], medium, [specie]; θin = θ)
        R_low = reflection_coefficient_halfspace(medium, eff_medium; θin = θ)
        Rw = wienerhopf_reflection_coefficient(ω, medium, [specie];
                θin = θ,
                tol=tol,
                hankel_order = hankel_order,
                num_coefs = 20000 + 200*Int(round(100*θ))
        )
        return abs(2*R_low-Rw-R1)/abs(2*R_low)
    end

    @test Rerror(0.5) < 1e-3
    @test Rerror(1.0) < 1e-3
end


text:
function Base.get{S<:VehicleState,D,I,R}(::Feature_PosFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(rec[pastframe][vehicle_index].state.posF.t)
end
function Base.get{S<:VehicleState,D,I,R}(::Feature_PosFyaw, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(rec[pastframe][vehicle_index].state.posF.ϕ)
end
function Base.get{S<:VehicleState,D,I,R}(::Feature_Speed, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    FeatureValue(rec[pastframe][vehicle_index].state.v)
end
function Base.get{S<:VehicleState,D,I,R}(::Feature_VelFs, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    FeatureValue(veh.state.v*cos(veh.state.posF.ϕ))
end
function Base.get{S<:VehicleState,D,I,R}(::Feature_VelFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    FeatureValue(veh.state.v*sin(veh.state.posF.ϕ))
end

generate_feature_functions("TurnRateG", :turnrateG, Float64, "rad/s")
function Base.get{S<:VehicleState,D,I,R}(::Feature_TurnRateG, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0; frames_back::Int=1)

    id = rec[pastframe][vehicle_index].id

    retval = FeatureValue(0.0, FeatureState.INSUF_HIST)
    pastframe2 = pastframe - frames_back
    if pastframe_inbounds(rec, pastframe2)

        veh_index_curr = vehicle_index
        veh_index_prev = findfirst(rec[pastframe2], id)

        if veh_index_prev != 0
            curr = rec[pastframe][veh_index_curr].state.posG.θ
            past = rec[pastframe2][veh_index_prev].state.posG.θ
            Δt = get_elapsed_time(rec, pastframe2, pastframe)
            retval = FeatureValue(deltaangle(past, curr) / Δt)
        end
    end

    retval
end
generate_feature_functions("TurnRateF", :turnrateF, Float64, "rad/s")
function Base.get{S<:VehicleState,D,I,R}(::Feature_TurnRateF, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    _get_feature_derivative_backwards(POSFYAW, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("AngularRateG", :angrateG, Float64, "rad/s²")
function Base.get{S<:VehicleState,D,I,R}(::Feature_AngularRateG, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    _get_feature_derivative_backwards(TURNRATEG, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("AngularRateF", :angrateF, Float64, "rad/s²")
function Base.get{S<:VehicleState,D,I,R}(::Feature_AngularRateF, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0)
    _get_feature_derivative_backwards(TURNRATEF, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("DesiredAngle", :desang, Float64, "rad")
function Base.get{S<:VehicleState,D,I,R}(::Feature_DesiredAngle, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0;
    kp_desired_angle::Float64 = 1.0,
    )

    retval = FeatureValue(0.0, FeatureState.INSUF_HIST)
    if pastframe_inbounds(rec, pastframe) && pastframe_inbounds(rec, pastframe-1)

        id = rec[pastframe][vehicle_index].id

        pastϕ = get_state(rec, id, pastframe-1).posF.ϕ
        currϕ = get_state(rec, id, pastframe).posF.ϕ

        Δt = rec.timestep
        expconst = exp(-kp_desired_angle*Δt)
        retval = FeatureValue((currϕ - pastϕ*expconst) / (1.0 - expconst))
    end
    retval
end

generate_feature_functions("MarkerDist_Left", :d_ml, Float64, "m")
Base.get(::Feature_MarkerDist_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0) =
    FeatureValue(get_markerdist_left(rec[pastframe][vehicle_index], roadway))

generate_feature_functions("MarkerDist_Right", :d_mr, Float64, "m")
Base.get(::Feature_MarkerDist_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0) =
    FeatureValue(get_markerdist_right(rec[pastframe][vehicle_index], roadway))

generate_feature_functions("MarkerDist_Left_Left", :d_mll, Float64, "m", can_be_missing=true)
function Base.get(::Feature_MarkerDist_Left_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    #=
    Distance to the left lane marker one lane to the left
    =#

    veh = rec[pastframe][vehicle_index]
    lane = roadway[veh.state.posF.roadind.tag]
    if n_lanes_left(lane, roadway) > 0
        offset = veh.state.posF.t
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        FeatureValue(lane.width/2 - offset + lane_left.width)
    else
        FeatureValue(NaN, FeatureState.MISSING) # there is no left lane
    end
end
generate_feature_functions("MarkerDist_Right_Right", :d_mrr, Float64, "m", can_be_missing=true)
function Base.get(::Feature_MarkerDist_Right_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    #=
    Distance to the right lane marker one lane to the right
    =#

    veh = rec[pastframe][vehicle_index]
    lane = roadway[veh.state.posF.roadind.tag]
    if n_lanes_right(lane, roadway) > 0
        offset = veh.state.posF.t
        lane_right = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        FeatureValue(lane.width/2 + offset + lane_right.width)
    else
        FeatureValue(NaN, FeatureState.MISSING) # there is no right lane
    end
end
generate_feature_functions("RoadEdgeDist_Left", :d_edgel, Float64, "m")
function Base.get(::Feature_RoadEdgeDist_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    offset = veh.state.posF.t
    footpoint = get_footpoint(veh)
    seg = roadway[veh.state.posF.roadind.tag.segment]
    lane = seg.lanes[end]
    roadproj = proj(footpoint, lane, roadway)
    curvept = roadway[RoadIndex(roadproj)]
    lane = roadway[roadproj.tag]
    FeatureValue(lane.width/2 + abs(curvept.pos - footpoint) - offset)
end
generate_feature_functions("RoadEdgeDist_Right", :d_edger, Float64, "m")
function Base.get(::Feature_RoadEdgeDist_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    offset = veh.state.posF.t
    footpoint = get_footpoint(veh)
    seg = roadway[veh.state.posF.roadind.tag.segment]
    lane = seg.lanes[1]
    roadproj = proj(footpoint, lane, roadway)
    curvept = roadway[RoadIndex(roadproj)]
    lane = roadway[roadproj.tag]
    FeatureValue(lane.width/2 + abs(curvept.pos - footpoint) + offset)
end
generate_feature_functions("LaneOffsetLeft", :posFtL, Float64, "m", can_be_missing=true)
function Base.get(::Feature_LaneOffsetLeft, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh_ego = rec[pastframe][vehicle_index]
    t = veh_ego.state.posF.t
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_left(lane, roadway) > 0
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        lane_offset = t - lane.width/2 - lane_left.width/2
        FeatureValue(lane_offset)
    else
        FeatureValue(NaN, FeatureState.MISSING)
    end
end
generate_feature_functions("LaneOffsetRight", :posFtR, Float64, "m", can_be_missing=true)
function Base.get(::Feature_LaneOffsetRight, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh_ego = rec[pastframe][vehicle_index]
    t = veh_ego.state.posF.t
    lane = roadway[veh_ego.state.posF.roadind.tag]
    if n_lanes_right(lane, roadway) > 0
        lane_right = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        lane_offset = t + lane.width/2 + lane_right.width/2
        FeatureValue(lane_offset)
    else
        FeatureValue(NaN, FeatureState.MISSING)
    end
end
generate_feature_functions("N_Lane_Right", :n_lane_right, Int, "-", lowerbound=0.0)
function Base.get(::Feature_N_Lane_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    nlr = rec[pastframe][vehicle_index].state.posF.roadind.tag.lane - 1
    FeatureValue(convert(Float64, nlr))
end
generate_feature_functions("N_Lane_Left", :n_lane_left, Int, "-", lowerbound=0.0)
function Base.get(::Feature_N_Lane_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    seg = roadway[veh.state.posF.roadind.tag.segment]
    nll = length(seg.lanes) - veh.state.posF.roadind.tag.lane
    FeatureValue(convert(Float64, nll))
end
generate_feature_functions("Has_Lane_Right", :has_lane_right, Bool, "-", lowerbound=0.0, upperbound=1.0)
function Base.get(::Feature_Has_Lane_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    val = get(N_LANE_RIGHT, rec, roadway, vehicle_index, pastframe).v > 0.0
    FeatureValue(convert(Float64, val))
end
generate_feature_functions("Has_Lane_Left", :has_lane_left, Bool, "-", lowerbound=0.0, upperbound=1.0)
function Base.get(::Feature_Has_Lane_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    val = get(N_LANE_LEFT, rec, roadway, vehicle_index, pastframe).v > 0.0
    FeatureValue(convert(Float64, val))
end
generate_feature_functions("LaneCurvature", :curvature, Float64, "1/m", can_be_missing=true)
function Base.get(::Feature_LaneCurvature, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    veh = rec[pastframe][vehicle_index]
    curvept = roadway[veh.state.posF.roadind]
    val = curvept.k
    if isnan(val)
        FeatureValue(0.0, FeatureState.MISSING)
    else
        FeatureValue(val)
    end
end

# Dist_Merge
# Dist_Split

generate_feature_functions("TimeToCrossing_Right", :ttcr_mr, Float64, "s", lowerbound=0.0, censor_hi=10.0)
function Base.get(::Feature_TimeToCrossing_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    d_mr = get(MARKERDIST_RIGHT, rec, roadway, vehicle_index, pastframe).v
    velFt = get(VELFT, rec, roadway, vehicle_index, pastframe).v

    if d_mr > 0.0 && velFt < 0.0
        FeatureValue(-d_mr / velFt)
    else
        FeatureValue(Inf)
    end
end
generate_feature_functions("TimeToCrossing_Left", :ttcr_ml, Float64, "s", lowerbound=0.0, censor_hi=10.0)
function Base.get(::Feature_TimeToCrossing_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    d_ml = get(MARKERDIST_RIGHT, rec, roadway, vehicle_index, pastframe).v
    velFt = get(VELFT, rec, roadway, vehicle_index, pastframe).v

    if d_ml > 0.0 && velFt < 0.0
        FeatureValue(-d_ml / velFs)
    else
        FeatureValue(Inf)
    end
end
generate_feature_functions("EstimatedTimeToLaneCrossing", :est_ttcr, Float64, "s", lowerbound=0.0, censor_hi=10.0)
function Base.get(::Feature_EstimatedTimeToLaneCrossing, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    ttcr_left = get(TIMETOCROSSING_LEFT, rec, roadway, vehicle_index, pastframe).v
    ttcr_right = get(TIMETOCROSSING_RIGHT, rec, roadway, vehicle_index, pastframe).v
    FeatureValue(min(ttcr_left, ttcr_right))
end
generate_feature_functions("A_REQ_StayInLane", :a_req_stayinlane, Float64, "m/s²", can_be_missing=true)
function Base.get(::Feature_A_REQ_StayInLane, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    velFt = get(VELFT, rec, roadway, vehicle_index, pastframe).v

    if velFt > 0.0
        d_mr = get(MARKERDIST_RIGHT, rec, roadway, vehicle_index, pastframe).v
        if d_mr > 0.0
            return FeatureValue(0.5velFt*velFt / d_mr)
        else
            return FeatureValue(NaN, FeatureState.MISSING)
        end
    else
        d_ml = get(MARKERDIST_LEFT, rec, roadway, vehicle_index, pastframe)
        if d_ml < 0.0
            return FeatureValue(-0.5velFt*velFt / d_ml)
        else
            return FeatureValue(NaN, FeatureState.MISSING)
        end
    end
end

generate_feature_functions("Time_Consecutive_Brake", :time_consec_brake, Float64, "s", lowerbound=0.0)
function Base.get(::Feature_Time_Consecutive_Brake, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)

    prev_accel = convert(Float64, get(ACC, rec, roadway, vehicle_index, pastframe))
    if prev_accel ≥ 0.0
        FeatureValue(0.0)
    else
        pastframe_orig = pastframe
        id = rec[pastframe][vehicle_index].id
        while pastframe_inbounds(rec, pastframe-1) &&
              get(ACC, rec, roadway, findfirst(rec[pastframe-1], id)) < 0.0

            pastframe -= 1
        end

        FeatureValue(get_elapsed_time(rec, pastframe, pastframe_orig))
    end
end
generate_feature_functions("Time_Consecutive_Accel", :time_consec_accel, Float64, "s", lowerbound=0.0)
function Base.get(::Feature_Time_Consecutive_Accel, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)

    prev_accel = convert(Float64, get(ACC, rec, roadway, vehicle_index, pastframe))
    if prev_accel ≤ 0.0
        FeatureValue(0.0)
    else

        pastframe_orig = pastframe
        id = rec[pastframe][vehicle_index].id
        while pastframe_inbounds(rec, pastframe-1) &&
              get(ACC, rec, roadway, findfirst(rec[pastframe-1], id)) > 0.0

            pastframe -= 1
        end

        FeatureValue(get_elapsed_time(rec, pastframe, pastframe_orig))
    end
end
generate_feature_functions("Time_Consecutive_Throttle", :time_consec_throttle, Float64, "s", lowerbound=0.0)
function Base.get(::Feature_Time_Consecutive_Throttle, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    tc_accel = get(TIME_CONSECUTIVE_ACCEL, rec, roadway, vehicle_index, pastframe).v
    tc_brake = get(TIME_CONSECUTIVE_BRAKE, rec, roadway, vehicle_index, pastframe).v
    FeatureValue(tc_accel ≥ tc_brake ? tc_accel : -tc_brake)
end

#############################################
#
# FRONT
#
#############################################

generate_feature_functions("Dist_Front", :d_front, Float64, "m", lowerbound=0.0, can_be_missing=true)
function Base.get(::Feature_Dist_Front, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_lane(rec[pastframe], vehicle_index, roadway),
    censor_hi::Float64=100.0,
    )

    if neighborfore.ind == 0
        FeatureValue(100.0, FeatureState.CENSORED_HI)
    else
        scene = rec[pastframe]
        len_ego = scene[vehicle_index].def.length
        len_oth = scene[neighborfore.ind].def.length
        FeatureValue(neighborfore.Δs - len_ego/2 - len_oth/2)
    end
end
generate_feature_functions("Speed_Front", :v_front, Float64, "m/s", can_be_missing=true)
function Base.get(::Feature_Speed_Front, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_lane(rec[pastframe], vehicle_index, roadway),
    )

    if neighborfore.ind == 0
        FeatureValue(0.0, FeatureState.MISSING)
    else
        FeatureValue(rec[pastframe][neighborfore.ind].state.v)
    end
end
generate_feature_functions("Timegap", :timegap, Float64, "s", can_be_missing=true)
function Base.get(::Feature_Timegap, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_lane(rec[pastframe], vehicle_index, roadway),
    censor_hi::Float64 = 10.0,
    )

    v = rec[pastframe][vehicle_index].state.v

    if v ≤ 0.0 || neighborfore.ind == 0
        FeatureValue(censor_hi, FeatureState.CENSORED_HI)
    else
        scene = rec[pastframe]
        len_ego = scene[vehicle_index].def.length
        len_oth = scene[neighborfore.ind].def.length
        Δs = neighborfore.Δs - len_ego/2 - len_oth/2

        if Δs > 0.0
            FeatureValue(Δs / v)
        else
            FeatureValue(0.0) # collision!
        end
    end
end

generate_feature_functions("Inv_TTC", :inv_ttc, Float64, "1/s", can_be_missing=true)
function Base.get(::Feature_Inv_TTC, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_lane(rec[pastframe], vehicle_index, roadway),
    censor_hi::Float64 = 10.0,
    )

    if neighborfore.ind == 0
        FeatureValue(0.0, FeatureState.MISSING)
    else
        scene = rec[pastframe]
        veh_fore = scene[neighborfore.ind]
        veh_rear = scene[vehicle_index]

        len_ego = veh_fore.def.length
        len_oth = veh_rear.def.length
        Δs = neighborfore.Δs - len_ego/2 - len_oth/2

        Δv = veh_fore.state.v - veh_rear.state.v

        if Δs < 0.0 # collision!
            FeatureValue(censor_hi, FeatureState.CENSORED_HI)
        elseif Δv > 0.0 # front car is pulling away
            FeatureValue(0.0)
        else
            f = -Δv/Δs
            if f > censor_hi
                FeatureValue(f, FeatureState.CENSORED_HI)
            else
                FeatureValue(f)
            end
        end
    end
end
generate_feature_functions("TTC", :ttc, Float64, "s", can_be_missing=true)
function Base.get(::Feature_TTC, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_lane(rec[pastframe], vehicle_index, roadway),
    censor_hi::Float64 = 10.0,
    inv_ttc::FeatureValue = get(INV_TTC, rec, roadway, vehicle_index, pastframe, neighborfore=neighborfore, censor_hi=censor_hi),
    )

    if inv_ttc.i == FeatureState.MISSING
        # if the value is missing then front car not found and set TTC to censor_hi
        return FeatureValue(censor_hi, FeatureState.MISSING)
    else
        @assert is_feature_valid(inv_ttc) || inv_ttc.i == FeatureState.CENSORED_HI
        return FeatureValue(min(1.0 / inv_ttc.v, censor_hi))
    end
end

#############################################
#
# FRONT LEFT
#
#############################################
generate_feature_functions("Dist_Front_Left", :d_front_left, Float64, "m", lowerbound=0.0, can_be_missing=true)
function Base.get(::Feature_Dist_Front_Left, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_left_lane(rec[pastframe], vehicle_index, roadway),
    )

    get(DIST_FRONT, rec, roadway, vehicle_index, pastframe, neighborfore=neighborfore)
end

#############################################
#
# FRONT RIGHT
#
#############################################

generate_feature_functions("Dist_Front_Right", :d_front_right, Float64, "m", lowerbound=0.0, can_be_missing=true)
function Base.get(::Feature_Dist_Front_Right, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborfore::NeighborLongitudinalResult = get_neighbor_fore_along_right_lane(rec[pastframe], vehicle_index, roadway),
    )

    get(DIST_FRONT, rec, roadway, vehicle_index, pastframe, neighborfore=neighborfore)
end

#############################################
#
# REAR
#
#############################################

generate_feature_functions("Dist_Rear", :d_rear, Float64, "m", lowerbound=0.0, can_be_missing=true)
function Base.get(::Feature_Dist_Rear, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborrear::NeighborLongitudinalResult = get_neighbor_rear_along_lane(rec[pastframe], vehicle_index, roadway),
    censor_hi::Float64=100.0,
    )

    if neighborrear.ind == 0
        FeatureValue(100.0, FeatureState.CENSORED_HI)
    else
        scene = rec[pastframe]
        len_ego = scene[vehicle_index].def.length
        len_oth = scene[neighborrear.ind].def.length
        FeatureValue(neighborrear.Δs - len_ego/2 - len_oth/2)
    end
end
generate_feature_functions("Speed_Rear", :v_rear, Float64, "m/s", can_be_missing=true)
function Base.get(::Feature_Speed_Rear, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0;
    neighborrear::NeighborLongitudinalResult = get_neighbor_rear_along_lane(rec[pastframe], vehicle_index, roadway),
    )

    if neighborrear.ind == 0
        FeatureValue(0.0, FeatureState.MISSING)
    else
        FeatureValue(rec[pastframe][neighborrear.ind].state.v)
    end
end

#############################################
#
# SCENE WISE
#
#############################################

generate_feature_functions("Is_Colliding", :is_colliding, Bool, "-", lowerbound=0.0, upperbound=1.0)
function Base.get{S,D,I,R}(::Feature_Is_Colliding, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0;
    mem::CPAMemory=CPAMemory(),
    )

    scene = rec[pastframe]
    is_colliding = convert(Float64, get_first_collision(scene, vehicle_index, mem).is_colliding)
    FeatureValue(is_colliding)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_AddressSpace.jl")
include("model_DhcpOptions.jl")
include("model_IPAddressAvailabilityResult.jl")
include("model_Resource.jl")
include("model_ResourceNavigationLinkFormat.jl")
include("model_SubResource.jl")
include("model_SubnetListResult.jl")
include("model_SubnetPropertiesFormat.jl")
include("model_VirtualNetworkListResult.jl")
include("model_VirtualNetworkListUsageResult.jl")
include("model_VirtualNetworkPeeringListResult.jl")
include("model_VirtualNetworkPeeringPropertiesFormat.jl")
include("model_VirtualNetworkPropertiesFormat.jl")
include("model_VirtualNetworkUsage.jl")
include("model_VirtualNetworkUsageName.jl")
include("model_IPConfiguration.jl")
include("model_NetworkSecurityGroup.jl")
include("model_ResourceNavigationLink.jl")
include("model_RouteTable.jl")
include("model_Subnet.jl")
include("model_VirtualNetwork.jl")
include("model_VirtualNetworkPeering.jl")


text:
# === PRIVATE FUNCTIONS THAT ARE NOT EXPORTED =============================================================== #
# Objective functions -
function _obj_function_logistics_regression(parameters::Array{Float64,1}, labels::Array{Int64,1},
    dataMatrix::Array{Float64,2}, bias::Float64 = 0.0)::Float64

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)
    prob_array = Array{Float64,1}(undef,number_of_rows)
    term_array = Array{Float64,2}(undef,number_of_rows,3)

    # augment the dataMatrix -
    X = [ones(number_of_rows) dataMatrix]

    # ok, so let's compute the log liklehood -
    # start w/the prob -
    for row_index = 1:number_of_rows
        prob_array[row_index] = _logistics_classifier_logic(parameters, X[row_index,:], bias)
    end

    # compute the term array -
    for row_index = 1:number_of_rows

        # compute term 1 and 2 -
        term_1_value = labels[row_index]*log(prob_array[row_index])
        term_2_value = (1 - labels[row_index])*log((1 - prob_array[row_index]))

        if (isnan(term_1_value) == true)
            term_1_value = 0.0
        end

        if (isnan(term_2_value) == true)
            term_2_value = 0.0
        end

        # package -
        term_array[row_index,1] = term_1_value
        term_array[row_index,2] = term_2_value
        term_array[row_index,3] = term_1_value + term_2_value
    end

    # check for inf -
    tmp = sum(term_array[:,3])
    if (isinf(tmp) == true)
        tmp = 0.0
    end

    # compute log liklehood -
    LL = -1*tmp

    # return -
    return LL
end

function _obj_function_linear_regression(parameters::Array{Float64,1}, outputVector::Array{Float64,1},
    dataMatrix::Array{Float64,2}, lambda::Float64)::Float64

    # compute the Y_model -
    results_tuple = _evaluate_ols_linear_model(outputVector, dataMatrix, parameters)

    # the results_tuple contains the residual_value already (yes!)
    # we add a "regularization" term (L2 ridge regression) to avoid the overfitting issue -
    rms_error = results_tuple.residual + (lambda*sum(parameters.^2))

    # return -
    return rms_error
end

function _obj_function_linear_regression(parameters::Array{Float64,1}, outputVector::Array{Float64,1},
    dataMatrix::Array{Float64,1})::Float64

    # compute the Y_model -
    results_tuple = _evaluate_ols_linear_model(outputVector, dataMatrix, parameters)

    # the results_tuple contains the residual_value already (yes!)
    # we add a "regularization" term (L2 ridge regression) to avoid the overfitting issue -
    rms_error = results_tuple.residual + (lambda*sum(parameters.^2))

    # return -
    return rms_error
end

# evaluate linear model -
function _evaluate_ols_linear_model(dataMatrix::Array{Float64,1}, paramaterArray::Array{Float64,1})::Float64

    # initialize -
    number_of_rows = length(dataMatrix)

    # augement the data array 1
    X = [1.0 transpose(dataMatrix)]

    # compute the Y_model -
    Y_model = X*paramaterArray

    return Y_model
end

function _evaluate_ols_linear_model(outputArray::Array{Float64,1}, dataMatrix::Array{Float64,2},
    paramaterArray::Array{Float64,1})::NamedTuple

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)
    Y_measured = outputArray

    # augement the data array 0
    ones_array = ones(number_of_rows)
    X = [ones_array dataMatrix]

    # compute the Y_model -
    Y_model = X*paramaterArray

    # what is the correlation?
    correlation_value = cor(Y_measured, Y_model)

    # what is the residual?
    residual_value = rmsd(Y_measured, Y_model)

    # package -
    results_tuple = (model_prediction=Y_model, residual=residual_value, correlation=correlation_value)

    # return -
    return results_tuple
end

function _evaluate_ols_linear_model(outputArray::Array{Float64,1}, dataMatrix::Array{Float64,1},
    paramaterArray::Array{Float64,1})::NamedTuple

    # initialize -
    number_of_rows = length(outputArray)
    Y_measured = outputArray

    # augement the data array 0
    ones_array = ones(number_of_rows)
    X = [ones_array dataMatrix]

    # compute the Y_model -
    Y_model = X*paramaterArray

    # what is the correlation?
    correlation_value = cor(Y_measured, Y_model)

    # what is the residual?
    residual_value = rmsd(Y_measured, Y_model)

    # package -
    results_tuple = (model_prediction=Y_model, residual=residual_value, correlation=correlation_value)

    # return -
    return results_tuple
end

# default: leave one out logic -
function _leave_one_out_logic(index::Int64, outputVector::Array{Float64,1}, dataMatrix::Array{Float64,2})::NamedTuple

    # ok, so need to impl leave one out -
    (number_of_rows, number_of_cols) = size(dataMatrix)

    # generate the "full" range -
    idx_full_index_array = range(1,stop=number_of_rows,step=1) |> collect

    # generate index array with a row = index missing - this is the training index array
    idx_missing_index_array = setdiff(idx_full_index_array, index)

    # what index's did we leave out?
    prediction_index_array = Array{Int64,1}()
    push!(prediction_index_array, index)

    # collect -
    Yhat = outputVector[idx_missing_index_array]
    Xhat = dataMatrix[idx_missing_index_array,:]

    # package and return -
    results_tuple = (output_vector=Yhat, input_matrix=Xhat,
        training_index_array=idx_missing_index_array,
        prediction_index_array=prediction_index_array)

    # return -
    return results_tuple
end
# =========================================================================================================== #

# === PUBLIC FUNCTIONS THAT ARE EXPORTED ==================================================================== #
"""
    mle_logistic_model_classifier_cross_validation(labelVector::Array{Int64,1}, dataMatrix::Array{Float64,2},
        numberOfGroups::Int64; selectionFunction::Union{Nothing, Function} = nothing)::VLResult
"""
function mle_logistic_model_classifier_cross_validation(labelVector::Array{Int64,1}, dataMatrix::Array{Float64,2},
    numberOfGroups::Int64; selectionFunction::Union{Nothing, Function} = nothing)::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)

    try

        ex = ErrorException("Ooops! mle_logistic_model_classifier_cross_validation not yet implemented.");
        throw(ex)

    catch error
        return VLResult(error)
    end
end

"""
    mle_fit_logistic_model_classifier(labelVector::Array{Int64,1}, dataMatrix::Array{Float64,2};
        initialParameterArray::Union{Nothing,Array{Float64,1}} = nothing, maxIterations::Int64=10000,
        showTrace::Bool = false, bias::Float64=0.0)::VLResult
"""
function mle_fit_logistic_model_classifier(labelVector::Array{Int64,1}, dataMatrix::Array{Float64,2};
    initialParameterArray::Union{Nothing,Array{Float64,1}} = nothing, maxIterations::Int64=10000,
    showTrace::Bool = false, bias::Float64=0.0)::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)

    try

        # setup the objective function -
        OF(p) = _obj_function_logistics_regression(p,labelVector,dataMatrix, bias)

        # setup initial guess -
        pinitial = 0.1*ones(number_of_cols+1)
        if (isnothing(initialParameterArray) == false)
            pinitial = initialParameterArray
        end

        # call the optimizer -
        opt_result = optimize(OF, pinitial, NelderMead(),
            Optim.Options(iterations=maxIterations, show_trace=showTrace))

        # get the optimal parameters -
        β = Optim.minimizer(opt_result)

        # return -
        return VLResult(β)
    catch error
        return VLResult(error)
    end
end

"""
    ols_fit_linear_model(outputVector::Array{Float64,1}, dataMatrix::Array{Float64,1})::VLResult
"""
function ols_fit_linear_model(outputVector::Array{Float64,1}, dataMatrix::Array{Float64,1};
    lambda::Float64 = 0.0, initialParameterArray::Union{Nothing,Array{Float64,1}} = nothing, maxIterations::Int64=10000,
    showTrace::Bool = false)::VLResult

    # initialize -
    number_of_cols = 1 # we have a single data col

    try

        # setup the obj function -
        OF(p) = _obj_function_linear_regression(p,outputVector,dataMatrix, lambda)

        # setup initial guess -
        pinitial = 0.1*ones(number_of_cols+1)
        if (isnothing(initialParameterArray) == false)
            pinitial = initialParameterArray
        end

        # call the optimizer -
        opt_result = optimize(OF, pinitial, NelderMead(),
            Optim.Options(iterations=maxIterations, show_trace=showTrace))

        # get the optimal parameters -
        β = Optim.minimizer(opt_result)

        # compute some performance stuff -
        performance_tuple = _evaluate_ols_linear_model(outputVector,dataMatrix,β)

        # setup results tuple -
        results_tuple = (model_prediction=performance_tuple.model_prediction,
            correlation=performance_tuple.correlation, residual=performance_tuple.residual,
            parameters=β)

        # return -
        return VLResult(results_tuple)
    catch error
        return VLResult(error)
    end
end

"""
    ols_fit_linear_model(outputVector::Array{Float64,1}, dataMatrix::Array{Float64,2})::VLResult
"""
function ols_fit_linear_model(outputVector::Array{Float64,1}, dataMatrix::Array{Float64,2}; lambda::Float64 = 0.0,
    initialParameterArray::Union{Nothing,Array{Float64,1}} = nothing, maxIterations::Int64=10000,
    showTrace::Bool = false)::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)

    try

        # setup the obj function -
        OF(p) = _obj_function_linear_regression(p,outputVector,dataMatrix)

        # setup initial guess -
        pinitial = 0.1*ones(number_of_cols+1)
        if (isnothing(initialParameterArray) == false)
            pinitial = initialParameterArray
        end

        # call the optimizer -
        opt_result = optimize(OF, pinitial, NelderMead(),
            Optim.Options(iterations=maxIterations, show_trace=showTrace))

        # get the optimal parameters -
        β = Optim.minimizer(opt_result)

        # compute some performance stuff -
        performance_tuple = _evaluate_ols_linear_model(outputVector,dataMatrix,β)

        # setup results tuple -
        results_tuple = (model_prediction=performance_tuple.model_prediction,
            correlation=performance_tuple.correlation, residual=performance_tuple.residual,
            parameters=β)

        # return -
        return VLResult(results_tuple)
    catch error
        return VLResult(error)
    end
end

"""
    ols_fit_linear_model_cross_validation(outputVector::Array{Float64,1},
        dataMatrix::Array{Float64,2}; numberOfGroups::Int64 = 0, selectionFunction::Union{Nothing, Function} = nothing)::VLResult
"""
function ols_fit_linear_model_cross_validation(outputVector::Array{Float64,1},
    dataMatrix::Array{Float64,2}; numberOfGroups::Int64 = 0, numberOfElementsPerGroup::Int64 = 0, lambda::Float64 = 0.0,
    selectionFunction::Union{Nothing, Function} = nothing)::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataMatrix)

    try

        # check: if numberOfGroups = 0, then set to the number of rows -
        if (numberOfGroups == 0)
            numberOfGroups = number_of_rows
        end

        if (numberOfElementsPerGroup == 0)
            numberOfElementsPerGroup = 1
        end

        # check: if selectionFunction == nothing, then use _leave_one_out_logic -
        mySelectionFunction = _leave_one_out_logic
        if (isnothing(selectionFunction) == false)
            mySelectionFunction = selectionFunction
        end

        # initialize storage -
        parameter_storage_array = zeros((number_of_cols + 1), numberOfGroups)
        total_residual_array = Array{Float64,1}()
        total_correlation_array = Array{Float64,1}()
        model_training_array = zeros(number_of_rows, numberOfGroups)
        model_prediction_archive = Array{Float64,1}()
        measured_output_array = zeros(number_of_rows, numberOfGroups)
        selection_index_archive = Array{Array{Int64,1},1}()

        # ok, lets go ...
        for group_index = 1:numberOfGroups

            # pass the full data set, and the group index to the selection function
            # the selectionFunction gives back the output and input arrsy -
            selection_tuple = mySelectionFunction(group_index,outputVector,dataMatrix)
            Yhat = selection_tuple.output_vector
            Xhat = selection_tuple.input_matrix
            selection_index_array = selection_tuple.training_index_array
            prediction_index_array = selection_tuple.prediction_index_array

            # cache the selection index buffer -
            push!(selection_index_archive, selection_index_array)

            # YHat -
            # scale_result = z_score_transform_vector(Yhat)
            # if (isa(scale_result.value,Exception) == true)
            #     throw(scale_result.value)
            # end
            # Yhat_z_scaled = scale_result.value

            # # Input array -
            # scale_result = z_score_transform_array(Xhat)
            # if (isa(scale_result.value,Exception) == true)
            #     throw(scale_result.value)
            # end
            # Xhat_z_scaled = scale_result.value

            # fit the model -
            fit_model_result = ols_fit_linear_model(Yhat, Xhat; lambda)
            if (isa(fit_model_result.value,Exception) == true)
                throw(fit_model_result.value)
            end
            performance_tuple = fit_model_result.value
            theta_parameters = performance_tuple.parameters
            model_output = performance_tuple.model_prediction

            # compute the prediction - leave 1 out for now ...
            Y_prediction = _evaluate_ols_linear_model(dataMatrix[group_index,:], theta_parameters)

            # capture the model output -
            for (output_index, output_value) in enumerate(model_output)
                real_index = selection_index_array[output_index]
                model_training_array[real_index, group_index] = output_value
            end

            # capture the predicted output -
            push!(model_prediction_archive,Y_prediction)

            # capture the measured output -
            for (output_index, output_value) in enumerate(Yhat)
                real_index = selection_index_array[output_index]
                measured_output_array[real_index, group_index] = output_value
            end

            # put the parameters in storage -
            for (parameter_index,parameter_value) in enumerate(theta_parameters)
                parameter_storage_array[parameter_index, group_index] = parameter_value
            end

            # grab the residual and other stuff for this group -
            push!(total_residual_array, performance_tuple.residual)
            push!(total_correlation_array, performance_tuple.correlation)
        end

        # return the results in a NamedTuple -
        results_tuple = (correlation=total_correlation_array,
            residual=total_residual_array,parameters=parameter_storage_array,
            model_training_array=model_training_array,
            model_prediction_array=model_prediction_archive,
            measured_output_array=measured_output_array,
            selection_index_archive=selection_index_archive)

        # return -
        return VLResult(results_tuple)
    catch error
        return VLResult(error)
    end
end
# =========================================================================================================== #

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
module ForeignCallbacks

import Base: Libc

struct Node{T}
    next::Ptr{Node{T}} # needs to be first field
    data::T
end
Node(data::T) where T = Node{T}(C_NULL, data)

function calloc(::Type{T}) where T
    ptr = Libc.malloc(sizeof(T))
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), ptr, 0, sizeof(T))
    return convert(Ptr{T}, ptr)
end

mutable struct LockfreeQueue{T}
    @atomic head::Ptr{Node{T}}
    @atomic tail::Ptr{Node{T}}
    function LockfreeQueue{T}() where T
        @assert !Base.ismutabletype(T) && Base.isconcretetype(T) && Base.datatype_pointerfree(T)
        tmp = calloc(Node{T})
        new{T}(tmp, tmp)
    end
end

# Notes:
# We require at-least one node in the queue.
# When we dequeue, we discard the head node and return the data
# from the new head.

function enqueue!(q::LockfreeQueue{T}, data::T) where T
    node = Node(data)
    p_node = convert(Ptr{Node{T}}, Libc.malloc(sizeof(Node{T})))
    Base.unsafe_store!(p_node, node)

    # Update the tail node in queue
    p_tail = @atomicswap :acquire_release q.tail = p_node

    # Link former tail to new tail
    Core.Intrinsics.atomic_pointerset(convert(Ptr{Ptr{Node{T}}}, p_tail), p_node, :release)
    return nothing
end

# Manual implementation of `enqueue!` since `unsafe_pointer_to_objref(ptr)::T` creates a gcframe
function unsafe_enqueue!(ptr::Ptr{Cvoid}, data::T) where T
    node = Node(data)
    p_node = convert(Ptr{Node{T}}, Libc.malloc(sizeof(Node{T})))
    Base.unsafe_store!(p_node, node)

    # Update the tail node in queue
    ptr += fieldoffset(ForeignCallbacks.LockfreeQueue{T}, 2)
    p_tail = Core.Intrinsics.atomic_pointerswap(convert(Ptr{Ptr{Node{T}}}, ptr), p_node, :acquire_release)

    # Link former tail to new tail
    Core.Intrinsics.atomic_pointerset(convert(Ptr{Ptr{Node{T}}}, p_tail), p_node, :release)
    return nothing
end

function dequeue!(q::LockfreeQueue{T}) where T
    p_head = @atomic :acquire q.head

    success = false
    p_new_head = convert(Ptr{Node{T}}, C_NULL)
    while !success
        # Load new head
        p_new_head = Core.Intrinsics.atomic_pointerref(convert(Ptr{Ptr{Node{T}}}, p_head), :acquire)
        if p_new_head == convert(Ptr{Node{T}}, C_NULL)
            return nothing # never remove the last node, queue is empty
        end
        # Attempt replacement of current head with new head
        p_head, success = @atomicreplace :acquire_release :monotonic q.head p_head => p_new_head
    end

    # We have atomically advanced head and claimed a node.
    # We return the data from the new head
    # The lists starts of with a temporary node, which we will now free.
    head = unsafe_load(p_new_head) # p_head is now valid to free
    # TODO: Is there a potential race between `free(p_head)` and `unsafe_load(p_head)`
    #       on the previous `dequeue!`?
    #       As long as we only have one consumer this is fine.
    Libc.free(p_head)
    return Some(head.data)
end

mutable struct ForeignCallback{T}
    queue::LockfreeQueue{T}
    cond::Base.AsyncCondition

    function ForeignCallback{T}(callback) where T
        queue = LockfreeQueue{T}()

        cond = Base.AsyncCondition() do _
            data = dequeue!(queue)
            while data !== nothing
                Base.errormonitor(Threads.@spawn callback(something($data)))
                data = dequeue!(queue)
            end
            return
        end
        this = new{T}(queue, cond)
        finalizer(this) do this
            close(this.cond)
            # TODO: free queue we are leaking at least one node here
        end
    end
end

struct ForeignToken
    handle::Ptr{Cvoid}
    queue::Ptr{Cvoid}
end

"""
    Create a foreign token for a callback.

Note that the passed `ForeignCallback` must be GC preserved while this token is being used by
the foreign code. The user can pass the token around and use it to invoke the callback with the data
by calling `notify!`.
"""
ForeignToken(fc::ForeignCallback) = ForeignToken(fc.cond.handle, Base.pointer_from_objref(fc.queue))

function notify!(token::ForeignToken, data::T) where T
    unsafe_enqueue!(token.queue, data)
    ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), token.handle)
    return
end

end # module


text:
module Trapz
    export trapz, @trapz

    include("kernels.jl")

    function trapz(x, y)
        return integrate(x,y)
    end

    function trapz(x, y, val_i::Val{i}) where i
        return integrate(x,y,val_i)
    end

    """
# Syntax for integration over single axis
        trapz(x,y,axis=Val(End))
        Calculates ∫y[..., i (axis) ,...] dx[i]
        for performance reasons x must have the same float type as y.
        1-D Example:
```julia
    vx=0:0.01:1
    vy=(x->x^2).(vx)
    trapz(vx,vy)
```
            Result ≈ 1/3
# Multidimensional Syntax
        2-D Example:
```julia
    vx=0:0.01:1
    vy=0:0.01:2
    z=[x^2+y^2 for x=vx, y=vy]
    trapz((vx,vy),z) # equivalent to trapz(vx, trapz(vy, z))
```
            Result ≈ 4/3
## Partial integration Syntax
        2-D Example: partial integration
```julia
    vx=0:0.01:1
    vy=0:0.01:2
    z=[x^2+y^2 for x=vx, y=vy]
    trapz((vx,:),z) # equivalent to trapz(vx, z, Val(1))
```
            Result ≈ 201-element Array{Float64,1}:
                     0.33335
                     0.3334500000000001
                     0.33375000000000005
                     0.3342499999999999
                     ⋮
                     4.253749999999999
                     4.293449999999999
                     4.333350000000002

        2-D Example: partial integration
```julia
    vx=0:0.01:1
    vy=0:0.01:2
    z=[x^2+y^2 for x=vx, y=vy]
    trapz((:,vy),z) # equivalent to trapz(vy, z, Val(2))
```
            Result ≈ 101-element Array{Float64,1}:
                     2.666700000000002
                     2.6669
                     2.6674999999999995
                     2.6685000000000008
                     ⋮
                     4.587500000000002
                     4.626900000000002
                     4.6667000000000005

    """
    trapz
end # module


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

# -----------------------------------------------
# Basic Types
# -----------------------------------------------

"""
    Args(args...; kwargs...)

A structure for updating arguments to an actor's behavior.
"""
struct Args{A,B}
    args::A
    kwargs::B

    Args(args...; kwargs...) = new{typeof(args),typeof(kwargs)}(args, kwargs)
end

"""
    Bhv(func, a...; kw...)(c...)

A callable struct to represent actor behavior. It is executed
with parameters from the incoming communication.

# Parameters

- `f`: a callable object,
- `a...`: stored acquaintance parameters to `f`,
- `kw...`: stored keyword arguments,
- `c...`: parameters from the incoming communication.
"""
struct Bhv{F}
    f
    a::Tuple
    kw::Base.Iterators.Pairs
    ϕ::F

    function Bhv(f, a...; kw...)
        ϕ = (c...) -> f(a..., c...; kw...)
        new{typeof(ϕ)}(f, a, kw, ϕ)
    end
end
(p::Bhv)(c...) = p.ϕ(c...)

#
# Since Bhv contains an anonymous function, the following
# is needed to make it executable in another thread or worker.
# It returns a Bhv for the current world age.
#
_current(p::Bhv) = Bhv(p.f, p.a...; p.kw...)
_current(x) = x

"""
```
Link{C} <: ActorInterfaces.Classic.Addr
Link(chn::C, pid::Int, type::Symbol) where C
```

A mailbox for communicating with actors. A concrete type of
this must be returned by an actor on creation with [`spawn`](@ref).

# Fields/Parameters
- `chn::C`: C can be any type and characterizes the interface
    to an actor,
- `pid::Int`: the pid (worker process identifier) of the actor,
- `mode::Symbol`: a symbol characterizing the actor mode.
"""
mutable struct Link{C} <: Addr
    chn::C
    pid::Int
    mode::Symbol
end

"Abstract type for connections between actors."
abstract type Connection end

"""
```
_ACT
```
Internal actor status variable.

# Fields

1. `mode::Symbol`: the actor mode,
2. `bhv::Any`:  behavior - a callable object,
3. `init::Any`: initialization - a callable object,
4. `term::Any`: termination - a callable object,
5. `self::Link`: the actor's address,
6. `name::Union{Nothing,Symbol}`: the actor's registered name.
7. `res::Any`: the result of the last behavior execution,
8. `sta::Any`: a variable for representing state,
9. `usr::Any`: user variable for plugging in something,
10. `conn::Array{Connection,1}`: connected actors.

see also: [`Bhv`](@ref), [`Link`](@ref)
"""
mutable struct _ACT
    mode::Symbol
    bhv::Any
    init::Any
    term::Any
    self::Union{Nothing,Link}
    name::Union{Nothing,Symbol}
    res::Any
    sta::Any
    usr::Any
    conn::Array{Connection,1}
end

"""
    _ACT(mode=:default)

Return a actor variable `_ACT`.
"""
_ACT(mode=:default) = _ACT(mode, Bhv(+), fill(nothing, 7)..., Connection[])

"""
## Actor information
- `mode::Symbol`: actor mode,
- `bhvf::Any`: behavior function,
- `pid::Int`: process identifier,
- `thrd::Int`: thread,
- `task::Task`: actor task address,
- `tid::String`: proquint identifier based on task address,
- `name::Union{Nothing,Symbol}`: name under which the actor is
    registered, `nothing` if not registered.
"""
struct Info
    mode::Symbol
    bhvf::Any
    pid::Int
    thrd::Int
    task::UInt
    tid::String
    name::Union{Nothing,Symbol}
end

# -----------------------------------------------
# Public message types
# -----------------------------------------------
"Abstract type for messages to actors."
abstract type Msg end

"""
    Request(x, from::Link)

A generic [`Msg`](@ref) for user requests.
"""
struct Request <: Msg
    x
    from::Link
end

"""
    Response(y, from::Link=self())

A [`Msg`](@ref) representing a response to requests.

# Fields
- `y`: response content,
- `from::Link`: sender link.
"""
struct Response <: Msg
    y
    from::Link
end
Response(y) = Response(y, self())


text:
using PETSc, MPI, Printf
MPI.Initialized() || MPI.Init()
PETSc.initialize()

# Set up the RHS for a cell centered scheme for PDE
#
#  Δu = -cos(π * x) * cos(π * y)
#
# with zeros Neumann boundary conditions
function rhs!(
    ksp::PETSc.AbstractKSP{PetscScalar},
    b_vec::PETSc.AbstractVec{PetscScalar},
) where {PetscScalar}
    dm = PETSc.KSPGetDM(ksp)
    comm = PETSc.PetscObjectGetComm(ksp)
    corners = PETSc.DMDAGetCorners(dm)
    global_size = PETSc.DMDAGetInfo(dm).global_size[1:2]

    # Grid spacing in each direction
    h = PetscScalar(1) ./ global_size

    # Build the RHS vector for the cell-centered grid. Since the b_vec is not a
    # julia vector but a PETSc vector we need to convert it before operating on it
    PETSc.map_unsafe_localarray!(b_vec; read = false) do b
        b = reshape(b, Int64(corners.size[1]), Int64(corners.size[2]))
        for (iy, y) in enumerate(corners.lower[2]:corners.upper[2])
            y = (2y - 1) * h[2] / 2
            for (ix, x) in enumerate(corners.lower[1]:corners.upper[1])
                x = (2x - 1) * h[1] / 2
                b[ix, iy] = -cos(π * x) * cos(π * y) * h[1] * h[2]
            end
        end
    end

    # Assemble the PETSc vector
    PETSc.assemblybegin(b_vec)
    PETSc.assemblyend(b_vec)

    # Hack to remove constant from the vector since we are using all Neumann
    # boundary conditions
    nullspace = PETSc.MatNullSpace{PetscScalar}(comm, PETSc.PETSC_TRUE)
    PETSc.MatNullSpaceRemove!(nullspace, b_vec)
    PETSc.destroy(nullspace)

    # zero is the PETSc sucess code
    return 0
end

function jacobian!(
    ksp::PETSc.AbstractKSP{PetscScalar},
    J::PETSc.AbstractMat{PetscScalar},
    jac::PETSc.AbstractMat{PetscScalar},
) where {PetscScalar}
    dm = PETSc.KSPGetDM(ksp)
    corners = PETSc.DMDAGetCorners(dm)
    PetscInt = eltype(corners.size)

    global_size = PETSc.DMDAGetInfo(dm).global_size[1:2]

    # Grid spacing in each direction
    h = PetscScalar(1) ./ global_size

    # XXX: Would Mvectors be better?
    Sten = PETSc.MatStencil{PetscInt}
    col = Vector{Sten}(undef, 5)
    row = Vector{Sten}(undef, 1)
    val = Vector{PetscScalar}(undef, 5)

    for j in corners.lower[2]:corners.upper[2]
        for i in corners.lower[1]:corners.upper[1]
            row[1] = Sten(i = i, j = j)
            num = 1
            fill!(val, 0)
            if i > 1
                val[num] = -h[2] / h[1]
                col[num] = Sten(i = i - 1, j = j)
                num += 1
            end
            if i < global_size[1]
                val[num] = -h[2] / h[1]
                col[num] = Sten(i = i + 1, j = j)
                num += 1
            end
            if j > 1
                val[num] = -h[1] / h[2]
                col[num] = Sten(i = i, j = j - 1)
                num += 1
            end
            if j < global_size[2]
                val[num] = -h[1] / h[2]
                col[num] = Sten(i = i, j = j + 1)
                num += 1
            end
            val[num] = -sum(val)
            col[num] = Sten(i = i, j = j)
            PETSc.MatSetValuesStencil!(
                jac,
                row,
                col,
                val,
                PETSc.INSERT_VALUES;
                num_cols = num,
            )
        end
    end

    PETSc.assemblybegin(jac)
    PETSc.assemblyend(jac)

    # Since we have all Neumann there is constant nullspace
    comm = PETSc.PetscObjectGetComm(ksp)
    nullspace = PETSc.MatNullSpace{PetscScalar}(comm, PETSc.PETSC_TRUE)
    PETSc.MatSetNullSpace!(J, nullspace)
    PETSc.destroy(nullspace)
    return 0
end

function main(PetscScalar; comm = MPI.COMM_WORLD, options...)
    boundary_type = (PETSc.DM_BOUNDARY_NONE, PETSc.DM_BOUNDARY_NONE)
    stencil_type = PETSc.DMDA_STENCIL_STAR
    global_size = (11, 11)
    procs = (PETSc.PETSC_DECIDE, PETSc.PETSC_DECIDE)
    dof_per_node = 1
    stencil_width = 1
    points_per_proc = (nothing, nothing)

    opts =
        isinteractive() ? (ksp_monitor = nothing, ksp_view = nothing) :
        PETSc.parse_options(ARGS)

    da = PETSc.DMDACreate2d(
        PetscScalar,
        comm,
        boundary_type...,
        stencil_type,
        global_size...,
        procs...,
        dof_per_node,
        stencil_width,
        points_per_proc...;
        opts...,
    )

    PETSc.DMSetUp!(da)

    ksp = PETSc.KSP(da; opts...)

    PETSc.KSPSetComputeRHS!(ksp, rhs!)
    PETSc.KSPSetComputeOperators!(ksp, jacobian!)
    PETSc.setfromoptions!(ksp)
    PETSc.solve!(ksp)
end

main(Float64)
# main(Float32)


text:
"""
    location_is_out(u::AbstractArray{T,1},grid::gcmgrid)

Test whether location (x,y,fIndex) is out of domain. If true then
one typically needs to call update_location_cs! or update_location_dpdo!
"""

function location_is_out(u::AbstractArray{T,1},grid::gcmgrid) where T
    u[1]<0|| u[1]> grid.fSize[Int(u[end])][1]||
    u[2]<0|| u[2]> grid.fSize[Int(u[end])][2]
end

"""
    NeighborTileIndices_dpdo(ni::Int,nj::Int)

List of W, E, S, N neighbor tile IDs in the case of a doubly
periodic domain with ni x nj tiles.
"""
function NeighborTileIndices_dpdo(ni::Int,nj::Int)
    tmp=fill(0,ni*nj,4)
    for i=1:ni
        for j=1:nj
            k=i+ni*(j-1)
            kS=j-1; kS==0 ? kS=nj : nothing; kS=i+ni*(kS-1)
            kN=j+1; kN==nj+1 ? kN=1 : nothing; kN=i+ni*(kN-1)
            kW=i-1; kW==0 ? kW=ni : nothing; kW=kW+ni*(j-1)
            kE=i+1; kE==ni+1 ? kE=1 : nothing; kE=kE+ni*(j-1)
            tmp[k,1]=kW
            tmp[k,2]=kE
            tmp[k,3]=kS
            tmp[k,4]=kN
        end
    end
    return tmp
end

"""
    update_location_cs!

Update location (x,y,fIndex) when out of domain for cube sphere (cs) grid
as implemented by `MeshArrays.jl` (and MITgcm)
"""
function update_location_cs!(u::Array{Float64,1},𝑃::NamedTuple)
    x,y = u[1:2]
    fIndex = Int(u[end])
    nx,ny=𝑃.XC.fSize[fIndex]
    if x<0||x>nx||y<0||y>ny
        j = 0
        x<0 ? j=𝑃.aW[fIndex] : nothing
        x>nx ? j=𝑃.aE[fIndex] : nothing
        y<0 ? j=𝑃.aS[fIndex] : nothing
        y>ny ? j=𝑃.aN[fIndex] : nothing
        (x,y)=𝑃.RelocFunctions[j,fIndex](x,y)
        u[1]=x
        u[2]=y
        u[end]=j
    end
    #
    return u
end

"""
    update_location_llc!

Update location (x,y,fIndex) when out of domain for lat-lon-cap (llc) grid
as implemented by `MeshArrays.jl` (and MITgcm)
"""
function update_location_llc!(u::Array{Float64,1},𝑃::NamedTuple)
    x,y = u[1:2]
    fIndex = Int(u[end])
    nx,ny=𝑃.XC.fSize[fIndex]
    if y<0&&(fIndex==1||fIndex==2)
        u[2]=eps(y)
    elseif x>nx&&(fIndex==4||fIndex==5)
        u[1]=nx-eps(x)
    elseif x<0||x>nx||y<0||y>ny
        j = 0
        x<0 ? j=𝑃.aW[fIndex] : nothing
        x>nx ? j=𝑃.aE[fIndex] : nothing
        y<0 ? j=𝑃.aS[fIndex] : nothing
        y>ny ? j=𝑃.aN[fIndex] : nothing
        (x,y)=𝑃.RelocFunctions[j,fIndex](x,y)
        u[1]=x
        u[2]=y
        u[end]=j
    end
    #
    return u
end

"""
    update_location_dpdo!

Update location (x,y,fIndex) when out of domain. Note: initially, this
only works for the `dpdo` grid type provided by `MeshArrays.jl`.
"""
function update_location_dpdo!(u::AbstractArray{T,1},grid::gcmgrid) where T
    x,y = u[1:2]
    fIndex = Int(u[3])
    #
    nx,ny=grid.fSize[fIndex]
    ni,nj=Int.(transpose(grid.ioSize)./grid.fSize[1])
    WESN=NeighborTileIndices_dpdo(ni,nj)
    #
    if x<0
        x=x+nx
        u[1]=x
        fIndex=WESN[fIndex,1]
        u[3]=fIndex
    elseif x>=nx
        x=x-nx
        u[1]=x
        fIndex=WESN[fIndex,2]
        u[3]=fIndex
    end
    #
    if y<0
        y=y+ny
        u[2]=y
        fIndex=WESN[fIndex,3]
        u[3]=fIndex
    elseif y>=ny
        y=y-ny
        u[2]=y
        fIndex=WESN[fIndex,4]
        u[3]=fIndex
    end
    #
    return u
end

"""
    NeighborTileIndices_cs(grid::Dict)

Derive list of neighboring tile indices for a cs or llc grid + functions that
convert indices from one tile to another. Returns a Dict to merge later.
"""
function NeighborTileIndices_cs(grid::Dict)
    s = grid["XC"].fSize
    nFaces = length(s)
    nFaces == 5 ? s = vcat(s, s[3]) : nothing
    aW=Array{Int,1}(undef,nFaces)
    aE=similar(aW); aS=similar(aW); aN=similar(aW);
    for i = 1:nFaces
        (aW[i], aE[i], aS[i], aN[i], _, _, _, _) = MeshArrays.exch_cs_sources(i, s, 1)
    end
    RelocFunctions=RelocationFunctions_cs(grid["XC"])
    return Dict("aW" => aW, "aE" => aE, "aS" => aS, "aN" => aN, "RelocFunctions" => RelocFunctions)
end

"""
    RelocationFunctions_cs(xmpl)

Define matrix of functions to convert indices across neighboring tiles
"""
function RelocationFunctions_cs(xmpl::MeshArray)

# f1 : 0-n,0-n => -n-0,0-n     for 1->2, 3->4, 5->6
# f2 : 0-n,0-n => n-0,-n-0     for 2->4, 4->6, 6->2
# f3 : 0-n,0-n => 0-n,-n-0     for 2->3, 4->5, 6->1
# f4 : 0-n,0-n => -n-0,n-0     for 1->3, 3->5, 5->1
# g1, g2, g3, g4 : the reverse connections

    f1(x, y, nx, ny) = (x .- Float64(nx), y)
    f2(x, y, nx, ny) = (Float64(ny) .- y , x .- Float64(nx))
    f3(x, y, nx, ny) = (x, y .- Float64(ny))
    f4(x, y, nx, ny) = (y .- Float64(ny), Float64(nx) .- x )

    g1(x, y, nx, ny) = (x .+ Float64(nx), y)
    g2(x, y, nx, ny) = (y .+ Float64(ny), Float64(nx) .- x )
    g3(x, y, nx, ny) = (x, y .+ Float64(ny))
    g4(x, y, nx, ny) = (Float64(ny) .- y , x .+ Float64(nx))

#

    s = size.(xmpl.f)
    nFaces = length(s)
    tmp = Array{Function,2}(undef, 6, 6)

# f1, f2, f3, f4 : always get nx & ny from the source tile

    tmp[2, 1] = (x, y) -> f1(x, y, s[1][1], s[1][2])
    tmp[4, 3] = (x, y) -> f1(x, y, s[3][1], s[3][2])
    tmp[6, 5] = (x, y) -> f1(x, y, s[5][1], s[5][2])

    tmp[4, 2] = (x, y) -> f2(x, y, s[2][1], s[2][2])
    tmp[6, 4] = (x, y) -> f2(x, y, s[4][1], s[4][2])
    tmp[2, 6] = (x, y) -> f2(x, y, s[6][1], s[6][2])

    tmp[3, 2] = (x, y) -> f3(x, y, s[2][1], s[2][2])
    tmp[5, 4] = (x, y) -> f3(x, y, s[4][1], s[4][2])
    tmp[1, 6] = (x, y) -> f3(x, y, s[6][1], s[6][2])

    tmp[3, 1] = (x, y) -> f4(x, y, s[1][1], s[1][2])
    tmp[5, 3] = (x, y) -> f4(x, y, s[3][1], s[3][2])
    tmp[1, 5] = (x, y) -> f4(x, y, s[5][1], s[5][2])

# g1, g2, g3, g4 : nx or ny can come from source or target + notice nx/ny flips

    tmp[1, 2] = (x, y) -> g1(x, y, s[1][1], s[2][2])
    tmp[3, 4] = (x, y) -> g1(x, y, s[3][1], s[4][2])
    tmp[5, 6] = (x, y) -> g1(x, y, s[5][1], s[6][2])

    tmp[2, 4] = (x, y) -> g2(x, y, s[4][1], s[2][1])
    tmp[4, 6] = (x, y) -> g2(x, y, s[6][1], s[4][1])
    tmp[6, 2] = (x, y) -> g2(x, y, s[2][1], s[6][1])

    tmp[2, 3] = (x, y) -> g3(x, y, s[3][1], s[2][2])
    tmp[4, 5] = (x, y) -> g3(x, y, s[5][1], s[4][2])
    tmp[6, 1] = (x, y) -> g3(x, y, s[1][1], s[6][2])

    tmp[1, 3] = (x, y) -> g4(x, y, s[1][2], s[3][2])
    tmp[3, 5] = (x, y) -> g4(x, y, s[3][2], s[5][2])
    tmp[5, 1] = (x, y) -> g4(x, y, s[5][2], s[1][2])

    return tmp

end


text:
module Re2oAPIClients

using HTTP
using Base.Filesystem
using Dates
using Serialization
using JSON

TOKENFILEMANE = joinpath(homedir(), ".re2otoken")

loadtoken(tokenfilename) = deserialize(tokenfilename)
savetoken(token, tokenfilename=TOKENFILEMANE) = serialize(tokenfilename, token)

mutable struct Re2oAPIClient
    hostname::String
    username::String
    password::String
    tokenfilename::String
    use_tls::Bool
    token::Union{Pair{String,DateTime}, Nothing}
    Re2oAPIClient(hostname, username, password, token_file, use_tls, token) = new(hostname, username, password, token_file, use_tls, token)
end

Re2oAPIClient(hostname, username, password; token_file=TOKENFILEMANE, use_tls=true, token=nothing) =
    Re2oAPIClient(hostname, username, password, token_file, use_tls, token)

get_url_for(client, endpoint) = begin
    proto = if client.use_tls "https" else "http" end
    "$(proto)://$(client.hostname)/api/$(endpoint)"
end

loadservertoken(client::Re2oAPIClient) = begin
    params = Dict(
        "username" => client.username,
        "password" => client.password
    )
    r = HTTP.post(
        get_url_for(client, "token-auth"),
        ["Content-Type" => "application/json"],
        JSON.json(params)
    )
    response_body = JSON.parse(String(r.body))
    token = response_body["token"]
    expiration = response_body["expiration"]
    # This is uggly, don't know how to do it otherwise...
    token => DateTime(expiration[1:23])
end
loadservertoken!(client::Re2oAPIClient) = begin
    client.token = loadservertoken(client)
    savetoken(client.token, client.tokenfilename)
end

token!(client) = begin
    rightnow = DateTime(now())
    token = if !isnothing(client.token)
        client.token
    elseif ispath(client.tokenfilename)
        client.token = loadtoken(client.tokenfilename)
    else
        loadservertoken!(client)
    end
    if last(token) > rightnow
        loadservertoken!(client)
    end
    first(token)
end

request(client, method::String, url, params=nothing, headers=[]; kwargs...) = begin
    token = token!(client)
    push!(headers, "Authorization"=>"Token $token")

    r = if isnothing(params)
        HTTP.request(
            method,
            url,
            headers;
            kwargs...
        )
    else
        push!(headers, "Content-Type" => "application/json")
        HTTP.request(
            method,
            url,
            headers,
            JSON.json(params);
            kwargs...
        )
    end
    JSON.parse(String(r.body))
end

get(client, url, params=nothing, headers=[]; kwargs...) = request(client, "GET", url, params, headers; kwargs...)
post(client, url, params=nothing, headers=[]; kwargs...) = request(client, "POST", url, params, headers; kwargs...)
delete(client, url, params=nothing, headers=[]; kwargs...) = request(client, "DELETE", url, params, headers; kwargs...)
head(client, url, params=nothing, headers=[]; kwargs...) = request(client, "HEAD", url, params, headers; kwargs...)
option(client, url, params=nothing, headers=[]; kwargs...) = request(client, "OPTION", url, params, headers; kwargs...)
patch(client, url, params=nothing, headers=[]; kwargs...) = request(client, "PATCH", url, params, headers; kwargs...)
put(client, url, params=nothing, headers=[]; kwargs...) = request(client, "PUT", url, params, headers; kwargs...)

list(client, endpoint, max_results=nothing, params=Dict()) = begin
    url = get_url_for(client, endpoint)
    if !("page_size" in keys(params))
        push!(params, "page_size"=>if isnothing(max_results) "all" else max_results end)
    end
    response = get(client, url, params)
    results = response["results"]
    while !isnothing(response["next"]) && (isnothing(max_results) || length(results) < max_results)
        response = get(client, response["next"])
        append!(results, response["results"])
    end
    if !isnothing(max_results) results[1:max_results] else results end
end

export Re2oAPIClient, list

end # module


text:
module FeedbackChains

using Flux
import Flux: children, mapchildren

using ..Splitters
using ..Mergers
export FeedbackChain

struct FeedbackChain{T<:Tuple}
    layers::T
    FeedbackChain(xs...) = new{typeof(xs)}(xs)
end

# empty feedback chain just returns the input
function (c::FeedbackChain{Tuple{}})(h, x)
    return h, x
end # function (c::FeedbackChain{Tuple{}})

"""
    (c::FeedbackChain)(h, x)

Apply a `FeedbackChain` to input `x` with hidden state `h`. `h` should take the
form of a dictionary mapping `Splitter` names to states.
"""
function (c::FeedbackChain)(h, x)
    newh = Dict{String, Any}()
    for layer ∈ c.layers
        if layer isa Splitter
            newh[layer.name] = x
        elseif layer isa Merger
            x = layer(x, h[layer.splitname])
        else
            x = layer(x)
        end
    end
    return newh, x
end # function (c::FeedbackChain)

children(c::FeedbackChain) = c.layers
mapchildren(f, c::FeedbackChain) = FeedbackChain(f.(c.layers)...)

end # module FeedbackChains


text:
"""
    filter_w(w::FT, w_min::FT) where {FT}

Return velocity such that
`abs(filter_w(w, w_min)) >= abs(w_min)`
while preserving the sign of `w`.
"""
filter_w(w::FT, w_min::FT) where {FT} =
    max(abs(w), abs(w_min)) * (w < 0 ? sign(w) : 1)

"""
    enforce_unit_bounds(a_up_i::FT, a_min::FT, a_max::FT) where {FT}

Enforce variable to be positive.

Ideally, this safety net will be removed once we
have robust positivity preserving methods. For now,
we need this to avoid domain error in certain
circumstances.
"""
enforce_unit_bounds(a_up_i::FT, a_min::FT, a_max::FT) where {FT} =
    clamp(a_up_i, a_min, a_max)

"""
    enforce_positivity(x::FT) where {FT}

Enforce variable to be positive.

Ideally, this safety net will be removed once we
have robust positivity preserving methods. For now,
we need this to avoid domain error in certain
circumstances.
"""
enforce_positivity(x::FT) where {FT} = max(x, FT(0))

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct Asset
    package::Union{Nothing,String}
    path::String
end

function Asset(s::String)
    parts = strip.(split(s, ":"))
    if length(parts) == 1
        Asset(nothing, parts[1])
    elseif length(parts) == 2
        Asset(parts[1], parts[2])
    else
        error("invalid syntax for asset string: $s")
    end
end

macro asset_str(s::String)
    return Asset(s)
end

function Base.show(io::IO, x::Asset)
    print(io, "asset\"")
    if x.package !== nothing
        printstyled(io, x.package, ": "; color=:green)
    end
    printstyled(io, x.path, "\""; color=:cyan)
end

Base.convert(::Type{Asset}, s::String) = Asset(s)

"""
    Install

Installation configurations.

## Keywords

- `path`: installation path.
- `completion`: set to `true` to install shell auto-completion scripts.
- `quiet`: print logs or not, default is `false`.
- `compile`: julia compiler option for CLIs if not built as standalone application, default is "min".
- `optimize`: julia compiler option for CLIs if not built as standalone application, default is `2`.
"""
@option struct Install
    path::String = "~/.julia"
    completion::Bool = true
    quiet::Bool = false
    compile::String = "yes"
    optimize::Int = 2
end

"""
    Precompile

Precompilation files for `PackageCompiler`.

## Keywords

- `execution_file`: precompile execution file.
- `statements_file`: precompile statements file.
"""
@option struct Precompile
    execution_file::Vector{String} = String[]
    statements_file::Vector{String} = String[]
end

# See https://github.com/JuliaCI/julia-buildbot/blob/489ad6dee5f1e8f2ad341397dc15bb4fce436b26/master/inventory.py
function default_app_cpu_target()
    if Sys.ARCH === :i686
        return "pentium4;sandybridge,-xsaveopt,clone_all"
    elseif Sys.ARCH === :x86_64
        return "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)"
    elseif Sys.ARCH === :arm
        return "armv7-a;armv7-a,neon;armv7-a,neon,vfp4"
    elseif Sys.ARCH === :aarch64
        return "generic" # is this really the best here?
    elseif Sys.ARCH === :powerpc64le
        return "pwr8"
    else
        return "generic"
    end
end

"""
    SysImg

System image build configurations.

## Keywords

- `path`: system image path to generate into, default is "deps/lib".
- `incremental`: set to `true` to build incrementally, default is `true`.
- `filter_stdlibs`: set to `true` to filter out unused stdlibs, default is `false`.
- `cpu_target`: cpu target to build, default is `PackageCompiler.default_app_cpu_target()`.
- `precompile`: precompile configurations, see [`Precompile`](@ref), default is `Precompile()`.
"""
@option struct SysImg
    path::String = "deps"
    incremental::Bool = true
    filter_stdlibs::Bool = false
    cpu_target::String = default_app_cpu_target()
    precompile::Precompile = Precompile()
end

"""
    Download

Download information.

## Keywords

- `host`: where are the tarballs hosted, default is "github.com"
- `user`: required, user name on the host.
- `repo`: required, repo name on the host.

!!! note
    Currently this only supports github, and this is considered experimental.
"""
@option struct Download
    host::String = "github.com"
    user::String
    repo::String
end

"""
    Application

Application build configurations.

## Keywords

- `path`: application build path, default is "build".
- `incremental`: set to `true` to build incrementally, default is `true`.
- `filter_stdlibs`: set to `true` to filter out unused stdlibs, default is `false`.
- `cpu_target`: cpu target to build, default is `PackageCompiler.default_app_cpu_target()`.
- `precompile`: precompile configurations, see [`Precompile`](@ref), default is `Precompile()`.
- `c_driver_program`: driver program.
"""
@option struct Application
    path::String = "build"
    assets::Vector{Asset} = Asset[]
    incremental::Bool = false
    filter_stdlibs::Bool = true
    cpu_target::String = default_app_cpu_target()
    precompile::Precompile = Precompile()
    c_driver_program::Union{String,Nothing} = nothing
end

"""
    Comonicon <: AbstractConfiguration

Build configurations for Comonicon. One can set this option
via `Comonicon.toml` under the root path of a Julia
project directory and read in using [`read_configs`](@ref).

## Keywords

- `name`: required, the name of CLI file to install.
- `install`: installation options, see also [`Install`](@ref).
- `sysimg`: system image build options, see also [`SysImg`](@ref).
- `download`: download options, see also [`Download`](@ref).
- `application`: application build options, see also [`Application`](@ref).
"""
@option struct Comonicon
    name::String

    install::Install = Install()
    sysimg::Maybe{SysImg} = nothing
    download::Union{Download,Nothing} = nothing
    application::Union{Application,Nothing} = nothing
end


text:
mutable struct robot_state_t <: LCMType
    timestamp::Int64

    num_robots::Int32
    robot_name::Vector{String}

    num_joints::Int32
    joint_robot::Vector{Int32}
    joint_name::Vector{String}

    joint_position::Vector{Float32}
    joint_velocity::Vector{Float32}
end

@lcmtypesetup(robot_state_t,
    robot_name => (num_robots,),
    joint_robot => (num_joints,),
    joint_name => (num_joints,),
    joint_position => (num_joints,),
    joint_velocity => (num_joints,),
)


text:
# 2D SLAM with velocity states

"""
$(TYPEDEF)
"""
mutable struct DynPoint2VelocityPrior{T} <: IncrementalInference.FunctorSingleton where {T <: Distribution}
  z::T
  DynPoint2VelocityPrior{T}() where {T <: Distribution} = new{T}()
  DynPoint2VelocityPrior{T}(z1::T) where {T <: Distribution} = new{T}(z1)
end
DynPoint2VelocityPrior(z1::T) where {T <: Distribution} = DynPoint2VelocityPrior{T}(z1)
getSample(dp2v::DynPoint2VelocityPrior, N::Int=1) = (rand(dp2v.z,N), )

"""
$(TYPEDEF)
"""
mutable struct DynPoint2DynPoint2{T} <: IncrementalInference.FunctorPairwise where {T <: Distribution}
  z::T
  DynPoint2DynPoint2{T}() where {T <: Distribution} = new{T}()
  DynPoint2DynPoint2(z1::T) where {T <: Distribution} = new{T}(z1)
end
getSample(dp2dp2::DynPoint2DynPoint2, N::Int=1) = (rand(dp2dp2.z,N), )
function (dp2dp2::DynPoint2DynPoint2)(
            res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple,
            Xi::Array{Float64,2},
            Xj::Array{Float64,2}  )
  #
  z = meas[1][:,idx]
  xi, xj = Xi[:,idx], Xj[:,idx]
  dt = (userdata.variableuserdata[2].ut - userdata.variableuserdata[1].ut)*1e-6   # roughly the intended use of userdata
  res[1:2] = z[1:2] - (xj[1:2] - (xi[1:2]+dt*xi[3:4]))
  res[3:4] = z[3:4] - (xj[3:4] - xi[3:4])
  nothing
end

"""
$(TYPEDEF)
"""
mutable struct VelPoint2VelPoint2{T} <: IncrementalInference.FunctorPairwiseMinimize where {T <: Distribution}
  z::T
  VelPoint2VelPoint2{T}() where {T <: Distribution} = new{T}()
  VelPoint2VelPoint2{T}(z1::T) where {T <: Distribution} = new{T}(z1)
end
VelPoint2VelPoint2(z1::T) where {T <: Distribution} = VelPoint2VelPoint2{T}(z1)
getSample(vp2vp2::VelPoint2VelPoint2, N::Int=1) = (rand(vp2vp2.z,N), )
function (vp2vp2::VelPoint2VelPoint2{D})(
                res::Array{Float64},
                userdata,
                idx::Int,
                meas::Tuple,
                Xi::Array{Float64,2},
                Xj::Array{Float64,2}  ) where D
  #
  z = meas[1][:,idx]
  xi, xj = Xi[:,idx], Xj[:,idx]
  # change in time from microseconds with DynPoint2(ut=1_000_000) to seconds
  dt = (userdata.variableuserdata[2].ut - userdata.variableuserdata[1].ut)*1e-6   # roughly the intended use of userdata
  # change in psoition Xi \ Xj
  dp = (xj[1:2]-xi[1:2])
  # change in velocity Xi \ Xj
  dv = (xj[3:4]-xi[3:4])
  res[1] = 0.0
  res[1] += sum((z[1:2] - dp).^2) # (meas - predicted) change in position error term
  res[1] += sum((z[3:4] - dv).^2) # (meas - predicted) change in velocity error term

  ## now cross couple the change in position information, via timestamps to accompanying velocity
   # recompute integration of velocity influence
  # forward diff, "measured velocity"
  dp_dt = dp./dt
  # zeroth order integration
  res[1] += sum((dp_dt - xi[3:4]).^2) # (meas - predicted) velocity error term

    # first order integration
    # res[1] += sum((dp/dt - 0.5*(xj[3:4]+xi[3:4])).^2)

  # return objective cost
  return res[1]
end

"""
$(TYPEDEF)
"""
mutable struct Point2Point2Velocity{T} <: IncrementalInference.FunctorPairwiseMinimize where {T <: Distribution}
  z::T
  Point2Point2Velocity{T}() where {T <: Distribution} = new{T}()
  Point2Point2Velocity(z1::T) where {T <: Distribution} = new{T}(z1)
end
getSample(p2p2v::Point2Point2Velocity, N::Int=1) = (rand(p2p2v.z,N), )
function (p2p2v::Point2Point2Velocity)(
                res::Array{Float64},
                userdata,
                idx::Int,
                meas::Tuple,
                Xi::Array{Float64,2},
                Xj::Array{Float64,2}  )
  #
  z = meas[1][:,idx]
  xi, xj = Xi[:,idx], Xj[:,idx]
  dt = (userdata.variableuserdata[2].ut - userdata.variableuserdata[1].ut)*1e-6   # roughly the intended use of userdata
  dp = (xj[1:2]-xi[1:2])
  dv = (xj[3:4]-xi[3:4])
  res[1] = 0.0
  res[1] += sum((z[1:2] - dp).^2)
  res[1] += sum((dp/dt - 0.5*(xj[3:4]+xi[3:4])).^2)  # (dp/dt - 0.5*(xj[3:4]+xi[3:4])) # midpoint integration
  res[1]
end

## Packing Types================================================================

"""
$(TYPEDEF)
"""
mutable struct PackedDynPoint2VelocityPrior <: IncrementalInference.PackedInferenceType
  str::String
  PackedDynPoint2VelocityPrior() = new()
  PackedDynPoint2VelocityPrior(z1::String) = new(z1)
end

function convert(::Type{PackedDynPoint2VelocityPrior}, d::DynPoint2VelocityPrior)
  return PackedDynPoint2VelocityPrior(string(d.z))
end
function convert(::Type{DynPoint2VelocityPrior}, d::PackedDynPoint2VelocityPrior)
  distr = extractdistribution(d.str)
  return DynPoint2VelocityPrior(distr)
end

"""
$(TYPEDEF)
"""
mutable struct PackedVelPoint2VelPoint2 <: IncrementalInference.PackedInferenceType
  str::String
  PackedVelPoint2VelPoint2() = new()
  PackedVelPoint2VelPoint2(z1::String) = new(z1)
end

function convert(::Type{PackedVelPoint2VelPoint2}, d::VelPoint2VelPoint2)
  return PackedVelPoint2VelPoint2(string(d.z))
end
function convert(::Type{VelPoint2VelPoint2}, d::PackedVelPoint2VelPoint2)
  distr = extractdistribution(d.str)
  return VelPoint2VelPoint2(distr)
end

#