text:
function encodePolyline(coord::Array{Float64}, precision::Int64=5)
    #= Encodes an array of GPS coordinates to a polyline.

    Args:
        coord(Array{Float64, 1}(undef, 1)): GPS coordinates.
        precision(Int16): Exponent for rounding the GPS coordinates.
    Returns:
        String: Polyline encoded GPS coordinates.
    =#

    # Compute the rounding precision
    factor::Float64 = 10. ^precision
    coord = coord .* factor

    output = Array{Char, 1}(undef, 1)

    for c in range(1, stop=size(coord)[1])
        if c == 1
            writePolyline!(output, coordinate{Float64}(coord[c, 1], coord[c, 2]),
                           coordinate{Float64}(0., 0.))
        else
            writePolyline!(output, coordinate{Float64}(coord[c, 1], coord[c, 2]),
                           coordinate{Float64}(coord[c-1, 1], coord[c-1, 2]))
        end
    end

    return join(output[2:end])
end

function writePolyline!(output::Array{Char, 1}, currValue::coordinate{Float64},
                        prevValue::coordinate{Float64})
    #= Convert the given coordinate points in a polyline and mutate the output.

    Args:
        output(Array{Char, 1}): Holds the resultant polyline.
        currValue(coordinate{Float64}): Current GPS data point.
        prevValue(coordinate{Float64}): Previous GPS data point.

    Returns:
        output(Array{Char, 1}): Mutate output by adding the current addition to the
                                polyline.
    =#

    # Transform GPS coordinates to Integers and round
    roundCurrValue::coordinate{Int64} = roundCoordinate(currValue)
    roundPrevValue::coordinate{Int64} = roundCoordinate(prevValue)

    # Get the difference from the previous GPS coordinated
    diffCurrValue::coordinate{Int64} = diffCoordinate(roundCurrValue, roundPrevValue)

    # Left shift the data points
    leftShift::coordinate{Int64} = leftShiftCoordinate(diffCurrValue)

    # Transform into ASCII
    charCoordinate::coordinate{Array{Char, 1}} = convertToChar(leftShift)

    # Add the characters to the polyline
    append!(output, collect(charCoordinate.Lat))
    append!(output, collect(charCoordinate.Lng))
end

function encodeToChar(c::Int64)::Array{Char, 1}
    #= Perform the encoding of the character from a binary number to ASCII.

    Args:
        c(Int64): GPS coordinate.

    Returns:
        String: ASCII characters of the polyline.
    =#

    LatChars = Array{Char, 1}(undef, 1)

    # Invert a negative coordinate using two's complement
    if c < 0
        c = ~c
    end

    # Add a continuation bit at the LHS for non-last chunks using OR 0x20
    # (0x20 = 100000)
    while c >= 0x20
        # Get the last 5 bits (0x1f)
        # Add 63 (in order to get "better" looking polyline characters in ASCII)
        CharMod = (0x20 | (c & 0x1f)) + 63
        append!(LatChars, Char(CharMod))

        # Shift 5 bits
        c = c >> 5
    end

    # Modify the last chunk
    append!(LatChars, Char(c + 63))

    # The return string holds a beginning character at the start
    # skip it and return the rest
    return LatChars[2:end]
end

function roundCoordinate(currValue::coordinate{Float64})::coordinate{Int64}
    #= Convert the coordinate to an integer and round.

    Args:
        currValue (coordinate{Float64}): GPS data point as a real number.
    Returns:
        roundedValue (coordinate{Int64}): GPS data point as rounded integer.
    =#

    roundedValue::coordinate = coordinate{Int64}(copysign(floor(abs(currValue.Lat)), currValue.Lat),
                                                 copysign(floor(abs(currValue.Lng)), currValue.Lng))
    return roundedValue
end

function diffCoordinate(currValue::coordinate{Int64},
                        prevValue::coordinate{Int64})::coordinate{Int64}
    #= Polyline encoding only considers the difference between GPS data points
    in order to reduce memory. diffCoordinate obtains the difference between
    consecutive coordinate points.

    Args:
        currValue (coordinate{Int64}): The current GPS data point.
        prevValue (coordinate{Int64}): The previous GPS data point. The count
                                       starts from 0.
    Returns:
        coordinate{Int64}: The difference between the GPS data points.
    =#

    return coordinate{Int64}(currValue.Lat - prevValue.Lat,
                             currValue.Lng - prevValue.Lng)
end

function leftShiftCoordinate(currValue::coordinate{Int64})::coordinate{Int64}
    #= Left bitwise shift to leave space for a sign bit as the right most bit.

    Args:
        currValue(coordinate{Int64}): The difference between the last two
                                      consecutive GPS data points.
    Returns:
        coordinate{Int64}: Left bitwise shifted values.
    =#

    return coordinate{Int64}(currValue.Lat << 1,
                             currValue.Lng << 1)
end

function convertToChar(currValue::coordinate{Int64})::coordinate{Array{Char, 1}}
    #= Convert the coordinates into ascii symbols.

    Args:
        currValue(coordinate{Int64}): Integer GPS coordinates.

    Returns:
        coordinate{String}: GPS coordinates as ASCII characters.
    =#

    Lat::Int64 = currValue.Lat
    Lng::Int64 = currValue.Lng

    return coordinate{Array{Char, 1}}(encodeToChar(Lat), encodeToChar(Lng))
end


text:
import SoundPressureLevel

using Test
using Dates
using Libaudio
@everywhere using DeviceUnderTest

