text:
@diffusion_process LorenzAux begin
    :dimensions
    process --> 3
    wiener --> 3

    :parameters
    _ --> (3, Float64)
    (t, T) --> Float64
    (u, v) --> ℝ{3}
    σ --> Float64

    :additional
    constdiff --> true
end

B(t, P::LorenzAux) = @SMatrix [
    -P.p1         P.p1        0.0;
     P.p2-P.v[3]  -1.0    -P.v[1];
     P.v[2]      P.v[1]    -P.p3
]

β(t, P::LorenzAux) = ℝ{3}(
    0.0,
    P.v[1]*P.v[3],
    -P.v[1]*P.v[3]
)

σ(t, P::LorenzAux) = SDiagonal(P.σ, P.σ, P.σ)


text:
# Generator functions for measurements.
function core_events(events = default_events(), cores = default_cpu_fore_mask())
    return (
        timestamp = SystemSnoop.Timestamp(),
        counters = CoreMonitorWrapper(;events = events, cores = cores)
    )
end

function uncore_events(counters)
    return (
        timestamp = SystemSnoop.Timestamp(),
        counters = Uncore{2,2,6}(counters)
    )
end

ext: jl

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Groups
using Base.Test

const n = 1000
srand(100)

# Test unitary irreps for O3
@test begin
    for i in 1:n
        α, β, γ = 2π*rand(), π*rand(), 2π*rand()
        g = O3(α, β, γ)

        for i in 0:4
            unitary_irrep(g, i)
        end
    end

    return true
end

# Test orthogonal irreps for O3
@test begin
    for i in 1:n
        α, β, γ = 2π*rand(), π*rand(), 2π*rand()
        g = O3(α, β, γ)

        for i in 0:4
            orthogonal_irrep(g, i)
        end
    end

    return true
end


text:
using CUDA

# Each block is assigned a row.
# Each thread checks  the main row
# with one of the element of the columns pointed by the row.
# The columns aren't read at once but in strided fashion.
function countTri_elem_per_th_limit(csr_rows, col_indx, out_sum, max_mem_nnz, max_nnz)
#function countTri(out_sum)
    bid::Int32 = blockIdx().x + (blockIdx().y-1)*gridDim().x + (blockIdx().z-1)*gridDim().x*gridDim().y
    rows::Int32 = size(csr_rows)[1]
    nnz::Int32 = size(col_indx)[1]
    tid::Int32 = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride::Int32 = blockDim().x*blockDim().y
    if bid > rows
        return nothing
    end

    sh_row = @cuDynamicSharedMem(Int32, max_mem_nnz)
    sh_cols = @cuDynamicSharedMem(Int32, max_nnz + CUDA.max(stride, max_nnz),
                                sizeof(Int32)*max_mem_nnz)
    sh_sum = @cuDynamicSharedMem(Int32, blockDim().y,
            (max_nnz + CUDA.max(stride, max_nnz) + max_mem_nnz)*sizeof(Int32))

    @inbounds row_start::Int32 =  csr_rows[bid]
    @inbounds row_end::Int32 = (bid == rows) ?  nnz+1 : csr_rows[bid + 1]

    len::Int32 = row_end-row_start
    for s::Int32 = 0:stride:(len-1)
        if tid + s <= len
            @inbounds sh_row[tid+s] = col_indx[row_start + tid - 1 + s]
        end
    end
    sync_threads()
    #
    #Read all the other rows
    total_elems::Int32 = 0
    sum::Int32 = 0
    for i::Int32 = 1:len
        @inbounds col::Int32 = sh_row[i]
        @inbounds row_start = csr_rows[col]
        @inbounds row_end = (col == rows) ?  nnz+1 : csr_rows[col + 1]
        templen::Int32 = row_end-row_start
        sync_threads()
        for s::Int32 = 0:stride:(templen-1)
            if s + tid <= templen
                @inbounds sh_cols[total_elems + s + tid] = col_indx[row_start + s + tid - 1]
            end
        end
        total_elems  += templen
        sync_threads()
        if total_elems >= stride || i == len
            for s::Int32 = 0:stride:(total_elems-1)
                if tid + s <= total_elems
                    for l = 1:len
                        # sum += sh_row[(l%len)+1] == sh_cols[s + tid]
                        @inbounds sum += sh_row[((tid+l-2)%len)+1] == sh_cols[s + tid]
                    end
                end
            end
            total_elems = 0
        end
    end

    w::Int32 = 16
    while w > 0
        sum += shfl_down_sync(0xffffffff, sum, w)
        w = floor(w/2)
    end
    if threadIdx().x == 1
        @inbounds sh_sum[threadIdx().y] = sum
    end
    sync_threads()
    if tid == 1
        for i::Int32 = 2:blockDim().y
            @inbounds sum += sh_sum[i]
        end
        @inbounds out_sum[bid] = sum
    end
    return nothing
end

function CountTriangles_ElemPerThLimit(csr_rows, col_indx, max_per_row; th_groups=2)
    dev_csr_rows = CuArray(csr_rows)
    dev_col_indx = CuArray(col_indx)
    dev_out_sum = CUDA.zeros(Int32, size(csr_rows)[1])
    thread_groups = th_groups
    max_mem_nnz = 0
    while max_mem_nnz < max_per_row
        max_mem_nnz += 32
    end
    shmem = sizeof(Int32)*(thread_groups + max_mem_nnz + max_per_row + max(thread_groups*32,max_per_row))
    blocks = Int64(floor((size(csr_rows)[1])/(512*512)+ 1))
    println("3 Total blocks: ", 512*512*blocks, " vs rows: ", size(csr_rows)[1])
    println("\t Threads: 32*", thread_groups)
    @cuda(blocks=(512,512,blocks), threads=(32, thread_groups),
        shmem=shmem,
        countTri_elem_per_th_limit(dev_csr_rows, dev_col_indx, dev_out_sum,
                                max_mem_nnz, max_per_row))
    out_sum = Array(dev_out_sum)
    # println(dev_out_sum)
    return reduce(+, out_sum, init=0)
end


text:
using PrototypesCriticisms
using BenchmarkTools
using Clustering
using Distances
using KernelFunctions
using Random

function naiveprototypes(X::AbstractMatrix{<:Real}, n::Int, k::Kernel=RBFKernel())
    protoids = []
    while length(protoids) < n
        fs = []
        for i in setdiff(1:size(X, 2), protoids)
            P2 = view(X, :, [protoids; i])
            P1 = view(X, :, protoids)
            if length(protoids) > 0
                f = sqmmd(X, P2, k) - sqmmd(X, P1, k)
            else
                f = sqmmd(X, P2, k)
            end
            push!(fs, (f, i))
        end
        push!(protoids, fs[argmin(fs)][2])
    end
    return protoids
end

function naivecriticisms(X::AbstractMatrix{<:Real}, protoids::AbstractVector{Int}, n::Int, k::Kernel=RBFKernel())
    critids = []
    while length(critids) < n
        absws = []
        for i in setdiff(1:size(X, 2), union(critids, protoids))
            w = abs(witness(view(X, :, i), X, view(X, :, protoids), k))
            push!(absws, (abs(w), i))
        end
        push!(critids, absws[argmax(absws)][2])
    end
    return critids
end

function main()
    Random.seed!(42)

    n = 200
    p = 3
    c = 3
    k = 3
    D = rand(20, n)

    # k-medoids
    clustering = kmedoids(pairwise(Euclidean(), D), k)
    println("Prototypes implementation (k-medoids)")
    @btime prototypes($clustering, $p)
    println("Criticisms implementation (k-medoids)")
    @btime criticisms($clustering, $c)

    # k-means
    clustering = kmeans(D, k)
    println("Prototypes implementation (k-means)")
    @btime prototypes($clustering, $p)
    println("Criticisms implementation (k-means)")
    @btime criticisms($clustering, $c)

    # fuzzy c-means
    clustering = fuzzy_cmeans(D, k, 2)
    println("Prototypes implementation (fuzzy c-means)")
    @btime prototypes($clustering, $p)
    println("Criticisms implementation (fuzzy c-means)")
    @btime criticisms($clustering, $c)

    # affinity propagation
    S = -pairwise(Euclidean(), D, dims=2)
    S = S - diagm(0 => diag(S)) + median(S) * I
    clustering = affinityprop(S)
    println("Prototypes implementation (affinity propagation)")
    @btime prototypes($clustering, $D, $p)
    println("Criticisms implementation (affinity propagation)")
    @btime criticisms($clustering, $D, $c)

    # MMD-critic
    kernel = with_lengthscale(RBFKernel(), sqrt(20))
    println("Naive prototypes implementation (MMD-critic)")
    @btime naiveprototypes($D, $p, $kernel)
    println("Prototypes implementation (MMD-critic)")
    @btime prototypes($D, $kernel, $p)
    protoids = rand(1:n, p)
    println("Naive criticisms implementation (MMD-critic)")
    @btime naivecriticisms($D, $protoids, $c, $kernel)
    println("Criticisms implementation (MMD-critic)")
    @btime criticisms($D, $kernel, $protoids, $c)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
"""
    get_colormap(colormap)

Get the colormap as a ColorScheme if it is a named color or ColorScheme
"""
function get_colormap(colormap)
    if typeof(colormap) <: ColorScheme
        return colormap
    elseif typeof(colormap) <: Symbol
        return colorschemes[colormap]
    else
        error("colormap must be a ColorScheme")
    end
end

"""
    get_color(colormap, var <: AbstractArray, range_var)
    get_color(colormap, var, range_var)

Map value(s) to colors from a colormap based on a range of values
"""
function get_color(colormap, var::T, range_var) where {T<:AbstractArray}
    x2 = (range_var[2] - range_var[1])
    # get the color based on a colormap and the normalized attribute value
    [get(colormap, (i - range_var[1]) / x2) for i in var]
end

function get_color(colormap, var, range_var)
    # get the color based on a colormap and the normalized attribute value
    get(colormap, (var - range_var[1]) / (range_var[2] - range_var[1]))
end


text:
module LOBPCG
    using LinearAlgebra

    function eigenvalues_LOPBCG!(n,md,x,H,λ,eps,x_temp,itemax)
        gram_real!(n,md,x,x_temp)
        r = zero(x)
        p = zero(x)
        z = zeros(Float64,n,3*md)
        ztemp = zero(z)

        #for ii=1:md
        mul!(x,H,x_temp)
        #end
        xhx = x_temp'*x
        λ,v = eigen(xhx)
        mul!(x,x_temp,v)

        for ii=1:md
            for k=1:n
                r[k,ii] = x[k,ii]*λ[ii]
            end
        end

        #for ii=1:md
        mul!(x_temp,H,x)
        #end
        r .= x_temp .- r

        for ite=1:itemax
            for ii=1:md
                for k=1:n
                    z[k,ii] = x[k,ii]
                    z[k,md+ii] = r[k,ii]
                    x[k,2*md+ii] = p[k,ii]
                end
            end
            ztemp .= z

            nz = ifelse(ite == 1,2md,3md)
            gram_real!(n,nz,x,x_temp)

        end
    end

    function gram_real!(m,n,mat_v,mat_v_out)
        mat_v_out .= mat_v
        viold = zeros(Float64,m)
        vi = zeros(Float64,m)
        for i=1:n
            for k=1:m
                viold[k] = mat_v_out[k,i]
            end
            for j=1:(i-1)
                nai = 0
                for k=1:m
                    nai += mat_v_out[k,j]*viold[m]
                end
                for k=1:m
                    vi[k] = viold[k] - nai*mat_v_out[k,j]
                end
                viold .= vi
            end
            nor = sqrt(dot(viold,viold))
            for k=1:m
                mat_v_out[k,j] = viold[k]/nor
            end
        end
        return
    end
end


text:
# Autogenerated wrapper script for cddlib_jll for armv6l-linux-gnueabihf
export libcdd, libcddgmp

using GMP_jll
JLLWrappers.@generate_wrapper_header("cddlib")
JLLWrappers.@declare_library_product(libcdd, "libcdd.so.0")
JLLWrappers.@declare_library_product(libcddgmp, "libcddgmp.so.0")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libcdd,
        "lib/libcdd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcddgmp,
        "lib/libcddgmp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using repeat
using Test

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


text:
function show(io::IO, mC::MapUnitGrp)
  println(io, "UnitGroup map of $(codomain(mC))")
end

function unit_group_disc_exp(x::GrpAbFinGenElem, U::UnitGrpCtx)
  K = nf(order(U))
  y = FacElem([K(U.torsion_units_gen)], [x.coeff[1,1]])
  for i=1:length(U.units)
    mul!(y, y, U.units[i]^x.coeff[1,i+1])
  end
  return y
end

function unit_group_disc_log(x::FacElem{nf_elem, AnticNumberField} , U::UnitGrpCtx, G::GrpAbFinGen) #snf

  if length(U.units) == 0
    r = [-1]
  else
    r = _add_dependent_unit!(U, x, Val{true})
  end
  @assert r[end] == -1

  y = deepcopy(x)
  for i=1:length(r)-1
    y *= U.units[i]^-r[i]
  end

  K = nf(order(U))

  p = next_prime(2^30)
  while (p-1) % U.torsion_units_order != 0
    p = next_prime(p)
  end
  P = prime_decomposition(order(U), p)[1][1]
  F, mF = ResidueField(order(U), P)
  mK = extend(mF, K)

  yp = F(1)
  for (k,v) = y.fac
    yp *= mK(k)^v
  end

  zp = mF(U.torsion_units_gen)
  res = fmpz[]
  for i=0:U.torsion_units_order-1
    if zp^i == yp
      push!(res, i)
    end
  end
  @assert length(res) == 1

  for i = 1:length(r)-1
    push!(res, r[i])
  end
  return G(res)
end

function unit_group_fac_elem(c::ClassGrpCtx; redo::Bool = false)
  u = unit_group_ctx(c, redo = redo)
  if isdefined(u, :unit_map)
    mU = u.unit_map
    U = domain(mU)
    return U, mU
  end
  return unit_group_fac_elem(u)
end

function unit_group_fac_elem(u::UnitGrpCtx)
  O = order(u)

  zo = u.torsion_units_order
  if zo == -1
    u.torsion_units_gen, u.torsion_units_order = torsion_units_gen_order(O)
    zo = u.torsion_units_order
  end
  r = unit_group_rank(O)
  d = fmpz[zo]
  for i=1:r
    push!(d, fmpz(0))
  end
  U = abelian_group(d)

  r = MapUnitGrp{typeof(FacElemMon(nf(O)))}()

  r.header = MapHeader(U, FacElemMon(nf(O)),
    x->unit_group_disc_exp(x, u),
    x->unit_group_disc_log(x, u, U))

  u.unit_map = r

  return U, r
end


text:
function imap(I::sideal, R::PolyRing)
    # @info imap [g for g in I]
    basis = [Singular.imap(g, R) for g in I]
    if isempty(basis)
        return Singular.Ideal(R)
    end
    Singular.Ideal(R, basis)
end

function fetch(I::sideal, R::PolyRing)
    basis = [Singular.fetch(g, R) for g in I]
    if isempty(basis)
        return Singular.Ideal(R)
    end
    Singular.Ideal(R, basis)
end

Base.length(ideal::Singular.sideal) = ngens(ideal)
Base.iterate(ideal::Singular.sideal) = ngens(ideal) > 0 ? (ideal[1], 1) : nothing
Base.iterate(ideal::Singular.sideal, state) = state < ngens(ideal) ? (ideal[state+1], state+1) : nothing
Base.eltype(::Singular.sideal) = Singular.spoly

import Base.==

function ==(I::Singular.sideal, J::Singular.sideal)
    # TODO: is this the right kind of equality?
    for (a, b) in zip(I, J)
        if a != b
            return false
        end
    end
    return true
end


text:
@doc """
A RollCall object represents a complete roll-call data set along with metadata
about the encoding used by the dataset. This encoding is stored in the form
of a RollCallEncoding object.
""" ->
immutable RollCall{T1 <: Tuple, T2 <: Tuple, T3 <: Tuple, T4 <: Tuple}
    votes::Matrix{Int}
    encoding::RollCallEncoding{T1, T2, T3, T4}
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
export
    AugmentedLagrangianSolver,
    AugmentedLagrangianSolverOptions,
    get_constraints

@with_kw mutable struct ALStats{T}
    iterations::Int = 0
    iterations_total::Int = 0
    iterations_inner::Vector{Int} = zeros(Int,0)
    cost::Vector{T} = zeros(0)
    c_max::Vector{T} = zeros(0)
    penalty_max::Vector{T} = zeros(0)
end

function reset!(stats::ALStats, L=0)
    stats.iterations = 0
    stats.iterations_total = 0
    stats.iterations_inner = zeros(Int,L)
    stats.cost = zeros(L)*NaN
    stats.c_max = zeros(L)*NaN
    stats.penalty_max = zeros(L)*NaN
end

"""$(TYPEDEF)
Solver options for the augmented Lagrangian solver.
$(FIELDS)
"""
@with_kw mutable struct AugmentedLagrangianSolverOptions{T} <: AbstractSolverOptions{T}
    "Print summary at each iteration."
    verbose::Bool=false

    "unconstrained solver options."
    opts_uncon::UnconstrainedSolverOptions{T} = UnconstrainedSolverOptions{Float64}()

    "dJ < ϵ, cost convergence criteria for unconstrained solve or to enter outerloop for constrained solve."
    cost_tolerance::T = 1.0e-4

    "dJ < ϵ_int, intermediate cost convergence criteria to enter outerloop of constrained solve."
    cost_tolerance_intermediate::T = 1.0e-3

    "gradient_norm < ϵ, gradient norm convergence criteria."
    gradient_norm_tolerance::T = 1.0e-5

    "gradient_norm_int < ϵ, gradient norm intermediate convergence criteria."
    gradient_norm_tolerance_intermediate::T = 1.0e-5

    "max(constraint) < ϵ, constraint convergence criteria."
    constraint_tolerance::T = 1.0e-3

    "max(constraint) < ϵ_int, intermediate constraint convergence criteria."
    constraint_tolerance_intermediate::T = 1.0e-3

    "maximum outerloop updates."
    iterations::Int = 30

    "global maximum Lagrange multiplier. If NaN, use value from constraint"
    dual_max::T = NaN

    "global maximum penalty term. If NaN, use value from constraint"
    penalty_max::T = NaN

    "global initial penalty term. If NaN, use value from constraint"
    penalty_initial::T = NaN

    "global penalty update multiplier; penalty_scaling > 1. If NaN, use value from constraint"
    penalty_scaling::T = NaN

    "penalty update multiplier when μ should not be update, typically 1.0 (or 1.0 + ϵ)."
    penalty_scaling_no::T = 1.0

    "ratio of current constraint to previous constraint violation; 0 < constraint_decrease_ratio < 1."
    constraint_decrease_ratio::T = 0.25

    "type of outer loop update (default, feedback)."
    outer_loop_update_type::Symbol = :default

    "numerical tolerance for constraint violation."
    active_constraint_tolerance::T = 0.0

    "terminal solve when maximum penalty is reached."
    kickout_max_penalty::Bool = false

    reset_duals::Bool = true

    reset_penalties::Bool = true

    log_level::Base.CoreLogging.LogLevel = OuterLoop
end

function AugmentedLagrangianSolverOptions(opts::SolverOptions)
    opts_uncon = UnconstrainedSolverOptions(opts)
    AugmentedLagrangianSolverOptions(
        opts_uncon=opts_uncon,
        cost_tolerance=opts.cost_tolerance,
        cost_tolerance_intermediate=opts.cost_tolerance_intermediate,
        iterations=opts.iterations,
        penalty_initial=opts.penalty_initial,
        penalty_scaling=opts.penalty_scaling,
        verbose=opts.verbose
    )
end

function reset!(conSet::ConstraintSet{T}, opts::AugmentedLagrangianSolverOptions{T}) where T
    if !isnan(opts.dual_max)
        for con in conSet.constraints
            params = get_params(con)::ConstraintParams{T}
            params.λ_max = opts.dual_max
        end
    end
    if !isnan(opts.penalty_max)
        for con in conSet.constraints
            params = get_params(con)::ConstraintParams{T}
            params.μ_max = opts.penalty_max
        end
    end
    if !isnan(opts.penalty_initial)
        for con in conSet.constraints
            params = get_params(con)::ConstraintParams{T}
            params.μ0 = opts.penalty_initial
        end
    end
    if !isnan(opts.penalty_scaling)
        for con in conSet.constraints
            params = get_params(con)::ConstraintParams{T}
            params.ϕ = opts.penalty_scaling
        end
    end
    if opts.reset_duals
        reset_duals!(conSet)
    end
    if opts.reset_penalties
        reset_penalties!(conSet)
    end
end

function set_verbosity!(opts::AugmentedLagrangianSolverOptions)
    log_level = opts.log_level
    if opts.verbose
        set_logger()
        Logging.disable_logging(LogLevel(log_level.level-1))
        logger = global_logger()
        if opts.opts_uncon.verbose
            freq = 1
        else
            freq = 5
        end
        logger.leveldata[log_level].freq = freq
    else
        Logging.disable_logging(log_level)
    end
end

@doc raw""" ```julia
struct AugmentedLagrangianSolver <: TrajectoryOptimization.AbstractSolver{T}
```
Augmented Lagrangian (AL) is a standard tool for constrained optimization. For a trajectory optimization problem of the form:
```math
\begin{aligned}
  \min_{x_{0:N},u_{0:N-1}} \quad & \ell_f(x_N) + \sum_{k=0}^{N-1} \ell_k(x_k, u_k, dt) \\
  \textrm{s.t.}            \quad & x_{k+1} = f(x_k, u_k), \\
                                 & g_k(x_k,u_k) \leq 0, \\
                                 & h_k(x_k,u_k) = 0.
\end{aligned}
```
AL methods form the following augmented Lagrangian function:
```math
\begin{aligned}
    \ell_f(x_N) + &λ_N^T c_N(x_N) + c_N(x_N)^T I_{\mu_N} c_N(x_N) \\
           & + \sum_{k=0}^{N-1} \ell_k(x_k,u_k,dt) + λ_k^T c_k(x_k,u_k) + c_k(x_k,u_k)^T I_{\mu_k} c_k(x_k,u_k)
\end{aligned}
```
This function is then minimized with respect to the primal variables using any unconstrained minimization solver (e.g. iLQR).
    After a local minima is found, the AL method updates the Lagrange multipliers λ and the penalty terms μ and repeats the unconstrained minimization.
    AL methods have superlinear convergence as long as the penalty term μ is updated each iteration.
"""
struct AugmentedLagrangianSolver{T,S<:AbstractSolver} <: ConstrainedSolver{T}
    opts::AugmentedLagrangianSolverOptions{T}
    stats::ALStats{T}
    stats_uncon::Vector{STATS} where STATS
    solver_uncon::S