function simple_test()

    r = DeviceUnderTest.register()
    ms = zeros(1,8)
    ms[1,1] = 1
    mm = zeros(8,1)
    mm[1,1] = 1
    s,fs = Libaudio.wavread("test/acqua_ieee_male_250ms_10450ms.wav", Float64)
    SoundPressureLevel.setdba(
        r[0x9fefe994b7e95bf1],
        s[:,1],
        1,
        -6,
        ms,
        mm,
        48000,
        57,
        "D:/Depot\\Git\\SoundPressureLevel.jl\\foo\\",
        SoundPressureLevel.Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        SoundPressureLevel.Instrument("42AB",114,NaN,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        false)

    SoundPressureLevel.setdba(
        r[0x9fefe994b7e95bf1],
        s[:,2],
        1,
        -6,
        ms,
        mm,
        48000,
        57,
        "D:/Depot\\Git\\SoundPressureLevel.jl\\foo\\",
        SoundPressureLevel.Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        SoundPressureLevel.Instrument("42AB",114,NaN,Date("2018-07-24"),"26XX","12AA",0,"UFS"))

    SoundPressureLevel.setdba(
        r[0x9fefe994b7e95bf1],
        s[:,2:2],
        -6,
        ms,
        mm,
        48000,
        60,
        "D:/Depot\\Git\\SoundPressureLevel.jl\\foo\\",
        SoundPressureLevel.Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        SoundPressureLevel.Instrument("42AB",114,NaN,Date("2018-07-24"),"26XX","12AA",0,"UFS"))

    SoundPressureLevel.setdba(
        r[0x9fefe994b7e95bf1],
        s[:,1:1],
        -6,
        ms,
        mm,
        48000,
        60,
        "D:/Depot\\Git\\SoundPressureLevel.jl\\foo\\",
        SoundPressureLevel.Instrument("42AA",114,105.4,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        SoundPressureLevel.Instrument("42AB",114,NaN,Date("2018-07-24"),"26XX","12AA",0,"UFS"),
        false,
        -10,
        3.0,
        2.0,
        1,
        180,
        0.0,
        47999.6)
end

ext: jl


text:
using PlotlyJS

function random_line()
    n = 400
    rw() = cumsum(randn(n))
    trace1 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#1f77b4", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#1f77b4", width=1))
    trace2 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#9467bd", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="rgb(44, 160, 44)", width=1))
    trace3 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#bcbd22", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#bcbd22", width=1))
    layout = Layout(autosize=false, width=500, height=500,
                    margin=attr(l=0, r=0, b=0, t=65))
    PlotlyJS.plot([trace1, trace2, trace3], layout)
end
random_line()


text:
function merge(A,B)
    C = promote_type(eltype(A), eltype(B))[]
    i, j = 1, 1
    n, m = 1, 1
    while i <= length(A)  || j <= length(B)
        while i <= length(A) && ( A[i] <= B[m] || j > length(B))
            push!(C, A[i])
            i += 1
            n += 1
        end
        if n > length(A)
            n -= 1
        end
        while j <= length(B) && (B[j] <= A[n] || i > length(A))
            push!(C,B[j])
            j += 1
            m += 1
        end
        if m > length(B)
            m -= 1
        end
    end

    return C
end

function mergesort!(a)
    if length(a) == 1
        return a
    end
    a1 = mergesort!(a[1:end ÷ 2])
    a2 = mergesort!(a[end ÷ 2 + 1:end])
    merge(a1, a2)
end


text:
"""
    Laplace <: Kernel

TODO
"""
struct Laplace <: Kernel
    # TODO
end


text:
using IteratorAsArrayDisguises
using Test

needlessly_restricted_maximum(x::AbstractArray) = maximum(x)

@testset "needlessly restricted maximum" begin
    generator = (rand() for _ = 1:10)
    disguised = IteratorAsArrayDisguises.IteratorAsArrayDisguise(generator)
    @test needlessly_restricted_maximum(disguised) ≥ 0
end


text:
include("test.jl")
include("test_showvalues.jl")
include("test_map.jl")
include("test_float.jl")

println("")
println("All tests complete")


text:
using Pkg
Pkg.activate("./")
using RK

using Plots

function OutputSolution()
    # Output numerical solution in the end
    nz = length(history(tspan[1]))
    if nz == 1
        println("y(", tspan[2], ") = ", sol[2][end])
    else
        for i = 1 : nz
            println("y", i, "(", tspan[2], ") = ", (sol[2]')[end - nz + i])
        end
    end
end

# Example (Yukihiko)
idefun(t,y,z,i) = [ -2.5*i[1];
                    -7.5*i[2];
                    -15*i[3] ]
K(t,s,y)        = [ sin(y[1]);
                    sin(y[2]);
                    sin(y[3]) ]
delays_int(t)   = [ t-1;
                    t-1;
                    t-1 ]
history(t)      = [ 1.5;
                    1.5;
                    1.5 ]

tspan = [0 25]
stepsize = 1e-2

sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize)

OutputSolution()

display(plot(sol,
            linewidth=3,
            xlabel="TIME",
            ylabel="SOLUTION",
            xticks=trunc(Int,tspan[1]):(tspan[2]-trunc(Int,tspan[1]))/10:tspan[2],
            title="Example (Yukihiko)",
            legend=:outertopright))

max_stars_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_issues_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_forks_repo_name: vitesempl/IDE-Delay-Runge-Kutta


text:
using OrdinaryDiffEq, Flux, Plots
using Sundials
using DiffEqSensitivity
using ForwardDiff
using ForwardDiff: jacobian, jacobian!
using LinearAlgebra
using BandedMatrices
using Random
using Statistics
using ProgressBars, Printf
using Flux.Optimise:update!
using Flux.Losses:mae
using BSON: @save, @load
using Arrhenius
using YAML
using Base.Threads
using LatinHypercubeSampling
using DelimitedFiles
using SparseArrays

ENV["GKSwstype"] = "100"

runtime = YAML.load_file("./input.yaml")
expr_name = runtime["expr_name"]
is_restart = Bool(runtime["is_restart"])

conf = YAML.load_file("$expr_name/config.yaml")

mech = conf["mech"]
fuel = conf["fuel"]
fuel2air = Float64(conf["fuel2air"])
oxygen = conf["oxygen"]
inert = conf["inert"]
dT = Float64(conf["dTign"])
n_plot = Int64(conf["n_plot"])

if is_restart
    println("Continue to run $expr_name ...\n")
else
    println("Runing $expr_name ...\n")
end

fig_path = string(expr_name, "/figs")
ckpt_path = string(expr_name, "/checkpoint")

if !is_restart
    if ispath(fig_path)
        rm(fig_path, recursive=true)
    end
    if ispath(ckpt_path)
        rm(ckpt_path, recursive=true)
    end
end

if ispath(fig_path) == false
    mkdir(fig_path)
    mkdir(string(fig_path, "/conditions"))
end

if ispath(ckpt_path) == false
    mkdir(ckpt_path)
end

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
using StaticModules: StaticModule
using MLStyle: @match
using MacroTools: rmlines

abstract type AbstractFootingsElement end

Base.@kwdef struct Parameter <: AbstractFootingsElement
    name::Symbol
    type::Type
    description::String
end

Base.@kwdef struct Temporary <: AbstractFootingsElement
    name::Symbol
    type::Type
    description::String
end

Base.@kwdef struct Result <: AbstractFootingsElement
    name::Symbol
    type::Type
    description::String
end

Base.@kwdef struct Steps <: AbstractFootingsElement
    steps::Vector
end

Base.@kwdef struct Step
    fsym::Symbol
    assign::Symbol
    func::Function
    args::Tuple
end

Base.@kwdef struct Model
    name::String
    description::String = ""
    parameters::Vector = []
    temps::Vector = []
    results::Vector = []
    steps::Vector = []
end

function add_field!(model::Model, field::Parameter)
    push!(model.parameters, field)
    return nothing
end

function add_field!(model::Model, field::Temporary)
    push!(model.temps, field)
    return nothing
end

function add_field!(model::Model, field::Result)
    push!(model.results, field)
    return nothing
end

function add_field!(model::Model, field::Step)
    push!(model.steps, field)
    return nothing
end

function add_parameter!(model::Model; name::Symbol, type::Any, description::String)
    add_field!(model, Parameter(name = name, type = type, description = description))
    return nothing
end

function add_temporary!(model::Model; name::Symbol, type::Any, description::String)
    add_field!(model, Temporary(name = name, type = type, description = description))
    return nothing
end

function add_result!(model::Model; name::Symbol, type::Any, description::String)
    add_field!(model, Result(name = name, type = type, description = description))
    return nothing
end

function add_step!(
    model::Model,
    ns::StaticModule;
    assign::Symbol,
    func::Symbol,
    args::Tuple,
)
    func2 = getproperty(ns, func)
    add_field!(model, Step(fsym = func, assign = assign, func = func2, args = args))
    return nothing
end

function make_expr(step)
    return Expr(:(=), step.assign, Expr(:call, step.func, step.args...))
end

function create_func_parameters_expr(model::Model)
    p = [Expr(:(::), Symbol(p.name), p.type) for p in values(model.parameters)]
    return Expr(:parameters, p...)
end

function create_func_call_expr(model::Model, prefix::String)
    f_name = Symbol(prefix, "_", lowercase(model.name))
    f_params = create_func_parameters_expr(model)
    return Expr(:call, f_name, f_params)
end

function create_runner_expr(model::Model)
    steps = [make_expr(step) for step in model.steps]
    return Expr(:block, steps...)
end

function create_auditor_expr(model::Model)
    return Expr()
end

function build_documentation(model::Model)
    return "This is documentation."
end

function build_runner(model::Model, document::Bool = false)::Function
    f_call = create_func_call_expr(model, "run")
    f_body = create_runner_expr(model)
    f = eval(Expr(:function, f_call, f_body))
    if document
        docs = build_documentation(model)
        @doc docs f
    end
    return f
end

function build_auditor(model)::Function
    f_call = create_func_call_expr(model, "audit")
    f_body = create_auditor_expr(model)
    return eval(Expr(:function, f_call, f_body))
end

escape_args(args) = Tuple([esc(arg) for arg in args])

function _parse_expr(expr::Expr)
    function _parse(expr::Expr)
        @match expr begin
            Expr(:block, args) => _parse(args)
            Expr(:call, :Parameter, args...) =>
                :(add_parameter!($(:m), $(escape_args(args)...)))
            Expr(:call, :Temporary, args...) =>
                :(add_temporary!($(:m), $(escape_args(args)...)))
            Expr(:call, :Result, args...) => :(add_result!($(:m), $(escape_args(args)...)))
            Expr(:call, :Steps, args...) => map(_parse, args)
            Expr(:call, :(=>), assign, call) => :(add_step!(
                $(:m),
                $(:ns),
                assign = $(esc(assign)),
                func = $(QuoteNode(call.args[1])),
                args = $(Tuple(call.args[2:end])),
            ))
            Expr(:(=), assign, func, args...) => :(add_step!(
                $(:m),
                $(:ns),
                assign = $(esc(assign)),
                func = $(esc(func)),
                args = $args,
            ))
            _ => nothing
        end
    end
    return [_parse(e) for e in expr.args]
end

function _flatten(vec::Vector)
    v2 = []
    function _inner(x)
        if !isnothing(x)
            if isa(x, Expr)
                push!(v2, x)
            elseif isa(x, Vector)
                for xx in x
                    _inner(xx)
                end
            end
        end
    end
    for v in vec
        _inner(v)
    end
    return v2
end

parse_expr(expr) = (_flatten ∘ _parse_expr ∘ rmlines)(expr)

macro model(name::String, ns, expr::Expr)
    add_calls = parse_expr(expr)
    return (quote
        let ns = $(esc(ns))
            m = Model(name = $name)
            $(add_calls...)
            m
        end
    end |> rmlines)
end


text:
abstract type AbstractArc{Rule} end

# Mutable in order to make objectid fast
mutable struct Arc{Rule} <: AbstractArc{Rule}
    start::Int
    stop::Int
    rule::Rule
    constituents::LayeredVector{Union{Arc{Rule}, String}}
    score::Float64
end

Arc(start, stop, rule::R, constituents, score) where {R} = Arc{R}(start, stop, rule, constituents, score)

start(arc::Arc) = arc.start
stop(arc::Arc) = arc.stop
rule(arc::Arc) = arc.rule
score(arc::Arc) = arc.score
constituents(arc::Arc) = arc.constituents

num_constituents(arc::AbstractArc) = length(constituents(arc))
head(arc::AbstractArc) = lhs(rule(arc))

abstract type WrappedArc{Rule} <: AbstractArc{Rule} end

struct PassiveArc{Rule} <: WrappedArc{Rule}
    inner::Arc{Rule}
end

struct ActiveArc{Rule} <: WrappedArc{Rule}
    inner::Arc{Rule}
end

inner(arc::WrappedArc) = arc.inner

start(arc::WrappedArc) = start(inner(arc))
stop(arc::WrappedArc) = stop(inner(arc))
rule(arc::WrappedArc) = rule(inner(arc))
score(arc::WrappedArc) = score(inner(arc))
constituents(arc::WrappedArc) = constituents(inner(arc))

is_finished(arc::ActiveArc) = num_constituents(arc) == length(rhs(rule(arc)))
next_needed(arc::ActiveArc) = rhs(rule(arc))[num_constituents(arc) + 1]

function passive(arc::ActiveArc, scoring_function = arc -> 1.0)
    s = scoring_function(inner(arc))
    PassiveArc(
        Arc(start(arc),
            stop(arc),
            rule(arc),
            constituents(arc),
            score(arc) * s))
end
    # PassiveArc(inner(arc))

_show(io::IO, arc::AbstractArc) = print(io, arc)
_show(io::IO, s::AbstractString) = print(io, '"', s, '"')

function Base.show(io::IO, arc::AbstractArc)
    print(io, "($(start(arc)), $(stop(arc)), $(lhs(rule(arc))) -> ")
    for i in 1:length(rhs(rule(arc)))
        if i <= length(constituents(arc))
            _show(io, constituents(arc)[i])
            print(io, ' ')
        elseif i == length(constituents(arc)) + 1
            print(io, " . ")
        end
        if i > length(constituents(arc))
            print(io, rhs(rule(arc))[i])
            print(io, ' ')
        end
    end
    print(io, " ($(score(arc))))")
end

function geometric_mean(x)
    prod(x) ^ (1 / length(x))
end

"""
    combine(active::ActiveArc, passive::PassiveArc)

Combine two arcs according to the Fundamental Rule.
"""
function combine(a1::ActiveArc, a2::PassiveArc)
    new_constituents = push(constituents(a1), inner(a2))
    # new_score = score(rule(a1)) * geometric_mean(score.(new_constituents))
    # new_score = geometric_mean(vcat(score.(new_constituents), score(rule(a1))))
    # new_score = score(rule(a1)) * prod(score.(new_constituents))
    # new_score = score(rule(a1)) * geometric_mean(score.(new_constituents))

    # Geometric mean
    # constituent_score = reduce((s, arc) -> s * score(arc), new_constituents, init=1.0) ^ (1 / length(new_constituents))
    # new_score = score(rule(a1)) * constituent_score

    new_score = score(rule(a1)) * minimum(score, new_constituents)

    ActiveArc(Arc(start(a1), stop(a2), rule(a1), new_constituents, new_score))
end

combine(a1::PassiveArc, a2::ActiveArc) = combine(a2, a1)


text:
# # Empty result handling
#
# Transducible processes such as [`foldl`](@ref) try to
# do the right thing even when `init` is not given, _if_ the given
# binary operation `step` is supported by
# [InitialValues.jl](https://tkf.github.io/InitialValues.jl/dev/)
# (for example, `+`, `*`, `&`, and `|` are supported).  However, those
# functions _throw_ an exception if the given collection is empty or
# filtered out by the transducers:

using Test                                                             #src

using Transducers

add1(x) = x + 1

err = try                                                            # hide
foldl(*, Map(add1), [])
catch err; err; end                                                  # hide
@test err isa Transducers.EmptyResultError                           #src
#md showerror(stdout, err)                                           # hide

# To write robust code, it is recommended to use `init` if there is a
# reasonable default.  However, it may be useful to postpone
# "materializing" the result.  In such case, `Init` from InitialValues.jl
# can be used as a placeholder.

using InitialValues

result = foldl(*, Map(add1), [], init=Init)
nothing                                                              # hide

# `init=Init` is a short-hand notation of `init=Init(*)` (so that `*` does
# not have to be repeated):

@assert result === foldl(*, Map(add1), [], init=Init(*))

# Note also that [`transduce`](@ref) can be used for passing `init` as
# a positional argument:

@assert result === transduce(Map(add1), Completing(*), Init, [])

# Since the input collection `[]` is empty, `result` is `Init(*)` (which
# is an `InitialValues.InitialValue`):

using InitialValues: InitialValue
@assert result::InitialValue === Init(*)

# `Init(*)` is the left identity of `*`.  Multiplying it with any `x`
# from right returns `x` as-is.  This property may be useful, e.g., if
# `result` is known to be a scalar that is multiplied by a matrix just
# after the `foldl`:

@test ones(2, 2) ==                                                    #src
result * ones(2, 2)

# The identities `Init(*)` and `Init(+)` can be `convert`ed to numbers:

@test 1 ===                                                            #src
convert(Int, Init(*))

# `Init(*)` can also be `convert`ed to a `String`:

@test "" ===                                                           #src
convert(String, Init(*))

# This means that no special code is required if the result is going
# to be stored into, e.g., an `Array` or a `struct`:

@test begin                                                            #src
xs = [true, true]
xs[1] = Init(+)
xs
end == [false, true]                                                   #src

# They can be converted into numbers also by using `Integer`:

@test 0 ===                                                            #src
Integer(Init(+))

# or `float`:

@test 1.0 ===                                                          #src
float(Init(*))
#-


text:
# ---
# title: 948. Bag of Tokens
# id: problem948
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers, Greedy, Sort
# link: <https://leetcode.com/problems/bag-of-tokens/description/>
# hidden: true
# ---
#
# You have an initial **power** of `P`, an initial **score** of `0`, and a bag
# of `tokens` where `tokens[i]` is the value of the `ith` token (0-indexed).
#
# Your goal is to maximize your total **score** by potentially playing each
# token in one of two ways:
#
#   * If your current **power** is at least `tokens[i]`, you may play the `ith` token face up, losing `tokens[i]` **power** and gaining `1` **score**.
#   * If your current **score** is at least `1`, you may play the `ith` token face down, gaining `tokens[i]` **power** and losing `1` **score**.
#
# Each token may be played **at most** once and **in any order**. You do **not**
# have to play all the tokens.
#
# Return _the largest possible **score** you can achieve after playing any
# number of tokens_.
#
#
#
# **Example 1:**
#
#
#
#     Input: tokens = [100], P = 50
#     Output: 0
#     **Explanation** **:** Playing the only token in the bag is impossible because you either have too little power or too little score.
#
#
# **Example 2:**
#
#
#
#     Input: tokens = [100,200], P = 150
#     Output: 1
#     Explanation: Play the 0th token (100) face up, your power becomes 50 and score becomes 1.
#     There is no need to play the 1st token since you cannot play it face up to add to your score.
#
#
# **Example 3:**
#
#
#
#     Input: tokens = [100,200,300,400], P = 200
#     Output: 2
#     Explanation: Play the tokens in this order to get a score of 2:
#     1. Play the 0th token (100) face up, your power becomes 100 and score becomes 1.
#     2. Play the 3rd token (400) face down, your power becomes 500 and score becomes 0.
#     3. Play the 1st token (200) face up, your power becomes 300 and score becomes 1.
#     4. Play the 2nd token (300) face up, your power becomes 0 and score becomes 2.
#
#
#
#
# **Constraints:**
#
#   * `0 <= tokens.length <= 1000`
#   * `0 <= tokens[i], P < 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Generative
using Test

@testset "Generative.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: HamletWantToCode/Generative

max_issues_repo_name: HamletWantToCode/Generative

max_forks_repo_name: HamletWantToCode/Generative


text:
using Miletus.TermStructure
using Miletus.DayCounts
using Dates
using Test

@test 1.7    ≈ compound_factor(0.7, :Compounded, Miletus.TermStructure.Annual, Actual360(), today(), today()+Day(360))
@test 1.35^2 ≈ compound_factor(0.7, :Compounded, Miletus.TermStructure.Semiannual, Actual360(), today(), today()+Day(360))

@test 1/1.7  ≈ discount_factor(0.7, :Compounded, Miletus.TermStructure.Annual, Actual360(), today(), today()+Day(360))

d=ConstantYieldCurve(Actual360(), 0.2, :Continuous, :NoFrequency, today())

d=ConstantYieldCurve(Actual360(), 0.2, :Continuous, :NoFrequency, today())
@test 1/exp(0.2) ≈ d[today()+Day(360)]

@test forward_rate(d, today()+Day(360), today()+Day(360)) - 0.2 < 10e-10
forward_rate(d, Date(2016, Dec, 1), Date(2016, Dec, 31))


text:
@inline function ExaTronKernels.dscal(n::Int,da::Float64,dx::CuDeviceArray{Float64,1},incx::Int)
    tx = threadIdx().x

    # Ignore incx for now.
    if tx <= n
        @inbounds dx[tx] = da*dx[tx]
    end
    CUDA.sync_threads()

    return
end


text:
## Previous direct examples

pts1a = 180.0.+10.0*randn(100)
pts2a = 200.0.+10.0*randn(100)

pts1b = TU.wrapRad.(0.1*randn(100).-pi.+0.5)
pts2b = TU.wrapRad.(0.1*randn(100).+pi.-0.5)

pts1 = [pts1a';pts1b']
pts2 = [pts2a';pts2b']

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [3.0; 0.1], (+,AMP.addtheta), (-,AMP.difftheta))
pc2 = kde!(pts2, [3.0; 0.1], (+,AMP.addtheta), (-,AMP.difftheta))

dummy = kde!(rand(2,100),[1.0;], (+,AMP.addtheta), (-,AMP.difftheta));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1,
                          addop=(+,AMP.addtheta), diffop=(-,AMP.difftheta), getMu=(KDE.getEuclidMu,AMP.getCircMu));
#

lin1 = getBW(kde!(pGM[1,:]))[:,1][1]
cir1 = getBW(kde!_CircularNaiveCV(pGM[2,:]))[:,1][1]

pc12 = kde!(pGM, [lin1;cir1], (+,AMP.addtheta), (-,AMP.difftheta))

##

R_ = TU.R(10*pi/180.0)

pts1 = rand(MvNormal([100.0;0.0], R_*[40.0 0.0; 0.0 0.2]*(R_')), 100)
pts1[2,:] = TU.wrapRad.(pts1[2,:])
# plot(x=pts1[1,:],y=pts1[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

R_ = TU.R(-10*pi/180.0)

pts2 = rand(MvNormal([100.0;0.0], R_*[40.0 0.0; 0.0 0.2]*(R_')), 100)
pts2[2,:] = TU.wrapRad.(pts2[2,:])
plot(x=pts2[1,:],y=pts2[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

lin1 = getBW(kde!(pts1[1,:]))[:,1][1]
cir1 = getBW(kde!_CircularNaiveCV(pts1[2,:]))[:,1][1]

pc1 = kde!(pts1, [lin1; cir1], (+,AMP.addtheta), (-,AMP.difftheta))
pc2 = kde!(pts2, [lin1; cir1], (+,AMP.addtheta), (-,AMP.difftheta))

dummy = kde!(rand(2,100),[1.0;], (+,AMP.addtheta), (-,AMP.difftheta));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1,
                          addop=(+,AMP.addtheta), diffop=(-,AMP.difftheta), getMu=(KDE.getEuclidMu,AMP.getCircMu));
#

lin1 = getBW(kde!(pGM[1,:]))[:,1][1]
cir1 = getBW(kde!_CircularNaiveCV(pGM[2,:]))[:,1][1]

pc12 = kde!(pGM, [lin1;cir1], (+,AMP.addtheta), (-,AMP.difftheta))

plotKDE([pc1;pc2; pc12], dims=[1],c=["red";"green";"magenta"])

pl = plotKDECircular([marginal(pc1,[2]);marginal(pc2,[2]); marginal(pc12,[2])])

pl = Gadfly.plot(y=pGM[1,:],x=pGM[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

##  Do crazier single and multiple

R_ = TU.R(20*pi/180.0)

pts1 = rand(MvNormal([100.0;0.0], R_*[30.0 0.0; 0.0 0.2]*(R_')), 100)
pts1[2,:] = TU.wrapRad.(pts1[2,:])
# plot(x=pts1[1,:],y=pts1[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

pts2a = [5.0*randn(50).+85; 5.0*randn(50).+115]
pts2b = [TU.wrapRad.(0.1*randn(50).+2.6); TU.wrapRad.(0.1*randn(50).-2.6)]
pts2 = [pts2a';pts2b']

pts2 = pts2[:,randperm(100)]

# plot(x=pts2[1,:],y=pts2[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

lin1 = getBW(kde!(pts1[1,:]))[:,1][1]
cir1 = getBW(kde!_CircularNaiveCV(pts1[2,:]))[:,1][1]
pc1 = kde!(pts1, [lin1; cir1], (+,AMP.addtheta), (-,AMP.difftheta))

lin2 = getBW(kde!(pts2[1,:]))[:,1][1]
cir2 = getBW(kde!_CircularNaiveCV(pts2[2,:]))[:,1][1]
pc2 = kde!(pts2, [lin2; cir2], (+,AMP.addtheta), (-,AMP.difftheta))

dummy = kde!(rand(2,100),[1.0;], (+,AMP.addtheta), (-,AMP.difftheta));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1,
                          addop=(+,AMP.addtheta), diffop=(-,AMP.difftheta), getMu=(KDE.getEuclidMu,AMP.getCircMu));
#

lin1 = getBW(kde!(pGM[1,:]))[:,1][1]
cir1 = getBW(kde!_CircularNaiveCV(pGM[2,:]))[:,1][1]

pc12 = kde!(pGM, [lin1;cir1], (+,AMP.addtheta), (-,AMP.difftheta))

pl = Gadfly.plot(x=pGM[1,:],y=pGM[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

# plotKDE([pc1;pc2; pc12], dims=[1],c=["red";"green";"magenta"])

## testing Euclidean dimension separately

pc1a = kde!(pts1a)
pc2a = kde!(pts2a)

dummy = kde!(rand(100),[1.0;]);
pGM, = prodAppxMSGibbsS(dummy, [pc1a; pc2a], nothing, nothing, Niter=1, addop=(+,), diffop=(-,), getMu=(KDE.getEuclidMu,));

pc12a = kde!(pGM)
plotKDE([pc1a;pc2a; pc12a])

# Calculate the approximate product between the hybrid densities

pc1 = kde!(pts1, (+, addtheta), (-,difftheta))
pc2 = kde!(pts2, (+, addtheta), (-,difftheta))

# TODO: make circular KDE
dummy = kde!(rand(2,100),[1.0;]);
pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1, addop=(+,AMP.addtheta), diffop=(-,AMP.difftheta), getMu=(KDE.getEuclidMu, AMP.getCircMu));

pc12a = kde!(pGM[1,:])
getBW(pc12a)[:,1][1]
pc12b = kde!_CircularNaiveCV(pGM[2,:])
getBW(pc12b)[:,1][1]

pc12 = kde!(pGM, [1.4;0.01], (+,addtheta), (-,difftheta))

plotKDE([marginal(pc1,[1;]);marginal(pc2,[1;]); marginal(pc12,[1])])

pl = plotKDECircular([marginal(pc1,[2;]);marginal(pc2,[2;]); marginal(pc12,[2])])

Gadfly.plot((x)->marginal(pc12,[2])([x;])[1], -pi, pi)

## Doing a Pose2 example  var = [x;y;θ]

R_ = Matrix{Float64}(LinAlg.I, 3,3)
R_[1:2,1:2] = TU.R(-30*pi/180.0)

pts1 = rand(MvNormal([-50;100.0;2.5], R_*[20.0 0 0; 0.0 40.0 0.0; 0.0 0.0 0.05]*(R_')), 100)
pts1[3,:] = TU.wrapRad.(pts1[3,:])

plot(x=pts1[1,:],y=pts1[2,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pts1[2,:],y=pts1[3,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pts1[1,:],y=pts1[3,:], Geom.histogram2d(xbincount=50,ybincount=50))

lin1 = getBW(kde!(pts1[1:2,:]))[:,1]
pc1c = kde!_CircularNaiveCV(pts1[3,:])
cir1 = getBW(pc1c)[:,1][1]

pc1 = kde!(pts1, [lin1; cir1], (+,+,AMP.addtheta), (-,-,AMP.difftheta))

getBW(pc1c)[:,1]
plotKDECircular(pc1c)

pts2a = [5.0*randn(50).-85; 5.0*randn(50).-15]
pts2b = 20.0*randn(100).+130.0
pts2c = TU.wrapRad.(0.1*randn(100).-2.5)
pts2  = [pts2a';pts2b';pts2c']

pts2 = pts2[:,randperm(100)]

lin2 = getBW(kde!(pts2[1:2,:]))[:,1]
pc2c =kde!_CircularNaiveCV(pts2[3,:])
cir2 = getBW(pc2c)[:,1][1]

pc2 = kde!(pts2, [lin2; cir2], (+,+, AMP.addtheta), (-,-, AMP.difftheta))

plot(x=pts2[1,:],y=pts2[2,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pts2[2,:],y=pts2[3,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pts2[1,:],y=pts2[3,:], Geom.histogram2d(xbincount=50,ybincount=50))

dummy = kde!(rand(3,100),[1.0;], (+,+,AMP.addtheta), (-,-,AMP.difftheta));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1,
                          addop=(+,+,AMP.addtheta), diffop=(-,-,AMP.difftheta), getMu=(KDE.getEuclidMu,KDE.getEuclidMu,AMP.getCircMu));
#

lin12 = getBW(kde!(pGM[1:2,:]))[:,1]
pc12c= kde!_CircularNaiveCV(pGM[3,:])
cir12 = getBW(pc12c)[:,1]

pc12 = kde!(pGM, [lin12;cir12], (+,+,AMP.addtheta), (-,-,AMP.difftheta))

plot(x=pGM[1,:],y=pGM[2,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pGM[2,:],y=pGM[3,:], Geom.histogram2d(xbincount=50,ybincount=50))
plot(x=pGM[1,:],y=pGM[3,:], Geom.histogram2d(xbincount=50,ybincount=50))

plotKDE([pc1;pc2; pc12], dims=[1], c=["red";"green";"magenta"])
plotKDE([pc1;pc2; pc12], dims=[2], c=["red";"green";"magenta"])

bw12 = getBW(pc12)[:,1]

pc12ab = marginal(pc12,[1;2])
# getBW(pc12ab)[:,1]

plotKDE(pc12ab, levels=3)

pl = plotKDECircular([pc1c;pc2c;pc12c])

0

##

pl |> PDF("/tmp/test.pdf",20cm,15cm)
@async run(`evince /tmp/test.pdf`)

#


text:
export MicMac

struct MicMac

    data :: DataSet

end

function run(micmac :: MicMac, dt, ntau, schema_micmac)

    Tfinal    = micmac.data.Tfinal
    N         = micmac.data.N
    T         = micmac.data.T
    k         = transpose(micmac.data.k)

    u         = collect(transpose(micmac.data.u))
    v         = collect(transpose(micmac.data.v))

    epsilon   = micmac.data.epsilon

    dx        = micmac.data.dx

    llambda   = micmac.data.llambda
    sigma     = micmac.data.sigma

    tau       = zeros(Float64, ntau)
    tau      .= T * collect(0:ntau-1) / ntau
    ktau      = similar(tau)
    ktau     .= 2 * pi / T * vcat(0:ntau÷2-1,-ntau÷2:-1)
    ktau[1]   = 1.0

    matr      = zeros(ComplexF64,ntau)
    conjmatr  = zeros(ComplexF64,ntau)
    matr     .= exp.( 1im * tau)
    conjmatr .= exp.(-1im * tau)

    A1 = zeros(Float64, (1,N))
    A2 = zeros(Float64, (1,N))

    if epsilon > 0
        A1 .= (sqrt.(1 .+ epsilon * k.^2) .- 1) / epsilon
        A2 .= (1 .+ epsilon * k.^2) .^ (-1/2)
    else
        A1 .= 0.5 * k .^ 2
        A2 .= 1
    end

    t    = 0
    iter = 0
    tabt = [t]

    fft_u = zeros(ComplexF64,(ntau,N))
    fft_v = zeros(ComplexF64,(ntau,N))

    if schema_micmac == 0

        fft_u = fft(u)
        fft_v = fft(v)

    elseif schema_micmac == 1

        ubar = u
        vbar = v
        ug = 0 * u
        vg = 0 * v
        fft_ubar = fft(ubar)
        fft_vbar = fft(vbar)
        fft_ug = fft(ug)
        fft_vg = fft(vg)

    elseif (schema_micmac == 2 || schema_micmac == 25)

        fft_u0 = fft(u)
        fft_v0 = fft(v)

        println("fft_u0:", size(fft_u0))

        fft_ubar, fft_vbar, fft_ug, fft_vg = init_2(0, fft_u0, fft_v0, A1, A2,
                                                    matr, conjmatr,
                                                    sigma,
                                                    llambda,
                                                    ktau, epsilon, ntau)

    elseif (schema_micmac == 3 || schema_micmac == 35)

        fft_u0 = fft(u)
        fft_v0 = fft(v)
        fft_ubar, fft_vbar, fft_ug, fft_vg = init_3(0, fft_u0, fft_v0, A1, A2,
                                                    matr, conjmatr,
                                                    sigma,
                                                    llambda,
                                                    ktau, epsilon, ntau)

    elseif (schema_micmac == 4)

        fft_u0 = fft(u)
        fft_v0 = fft(v)

        fft_ubar, fft_vbar, fft_ug, fft_vg = init_4(0, fft_u0, fft_v0, A1, A2,
                                                    matr, conjmatr,
                                                    sigma,
                                                    llambda,
                                                    ktau, epsilon, ntau)
    end

    while t < Tfinal

        iter = iter + 1
        dt, lstat = adjust_step(Tfinal, t, dt)

        if schema_micmac == 0

            champu, champv, champumoy, champvmoy = ichampf(t, fft_u, fft_v, A1, A2,
                                                           ntau, matr, conjmatr,
                                                           ktau, llambda,
                                                           sigma)
            fft_u = fft_u + epsilon * reconstr(champu, (t + dt) / epsilon, T, ntau) \
                    - epsilon * reconstr(champu, t / epsilon, T, ntau) + dt * champumoy
            fft_v = fft_v + epsilon * reconstr(champv, (t + dt) / epsilon, T, ntau) \
                    - epsilon * reconstr(champv, t / epsilon, T, ntau) + dt * champvmoy
            t = t + dt

        elseif schema_micmac == 2

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_2(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr, conjmatr,
                         sigma, llambda,
                         ktau, epsilon, ntau)
            fft_ubar12 = fft_ubar + dt / 2 * champubaru
            fft_vbar12 = fft_vbar + dt / 2 * champubarv
            fft_ug12 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg12 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_2(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12, A1, A2,
                         matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar = fft_ubar + dt * champubaru
            fft_vbar = fft_vbar + dt * champubarv
            fft_ug = fft_ug + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                     - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt * champmoyu
            fft_vg = fft_vg + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                     - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt * champmoyv

            t = t + dt
            C1u, C1v = C1(t, fft_ubar, fft_vbar, A1, A2, matr, conjmatr,
                          sigma, llambda,
                          ktau, epsilon, ntau)
            uC1eval = reconstr(C1u, t / epsilon, T, ntau)
            vC1eval = reconstr(C1v, t / epsilon, T, ntau)
            fft_u = uC1eval + fft_ug
            fft_v = vC1eval + fft_vg

        elseif schema_micmac == 25

            fft_ubartmp = fft_ubar
            fft_vbartmp = fft_vbar
            fft_ugtmp   = fft_ug
            fft_vgtmp   = fft_vg
            dw          = 1
            iter2       = 0

            while (dw > 1.e-12 && iter2 < 100)

                iter2        = iter2 + 1
                fft_ubartmp2 = fft_ubartmp
                fft_vbartmp2 = fft_vbartmp
                fft_ugtmp2   = fft_ugtmp
                fft_vgtmp2   = fft_vgtmp

                champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                    champs_2(t + dt / 2, (fft_ubar + fft_ubartmp) / 2,
                             (fft_vbar + fft_vbartmp) / 2,
                             (fft_ug + fft_ugtmp) / 2,
                             (fft_vg + fft_vgtmp) / 2, A1, A2, matr, conjmatr,
                             sigma, llambda, ktau, epsilon, ntau)
                fft_ubartmp = fft_ubar + dt * champubaru
                fft_vbartmp = fft_vbar + dt * champubarv
                fft_ugtmp = fft_ug + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                            - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt * champmoyu
                fft_vgtmp = fft_vg + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                            - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt * champmoyv
                dw = linalg.norm(fft_ubartmp - fft_ubartmp2) \
                     + linalg.norm(fft_vbartmp - fft_vbartmp2) \
                     + linalg.norm(fft_ugtmp2 - fft_ugtmp) \
                     + linalg.norm(fft_vgtmp2 - fft_vgtmp)
            end

            fft_ubar = fft_ubartmp
            fft_vbar = fft_vbartmp
            fft_ug = fft_ugtmp
            fft_vg = fft_vgtmp
            t = t + dt

            C1u, C1v = C1(t, fft_ubar, fft_vbar, A1, A2, matr, conjmatr, sigma,
                          llambda, ktau, epsilon, ntau)
            uC1eval = reconstr(C1u, t / epsilon, T, ntau)
            vC1eval = reconstr(C1v, t / epsilon, T, ntau)
            fft_u = uC1eval + fft_ug
            fft_v = vC1eval + fft_vg

        elseif schema_micmac == 3

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr, conjmatr,
                         sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar14 = fft_ubar + dt / 4 * champubaru
            fft_vbar14 = fft_vbar + dt / 4 * champubarv
            fft_ug14 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 4 * champmoyu
            fft_vg14 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 4 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t + dt / 4, fft_ubar14, fft_vbar14, fft_ug14, fft_vg14, A1, A2,
                         matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar12 = fft_ubar + dt / 2 * champubaru
            fft_vbar12 = fft_vbar + dt / 2 * champubarv
            fft_ug12 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg12 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12, A1, A2, matr, conjmatr,
                         sigma, llambda, ktau, epsilon, ntau)

            fft_ubar34 = fft_ubar12 + dt / 4 * champubaru
            fft_vbar34 = fft_vbar12 + dt / 4 * champubarv
            fft_ug34 = fft_ug12 + epsilon * reconstr(ichampgu, (t + 3 * dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) + dt / 4 * champmoyu
            fft_vg34 = fft_vg12 + epsilon * reconstr(ichampgv, (t + 3 * dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) + dt / 4 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t + 3 * dt / 4, fft_ubar34, fft_vbar34, fft_ug34, fft_vg34,
                         A1, A2, matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar1 = fft_ubar12 + dt / 2 * champubaru
            fft_vbar1 = fft_vbar12 + dt / 2 * champubarv
            fft_ug1 = fft_ug12 + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg1 = fft_vg12 + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr, conjmatr,
                         sigma, llambda, ktau, epsilon, ntau)

            fft_ubar12 = fft_ubar + dt / 2 * champubaru
            fft_vbar12 = fft_vbar + dt / 2 * champubarv
            fft_ug12 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg12 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_3(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12, A1, A2,
                         matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar1bis = fft_ubar + dt * champubaru
            fft_vbar1bis = fft_vbar + dt * champubarv
            fft_ug1bis = fft_ug + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                         - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt * champmoyu
            fft_vg1bis = fft_vg + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                         - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt * champmoyv

            t = t + dt
            fft_ubar = (4 * fft_ubar1 - fft_ubar1bis) / 3
            fft_vbar = (4 * fft_vbar1 - fft_vbar1bis) / 3
            fft_ug = (4 * fft_ug1 - fft_ug1bis) / 3
            fft_vg = (4 * fft_vg1 - fft_vg1bis) / 3

            C2u, C2v = C2(t, fft_ubar, fft_vbar, A1, A2, matr, conjmatr,
                          sigma, llambda, ktau, epsilon, ntau)
            uC2eval = reconstr(C2u, t / epsilon, T, ntau)
            vC2eval = reconstr(C2v, t / epsilon, T, ntau)
            fft_u = uC2eval + fft_ug
            fft_v = vC2eval + fft_vg

        elseif schema_micmac == 4

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr, conjmatr,
                         sigma, llambda, ktau, epsilon, ntau)

            fft_ubar16 = fft_ubar + dt / 6 * champubaru
            fft_vbar16 = fft_vbar + dt / 6 * champubarv
            fft_ug16 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 6) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 6 * champmoyu
            fft_vg16 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 6) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 6 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 6, fft_ubar16, fft_vbar16, fft_ug16, fft_vg16,
                         A1, A2, matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar13 = fft_ubar + dt / 3 * champubaru
            fft_vbar13 = fft_vbar + dt / 3 * champubarv
            fft_ug13 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 3) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 3 * champmoyu
            fft_vg13 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 3) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 3 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 3, fft_ubar13, fft_vbar13, fft_ug13, fft_vg13,
                         A1, A2, matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar12 = fft_ubar13 + dt / 6 * champubaru
            fft_vbar12 = fft_vbar13 + dt / 6 * champubarv
            fft_ug12 = fft_ug13 + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, (t + dt / 3) / epsilon, T, ntau) + dt / 6 * champmoyu
            fft_vg12 = fft_vg13 + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, (t + dt / 3) / epsilon, T, ntau) + dt / 6 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12,
                         A1, A2, matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar23 = fft_ubar13 + dt / 3 * champubaru
            fft_vbar23 = fft_vbar13 + dt / 3 * champubarv
            fft_ug23 = fft_ug13 + epsilon * reconstr(ichampgu, (t + 2 * dt / 3) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, (t + dt / 3) / epsilon, T, ntau) + dt / 3 * champmoyu
            fft_vg23 = fft_vg13 + epsilon * reconstr(ichampgv, (t + 2 * dt / 3) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, (t + dt / 3) / epsilon, T, ntau) + dt / 3 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + 2 * dt / 3, fft_ubar23, fft_vbar23, fft_ug23, fft_vg23,
                         A1, A2, matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar56 = fft_ubar23 + dt / 6 * champubaru
            fft_vbar56 = fft_vbar23 + dt / 6 * champubarv
            fft_ug56 = fft_ug23 + epsilon * reconstr(ichampgu, (t + 5 * dt / 6) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, (t + 2 * dt / 3) / epsilon, T, ntau) + dt / 6 * champmoyu
            fft_vg56 = fft_vg23 + epsilon * reconstr(ichampgv, (t + 5 * dt / 6) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, (t + 2 * dt / 3) / epsilon, T, ntau) + dt / 6 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + 5 * dt / 6, fft_ubar56, fft_vbar56, fft_ug56, fft_vg56,
                         A1, A2, matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar0 = fft_ubar23 + dt / 3 * champubaru
            fft_vbar0 = fft_vbar23 + dt / 3 * champubarv
            fft_ug0 = fft_ug23 + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgu, (t + 2 * dt / 3) / epsilon, T, ntau) + dt / 3 * champmoyu
            fft_vg0 = fft_vg23 + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgv, (t + 2 * dt / 3) / epsilon, T, ntau) + dt / 3 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2, matr, conjmatr,
                         sigma, llambda, ktau, epsilon, ntau)

            fft_ubar14 = fft_ubar + dt / 4 * champubaru
            fft_vbar14 = fft_vbar + dt / 4 * champubarv
            fft_ug14 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 4 * champmoyu
            fft_vg14 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 4 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 4, fft_ubar14, fft_vbar14, fft_ug14, fft_vg14, A1, A2,
                         matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar12 = fft_ubar + dt / 2 * champubaru
            fft_vbar12 = fft_vbar + dt / 2 * champubarv
            fft_ug12 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg12 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12, A1, A2,
                         matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar34 = fft_ubar12 + dt / 4 * champubaru
            fft_vbar34 = fft_vbar12 + dt / 4 * champubarv
            fft_ug34 = fft_ug12 + epsilon * reconstr(ichampgu, (t + 3 * dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) + dt / 4 * champmoyu
            fft_vg34 = fft_vg12 + epsilon * reconstr(ichampgv, (t + 3 * dt / 4) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) + dt / 4 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + 3 * dt / 4, fft_ubar34, fft_vbar34, fft_ug34, fft_vg34,
                         A1, A2, matr, conjmatr, sigma, llambda, ktau, epsilon, ntau)

            fft_ubar1 = fft_ubar12 + dt / 2 * champubaru
            fft_vbar1 = fft_vbar12 + dt / 2 * champubarv
            fft_ug1 = fft_ug12 + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg1 = fft_vg12 + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t, fft_ubar, fft_vbar, fft_ug, fft_vg, A1, A2,
                         matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar12 = fft_ubar + dt / 2 * champubaru
            fft_vbar12 = fft_vbar + dt / 2 * champubarv
            fft_ug12 = fft_ug + epsilon * reconstr(ichampgu, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt / 2 * champmoyu
            fft_vg12 = fft_vg + epsilon * reconstr(ichampgv, (t + dt / 2) / epsilon, T, ntau) \
                       - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt / 2 * champmoyv

            champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv = \
                champs_4(t + dt / 2, fft_ubar12, fft_vbar12, fft_ug12, fft_vg12,
                         A1, A2, matr, conjmatr, sigma, llambda,
                         ktau, epsilon, ntau)

            fft_ubar2 = fft_ubar + dt * champubaru
            fft_vbar2 = fft_vbar + dt * champubarv
            fft_ug2 = fft_ug + epsilon * reconstr(ichampgu, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgu, t / epsilon, T, ntau) + dt * champmoyu
            fft_vg2 = fft_vg + epsilon * reconstr(ichampgv, (t + dt) / epsilon, T, ntau) \
                      - epsilon * reconstr(ichampgv, t / epsilon, T, ntau) + dt * champmoyv

            t = t + dt
            fft_ubar = (27 * fft_ubar0 - 16 * fft_ubar1 + fft_ubar2) / 12
            fft_vbar = (27 * fft_vbar0 - 16 * fft_vbar1 + fft_vbar2) / 12
            fft_ug = (27 * fft_ug0 - 16 * fft_ug1 + fft_ug2) / 12
            fft_vg = (27 * fft_vg0 - 16 * fft_vg1 + fft_vg2) / 12

            (C3u, C3v) = C3(t, fft_ubar, fft_vbar, A1, A2, matr, conjmatr,
                            sigma, llambda, ktau, epsilon, ntau)
            uC3eval = reconstr(C3u, t / epsilon, T, ntau)
            vC3eval = reconstr(C3v, t / epsilon, T, ntau)
            fft_u = uC3eval + fft_ug
            fft_v = vC3eval + fft_vg

        end

    end

    fft_u = exp(1im * sqrt(1 + epsilon * k^2) * t / epsilon) * fft_u
    fft_v = exp(1im * sqrt(1 + epsilon * k^2) * t / epsilon) * fft_v

    u = ifft(fft_u)
    v = ifft(fft_v)

    u, v

end


text:
abstract type AbstractGeometricalAxisWeights{T} <: AbstractArray{T, 2} end

function sizeof(gw::AbstractGeometricalAxisWeights{T})::Int where {T}
    return sizeof(gw.weights)
end
@inline function size(gw::AbstractGeometricalAxisWeights{T}) where {T}
    return size(gw.weights)
end

@inline getindex(gw::AbstractGeometricalAxisWeights{T}, I::Vararg{Int, N}) where {T, N} = getindex(gw.weights, I...)

"""
    struct GeometricalCartesianAxisWeights{T} <: AbstractGeometricalAxisWeights{T}

`GeometricalCartesianAxisWeights` stores certain precalculated (in its construction) fixed values
of one of the cartesian (linear) axes, `ax`, of the grid on which the field calculation is performed.
These precalculated values are further used in the field calculation
to calculate the six weights for the neighboring grid points in the SOR.

`GeometricalCartesianAxisWeights` has only one field: `weights::Array{T}, 2`, which is
of size `(4, length(ax.ticks))`. Each row holding the 4 precalculated values for one axis tick.

Axis ticks, `t = ax.ticks`, and midpoints, `mp = midpoints(get_extended_ticks(ax))`,\\
(the middle points between to axis ticks) inbetween:\\
`... t[i-1] --- mp[i] --- t[i] --- mp[i+1] --- t[i+1] ...`\\
are required for the understanding of the precalculated values.

The columns store the following quantities:
* `weights[1, i]`: `(mp[i+1] - t[i])` / `weights[3, i]`
* `weights[2, i]`: `(t[i] - mp[i])` / `weights[3, i]`
* `weights[3, i]`: `mp[i+1] - mp[i]`
* `weights[4, i]`: `inv(t[i] - t[i-1])`
"""
struct GeometricalCartesianAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end

"""
    struct GeometricalAzimutalAxisWeights{T} <: AbstractGeometricalAxisWeights{T}

`GeometricalAzimutalAxisWeights` stores certain precalculated (in its construction) fixed values
of the azimulal (or polar) axis, `ax`, of an cylindrical grid on which the field calculation is performed.
These precalculated values are further used in the field calculation
to calculate the six weights for the neighboring grid points in the SOR.

`GeometricalAzimutalAxisWeights` has only one field: `weights::Array{T}, 2`, which is
of size `(4, length(ax.ticks))`. Each row holding the 4 precalculated values for one axis tick.

Axis ticks, `t = ax.ticks`, and midpoints, `mp = midpoints(get_extended_ticks(ax))`,\\
(the middle points between to axis ticks) inbetween:\\
`... t[i-1] --- mp[i] --- t[i] --- mp[i+1] --- t[i+1] ...`\\
are required for the understanding of the precalculated values.

The columns store the following quantities:
* `weights[1, i]`: `(mp[i+1] - t[i])` / `weights[3, i]`
* `weights[2, i]`: `(t[i] - mp[i])` / `weights[3, i]`
* `weights[3, i]`: `mp[i+1] - mp[i]`
* `weights[4, i]`: `inv(t[i] - t[i-1])`

Developer note: This actually seems to be the same as `GeometricalCartesianAxisWeights`.
Thus, we actually could remove that. Or change this into an alias.
"""
struct GeometricalAzimutalAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end

"""
    struct GeometricalRadialAxisWeights{T} <: AbstractGeometricalAxisWeights{T}

`GeometricalRadialAxisWeights` stores certain precalculated (in its construction) fixed values
of the radial axis, `ax`, of an cylindrical grid on which the field calculation is performed.
These precalculated values are further used in the field calculation
to calculate the six weights for the neighboring grid points in the SOR.

`GeometricalRadialAxisWeights` has only one field: `weights::Array{T}, 2`, which is
of size `(6, length(ax.ticks))`. Each row holding the 6 precalculated values for one axis tick.

Axis ticks, `t = ax.ticks`, and midpoints, `mp = midpoints(get_extended_ticks(ax))`,\\
(the middle points between to axis ticks) inbetween:\\
`... t[i-1] --- mp[i] --- t[i] --- mp[i+1] --- t[i+1] ...`\\
are required for the understanding of the precalculated values.

The columns store the following quantities:
* `weights[1, i]`: `(mp[i+1] - t[i])` / `weights[3, i]`
* `weights[2, i]`: `(t[i] - mp[i])` / `weights[3, i]`
* `weights[3, i]`: `(mp[i+1] - mp[i]) / t[i]`
* `weights[4, i]`: `mp[i+1] / (t[i+1] - t[i])`
* `weights[5, i]`: `mp[i] / (t[i] - t[i-1])`\\
  # Developer note: maybe this one could be removed as it basically stores the same as `weights[4, i-1]`
* `weights[6, i]`: `(mp[i+1]^2 - mp[i]^2)/2`

Developer note: Some of the weights for `i == 1` are set manually as `r = 0` requires some special treatment.
"""
struct GeometricalRadialAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end


text:
using Test
using QuantumLattices.Prerequisites.TypeTraits

struct WithTrait{F1, F2, F3}
    f1::F1
    f2::F2
    f3::F3
end
Base.:(==)(fc1::WithTrait, fc2::WithTrait) = ==(efficientoperations, fc1, fc2)
Base.isequal(fc1::WithTrait, fc2::WithTrait) = isequal(efficientoperations, fc1, fc2)
Base.:<(fc1::WithTrait, fc2::WithTrait) = <(efficientoperations, fc1, fc2)
Base.isless(fc1::WithTrait, fc2::WithTrait) = isless(efficientoperations, fc1, fc2)
Base.isapprox(fc1::WithTrait, fc2::WithTrait; atol::Real=10^-5, rtol::Real=10^-5) = isapprox(efficientoperations, Val((:f1, :f2)), fc1, fc2; atol=atol, rtol=rtol)
Base.replace(fc::WithTrait; kwargs...) = replace(efficientoperations, fc; kwargs...)

@testset "efficientoperations" begin
    @test ==(efficientoperations, (), ())
    @test isequal(efficientoperations, (), ())
    @test ==(efficientoperations, (1, 2), (1, 2, 3)) == false
    @test isequal(efficientoperations, (1, 2), (1, 2, 3)) == false

    fc1, fc2 = WithTrait(1.0, 2, 3), WithTrait(1, 2, 3)
    @test fc1 == fc2
    @test isequal(fc1, fc2)
    fc1, fc2 = WithTrait(1.0, 2, 3), WithTrait(1, 2, 4.0)
    @test fc1 < fc2
    @test isless(fc1, fc2)

    fc1, fc2 = WithTrait(1.0+10^-6, 2, 3), WithTrait(1, 2, 3)
    @test fc1 ≈ fc2
    fc1, fc2 = WithTrait(1.0, 2-10^-6, 3), WithTrait(1, 2, 3)
    @test fc1 ≈ fc2
    fc1, fc2 = WithTrait(1.0, 2, 3-10^-6), WithTrait(1, 2, 3)
    @test fc1 ≉  fc2

    @test replace(WithTrait(1, 2, 3), f1='c') == WithTrait('c', 2, 3)
end

@testset "forder/corder" begin
    dims = (2, 2, 2)
    finds = [(1, 1, 1), (2, 1, 1), (1, 2, 1), (2, 2, 1), (1, 1, 2), (2, 1, 2), (1, 2, 2), (2, 2, 2)]
    cinds = [(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2), (2, 1, 1), (2, 1, 2), (2, 2, 1), (2, 2, 2)]
    @test collect(indtosub(dims, i, forder) for i = 1:prod(dims)) == finds
    @test collect(indtosub(dims, i, corder) for i = 1:prod(dims)) == cinds
    @test collect(subtoind(dims, inds, forder) for inds in finds) == collect(1:prod(dims))
    @test collect(subtoind(dims, inds, corder) for inds in cinds) == collect(1:prod(dims))
end


text:
module FedeDTW
using Distances: PreMetric
using Statistics: median, quantile
using PyCall: pyimport

"""
Computes the 1-1 average difference at various shift and penalizes each difference by
elevating to `1/shift`.

Returns the meainum
"""
function win_dist(
    arr1::Vector{Vector{Int64}},
    arr2::Vector{Vector{Int64}},
    dist_fn::PreMetric,
)
    # compute distance matrix
    # for j in 1:length(arr2)
    #     @simd for i in 1:length(arr1)
    #         diff[i, j] = dist_fn(arr1[i], arr2[j])
    #     end
    # end

    # # apply hungarian algorithm and find all indices
    # indices = findall(==(Hungarian.STAR), Hungarian.munkres(diff))

    # # compute cost of the matching
    # return median(diff[indices])
    # return quantile(diff[indices], 0.25)
    #
    out = Inf
    K = length(arr1)
    K1 = K + 1
    @simd for shift = 1:K
        v = mean(dist_fn.(arr1[shift:K], arr2[1:K1-shift]))^(1 / shift)
        if v < out
            out = v
        end
    end
    return out
end

"""
given an idx, a radius and a maximum length, returns starting and ending
indices of a a window centered at that idx and having that radius, without
indices > length nor < 0
"""
function idx_range(idx, radius, length)
    return max(1, idx - radius), min(length, idx + radius)
end

struct FedeWindowParam
    dist::PreMetric
    min_radius::Int64
    α::Int64
    β::Int64
end

"""
A windowing function which computes a different slanted-band at each point
based on the local difference of the main slanted diagonal; the local
radius is computed as:

`max(
   min_radius,
   floor(
       α * avg_dist_fn(
               x[i - β : i + β],
               y[j - β : j + β]
           )
   )
)`

where:
* N is the length of x
* M is the length of y
* avg_dist_fn is the average of dist_fn on each corresponding sample
* j = floor(i * M / N)

By words, `β` is half the length of a sliding window used to compute
distances between excerpts of `x` and `y` taken along the slanted diagonal.
The distance is multiplied by `α` to get the local radius length.

`x` and `y` are sequences with shape ([M, N], features)
"""
struct FedeWindowFunc <: Function
    mask::Matrix{Bool}
end

# adding the callable action
(self::FedeWindowFunc)(i::Int64, j::Int64)::Bool = self.mask[i, j]

function get_fede_win(
    params::FedeWindowParam,
    x::Vector{Vector{T}},
    y::Vector{Vector{T}},
) where {T<:Number}
    # take the distance function
    N = length(x)
    M = length(y)
    _transpose = false
    if M > N
        # x should always be longer than y
        x, y = y, x
        N, M = M, N
        # if we swap x and y, we need to swap the mask too
        _transpose = true
    end

    # a matrix where distances between windows are computed
    # R = params.β * 2 + 1
    # diff = fill(Inf, R, R)

    # a mask to remember points
    mask = falses(N, M)

    # for each point in x
    for i = 1:N
        # compute the point in y along the diagonal
        j = floor(Int64, (i - 1) * M / N) + 1

        # compute the sliding windows
        start_x, end_x = idx_range(i, params.β, N)
        start_y, end_y = idx_range(j, params.β, M)
        _x = x[start_x:end_x]
        _y = y[start_y:end_y]

        # pad the windows
        if start_x == 1
            pushfirst!(_x, [[0] for k = 1:(params.β-i+1)]...)
        end
        if end_x == N
            push!(_x, [[0] for k = 1:(i+params.β-N)]...)
        end
        if start_y == 1
            pushfirst!(_y, [[0] for k = 1:(params.β-j+1)]...)
        end
        if end_y == M
            push!(_y, [[0] for k = 1:(j+params.β-M)]...)
        end

        # compute the local radius
        lr = max(
            params.min_radius,
            floor(Int64, params.α * quantile(params.dist.(_x, _y), 0.25)),
        )
        # floor(Int64, params.α * win_dist(_x, _y, params.dist)))

        # set the points inside the local radius to True
        idx_i = idx_range(i, lr, N)
        idx_j = idx_range(j, lr, M)
        mask[idx_i[1]:idx_i[2], idx_j[1]:idx_j[2]] .= true
    end

    if _transpose
        mask = transpose(mask)
    end

    return FedeWindowFunc(mask)
end

"""
Returns a Vector of Tuples. Each tuple contains a value from `x` and a Vector
with the indices where that value is in `x`. Only non-unique values are
returned.
"""
function my_nonunique(x::AbstractArray{T}) where {T}
    # counting each value and remembering indices
    indices = Dict{T,Vector{Int64}}()
    for (idx, val) in enumerate(x)
        if haskey(indices, val)
            # this is a duplicated value, push its index
            push!(indices[val], idx)
        else
            # this is the first time we meet this value
            indices[val] = [idx]
        end
    end
    # now filter out unique values
    out = Tuple{T,Vector{Int64}}[]
    for (key, val) in indices
        if length(val) > 1
            push!(out, (key, val))
        end
    end
    return out
end

"""
1. look for repeated values in `arr_x` or `arr_y`, depending on `target`
2. look for the maximum value in `graph_matrix[2]`, at the indices in
`arr_x` and `arr_y` relative to the repeated values
3. among the repeated values in the target, chose the ones corresponding to
the maximum in `graps_matrix[2]`
4. return `arr_x` and `arr_y` without the removed indices
"""
@views function _remove_conflicting_match(
    # TODO: rewrite in the light of weights...
    arr_x,
    arr_y,
    graph_matrix::AbstractArray{Int64,2},
    target::Int64,
)

    if target == 0
        _target = arr_x
    elseif target == 1
        _target = arr_y
    end

    arr_mask = ones(Bool, size(_target, 1))
    nonunique = my_nonunique(_target)
    for (val, conflicting_idx) in nonunique
        to_keep_idx_of_idx =
            argmax(graph_matrix[arr_x[conflicting_idx], arr_y[conflicting_idx]])[2]
        arr_mask[conflicting_idx] .= false
        arr_mask[conflicting_idx[to_keep_idx_of_idx]] = true
    end

    return arr_x[arr_mask], arr_y[arr_mask]
end

"""
OLD

Takes a list of mapping indices, fills the graph matrix counting the number
of times a match happens in the mappings. Then start taking matching from
the most matched and iteratively adding new matching. If two conflicting
matching have the same number of counts, takes the matching which appears
in the longest mapping; in case of parity the first one is taken
"""
@views function merge_matching_indices(args::Vector{Matrix{Int64}}, weights::AbstractArray{Float64})::Matrix{Int64}

    # voting
    votes::Dict{Vector{Int64}, Float64} = Dict()
    for i in 1:length(args)
        weight = weights[i]
        for match in eachcol(args[i]')
            if haskey(votes, match)
                votes[match] += weight
            else
                votes[match] = weight
            end
        end
    end
    if length(votes) == 0
        # no matches found!
        return Matrix{Float64}(undef, 0, 0)
    end

    # sorting by vote
    sorted_votes::Vector{Pair{Vector{Int64}, Float64}} = sort(
        collect(votes); by=x->x[2], rev=true)

    # number of notes
    matches = getindex.(sorted_votes, 1)
    N = maximum(getindex.(matches, 1))
    M = maximum(getindex.(matches, 2))

    # taking unique matches
    rows = falses(N)
    cols = falses(M)
    merged = Vector{Vector{Int64}}(undef, 0)
    for (match, vote) in sorted_votes
        row, col = match
        if !rows[row] && !cols[col]
            rows[row] = true
            cols[col] = true
            push!(merged, match)
        end
    end
    # re-sort everything and return
    out = vcat(merged'...)
    # print(f"Added notes from merging: {len(ref) - L}")
    return out[sortperm(out[:, 1]), :]
end
end # module


text:
@testset "Parameterized Models" begin

    @testset "Makeham" begin

        g = Gompertz(a=0.0002, b=.13)

        @test survival(g, 45) ≈ 0.5870365016720939
        @test survival(g, 45, 46) ≈ 0.9285202788707242
        @test decrement(g, 45, 46) ≈ 1 - 0.9285202788707242
        @test hazard(g, 45) ≈ 0.06944687609574696

        @testset "AMLCR" begin
            m = Makeham(a=2.7e-6, b=log(1.124), c=0.00022)

            @test MortalityTables.μ(m, 20) == 0.00022 + 2.7e-6 * 1.124^20
            @test m[20] == MortalityTables.μ(m, 20)
            @test m(20) == MortalityTables.μ(m, 20)

            # vs manually calculated (via QuadGK) integrals
            @test decrement(m, 20, 25) ≈ 0.0012891622754368504
            @test survival(m, 20, 25) ≈ 1 - 0.0012891622754368504
            @test decrement(m, 25) ≈ 0.005888764668801838
            @test survival(m, 25) ≈ 1 - 0.005888764668801838

            # these values come from the 'Standard Select and Ultimate Survival Model'
            # from Actuarial Mathematics for Life Contingent Risks, 2nd end

            ℓ = 100_000
            ℓs = [survival(m, 20, age) for age in 21:100] .* ℓ

            ℓ_age(x) = ℓs[x - 20]
            @test isapprox(ℓ_age(21),  99975.04, atol=0.01)
            @test isapprox(ℓ_age(31),  99695.83, atol=0.01)
            @test isapprox(ℓ_age(82),  70507.19, atol=0.01)
            @test isapprox(ℓ_age(100),  6248.17, atol=0.01)
        end

    end

    @testset "Gompertz and Makeham equality" begin

        # Gompertz is Makeham's where c = 0
        m = Makeham(a=2.7e-6, b=1.124, c=0.0)
        g = Gompertz(a=2.7e-6, b=1.124)

        for age ∈ 20:100
            @test survival(m, age) == survival(g, age)
            @test survival(m, age, 1) == survival(g, age, 1)
        end
    end

    @testset "MortalityLaws R package" begin
        model_tests = [(rmodel = "gompertz", juliamodel = MortalityTables.Gompertz()),
                        # (rmodel="gompertz0",juliamodel=MortalityTables.nothing),
                        (rmodel = "invgompertz", juliamodel = MortalityTables.InverseGompertz()),
                        (rmodel = "makeham", juliamodel = MortalityTables.Makeham()),
                        # (rmodel="makeham0",juliamodel=MortalityTables.nothing),
                        (rmodel = "opperman", juliamodel = MortalityTables.Opperman()),
                        (rmodel = "thiele", juliamodel = MortalityTables.Thiele()),
                        (rmodel = "wittstein", juliamodel = MortalityTables.Wittstein()),
                        (rmodel = "perks", juliamodel = MortalityTables.Perks()),
                        (rmodel = "weibull", juliamodel = MortalityTables.Weibull()),
                        (rmodel = "invweibull", juliamodel = MortalityTables.InverseWeibull()),
                        (rmodel = "vandermaen", juliamodel = MortalityTables.VanderMaen()),
                        (rmodel = "vandermaen2", juliamodel = MortalityTables.VanderMaen2()),
                        (rmodel = "strehler_mildvan", juliamodel = MortalityTables.StrehlerMildvan()),
                        (rmodel = "quadratic", juliamodel = MortalityTables.Quadratic()),
                        (rmodel = "beard", juliamodel = MortalityTables.Beard()),
                        (rmodel = "beard_makeham", juliamodel = MortalityTables.MakehamBeard()),
                        (rmodel = "ggompertz", juliamodel = MortalityTables.GammaGompertz()),
                        (rmodel = "siler", juliamodel = MortalityTables.Siler()),
                        (rmodel = "HP", juliamodel = MortalityTables.HeligmanPollard()),
                        (rmodel = "HP2", juliamodel = MortalityTables.HeligmanPollard2()),
                        (rmodel = "HP3", juliamodel = MortalityTables.HeligmanPollard3()),
                        (rmodel = "HP4", juliamodel = MortalityTables.HeligmanPollard4()),
                        (rmodel = "rogersplanck", juliamodel = MortalityTables.RogersPlanck()),
                        (rmodel = "martinelle", juliamodel = MortalityTables.Martinelle()),
                        (rmodel = "kostaki", juliamodel = MortalityTables.Kostaki()),
                        (rmodel = "kannisto", juliamodel = MortalityTables.Kannisto()),
                        (rmodel = "kannisto_makeham", juliamodel = MortalityTables.KannistoMakeham())
                        # the next two requre adding an autodiff dependency:
                        # (rmodel="carriere1",juliamodel=MortalityTables.Carriere()),
                        # (rmodel="carriere2",juliamodel=MortalityTables.Carriere2()),
                    ]

        # load test targets from data
        dir = joinpath(pwd(), "data", "parametric")

        rmodels = Dict()
        for (root, dirs, files) in walkdir(dir)
            for file in files
                if file[end - 3:end] == "json"
                    json = JSON.parse(MortalityTables.open_and_read(joinpath(root, file)))
                    rmodels[json["modelname"][1]] = json
                end
            end
        end

        # compare values
        ages = 20:100
        @testset "Model: $(model.rmodel)" for model in model_tests
            rmodel = rmodels[model.rmodel]
            # test hazard
            if "hx" in keys(rmodel)
                for (i, age) in enumerate(ages)
                    @test hazard(model.juliamodel, age) ≈ rmodel["hx"][i]
                end
            end

            # test cumulative hazard

            if "Hx" in keys(rmodel)
                for (i, age) in enumerate(ages)
                    @test cumhazard(model.juliamodel, age) ≈ rmodel["Hx"][i]
                end
            end

            # test Survival
            if "Sx" in keys(rmodel)
                for (i, age) in enumerate(ages)
                    @test survival(model.juliamodel, age) ≈ rmodel["Sx"][i]
                end
            end

            # test other characteristics

            @test survival(model.juliamodel, 20, 20) == 1.0
            if model.rmodel in ["quadratic","perks","vandermaen","vandermaen2"]
                # the default params create a crazy hazard function
                @test_broken survival(model.juliamodel, 50, 51) < 1.0
            else
                @test survival(model.juliamodel, 50, 51) < 1.0
            end

            @test model.juliamodel[20] >= 0

        end

    end

end


text:
int(x) = (Int ∘ round)(x)

"""
  loop_last(v::Vector)

  Returns an iterable yielding tuples (is_last, value).
"""
function loop_last(v::Vector)
    is_last = [i == length(v) for i in 1:length(v)]
    return zip(is_last, v)
end


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG

const R123Array1x{T<:Union{UInt128}}        = NTuple{1, T}
const R123Array2x{T<:Union{UInt32, UInt64}} = NTuple{2, T}
const R123Array4x{T<:Union{UInt32, UInt64}} = NTuple{4, T}

"The base abstract type for RNGs in [Random123 Family](@ref)."
abstract type AbstractR123{T<:Union{UInt32, UInt64, UInt128}} <: AbstractRNG{T} end

"RNG that generates one number at a time."
abstract type R123Generator1x{T} <: AbstractR123{T} end
"RNG that generates two numbers at a time."
abstract type R123Generator2x{T} <: AbstractR123{T} end
"RNG that generates four numbers at a time."
abstract type R123Generator4x{T} <: AbstractR123{T} end

_value(r::AbstractR123{T}, i = 1, ::Type{T2} = T) where {T, T2} = unsafe_load(Ptr{T2}(pointer_from_objref(r)), i)

"Set the counter of a Random123 RNG."
@inline function set_counter!(r::R123Generator1x{T}, ctr::Integer) where T <: UInt128
    r.ctr = ctr % T
    random123_r(r)
    r
end

@inline function set_counter!(
    r::R123Generator2x{T},
    (ctr1, ctr2)::NTuple{2, Integer}
) where T <: Union{UInt32, UInt64}
    r.p = 0
    r.ctr1 = ctr1 % T
    r.ctr2 = ctr2 % T
    random123_r(r)
    r
end

@inline function set_counter!(
    r::R123Generator4x{T},
    (ctr1, ctr2, ctr3, ctr4)::NTuple{4, Integer}
) where T <: Union{UInt32, UInt64}
    r.p = 0
    r.ctr1 = ctr1 % T
    r.ctr2 = ctr2 % T
    r.ctr3 = ctr3 % T
    r.ctr4 = ctr4 % T
    random123_r(r)
    r
end
@inline set_counter!(r::R123Generator2x, ctr::Integer) = set_counter!(r, (ctr, 0))
@inline set_counter!(r::R123Generator4x, ctr::Integer) = set_counter!(r, (ctr, 0, 0, 0))

@inline inc_counter!(r::R123Generator1x{T}) where T = (r.ctr += one(T); r)
@inline function inc_counter!(r::R123Generator2x{T}) where T
    r.ctr1 += one(T)
    r.ctr1 == zero(T) && (r.ctr2 += one(T))
    r
end
@inline function inc_counter!(r::R123Generator4x{T}) where T
    r.ctr1 += one(T)
    if r.ctr1 == zero(T)
        r.ctr2 += one(T)
        if r.ctr2 == zero(T)
            r.ctr3 += one(T)
            if r.ctr3 == zero(T)
                r.ctr4 += one(T)
            end
        end
    end
    r
end

@inline function rand(r::R123Generator1x{T}, ::Type{T}) where T <: UInt128
    r.ctr += one(T)
    random123_r(r)
    _value(r)
end

@inline function rand(r::R123Generator2x{T}, ::Type{T}) where T <: Union{UInt32, UInt64}
    to = one(T)
    if r.p == 1
        inc_counter!(r)
        random123_r(r)
        r.p = 0
        return r.x2
    end
    r.p = 1
    r.x1
end

@inline function rand(r::R123Generator4x{T}, ::Type{T}) where T <: Union{UInt32, UInt64}
    if r.p == 4
        inc_counter!(r)
        random123_r(r)
        r.p = 0
    end
    r.p += 1
    _value(r, r.p)
end

@inline function rand(r::R123Generator1x{T}, ::Type{R123Array1x{T}}) where T <: UInt128
    inc_counter!(r)
    random123_r(r)
end

@inline function rand(r::R123Generator2x{T}, ::Type{R123Array2x{T}}) where T <: Union{UInt32, UInt64}
    if r.p > 0
        inc_counter!(r)
    end
    ret = random123_r(r) # which returns a Tuple{T, T}
    r.p = 1
    ret
end

@inline function rand(r::R123Generator4x{T}, ::Type{R123Array4x{T}}) where T <: Union{UInt32, UInt64}
    if r.p > 0
        inc_counter!(r)
    end
    ret = random123_r(r)
    r.p = 4
    ret
end

for (T, DT) in ((UInt32, UInt64), (UInt64, UInt128))
    @eval @inline function rand(r::R123Generator2x{$T}, ::Type{$DT})
        if r.p == 1
            inc_counter!(r)
            random123_r(r)
        end
        r.p = 1
        _value(r, 1, $DT)
    end
end

@inline function rand(r::R123Generator4x{UInt32}, ::Type{UInt128})
    if r.p > 0
        inc_counter!(r)
        random123_r(r)
    end
    r.p = 4
    _value(r, 1, UInt128)
end

"Do one iteration and
return the result tuple of a Random123 RNG object."
random123_r


text:
using AtomicLevels: CSF
using HalfIntegers

angularmomentum(o::RelativisticOrbital) = AtomicLevels.κ2j(o.κ)
angularmomentum(csf::CSF{O,<:HalfInteger}) where O = last(csf.terms)

# Relativistic CSLs and parsing of GRASP CSL files

function parse_rcsf(filename)
    open(filename, "r") do io
        # First line should be "Core subshells:"
        line = readline(io); @assert strip(line) == "Core subshells:"
        line_cores = strip(readline(io), ['\n', '\r'])
        line = readline(io); @assert strip(line) == "Peel subshells:"
        line_peels = readline(io)
        line = readline(io); @assert strip(line) == "CSF(s):"

        core_orbitals = parse_cores(line_cores)
        core_couplings = zeros(HalfInt, length(core_orbitals))
        core_occupations = map(degeneracy, core_orbitals)

        blockid, csfid = 1, 1
        csfs = CSF{RelativisticOrbital{Int},HalfInt,Seniority}[]
        while ! eof(io)
            line1 = readline(io)
            if startswith(line1, " *")
                blockid += 1
                csfid = 1
                line1 = readline(io)
            end
            line1 = strip(line1, ['\n', '\r'])
            line2 = strip(readline(io), ['\n', '\r'])
            line3 = strip(readline(io), ['\n', '\r'])

            orbitals, noccupations, orbcouplings, csfcouplings = parse_csflines(line1, line2, line3)

            @assert !any(isequal(0), noccupations) # we should never get 0-electron orbitals

            # Fix orbcouplings that are not explicitly written in the CSL file (represented
            # with nothings in the array). It is assumed that this is the case only if the
            # orbital is fully occupied.
            #
            # NOTE: we could, though, also omit values if there is only 1 electron or
            # maxelectrons(orb) - 1 electrons, this which case the angular momentum can
            # only have only one value too.
            for i = 1:length(orbitals)
                orbital, nelec = orbitals[i], noccupations[i]
                if orbcouplings[i] === nothing
                    nelec == degeneracy(orbital) || error("""
                    Unable to fix missing orbital coupling.
                      Orbital $i, orbital=$(repr(orbital)), nelec=$nelec
                    1: $(line1)
                    2: $(line2)
                    3: $(line3)
                    """)
                    orbcouplings[i] = 0
                elseif (nelec == 1 || nelec == degeneracy(orbital) - 1) && orbcouplings[i] != angularmomentum(orbital)
                    @warn "Bad orbital coupling" orbcouplings[i] nelec angularmomentum(orbital)
                elseif orbcouplings[i] > nelec * angularmomentum(orbital)
                    # If the orbital coupling is larger than (# particles) * (l of orbital),
                    # then that has to be an error.
                    @warn "Bad orbital coupling" orbcouplings[i] nelec angularmomentum(orbital)
                end
            end

            # Fix csfcouplings which are not explicitly written to the CSL file. This appears
            # to be the case for "trivial" zero couplings, if the previous CSF layer and
            # current orbital are both zeros.
            for i = 1:length(orbitals)
                oj = orbcouplings[i]
                cj = (i > 1) ? csfcouplings[i-1] : zero(HalfInt)
                Δupper, Δlower = oj+cj, abs(oj-cj)
                if csfcouplings[i] === nothing
                    Δupper == Δlower || error("""
                    Unable to fix missing CSF coupling.
                      Orbital $i, orbital=$(repr(orbitals[i])), oj=$(repr(oj)), cj=$(repr(cj))
                    1: $(line1)
                    2: $(line2)
                    3: $(line3)
                    """)
                    csfcouplings[i] = Δlower
                elseif !(Δlower <= csfcouplings[i] <= Δupper)
                    @warn "Invalid csfcoupling value?" csfcouplings[i] Δupper Δlower
                end
            end

            config = Configuration(vcat(core_orbitals, orbitals), vcat(core_occupations, noccupations), sorted=true)
            # I have no idea how to generate correct seniority numbers, even in the case where
            # it does not matter (no degeneracy). So instead, I'll try to match the couplings
            # to the intermediate terms generated by AtomicLevels.
            subshell_terms = let its = intermediate_terms.(vcat(core_orbitals, orbitals), vcat(core_occupations, noccupations))
                map(its, vcat(core_couplings, Vector{HalfInt}(orbcouplings))) do its, grasp_term
                    idxs = findall(it -> it.term == grasp_term, its)
                    @assert length(idxs) == 1 # make sure there is only one intermediate term with the coupling we're interested in
                    its[first(idxs)]
                end
            end
            terms = map(x -> convert(Rational{Int}, x),
                vcat(core_couplings, Vector{HalfInt}(csfcouplings)))
            csf = CSF(config, subshell_terms, terms)
            push!(csfs, csf)
        end
        return csfs
    end
end

function parse_csflines(line1, line2, line3)
    # Assuming that the CSF line consists of NNNLL(EE) blocks, each 9 characters long.
    @assert length(line1) % 9 == 0

    orbs = RelativisticOrbital{Int}[]
    orbs_nelectrons = Int[]
    orbs_orbcouplings = Union{HalfInt,Nothing}[]
    orbs_csfcouplings = Union{HalfInt,Nothing}[]
    norbitals = div(length(line1), 9) # number of orbitals in this group
    for i = 1:norbitals
        orb = line1[9*(i-1)+1:9*i]
        @assert orb[6]=='(' && orb[9]==')'
        orbital = parse(RelativisticOrbital, strip(orb[1:5]))
        # n = parse(Int, orb[1:3])
        # kappa = parse_j(orb[4:5])
        nelec = parse(Int, orb[7:8])

        # pick out coupled angular momentum from the second line:
        angmom =  if length(line2) < 9*i
            nothing
        else
            parse_angmom_string(line2[9*(i-1)+1:9*i])
        end

        # Pick the J-coupling from between the orbitals (line3).
        # The items in that line are shifted by three characters to the right for
        # some reason.. except for the last one, which defines the J^P values of
        # the CSF. That one is shifted by 1 character.. and then one after that
        # is the parity +/- symbol.
        c2J_idx_first, c2J_idx_last = if i < norbitals
            # So, for the non-last ones we assume a 9 character block that has been
            # shifted by 3 characters to the right.
            # +1 to go from 0-based to 1-based
            9*(i-1)+3+1, 9*i+3
        else # i == norbitals -- the last non-regular coupling
            9*(i-1)+3+1, 9*i+1
        end
        c2J_string = line3[c2J_idx_first:c2J_idx_last]
        coupled_angmom = try
            parse_angmom_string(c2J_string)
        catch e
            error("""
            Error parsing 2J value on line 3 (i=$i)
              range $(c2J_idx_first):$(c2J_idx_last) -> '$(c2J_string)'
            1: $(line1)
            2: $(line2)
            3: $(line3)
            $(' '^(c2J_idx_first+2))^$('-'^(c2J_idx_last-c2J_idx_first-1))^
            """)
        end
        push!(orbs_csfcouplings, coupled_angmom)

        push!(orbs, orbital)
        push!(orbs_nelectrons, nelec)
        push!(orbs_orbcouplings, angmom)
    end
    @assert length(orbs_csfcouplings) == norbitals

    return orbs, orbs_nelectrons, orbs_orbcouplings, orbs_csfcouplings
end

function parse_angmom_string(s)
    length(strip(s)) == 0 && return nothing
    parse(HalfInt, s)
end

function parse_cores(line)
    orbstrings = split(line)
    orbs = RelativisticOrbital{Int}[]
    for os in orbstrings
        push!(orbs, parse(RelativisticOrbital, strip(os)))
    end
    orbs
end


text:
# overload model function

mutable struct DiscreteVariable{T} <: AbstractVariable
    wdth::T         # gap width
    vpop::Vector{T} # current values of the variable over population
    bstv::T         # best value so far
    modl            # probability model
    lini            # initial lower bound
end

function DiscreteVariable(r)
    T = eltype(r)
    l = first(r)
    u = last(r)
    wdth = step(r)
    vpop = Vector{T}()
    bstv = initv(T)
    modl = Uniform(0, (u - l) / wdth)
    DiscreteVariable(wdth, vpop, bstv, modl, l)
end

fit_model(v::DiscreteVariable, s) = fit_mle(Normal, v.vpop[s] ./ v.wdth)

function Base.rand(v::DiscreteVariable)
    if isa(v.modl, Uniform)
        s = (v.modl |> rand |> round) * v.wdth + v.lini
    else
        s = (v.modl |> rand |> round) * v.wdth
    end
    convert(eltype(v.vpop), s)
end

mutable struct DiscreteBoundedVariable{T} <: AbstractVariable
    lbnd::T         # lower bound
    ubnd::T         # upper bound
    wdth::T         # gap width
    vpop::Vector{T} # current values of the variable over population
    bstv::T         # best value so far
    modl            # probability model
end

function DiscreteBoundedVariable(r)
    T = eltype(r)
    lbnd = first(r)
    ubnd = last(r)
    wdth = step(r)
    vpop = Vector{T}()
    bstv = initv(T)
    modl = Uniform(0, (ubnd - lbnd) / wdth)
    DiscreteBoundedVariable(lbnd, ubnd, wdth, vpop, bstv, modl)
end

fit_model(v::DiscreteBoundedVariable, s) =
    truncated(fit_mle(Normal, (v.vpop[s] .- v.lbnd) ./ v.wdth), 0, (v.ubnd - v.lbnd) / v.wdth)

function Base.rand(v::DiscreteBoundedVariable)
    s = (v.modl |> rand |> round) * v.wdth + v.lbnd
    convert(eltype(v.vpop), s)
end