end

AbstractSolver(prob::Problem{Q,T},
    opts::AugmentedLagrangianSolverOptions{T}=AugmentedLagrangianSolverOptions{T}()) where {Q,T} =
    AugmentedLagrangianSolver(prob,opts)

"""$(TYPEDSIGNATURES)
Form an augmented Lagrangian cost function from a Problem and AugmentedLagrangianSolver.
    Does not allocate new memory for the internal arrays, but points to the arrays in the solver.
"""
function AugmentedLagrangianSolver(prob::Problem{Q,T}, opts=SolverOptions{T}();
        solver_uncon=iLQRSolver) where {Q,T}
    # Init solver statistics
    stats = ALStats()
    stats_uncon = Vector{iLQRStats{T}}()

    # Convert problem to AL problem
    opts_al = AugmentedLagrangianSolverOptions(opts)
    alobj = ALObjective(prob.obj, prob.constraints)
    rollout!(prob)
    prob_al = Problem(prob.model, alobj, ConstraintSet(size(prob)...),
        prob.x0, prob.xf, prob.Z, prob.N, prob.t0, prob.tf)

    solver_uncon = solver_uncon(prob_al, opts)

    solver = AugmentedLagrangianSolver(opts_al, stats, stats_uncon, solver_uncon)
    reset!(solver)
    return solver
end

function reset!(solver::AugmentedLagrangianSolver{T}) where T
    reset!(solver.stats, solver.opts.iterations)
    reset!(solver.solver_uncon)
    reset!(get_constraints(solver), solver.opts)
end

Base.size(solver::AugmentedLagrangianSolver) = size(solver.solver_uncon)
@inline cost(solver::AugmentedLagrangianSolver) = cost(solver.solver_uncon)
@inline get_trajectory(solver::AugmentedLagrangianSolver) = get_trajectory(solver.solver_uncon)
@inline get_objective(solver::AugmentedLagrangianSolver) = get_objective(solver.solver_uncon)
@inline get_model(solver::AugmentedLagrangianSolver) = get_model(solver.solver_uncon)
@inline get_initial_state(solver::AugmentedLagrangianSolver) = get_initial_state(solver.solver_uncon)
@inline iterations(solver::AugmentedLagrangianSolver) = solver.stats.iterations_total

function get_constraints(solver::AugmentedLagrangianSolver{T}) where T
    obj = get_objective(solver)::ALObjective{T}
    obj.constraints
end

struct ALObjective{T,O<:Objective} <: AbstractObjective
    obj::O
    constraints::ConstraintSet{T}
end

get_J(obj::ALObjective) = obj.obj.J
Base.length(obj::ALObjective) = length(obj.obj)

# TrajectoryOptimization.num_constraints(prob::Problem{Q,T,<:ALObjective}) where {T,Q} = prob.obj.constraints.p

function Base.copy(obj::ALObjective)
    ALObjective(obj.obj, ConstraintSet(copy(obj.constraints.constraints), length(obj.obj)))
end

function cost!(obj::ALObjective, Z::Traj)
    # Calculate unconstrained cost
    cost!(obj.obj, Z)

    # Calculate constrained cost
    evaluate!(obj.constraints, Z)
    update_active_set!(obj.constraints, Z, Val(0.0))
    for con in obj.constraints.constraints
        cost!(obj.obj.J, con, Z)
    end
end

function cost_expansion!(E, G, obj::ALObjective, model::AbstractModel, Z::Traj)
    # Update constraint jacobians
    jacobian!(obj.constraints, Z)

    ix, iu = Z[1]._x, Z[1]._u

    # Calculate expansion of original objective
    cost_expansion!(E, G, obj.obj, model, Z)

    # Add in expansion of constraints
    for con in obj.constraints.constraints
        cost_expansion(E, G, con, Z)
    end
end

function cost_expansion!(E::Vector{<:AbstractExpansion}, obj::ALObjective,Z::Traj)
    # Update constraint jacobians
    jacobian!(obj.constraints, Z)

    ix, iu = Z[1]._x, Z[1]._u

    # Calculate expansion of original objective
    cost_expansion!(E, obj.obj, Z)

    # Add in expansion of constraints
    for con in obj.constraints.constraints
        cost_expansion!(E, con, Z)
    end
end


text:
# code from https://github.com/JuliaGraphics/FreeType.jl/tree/master/gen

using Clang.Generators

# https://github.com/EsotericSoftware/spine-runtimes/releases

const SPINE_VERSION = v"3.8.95"
const spine_runtimes_dir = normpath(@__DIR__, "spine-runtimes-$SPINE_VERSION")
!isdir(spine_runtimes_dir) && throw(string("need ", spine_runtimes_dir))
const spine_include_dir = normpath(spine_runtimes_dir, "spine-c/spine-c/include")
const spine_include_spine_dir = normpath(spine_include_dir, "spine")
((root, dirs, files),) = walkdir(spine_include_spine_dir)
const spine_headers = normpath.(root, filter(x -> endswith(x, ".h"), files))

const headers = spine_headers
const args = vcat(get_default_args(), "-I$spine_include_dir")
const options = load_options(joinpath(@__DIR__, "generator.toml"))

ctx = create_context(headers, args, options)

build!(ctx)


text:
using Luxor, Test, Colors

function strokescale()
    Drawing(10, 10, :png, strokescale=true)

    # was the scale set correctly when the Drawing was created?
    @test setstrokescale() == true

    background("white")
    setcolor(0, 0, 0)

    # scale the canvas to use coordinates [0, 1]
    scale(10)

    # set our stroke width to the width of the canvas
    setline(1)

    # a line the width of the canvas drawn at 0 should cover 50% of the canvas
    line(Point(0, 0), Point(0, 1), :stroke)

    mat = image_as_matrix()
    @test RGBA{Float32}(0, 0, 0, 1) == convert(RGBA{Float32}, mat[5, 3])

    background("white")

    setstrokescale(false)
    @test setstrokescale() == false

    # set our stroke width to size 1, independent of scaling
    setline(1)

    # An unscaled line of size 1 will not cover 50% of the canvas
    line(Point(0, 0), Point(0, 1), :stroke)
    mat = image_as_matrix()

    # assert that the same pixel is now white
    @test RGBA{Float32}(1, 1, 1, 1) == convert(RGBA{Float32}, mat[5, 3])
    @test finish() == true
end

strokescale()

println("...finished strokescaletest")


text:
using Test
using EvolutionaryModelingTools
using EvolutionaryModelingTools: parse_func, parse_eq, parse_args_body, collectaxes
using LoopVectorization

@testset "parse function" begin
    fn, args = parse_func(:(@inline f(rng, ind, v) = v[ind] = rand(rng)), (:rng, :ind))
    @test fn == :f
    @test args == [:rng, :ind, :(args.v)]
end
@testset "@quickloop" begin
    A = rand(4)
    B = rand(4)
    C = rand(4, 4)
    D = rand(4, 4)
    # those functions is sensitive to the order of arguments, which is undetermined
    f1 = @quickloop R1[i, j] := A[i] * B[j]
    f2 = @quickloop R2[j, i] := A[i] * B[j]
    f3 = @quickloop ifelse(i == j, zero(eltype(A)), A[i, j]) # make diagonal zero
    f4 = @quickloop R3[i, j, k] := A[i, j] * B[i, k]
    @test (f1(A, B) == A .* B' == transpose(f2(A, B)) ||
           f1(A, B) == B .* A' == transpose(f2(A, B)))
    F = f3(C)
    equality_f3 = true
    for i in 1:4, j in 1:4
        equality_f3 &= ifelse(i == j, F[i, j] == 0, F[i, j] == C[i, j])
    end
    @test equality_f3
    E = f4(C, D)
    equality_f4 = true
    for k in 1:4, j in 1:4, i in 1:4
        equality_f4 &= ((E[i, j, k] == C[i, j] * D[i, k]) ||
            E[i, j, k] == D[i, j] * C[i, k])
    end
    @test equality_f4
end
@testset "collectaxes" begin
    let
        inds, inds_ax = collectaxes(:(C[i, j]), :(A[i] * B[j]))
        @test inds == [:i, :j]
        @test inds_ax[:i] == Dict(:A => [1])
        @test inds_ax[:j] == Dict(:B => [1])
    end
    let
        inds, inds_ax = collectaxes(:(E[]), :(A[i] * C[i, j] * B[j]))
        @test inds == [:i, :j]
        @test inds_ax[:i] == Dict(:A => [1], :C => [1])
        @test inds_ax[:j] == Dict(:B => [1], :C => [2])
    end
    let
        inds, inds_ax = collectaxes(:(E[j, i]), :(A[i] * C[i, j] * B[j]))
        @test inds == [:j, :i]
        @test inds_ax[:i] == Dict(:A => [1], :C => [1])
        @test inds_ax[:j] == Dict(:B => [1], :C => [2])
    end
    let
        inds, inds_ax = collectaxes(:(E[i, j, k]), :(c * A[i] * B[j, k]))
        @test inds == [:i, :j, :k]
        @test inds_ax[:i] == Dict(:A => [1])
        @test inds_ax[:j] == Dict(:B => [1])
        @test inds_ax[:k] == Dict(:B => [2])
    end
    let
        inds, inds_ax = collectaxes(:(E[i, j]), :(A[i] * B[j, j]))
        @test inds == [:i, :j]
        @test inds_ax[:i] == Dict(:A => [1])
        @test inds_ax[:j] == Dict(:B => [1, 2])
    end
    @test_throws ArgumentError collectaxes(:(C[i, k]),  :(A[i] * B[j]))
end
@testset "parse_eq" begin
    @test parse_eq(:(S + I[i] --> 2*I[i])) == (
        Dict{Any, Int}(:S => 1, :(I[i]) => 1),
        Dict{Any, Int}(:S => -1, :(I[i]) => 1)
    )
    @test parse_eq(:(S --> ∅)) == (
        Dict{Any, Int}(:S => 1),
        Dict{Any, Int}(:S => -1)
    )
    @test parse_eq(:(0 --> S)) == (
        Dict{Any, Int}(),
        Dict{Any, Int}(:S => 1)
    )
    @test parse_eq(:(S[i] --> S[i] + S[i])) == (
        Dict{Any, Int}(:(S[i]) => 1),
        Dict{Any, Int}(:(S[i]) => 1)
    )
    @test parse_eq(:(S + I --> 2(S + I))) == (
        Dict{Any, Int}(:S => 1, :I => 1),
        Dict{Any, Int}(:S => 1, :I => 1)
    )
    @test_throws ArgumentError parse_eq(:(S -> S))
    @test_throws ArgumentError("unexpected expression: f(S)") parse_eq(:(S --> f(S)))
end
@testset "parse_args_body" begin
    @test parse_args_body(quote
        global exp
        iv = f1(arg) # internal variable
        mpd = map(exp, iv) # global variable
        A[1] = 1 # setindex
        A[begin+1:end] # getindex
        A[(i = 1; begin+i:end-i)]
        i # i is a local variable defined in above expression
        A[begin+j:end-j] # j is a arguments needed
        arg.x
        f2(x) = x + 1 # local function
        (x, y) -> x + y # anonymous function
        @inbounds iv[ind] # macro call
        rand(rng)::Float64 # type annotation
    end)[1] == Set([:A, :arg, :j, :ind, :rng])
    @test parse_args_body(quote
        iv = f1(arg) # internal variable
        mpd = map(exp, iv) # global variable
        A[1] = 1 # setindex
        A[begin+1:end] # getindex
        arg.x
        f2(x) = x + 1 # local function
        (x, y) -> x + y # anonymous function
        @inbounds iv[ind] # macro call
        rand(rng)::Float64 # type annotation
    end)[1] == Set([:A, :arg, :exp, :ind, :rng]) # if exp is not mark as global, it will be collected
    @test parse_args_body(quote
        x = x.v # redefine argument
        f2(x) = x + y # local argument x and argument y
        (x, y) -> x + y + z # anonymous function use argument z
    end)[1] == Set([:x, :y, :z])
    @test parse_args_body(:((x::Int, y) -> x * y)) == (Set([:x, :y]), :(x * y))
    @test parse_args_body(:((x, y) -> @ein C[i, j] := x[i] * y[j]))[1] == Set([:x, :y])
end


text:
module SimpleTraits
using MacroTools
const curmod = nameof(@__MODULE__)

import InteractiveUtils

# This is basically just adding a few convenience functions & macros
# around Holy Traits.

export Trait, istrait, @traitdef, @traitimpl, @traitfn, Not, @check_fast_traitdispatch

# All traits are concrete subtypes of this trait.  SUPER is not used
# but present to be compatible with Traits.jl.
## @doc """
## `abstract Trait{SUPER}`

abstract type Trait end #{SUPER}

"""
All Traits are subtypes of abstract type Trait.

A concrete Trait will look like:
```julia
struct Tr1{X,Y} <: Trait end
```
where X and Y are the types involved in the trait.

(SUPER is not used here but in Traits.jl, thus retained for possible
future compatibility.)
"""
Trait

abstract type Not{T<:Trait} <: Trait end

"""
The set of all types not belonging to a trait is encoded by wrapping
it with Not{}, e.g.  Not{Tr1{X,Y}}
"""
Not

# Helper to strip an even number of Not{}s off: Not{Not{T}}->T
stripNot(::Type{T}) where {T<:Trait} = T
stripNot(::Type{Not{T}}) where {T<:Trait} = Not{T}
stripNot(::Type{Not{Not{T}}}) where {T<:Trait} = stripNot(T)

"""
A trait is defined as full-filled if this function is the identity
function for that trait.  Otherwise it returns the trait wrapped in
`Not`.

Example:
```
trait(IsBits{Int}) # returns IsBits{Int}
trait(IsBits{Array}) # returns Not{IsBits{Array}}
```

Usually this function is defined when using the `@traitimpl` macro.

However, instead of using `@traitimpl` one can define a method for
`trait` to implement a trait, see the README.
"""
trait(::Type{T}) where {T<:Trait} = Not{T}
trait(::Type{Not{T}}) where {T<:Trait} = trait(T)

## Under the hood, a trait is then implemented for specific types by
## defining:
#   trait(::Type{Tr1{Int,Float64}}) = Tr1{Int,Float64}
# or
#   trait(::Type{Tr1{I,F}}) where {I<:Integer,F<:FloatingPoint} = Tr1{I,F}
#
# Note due to invariance, this does probably not the right thing:
#   trait(::Type{Tr1{Integer,FloatingPoint}}) = Tr1{Integer, FloatingPoint}

"""
This function checks whether a trait is fulfilled by a specific
set of types.
```
istrait(Tr1{Int,Float64}) => return true or false
```
"""
istrait(::Any) = error("Argument is not a Trait.")
istrait(tr::Type{T}) where {T<:Trait} = trait(tr)==stripNot(tr) ? true : false # Problem, this can run into issue #265
                                                                        # thus is redefine when traits are defined
"""

Used to define a trait.  Traits, like types, are camel cased.  I
suggest to start them with a verb, e.g. `IsImmutable`, to distinguish
them from actual types, which are usually nouns.

Traits need to have one or more (type-)parameters to specify the type
to which the trait is applied.  For instance `IsImmutable{Int}`
signifies that `Int` is part of `IsImmutable` (although whether that
is true needs to be checked with the `istrait` function).  Most traits
will be one-parameter traits, however, several parameters are useful
when there is a "contract" between several types.

Examples:
```julia
@traitdef IsFast{X}
@traitdef IsSlow{X,Y}
```
"""
macro traitdef(tr)
    :(struct $(esc(tr)) <: Trait end)
end

"""
Used to add a type or type-tuple to a trait.  By default a type does
not belong to a trait.

Example:
```julia
@traitdef IsFast{X}
@traitimpl IsFast{Array{Int,1}}
```

Often a trait is dependent on some check-function returning true or
false.  This can be done with:
```julia
@traitimpl IsFast{T} <- isfast(T)
```
where `isfast` is that check-function.

Note that traits implemented with the former of above methods will
override an implementation with the latter method.  Thus it can be
used to define exceptions to the rule.
"""
macro traitimpl(tr)
    if tr.head==:curly || (tr.head==:call && tr.args[1]==:!)
        # makes
        # trait{X1<:Int,X2<:Float64}(::Type{Tr1{X1,X2}}) = Tr1{X1,X2}
        if tr.args[1]==:Not || isnegated(tr)
            tr = tr.args[2]
            negated = true
        else
            negated = false
        end
        typs = tr.args[2:end]
        trname = esc(tr.args[1])
        curly = Any[]
        paras = Any[]
        for (ty,v) in zip(typs, GenerateTypeVars{:upcase}())
            push!(curly, Expr(:(<:), esc(v), esc(ty)))  #:($v<:$ty)
            push!(paras, esc(v))
        end
        arg = :(::Type{$trname{$(paras...)}})
        fnhead = :($curmod.trait($arg) where {$(curly...)})
        isfnhead = :($curmod.istrait($arg) where {$(curly...)})
        if !negated
            return quote
                $fnhead = $trname{$(paras...)}
                nothing
            end
        else
            return quote
                $fnhead = Not{$trname{$(paras...)}}
                nothing
            end
        end
    elseif tr.head==:call
        @assert tr.args[1]==:<
        negated,Tr,P1,fn,P2 = @match tr begin
            Not{Tr_{P1__}} <- fn_(P2__) => (true, Tr, P1, fn, P2)
            Tr_{P1__} <- fn_(P2__) => (false, Tr, P1, fn, P2)
        end
        if negated
            fn = Expr(:call, GlobalRef(SimpleTraits, :!), fn)
        end
        return esc(quote
                   function SimpleTraits.trait(::Type{$Tr{$(P1...)}}) where {$(P1...)}
                       return $fn($(P2...)) ? $Tr{$(P1...)} : Not{$Tr{$(P1...)}}
                   end
                   nothing
                   end)
    else
        error("Cannot parse $tr")
    end
end

# Defining a function dispatching on the trait (or not)
# @traitfn f(x::X,y::Y) where {X,Y;  Tr1{X,Y}} = ...
# @traitfn f(x::X,y::Y) where {X,Y; !Tr1{X,Y}} = ... # which is just sugar for:
# @traitfn f(x::X,y::Y) where {X,Y; Not{Tr1{X,Y}}} = ...

dispatch_cache = Dict()  # to ensure that the trait-dispatch function is defined only once per pair
let
    global traitfn
    function traitfn(tfn, macro_module)
        # Need
        # f(x::X,Y::Y) where {X,Y} = f(trait(Tr1{X,Y}), x, y)
        # f(::False, x, y)= ...
        if tfn.head==:macrocall
            hasmac = true
            mac = tfn.args[1]
            tfn = tfn.args[3]
        else
            hasmac = false
        end

        # Dissect AST into:
        # fbody: body
        # fname: symbol of name
        # args0: vector of all arguments except wags (without any gensym'ed symbols but stripped of Traitor-traits)
        # args1: like args0 but with gensym'ed symbols where necessary
        # kwargs: all kwargs
        # typs0: vector of all function parameters (without the trait-ones, no gensym'ed)
        # typs: vector of all function parameters (without the trait-ones, with gensym'ed for Traitor)
        # trait: expression of the trait
        # trait0: expression of the trait without any gensym'ed symbols
        #
        # oldfn_syntax: set to true if the old parametric syntax is used
        #
        # (The variables without gensym'ed symbols are mostly used for the key of the dispatch cache)

        fhead = tfn.args[1]
        fbody = tfn.args[2]

        # TODO 1.0: remove
        out = @match fhead begin
            f_{paras__}(args0__;kwargs__)       => (f,paras,args0,kwargs,true)
            f_(args0__; kwargs__)               => (f,[],args0,kwargs,false)
            f_{paras__}(args0__)                => (f,paras,args0,[],true)
            f_(args0__)                         => (f,[],args0,[],false)
            f_(args0__; kwargs__) where paras__ => (f,paras,args0,kwargs,false)
            f_(args0__) where paras__           => (f,paras,args0,[],false)
        end
        if out==nothing
            error("Could not parse function-head: $fhead. Note that several `where` are not supported.")
        end
        fname, paras, args0, kwargs, oldfn_syntax = out
        haskwargs = length(kwargs)>0
        # extract parameters and traits from paras and/or args0

        if length(paras)==0
            maybe_traitor = true
        else
            if paras[1] isa Expr && paras[1].head==:parameters
                maybe_traitor = false
                length(paras)<2 && error("Cannot parse function parameters: $para")
                typs = paras[2:end]
                trait = paras[1].args[1]
            elseif paras[1] isa Expr && paras[1].head==:bracescat
                maybe_traitor = false
                length(paras)!=1 && error("Cannot parse function parameters: $para")
                typs = paras[1].args[1:1]
                trait = paras[1].args[2]
            else
                maybe_traitor = true
                # the processing happens below
            end
        end
        if !maybe_traitor
            trait0 = trait # without gensym'ed types, here identical
            typs0 = typs # without gensym'ed types, here identical
            args1 = insertdummy(args0)
        else # This is a Traitor.jl style function (or invalid).
            # Change paras & args0 into a Traits.jl function.
            # Find the traitor:
            typs0 = deepcopy(paras) # without gensym'ed types
            typs = paras
            out = nothing
            i = 0 # index of function argument with Traitor trait
            vararg = false
            for outer i in eachindex(args0)
                a = args0[i]
                vararg = a.head==:...
                if vararg
                    a = a.args[1]
                end
                out = @match a begin
                    ::::Tr_           => (nothing,nothing,Tr)
                    ::T_::Tr_         => (nothing,T,Tr)
                    x_Symbol::::Tr_   => (x,nothing,Tr)
                    x_Symbol::T_::Tr_ => (x,T,Tr)
                end
                out!=nothing && break
            end
            out==nothing && error("No trait found in function signature")
            arg,typ,trait0 = out
            if typ==nothing
                typ = gensym()
                push!(typs, typ)
            end
            if isnegated(trait0)
                trait = :(!($(trait0.args[2]){$typ}))
            else
                trait = :($trait0{$typ})
            end
            args1 = deepcopy(args0)
            if vararg
                args0[i] = arg==nothing ? nothing : :($arg...,).args[1]
                args1[i] = arg==nothing ? :(::$typ...,).args[1] : :($arg::$typ...,).args[1]
            else
                args0[i] = arg==nothing ? nothing : :($arg)
                args1[i] = arg==nothing ? :(::$typ) : :($arg::$typ)
            end
            args1 = insertdummy(args1)
        end

        # Process dissected AST
        if isnegated(trait)
            trait0_opposite = trait0 # opposite of `trait` below as that gets stripped of !
            trait = trait.args[2]
            val = :(::Type{$curmod.Not{$trait}})
        else
            trait0_opposite = Expr(:call, :!, trait0)  # generate the opposite
            val = :(::Type{$trait})
        end
        # Get line info for better backtraces
        ln = findline(tfn)
        if isa(ln, Expr)
            pushloc = Expr(:meta, :push_loc, ln.args[2], fname, ln.args[1])
            poploc = Expr(:meta, :pop_loc)
        else
           pushloc = poploc = nothing
        end
        # create the function containing the logic.  Do it separately if old-school parametric functions
        # (TODO: delete the oldfn_syntax branches in Julia 1.0 (or 0.7 already?))
        retsym = gensym()
        if hasmac
            fn = if oldfn_syntax
                :(@dummy $fname{$(typs...)}($val, $(strip_kw(args1)...); $(kwargs...)) = ($pushloc; $retsym = $fbody; $poploc; $retsym))
            else
                :(@dummy $fname($val, $(strip_kw(args1)...); $(kwargs...)) where {$(typs...)} = ($pushloc; $retsym = $fbody; $poploc; $retsym))
            end
            fn.args[1] = mac # replace @dummy
        else
            fn = if oldfn_syntax
                :($fname{$(typs...)}($val, $(strip_kw(args1)...); $(kwargs...)) = ($pushloc; $retsym = $fbody; $poploc; $retsym))
            else
                :($fname($val, $(strip_kw(args1)...); $(kwargs...)) where {$(typs...)} = ($pushloc; $retsym = $fbody; $poploc; $retsym))
            end
        end
        # Create the trait dispatch function
        ex = fn
        key = (macro_module, fname, typs0, strip_kw(args0), trait0_opposite)
        if !(key ∈ keys(dispatch_cache)) # define trait dispatch function
            if !haskwargs
                ex = if oldfn_syntax
                    quote
                        $fname{$(typs...)}($(args1...)) = (Base.@_inline_meta(); $fname($curmod.trait($trait),
                                                                                $(strip_tpara(strip_kw(args1))...)
                                                                                )
                                                           )
                        $ex
                    end
                else
                    quote
                        $fname($(args1...)) where {$(typs...)} = (Base.@_inline_meta(); $fname($curmod.trait($trait),
                                                                                $(strip_tpara(strip_kw(args1))...)
                                                                                )
                                                                  )
                        $ex
                    end
                end
            else
                ex = if oldfn_syntax
                    quote
                        $fname{$(typs...)}($(args1...);kwargs...) = (Base.@_inline_meta(); $fname($curmod.trait($trait),
                                                                                $(strip_tpara(strip_kw(args1))...);
                                                                                kwargs...
                                                                                )
                                                                     )
                        $ex
                    end
                else
                    quote
                        $fname($(args1...);kwargs...) where {$(typs...)} = (Base.@_inline_meta(); $fname($curmod.trait($trait),
                                                                                $(strip_tpara(strip_kw(args1))...);
                                                                                kwargs...
                                                                                )
                                                                            )
                        $ex
                    end
                end
            end
            dispatch_cache[key] = (haskwargs, args0)
        else # trait dispatch function already defined
            if dispatch_cache[key][1]!=haskwargs
                ex = :(error("""
                             Trait-functions can have keyword arguments.
                             But if so, add the same to both `Tr` and `!Tr`, but they can have different default values.
                             """))
            end
            if dispatch_cache[key][2]!=args0
                ex = :(error("""
                             Trait-functions can have default arguments.
                             But if so, add them to both `Tr` and `!Tr`, and they both need identical values!
                             """))
            end
            delete!(dispatch_cache, key) # permits function redefinition if that's what we want
        end
        ex
    end
end

"""
Defines a function dispatching on a trait. Examples:

```julia
@traitfn f(x::X,y) where {X;  Tr1{X}} = ...
@traitfn f(x::X,y) where {X; !Tr1{X}} = ...

@traitfn f(x::X,y::Y) where {X,Y;  Tr2{X,Y}} = ...
@traitfn f(x::X,y::Y) where {X,Y; !Tr2{X,Y}} = ...
```

Note that the second example is just syntax sugar for `@traitfn f(x::X,y::Y) where {X,Y; Not{Tr1{X,Y}}} = ...`.
"""
macro traitfn(tfn)
    esc(traitfn(tfn, __module__))
end

######
## Helpers
######

# true if :(!(Tr{x}))
isnegated(t::Expr) = t.head==:call
isnegated(t::Symbol) = false

# [:(x::X)] -> [:x]
# also takes care of :...
strip_tpara(args::Vector) = Any[strip_tpara(a) for a in args]
strip_tpara(a::Symbol) = a
function strip_tpara(a::Expr)
    if a.head==:(::)
        return a.args[1]
    elseif a.head==:...
        return Expr(:..., strip_tpara(a.args[1]))
    elseif a.head==:kw
        @assert length(a.args)==2
        return Expr(:kw, strip_tpara(a.args[1]), a.args[2])
    else
        error("Cannot parse argument: $a")
    end
end

# [:(x::X=4)] -> [:x::X]
# i.e. strips defaults arguments
# also takes care of :...
strip_kw(args::Vector) = Any[strip_kw(a) for a in args]
strip_kw(a) = a
function strip_kw(a::Expr)
    if a.head==:(::) || a.head==:...
        return a
    elseif a.head==:kw
        @assert length(a.args)==2
        return a.args[1]
    else
        error("Cannot parse argument: $a")
    end
end

# insert dummy: ::X -> gensym()::X
# also takes care of :...
# not needed for kwargs
insertdummy(args::Vector) = Any[insertdummy(a) for a in args]
insertdummy(a::Symbol) = a
function insertdummy(a::Expr)
    if a.head==:(::) && length(a.args)==1
        return Expr(:(::), gensym(), a.args[1])
    elseif a.head==:...
        return Expr(:..., insertdummy(a.args[1]))
    else
        return a
    end
end

# generates: X1, X2,... or x1, x2.... (just symbols not actual TypeVar)
struct GenerateTypeVars{CASE} end
Base.iterate(::GenerateTypeVars{:upcase}, state=1) = (Symbol("X$state"), state+1) # X1,..
Base.iterate(::GenerateTypeVars{:lcase}, state=1) = (Symbol("x$state"), state+1)  # x1,...

####
# Annotating the source location
####

function findline(ex::Expr)
    ex.head == :line && return ex
    for a in ex.args
        ret = findline(a)
        isa(ret, Expr) && return ret
    end
    nothing
end
findline(arg) = nothing

####
# Extras
####

# # This does not work.  Errors on compilation with:
# # ERROR: LoadError: UndefVarError: Tr not defined
# function check_traitdispatch(Tr; nlines=5, Args=(Int,))
#     @traitfn fn_test(x::::Tr) = 1
#     @traitfn fn_test(x::::(!Tr)) = 2
#     @assert llvm_lines(fn_test, Args) == nlines
# end

"""
    check_fast_traitdispatch(Tr, Args=(Int,), nlines=6, verbose=false)

Macro to check whether a trait-dispatch is fast (i.e. as fast as an
ordinary function call) or whether dispatch is slow (dynamic).  Only
works with single parameters traits (so far).

Optional arguments are:
- Type parameter to the trait (default `Int`)
- Verbosity (default `false`)

Example:

    @check_fast_traitdispatch IsBits
    @check_fast_traitdispatch IsBits String true

NOTE: This only works when code-coverage is disabled!  Thus, do not
use this macro in tests (or disable `coverage=true` in your
`.travis.yml` script), or it will error.

NOTE: This does not seem to work all the time.  At least in Julia 0.7 this test is
negative for IsIndexLinear, however, it's fast.

TODO: This is rather ugly.  Ideally this would be a function but I ran
into problems, see source code.  Also the macro is ugly.  PRs welcome...
"""
macro check_fast_traitdispatch(Tr, Arg=:Int, verbose=false)
    if Base.JLOptions().code_coverage==1
        warn("The SimpleTraits.@check_fast_traitdispatch macro only works when running Julia without --code-coverage")
        return nothing
    end
    test_fn = gensym()
    test_fn_null = gensym()
    nl = gensym()
    nl_null = gensym()
    out = gensym()
    esc(quote
        $test_fn_null(x) = 1
        $nl_null = SimpleTraits.llvm_lines($test_fn_null, ($Arg,))
        @traitfn $test_fn(x::::$Tr) = 1
        @traitfn $test_fn(x::::(!$Tr)) = 2
        $nl = SimpleTraits.llvm_lines($test_fn, ($Arg,))
        $out = $nl == $nl_null
        if $verbose && !$out
            println("Number of llvm code lines $($nl) but should be $($nl_null).")
        end
        $out
    end)
end

"Returns number of llvm-IR lines for a call of function `fn` with argument types `args`"
function llvm_lines(fn, args)
    io = IOBuffer()
    InteractiveUtils.code_llvm(io, fn, args)
    #Base.code_native(io, fn, args)
    count(c->c=='\n', String(take!(copy(io))))
end

include("base-traits.jl")

end # module


text:
"""
    @exception_with_line(
        exception_name::Symbol,
        docstring::Union{Expr, String},
        error_message_bits::Tuple{Vararg{Union{Expr, String}}},
    ) -> Expr

Create an exception with two fields: a path to the file and a line number.

# Arguments
- `exception_name::Symbol`: name of the exception
- `docstring::Union{Expr, String}`: documentation string
- `error_message_bits::Tuple{Vararg{Union{Expr, String}}}`: strings and expressions which
  will be interpolated in the `showerror` output

# Returns
- `Expr`: an exception definition (struct + `showerror` overload)
"""
@exception exception_with_line file::String line::Int context = begin
    @fields(
        "`file::String`: absolute or relative path to a file that is associated with the
        exception",
        "`line::Int`: line number of the specified file to which this exception is
        associated",
    )
    @docstring exception_with_line
end


text:
__precompile__()

using Ensemblearn

# Based on https://github.com/eriklindernoren/ML-From-Scratch/blob/master/mlfromscratch/supervised_learning/random_forest
# - https://machinelearningmastery.com/implement-random-forest-scratch-python/
# - https://towardsdatascience.com/random-forests-and-decision-trees-from-scratch-in-python-3e4fa5ae4249

export
    RandomForest,
    fit!,
    predict,
    predict_probs

# TODO: type?
mutable struct RandomForest
    n_trees::Int64
    max_depth::Int64  # TODO: any of None or Int
    max_features::Union{Int64, Nothing}
    min_samples_split::Int64
    min_gain::Float64
    impurity::String
    #bootstrap::Bool
    #oob_score::Bool
    seed::Int64
    trees::Union{Array{DecisionTree, 1}, Nothing}
    feature_indices::Union{Array{Array{Int64,1},1}, Nothing}  # TODO: Array{Int64, 2} ?

    function RandomForest(n_trees=10::Int64, max_depth=3::Int64, max_features=nothing::Union{Int64, Nothing},
                          min_samples_split=2::Int64, min_gain=0.0::Float64, impurity="gini"::String, seed=123::Int64)
        trees = []

        for i=1:n_trees
            push!(trees, DecisionTree(min_samples_split, 1e-7, max_depth, impurity))
        end

        new(n_trees, max_depth, max_features, min_samples_split, min_gain, impurity, seed, trees, nothing)
    end

end

function fit!(model::RandomForest, x::Array{T,2}, y::Array{T2,1}) where {T <: Real, T2 <: Real}
    n_samples, n_features = size(x)

    # If max_features have not been defined => select it as
    # sqrt(n_features)
    if model.max_features == nothing
        model.max_features = round(Integer, sqrt(n_features))
    end

    if model.feature_indices == nothing
        model.feature_indices = []
    end

    # TODO assert dims of x and y matching
    x_y = hcat(x, y)

    # TODO: parallel -> https://docs.julialang.org/en/v1/manual/parallel-computing/
    for i in 1:model.n_trees
        # Choose one random subset of the data for each tree
        subset = sample(x_y, n_samples, seed=model.seed)  # TODO: replacement?
        x_subset, y_subset = subset[:, 1:end-1], subset[:, end]

        # Feature bagging (select random subsets of the features)
        idx = rand(1:n_features, model.max_features)

        # Save the indices of the features for prediction
        push!(model.feature_indices, idx)

        # Choose the features corresponding to the indices
        x_subset = x_subset[:, idx]

        # Fit the tree to the data
        fit!(model.trees[i], x_subset, y_subset)
    end
end

function predict(model::RandomForest, x::Array{T,2}) where T <: Real
    n_samples = size(x,1)
    y_preds = Array{Float64, 2}(undef, n_samples, model.n_trees)
    # Let each tree make a prediction on the data
    for i_tree in enumerate(model.trees)
        i, tree = i_tree
        # Indices of the features that the tree has trained on
        idx = model.feature_indices[i]

        # Make a prediction based on those features
        prediction = predict(tree, x[:, idx])

        y_preds[:, i] = prediction
    end

    y_pred = []
    # For each sample
    for i in 1:n_samples
        y = modes(y_preds[i, :])[1]
        push!(y_pred, y)
    end

    y_pred
end

function predict_probs(model::RandomForest, x::Array{T,2}; get_probs=false::Bool, index_class=true::Bool) where T <: Real
    n_samples = size(x,1)

    y_preds = Dict{Any, Array{Float64,1}}()

    # Let each tree make a prediction on the data
    for i_tree in enumerate(model.trees)
        i, tree = i_tree
        # Indices of the features that the tree has trained on
        idx = model.feature_indices[i]

        # Make a prediction based on those features
        prediction = predict_probs(tree, x[:, idx])

        #for kv in prediction:
        #    if !haskey(y_preds, )
        #end

        y_preds[:, i] = prediction
    end

    y_pred = []
    # For each sample
    for i in 1:n_samples
        if get_probs == true
            y = nothing
        else
            y = modes(y_preds[i, :])[1]
        end

        push!(y_pred, y)
    end

    y_pred
end


text:
module LinearKernelModels

using OffsetArrays, LinearAlgebra

export solve_for_kernel, compute_Cd, compute_r, compute_kerr

"""
    C, d = compute_Cd(S::AbstractVecOrMat, r::AbstractVector, trange::AbstractUnitRange)
    C, d = compute_Cd(S::AbstractVecOrMat, r::AbstractVector, nt::Integer)

Given one or more stimuli, where `S[t, i]` is the strength of the
`i`th stimulus at time `t`, and a response vs. time `r[t]`, calculate
`C` and `d` that can be used to solve for a kernel-matrix `k`, where
`k` convolved with `S` best reconstructs `r`. In particular,

    k = reshape(C \\ d, (nt, nstim))

`nt` (an integer) is the number of timepoints in the kernel.
If the kernel is allowed to look into `r`'s future, instead supply `trange = tpast:tfuture`
where `tpast` might typically be negative and `tfuture` typically positive.
(The meaning is the offset relative to the current moment `t=0`.)

If you have temporal gaps in your measured response, splice the blocks
together into a single vector, separating the blocks by stretches of
NaNs of length `nt` or larger. All entries in `S` should be
finite and represent the actual value of the stimulus at these
times. If you have periods of uncertainty about the stimulus history,
then for any time points within `nt` after these unknown
blocks, replace the corresponding entries in `r` with NaN. You can
then set the unknown entries in `S` to an arbitrary value (e.g.,
zero).

See also [`compute_r`](@ref), [`compute_kerrs`](@ref).
"""
function compute_Cd(S::AbstractVecOrMat, r::AbstractVector, trange::AbstractUnitRange)
    T, nstim = size(S, 1), size(S, 2)
    axes(S, 1) == 1:T || error("S must start indexing at t=1, got $(axes(S))")
    axes(r) == (1:T,) || error("r must agree with S's indexing, got $(axes(r)) versus $(axes(S))")
    kernelaxes = (trange, oftype(trange, Base.OneTo(nstim)))
    L = prod(length.(kernelaxes))
    C = zeros(L, L)
    LI = LinearIndices(OffsetArrays.IdentityUnitRange.(kernelaxes))
    w = (!).(isnan.(r))
    for i1 = 1:nstim, i2 = i1:nstim
        for τ1 in trange, τ2 in trange
            s = 0.0
            for t = max(1, 1-τ1, 1-τ2):min(T, T-τ1, T-τ2)
                s += w[t] * S[t+τ1, i1] * S[t+τ2, i2]
            end
            ind1, ind2 = LI[τ1, i1], LI[τ2, i2]
            C[ind1, ind2] = C[ind2, ind1] = s
        end
    end
    d = zeros(L)
    for i = 1:nstim
        for τ in trange
            s = 0.0
            for t = max(1, 1-τ):min(T, T-τ)
                p = r[t]*S[t+τ, i]
                s += ifelse(w[t], p, zero(p))
            end
            d[LI[τ, i]] = s
        end
    end
    return C, d
end
compute_Cd(S::AbstractVecOrMat, r::AbstractVector, nt::Integer) =
    compute_Cd(S, r, -(nt-1):0)

"""
    r = compute_r(S::AbstractVecOrMat, k::AbstractVecOrMat)

Given one or more stimuli, where `S[t, i]` is the strength of the
`i`th stimulus at time `t`, and a kernel `k[τ, i]` for each
stimulus, compute the summed-convolution between `S` and `k` to
predict the response `r`.
"""
function compute_r(S, k)
    T, nstim = size(S, 1), size(S, 2)
    trange = axes(k, 1)
    @assert(axes(k, 2) == 1:nstim)
    r = zeros(T)
    for i = 1:nstim
        for τ in trange
            kτi = k[τ,i]
            for t = max(1, 1-τ):min(T, T-τ)
                r[t] += S[t+τ,i]*kτi
            end
        end
    end
    return r
end

"""
    k = solve_for_kernel(S, r, isnz; rtol=sqrt(eps()))

Solve for a kernel, imposing constraints on entries that are allowed to be non-zero.
`S` is the stimulus and `r` the response vector. (See [`compute_Cd`](@ref) for more information.)
`isnz` must have the axes of the kernel, and be true in every slot where the solved kernel
may be non-zero. If you have no need for constraints, supply `trues(axs)` where `axs`
encodes the axes of the kernel.

# Example
```
julia> using LinearKernelModels, OffsetArrays

julia> S = rand(100,3);

julia> ktrue = OffsetArray(rand(5, 3), -4:0, 1:3)
OffsetArray(::Array{Float64,2}, -4:0, 1:3) with eltype Float64 with indices -4:0×1:3:
 0.155762  0.877222  0.457633
 0.88788   0.355702  0.0268748
 0.276609  0.914133  0.230344
 0.718295  0.158363  0.0402396
 0.195513  0.783832  0.5007

julia> ktrue[-4, 1] = 0
0

julia> isnz = ktrue .!= 0
OffsetArray(::BitArray{2}, -4:0, 1:3) with eltype Bool with indices -4:0×1:3:
 false  true  true
  true  true  true
  true  true  true
  true  true  true
  true  true  true

julia> r = compute_r(S, ktrue);

julia> k = solve_for_kernel(S, r, isnz)
OffsetArray(::Array{Float64,2}, -4:0, 1:3) with eltype Float64 with indices -4:0×1:3:
 6.99441e-15  0.877222  0.457633
 0.88788      0.355702  0.0268748
 0.276609     0.914133  0.230344
 0.718295     0.158363  0.0402396
 0.195513     0.783832  0.5007
```
"""
function solve_for_kernel(S::AbstractVecOrMat, r::AbstractVector, isnz::AbstractVecOrMat{Bool}; rtol=sqrt(eps()))
    trange = axes(isnz, 1)
    C, d = compute_Cd(S, r, trange)
    nconstrained = prod(length.(axes(isnz))) - sum(isnz)
    Q = zeros(size(C,1), nconstrained)
    col = 0
    for (i, isunconstrained) in enumerate(isnz)
        if !isunconstrained
            col += 1
            Q[i,col] = 1
        end
    end
    A = [C Q; Q' zeros(nconstrained, nconstrained)]
    dcat = [d; zeros(nconstrained)]
    F = svd(A)
    S = F.S
    Smax = maximum(S)
    flag = S .< rtol*Smax
    S[flag] .= Inf
    soln = F \ dcat
    reshape(soln[1:length(d)], axes(isnz))
end

"""
    kerr = compute_kerr(k, C, S, r)

Compute the estimated standard error in the parameters `k` of the
kernels, where `C` is determined from [`compute_Cd`](@ref), `S` is the
stimulus, and `r` the response.
"""
function compute_kerr(k, C, S, r)
    rbar = compute_r(S, k)
    N = 0
    Σ = 0.0
    for t = 1:length(r)
        rt = r[t]
        if !isnan(rt)
            Σ += (rt - rbar[t])^2
            N += 1
        end
    end
    lk = prod(length.(axes(k)))
    σ2 = Σ/(N-lk)  # replace with length(k)
    return reshape(sqrt(σ2) * sqrt.(diag(inv(C))), axes(k))
end

end # module


text:
using Documenter, DataUtils

makedocs(;
    modules=[DataUtils],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/DataUtils.jl/blob/{commit}{path}#L{line}",
    sitename="DataUtils.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/tkf/DataUtils.jl",
)


text:
# multiple_dispatch test

ENV["MPLBACKEND"]="tkagg"
using PyPlot
pygui(true)

mutable struct Grid
    nx::Int64
    dx::Float64
end

mutable struct Params
    K::Float64
    dt::Float64
    nsteps::Int64
    Tend::Float64
    Courant::Float64
end

function preprocess(method)
    if method == "advection_explicit"
        advection       = true
        advection_theta = false
        diffusion       = 0
        diffusion_theta = false
    elseif method == "advection_implicit"
        advection       = true
        advection_theta = 1.0
        diffusion       = 0
        diffusion_theta = false
    elseif method == "diffusion_explicit"
        advection       = 0
        advection_theta = false
        diffusion       = true
        diffusion_theta = false
    elseif method == "diffusion_implicit"
        advection_theta = false
        advection       = 0
        diffusion_theta = 1.0
        diffusion       = true
    elseif method == "advection_explicit_diffusion_explicit"
        advection       = true
        advection_theta = false
        diffusion       = true
        diffusion_theta = false
    elseif method == "advection_implicit_diffusion_explicit"
        advection       = true
        advection_theta = 1.0
        diffusion       = true
        diffusion_theta = false
    elseif method == "advection_explicit_diffusion_implicit"
        advection       = true
        advection_theta = false
        diffusion       = true
        diffusion_theta = 1.0
    elseif method == "advection_implicit_diffusion_implicit"
        advection       = true
        advection_theta = 1.0
        diffusion       = true
        diffusion_theta = 1.0
    end

    #model_components = [advection, advection_theta, diffusion, diffusion_theta]
    #return model_components
    return advection, advection_theta, diffusion, diffusion_theta

end

function initialize!(model_components)

    L      = 1000
    nx     = 1000
    dx     = L/nx
    K      = 0.1
    nsteps = 1000
    dt0    = 5.0
    Tend   = 100000.0

    # Courant number
    Courant = 0.98
    maxdx  = dx #max(dx,dy)

    # Unpack model_components
    advection       = model_components[1]
    advection_theta = model_components[2]
    diffusion       = model_components[3]
    diffusion_theta = model_components[4]

    # time step size based on diffusion
    if diffusion==true && diffusion_theta == false
        dx2   = maxdx*maxdx
        dt0    = min(dt0, Courant*0.5*dx2/K)
    end
    dt = dt0

    params = Params(K, dt, nsteps, Tend, Courant)

    h    = Array{Float64,1}(undef, nx)
    diag = 1.0*ones(Float64, nx)
    rhs  = Array{Float64,1}(undef, nx)

    # Save state
    state = [h, diag, rhs]

    # Set initial solution
    for i = 1:nx
        if nx*0.4 < i < nx*0.6
            h[i] = 0.4
        else
            h[i] = 0.2
        end
    end
    initial_condition = copy(h)

    grid   = Grid(nx, dx)

    return grid, state, initial_condition, params
end

function compdt!(params, model_components, grid, time, initial_condition)
    # time step size based on advection

    # Unpack model_components
    advection       = model_components[1]
    advection_theta = model_components[2]
    diffusion       = model_components[3]
    diffusion_theta = model_components[4]

    # Local variable dt
    #dt = 0.0

    if advection==true && advection_theta==false
        hmax = maximum(abs.(initial_condition))
        params.dt    = min(params.dt, params.Courant*grid.dx/hmax)

        # Don't exceed end time
        params.dt    = min(params.dt, params.Tend-time)
    end
    #println(dt)
    return params.dt
end

# Update function for explicit advection
function update!(state, grid, params, adv::Bool, diff::Int64, adv_theta::Bool, diff_theta::Bool)
    # Implementation corresponding to method:
    # explicit advection
    #println("Update state using explicit advection")
    state = advection_explicit!(state, grid, params)
    return state
end

# Update function for implicit advection
function update!(state, grid, params, adv::Bool, diff::Int64, adv_theta::Float64, diff_theta::Bool)
    # Implementation corresponding to method:
    # implicit advection
    #println("Update state using implicit advection")
    state = advection_implicit!(state, grid, params, adv_theta)
    return state
end

# Update function for explicit diffusion
function update!(state, grid, params, adv::Int64, diff::Bool, adv_theta::Bool, diff_theta::Bool)
    # Implementation corresponding to method:
    # explicit diffusion
    #println("Update state using explicit diffusion")
    state = diffusion_explicit!(state, grid, params)
    return state
end

# Update function for implicit diffusion
function update!(state, grid, params, adv::Int64, diff::Bool, adv_theta::Bool, diff_theta::Float64)
    # Implementation corresponding to method:
    # implicit diffusion
    #println("Update state using implicit diffusion")
    state = diffusion_implicit!(state, grid, params, diff_theta)
    return state
end

# Update function for explicit advection and explicit diffusion
function update!(state, grid, params, adv::Bool, diff::Bool, adv_theta::Bool, diff_theta::Bool)
    # Implementation corresponding to method:
    # explicit advection and explicit diffusion
    #println("Update state using explicit advection and explicit diffusion")
    state = advection_explicit!(state, grid, params)
    state = diffusion_explicit!(state, grid, params)
    return state
end

# Update function for implicit advection and explicit diffusion
function update!(state, grid, params, adv::Bool, diff::Bool, adv_theta::Float64, diff_theta::Bool)
    # Implementation corresponding to method:
    # implicit advection and explicit diffusion
    #println("Update state using implicit advection and explicit diffusion")
    state = advection_implicit!(state, grid, params, adv_theta)
    state = diffusion_explicit!(state, grid, params)
    return state
end

# Update function for explicit advection and implicit diffusion
function update!(state, grid, params, adv::Bool, diff::Bool, adv_theta::Bool, diff_theta::Float64)
    # Implementation corresponding to method:
    # explicit advection and implicit diffusion
    #println("Update state using explicit advection and implicit diffusion")
    state = advection_explicit!(state, grid, params)
    state = diffusion_implicit!(state, grid, params, diff_theta)
    return state
end

# Update function for implicit advection and implicit diffusion
function update!(state, grid, params, adv::Bool, diff::Bool, adv_theta::Float64, diff_theta::Float64)
    # Implementation corresponding to method:
    # implicit advection and implicit diffusion
    #println("Update state using implicit advection and implicit diffusion")
    state = advection_implicit!(state, grid, params, adv_theta)
    state = diffusion_implicit!(state, grid, params, diff_theta)
    return state
end

# Explicit advection computation
function advection_explicit!(state, grid, params)

    nx  = grid.nx
    dx  = grid.dx

    dt  = params.dt

    h    = state[1]
    diag = state[2]
    rhs  = state[3]

    for i=2:nx-1
        uxm = (h[i]-h[i-1])/dx
        uxp = (h[i+1]-h[i])/dx
        am  = min(h[i],0.0)
        ap  = max(h[i],0.0)
        adv = ap*uxm + am*uxp
        @inbounds rhs[i] +=  -dt * adv
    end

    state = [h, diag, rhs]

    return state
end

# Implicit advection computation
function advection_implicit!(state, grid, params, adv_theta)
    nx  = grid.nx
    dx  = grid.dx

    dt  = params.dt

    h    = state[1]
    diag = state[2]
    rhs  = state[3]

    for i=2:nx-1
        uxm = -h[i-1]/dx
        uxp = h[i+1]/dx
        am  = min(h[i],0.0)
        ap  = max(h[i],0.0)
        adv_expl = ap*uxm + am*uxp
        adv_impl = ap/dx - am/dx
        @inbounds rhs[i]  += -dt * adv_expl
        @inbounds diag[i] +=  dt * adv_impl
    end

    state = [h, diag, rhs]

    return state
end

# Explicit advection computation
function diffusion_explicit!(state, grid, params)

    nx  = grid.nx
    dx  = grid.dx
    dx2 = dx*dx

    K   = params.K
    dt  = params.dt
    dtK = dt*K

    h    = state[1]
    diag = state[2]
    rhs  = state[3]

    #hnew = copy(h)

    for i=2:nx-1
        @inbounds rhs[i] += dtK * (h[i+1] - 2.0*h[i] + h[i-1]) / dx2
        #@inbounds hnew[i,j] = h[i,j] + dt*K * ( (h[i+1,j] + h[i-1,j] - 2.0*h[i,j])/(dx*dx) + (h[i,j+1]+h[i,j-1] - 2.0*h[i,j])/(dy*dy) )
    end

    state = [h, diag, rhs]

    return state

end

# Implicit diffusion computation
function diffusion_implicit!(state, grid, params, diff_theta)
    nx  = grid.nx
    dx  = grid.dx
    dx2 = dx*dx

    K   = params.K
    dt  = params.dt
    dtK = dt*K

    h    = state[1]
    diag = state[2]
    rhs  = state[3]

    #hnew = copy(h)

    for i=2:nx-1
        @inbounds diag[i] += dtK * 2.0 / dx2
        @inbounds rhs[i] += dtK * (h[i+1] + h[i-1]) / dx2
        #@inbounds hnew[i,j] = h[i,j] + dt*K * ( (h[i+1,j] + h[i-1,j] - 2.0*h[i,j])/(dx*dx) + (h[i,j+1]+h[i,j-1] - 2.0*h[i,j])/(dy*dy) )
    end

    state = [h, diag, rhs]

    return state
end

# Functio for solving a diagonal system of equations
function solve!(state, grid)
    h    = state[1]
    diag = state[2]
    rhs  = state[3]
    nx   = grid.nx

    for i = 1:nx
        h[i] = rhs[i]/diag[i]
        diag[i] = 1.0
    end

    return state
end

function plot_solution(n, time, h)
    clf()
    display(plot(h, linestyle="-",marker="o"))
    title("Time = $time s. (timestep $n)")
    pause(0.1)
end

# Test multiple dispatch for update function
# based on the contents of "parameters"
# For now for a 1D model
# depending on the argument, 8 methods can be invoked
# (through multiple dispatch of UPDATE function):
#   - explicit advection
#   - explicit diffusion
#   - (diagonally-)implicit advection
#   - (diagonally-)implicit diffusion
#   - explicit advection, explicit diffusion
#   - (diagonally-)implicit advection, explicit diffusion
#   - (diagonally-)implicit advection, explicit diffusion
#   - (diagonally-)implicit advection, (diagonally-)implicit diffusion

function test_multiple_dispatch(;method="advection_explicit")

    # Preprocess function arguments to determine model components
    # needed for multiple dispatch of UPDATE function
    advection, advection_theta, diffusion, diffusion_theta = preprocess(method)

    # Pack model components
    model_components = [advection, advection_theta, diffusion, diffusion_theta]

    # Initialize model
    grid, state, initial_condition, params = initialize!(model_components)

    # Number of grid points in x-direction
    nx = grid.nx

    # Unpack state
    h    = state[1]
    diag = state[2]
    rhs  = state[3]

    time = 0.0
    for n = 1:params.nsteps

        # Determine time step (possible time step restriction from explicit model components)
        params.dt     = compdt!(params, model_components, grid, time, initial_condition)

        # Update time
        time  += params.dt

        # Update state
        rhs    = copy(h)
        state  = [h, diag, rhs]
        state  = update!(state, grid, params, advection, diffusion, advection_theta, diffusion_theta)

        # Solve system of equations (diagonal for now)
        state  = solve!(state, grid)

        # Plot solution
        if mod(n,50) == 0
            plot_solution(n, time, h)
        end

        # Exit at end time
        if time == params.Tend
            break
        end
    end

    close()
    println("\n")
end

#@time test_multiple_dispatch(;method="advection_explicit")
#@time test_multiple_dispatch(;method="advection_implicit")
#@time test_multiple_dispatch(;method="diffusion_explicit")
#@time test_multiple_dispatch(;method="diffusion_implicit")
#@time test_multiple_dispatch(;method="advection_explicit_diffusion_explicit")
#@time test_multiple_dispatch(;method="advection_implicit_diffusion_explicit")
#@time test_multiple_dispatch(;method="advection_explicit_diffusion_implicit")
@time test_multiple_dispatch(;method="advection_implicit_diffusion_implicit")


text:
using Documenter
using BioMart

makedocs(
    sitename = "BioMart",
    format = Documenter.HTML(),
    modules = [BioMart]
)


text:
# These functions are only conditionally loaded with ColorSchemes.jl
function _colordither(T, img, alg, cs::ColorSchemes.ColorScheme; kwargs...)
    return _colordither(T, img, alg, cs.colors)
end

function _colordither(T, img, alg, csname::Symbol; kwargs...)
    cs = ColorSchemes.colorschemes[csname]
    return _colordither(T, img, alg, cs.colors; kwargs...)
end


text:
using ValueOrientedRiskManagementInsurance
using DataFrames
using Test
using Distributions

println("start NonLife.jl...")

function logpar2statpar(m_s::Vector)
  Real[exp(m_s[1] + m_s[2]^2 / 2),
       √(exp(m_s[2]^2)-1) * exp(m_s[1] + m_s[2]^2 / 2)]
end

function statpar2logpar(μ_σ::Vector)
  vc = μ_σ[2]/μ_σ[1]
  Real[log(μ_σ[1]) - 0.5 * log(1+vc^2), √(log(1+vc^2))]
end

"""
   dr = ( θ(t)-ar ) dt + σ dW
"""
mutable struct HullWhite
  t₀::Real
  T::Real
  δt::Real
  θ::Vector{Real}
  a::Real
  σ::Real
end

function HullWhite(t₀::Real,
                   δt::Real,
                   θ::Vector{Real},
                   a::Real,
                   σ::Real)
  T = t₀ + length(θ) * δt
  HullWhite(t₀, T, δt, θ, a, σ)
end

periods_int(t::Real, hw::HullWhite) = int(div(t, hw.Δt))
periods_frac(t::Real, hw::HullWhite) = t/hw.Δt - period_int(t, hw)

function r(rstart::Real, k::Int, Δk::Int, hw::HullWhite)
  μ =
    exp(-hw.a*Δk*hw.δt)/hw.a *
    (hw.a * rstart  +
     sum([(exp(-hw.a*𝑗*hw.δt) - exp(-hw.a*(𝑗-1)*hw.δt)) *
           hw.θ[k+𝑗-1] for 𝑗 ∈ 1:Δk]))
  std = hw.σ * √( (1-2exp(-2hw.a*Δk*hw.δt)) / (2hw.a) )
  Normal(μ, std)
end

nllobs = [:fire, :liab, :theft]
n_nl = length(nllobs)
nl_names = [ucfirst(string(nllobs[𝑖])) for 𝑖 ∈ 1:n_nl]
df_claims = Vector(n_nl)
claims = Vector(n_nl)
res = Vector(n_nl)
β = Vector(n_nl)

claimpath = "test/NonLife_Input_Claims_"
for 𝑖 ∈ 1:n_nl
  df_claims[𝑖] =
    CSV.read(claimpath * nl_names[𝑖]  * ".csv",
              header = false)
  res[𝑖] = Mack(df_claims[𝑖])
  β[𝑖] = res[𝑖].futureclaims / sum(res[𝑖].futureclaims)
end

#
# mean, sde -> logmean, logsigma
# R[i] = LogNormal(log_mean[i], log_sigma[i])
# v⋅mack[i].β * R[i]
#
# v: stochastic discont
# β
# R[i] = stochastic total undiscounted future paiments,
#

println("...end NonLife.jl")


text:
module CpelTdm

## Dependencies
using XAM                       # For BAM/SAM files
using FASTX                     # For FASTA/FASTQ files
using BED                       # For BED files
using GenomicFeatures           # Find overlaps between coordinates
using Calculus                  # For gradient computation
using Combinatorics             # For permutation test
using Dates                     # For printing with date and hour
using Distributed               # For pmap
using Distributions             # For simulations
using Optim                     # For estimating parameters
using Random                    # For simulations
using LinearAlgebra             # For matrix operations
using DelimitedFiles            # For reading and writing delimited files
using MultipleTesting           # For BH correction

## CONSTANTS
const BG_BUFFER = 50000                     # Bytes of bedGraph records until write
const GFF_BUFFER = 500000                   # Bytes of GFF records until write
const THRESH_MAPQ = 39                      # MAPQ threshold (-10*log(p)) only true uni-reads
const FLAGS_ALLOWED = [0,16,83,99,147,163]  # Flags allowed in BAM recs
const ETA_MAX_ABS=5.0                       # Parameters space
const LOG2 = log(2)                         # Ln(2)

## Includes
include("Structures.jl")
include("IO.jl")
include("Bioinformatics.jl")
include("Inference.jl")
include("HypothesisTesting.jl")
include("Simulations.jl")

## Exports
export cpel_tdm

end # module


text:
function AbstractFrames(arg0::TimeScales, arg1::Supplier)
    return AbstractFrames((TimeScales, Supplier), arg0, arg1)
end

function get_cirf(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getCIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_ecliptic(obj::AbstractFrames, arg0::IERSConventions)
    return jcall(obj, "getEcliptic", Frame, (IERSConventions,), arg0)
end

function get_eme2000(obj::AbstractFrames)
    return jcall(obj, "getEME2000", FactoryManagedFrame, ())
end

function get_eop_history(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getEOPHistory", EOPHistory, (IERSConventions, jboolean), arg0, arg1)
end

function get_frame(obj::AbstractFrames, arg0::Predefined)
    return jcall(obj, "getFrame", Frame, (Predefined,), arg0)
end

function get_gcrf(obj::AbstractFrames)
    return jcall(obj, "getGCRF", Frame, ())
end

function get_gtod(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getGTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_gtod(obj::AbstractFrames, arg0::jboolean)
    return jcall(obj, "getGTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_icrf(obj::AbstractFrames)
    return jcall(obj, "getICRF", Frame, ())
end

function get_itrf(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getITRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_itrf(obj::AbstractFrames, arg0::ITRFVersion, arg1::IERSConventions, arg2::jboolean)
    return jcall(obj, "getITRF", VersionedITRF, (ITRFVersion, IERSConventions, jboolean), arg0, arg1, arg2)
end

function get_itrf_equinox(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getITRFEquinox", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_mod(obj::AbstractFrames, arg0::IERSConventions)
    return jcall(obj, "getMOD", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_mod(obj::AbstractFrames, arg0::jboolean)
    return jcall(obj, "getMOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_pz9011(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getPZ9011", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_teme(obj::AbstractFrames)
    return jcall(obj, "getTEME", FactoryManagedFrame, ())
end

function get_tirf(obj::AbstractFrames, arg0::IERSConventions)
    return jcall(obj, "getTIRF", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_tirf(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getTIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(obj::AbstractFrames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(obj::AbstractFrames, arg0::jboolean)
    return jcall(obj, "getTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_veis1950(obj::AbstractFrames)
    return jcall(obj, "getVeis1950", FactoryManagedFrame, ())
end


text:
function policy_opt(p ;sⁱ = false, sᶠ = false, sᴱ = false)
  poltype = [sⁱ, sᶠ, sᴱ]
  npol    = sum(poltype)
  polinit = ones(npol)*.05

# first solve without policy
  eqInit = [ 2.00, 0.55, 0.06,
             0.72, 0.83, 1.73]

  eq,res = solveBGP(p,eqInit)
  eqInit.= res.zero

  cevMax = [calCev(eq.g, eq.cactivtot,eq.g, eq.cactivtot,p)]
  @printf("OPTIMAL POLICY\n")
  @printf("----------------------------\n")
  @printf("Baseline Welfare is %3.4f \n", cevMax[1])

  # optimization
  # ------------------------------------------------------------------------------------------------
  opt = Opt(:LN_NELDERMEAD, npol)
  ftol_rel!(opt,1e-8)
  min_objective!(opt, (x, grad) -> policy_obj(x, poltype, p, eqInit, cevMax, eq.g, eq.cactivtot))
  minf,minpol,ret = NLopt.optimize(opt, polinit)
  nEvals = opt.numevals
  #-------------------------------------------------------------------------------------------------

  if ret == :FTOL_REACHED
  	@printf("Optimal policy found!\n\n")
  	optpolval = zeros(3)
    optpolval[poltype] .= minpol
	res = Dict(:sⁱ=> optpolval[1], :sᶠ=>optpolval[2], :sᴱ=>optpolval[3], :cev=>-minf, :ret=>ret, :opt=>opt, :nEvals=>nEvals)
  	return res
  else
	@printf("Optimal policy could not solved!\n\n")
	return ret
  end

end

function policy_obj(polguess, poltype, p, eqInit, cevMax, gBase, cactivtotBase)

  polval = zeros(3)
  polval[poltype] .= polguess

  pNew = Params(p, (sⁱ = polval[1], sᶠ = polval[2], sᴱ = polval[3]))

  eq, res = solveBGP(pNew, eqInit)

  if res.f_converged
  	eqInit .= res.zero
  	# welfare calculation
  	cev = calCev(gBase, cactivtotBase, eq.g, eq.cactivtot, pNew)
  else
    cev = -Inf
  end
  if cev>cevMax[1]
  	@printf("Welfare improved to %3.4f (", cev)
	poltype[1] ? @printf("sⁱ = %2.2f", pNew.sⁱ) : nothing
	poltype[1]&(poltype[2] | poltype[3]) ? @printf(",") : nothing
	poltype[2] ? @printf("sᶠ = %2.2f", pNew.sᶠ) : nothing
	poltype[2]&poltype[3] ? @printf(",") : nothing
	poltype[3] ? @printf("sᴱ = %2.2f", pNew.sᴱ) : nothing
	@printf(")\n")
	cevMax[1] = cev
  end
  return -cev
end

function calCev(gBase, cactivtotBase, gPol, cactivtotPol, p)
#	welf  = (1.0/(1.0-p.σ))*((eq.cactivtot^((1.0 - p.σ)/(p.ε-1.0))/(p.ρ-(1.0-p.σ)*eq.g))-(1.0/p.ρ))
	cev  = ((p.ρ-(1.0-p.σ)*gBase)/(p.ρ-(1.0-p.σ)*gPol))^(1.0/(1.0-p.σ))*((cactivtotPol/cactivtotBase)^(1.0/(p.ε-1.0)));
end