text:
module Contour

using StaticArrays

include("interpolate.jl")

export
    ContourLevel,
    Curve2,
    contour,
    contours,
    level,
    levels,
    lines,
    coordinates

import Base: push!, length, eltype, show

struct Curve2{T}
    vertices::Vector{SVector{2,T}}
end
Curve2(::Type{T}) where {T} = Curve2(SVector{2,T}[])
show(io::IO, ::MIME"text/plain", c2::Curve2) = write(io, "$(typeof(c2))\n  with $(length(c2.vertices)-1) vertices")
show(io::IO, ::MIME"text/plain", c2s::Vector{Curve2{T}}) where {T} = write(io, "$(typeof(c2s))\n  $(length(c2s)) contour line(s)")

struct ContourLevel{T}
    level::T
    lines::Vector{Curve2{T}}
end
ContourLevel(h::T) where {T <: AbstractFloat} = ContourLevel(h, Curve2{T}[])
ContourLevel(h::T) where {T} = ContourLevel(Float64(h))
show(io::IO, ::MIME"text/plain", cl::ContourLevel) = write(io, "$(typeof(cl))\n  at $(level(cl)) with $(length(lines(cl))) line(s)")
show(io::IO, ::MIME"text/plain", cls::Vector{ContourLevel{T}}) where {T} = write(io, "$(typeof(cls))\n  $(length(cls)) contour level(s)")
"""
`lines(c)` Extracts an iterable collection of isolines from a contour level.
Use [`coordinates`](@ref) to get the coordinates of a line.
"""
lines(cl::ContourLevel) = cl.lines
"""
`level(c)` Indicates the `z`-value at which the contour level `c` was traced.
"""
level(cl::ContourLevel) = cl.level

struct ContourCollection{Tlevel<:ContourLevel}
    contours::Vector{Tlevel}
end
ContourCollection() = ContourCollection(Float64)
ContourCollection(::Type{Tlevel}) where {Tlevel} = ContourCollection(ContourLevel{Tlevel}[])
show(io::IO, ::MIME"text/plain", cc::ContourCollection) = write(io, "$(typeof(cc))\n with $(length(levels(cc))) level(s).")

"""
Turns the output of [`contours`](@ref) into an iterable with each of the traced
contour levels. Each of the objects support [`level`](@ref) and
[`coordinates`](@ref).
"""
levels(cc::ContourCollection) = cc.contours

"""
`contour(x, y, z, level::Number)` Trace a single contour level, indicated by the
argument `level`.

You'll usually call [`lines`](@ref) on the output of `contour`, and then iterate
over the result.
"""
function contour(x, y, z, level::Number)
    if !(axes(x) == (axes(z,1),) && axes(y) == (axes(z,2),) || axes(x) == axes(y) == axes(z))
        throw(ArgumentError("Incompatible input axes in `Contour.contour`."))
    end
    trace_contour(x, y, z, level, get_level_cells(z, level))
end

"""
`contours` returns a set of isolines.

You'll usually call [`levels`](@ref) on the output of `contours`.
"""
contours(::Any...) = error("This method exists only for documentation purposes")

"""
`contours(x,y,z,levels)` Trace the contour levels indicated by the `levels`
argument.
"""
contours(x, y, z, levels) = ContourCollection([contour(x, y, z, l) for l in levels])

"""
`contours(x,y,z,Nlevels::Integer)` Trace `Nlevels` contour levels at heights
chosen by the library (using the  [`contourlevels`](@ref) function).
"""
function contours(x, y, z, Nlevels::Integer)
    contours(x, y, z, contourlevels(z, Nlevels))
end

"""
`contours(x,y,z)` Trace 10 automatically chosen contour levels.
"""
contours(x, y, z) = contours(x, y, z, 10)

"""
`contourlevels(z,n)` Examines the values of `z` and chooses `n` evenly spaced
levels to trace.
"""
function contourlevels(z, n)
    zmin, zmax = extrema(z)
    dz = (zmax - zmin) / (n + 1)
    range(zmin + dz; step = dz, length = n)
end

"""
`coordinates(c)` Returns the coordinates of the vertices of the contour line as
a tuple of lists.
"""
function coordinates(c::Curve2{T}) where {T}
    N = length(c.vertices)
    xlist = Vector{T}(undef, N)
    ylist = Vector{T}(undef, N)

    for (i, v) in enumerate(c.vertices)
        xlist[i] = v[1]
        ylist[i] = v[2]
    end
    xlist, ylist
end

# The marching squares algorithm defines 16 cell types
# based on the edges that a contour line enters and exits
# through. The edges of the cells are identified using
# compass directions, while the vertices are ordered as
# follows:
#
#      N
#  4 +---+ 3
# W  |   |  E
#  1 +---+ 2
#      S
#
# Each cell type is identified with 4 bits, with each
# bit corresponding to a vertex (MSB -> 4, LSB -> 1).
# A bit is set for vertex v_i is set if z(v_i) > h. So a cell
# where a contour line only enters from the W edge and exits
# through the N edge will have the cell type: 0b0111
# Note that there are two cases where there are two
# lines crossing through the same cell: 0b0101, 0b1010.
const N, S, E, W = (UInt8(1)), (UInt8(2)), (UInt8(4)), (UInt8(8))
const NS, NE, NW = N|S, N|E, N|W
const SN, SE, SW = S|N, S|E, S|W
const EN, ES, EW = E|N, E|S, E|W
const WN, WS, WE = W|N, W|S, W|E
const NWSE = NW | 0x10 # special (ambiguous case)
const NESW = NE | 0x10 # special (ambiguous case)

# Maps cell type to crossing types for non-ambiguous cells
const edge_LUT = (SW, SE, EW, NE, 0x0, NS, NW, NW, NS, 0x0, NE, EW, SE, SW)

# The way a contour crossing goes through a cell is labeled
# by combining compass directions (e.g. a NW crossing connects
# the N edge and W edges of the cell).  The Cell type records
# the type of crossing that a cell contains.  While most
# cells will have only one crossing, cell type 5 and 10 will
# have two crossings.
function get_next_edge!(cells::Dict, key, entry_edge::UInt8)
    cell = pop!(cells, key)
    if cell == NWSE
        if entry_edge == N || entry_edge == W
            cells[key] = SE
            cell = NW
        else #SE
            cells[key] = NW
            cell = SE
        end
    elseif cell == NESW
        if entry_edge == N || entry_edge == E
            cells[key] = SW
            cell = NE
        else #SW
            cells[key] = NE
            cell = SW
        end
    end
    return cell ⊻ entry_edge
end

# N, S, E, W
const next_map = ((0,1), (0,-1), (1,0), (-1,0))
const next_edge = (S,N,W,E)

@inline function advance_edge(ind, edge)
    n = trailing_zeros(edge) + 1
    nt = ind .+ next_map[n]
    return nt, next_edge[n]
end

@inline function get_first_crossing(cell)
    if cell == NWSE
        return NW
    elseif cell == NESW
        return NE
    else
        return cell
    end
end

@inline function _get_case(z, h)
    case = z[1] > h ? 0x01 : 0x00
    z[2] > h && (case |= 0x02)
    z[3] > h && (case |= 0x04)
    z[4] > h && (case |= 0x08)
    case
end

function get_level_cells(z, h::Number)
    cells = Dict{Tuple{Int,Int},UInt8}()
    x_ax, y_ax = axes(z)

    @inbounds for xi in first(x_ax):last(x_ax)-1
        for yi in first(y_ax):last(y_ax)-1
            elts = (z[xi, yi], z[xi + 1, yi], z[xi + 1, yi + 1], z[xi, yi + 1])
            case = _get_case(elts, h)

            # Contour does not go through these cells
            if iszero(case) || case == 0x0f
                continue
            end

            # Process ambiguous cells (case 5 and 10) using
            # a bilinear interpolation of the cell-center value.
            if case == 0x05
                cells[(xi, yi)] = 0.25*sum(elts) >= h ? NWSE : NESW
            elseif case == 0x0a
                cells[(xi, yi)] = 0.25*sum(elts) >= h ? NESW : NWSE
            else
                cells[(xi, yi)] = edge_LUT[case]
            end
        end
    end

    return cells
end

# Given a cell and a starting edge, we follow the contour line until we either
# hit the boundary of the input data, or we form a closed contour.
function chase!(cells, curve, x, y, z, h, start, entry_edge, xi_range, yi_range, ::Type{VT}) where VT

    ind = start

    # When the contour loops back to the starting cell, it is possible
    # for it to not intersect with itself.  This happens if the starting
    # cell contains a saddle-point. So a loop is only closed if the
    # contour returns to the starting edge of the starting cell
    loopback_edge = entry_edge

    @inbounds while true
        exit_edge = get_next_edge!(cells, ind, entry_edge)

        push!(curve, interpolate(x, y, z, h, ind, exit_edge, VT))

        ind, entry_edge = advance_edge(ind, exit_edge)

        !((ind[1], ind[2], entry_edge) != (start[1], start[2], loopback_edge) &&
           ind[2] ∈ yi_range && ind[1] ∈ xi_range) && break
    end

    return ind
end

function trace_contour(x, y, z, h::Number, cells::Dict)

    contours = ContourLevel(h)

    x_ax, y_ax = axes(z)
    xi_range = first(x_ax):last(x_ax)-1
    yi_range = first(y_ax):last(y_ax)-1

    VT = SVector{2,promote_type(map(eltype, (x, y, z))...)}

    # When tracing out contours, this algorithm picks an arbitrary
    # starting cell, then first follows the contour in one direction
    # until it either ends up where it started # or at one of the boundaries.
    # It then tries to trace the contour in the opposite direction.

    @inbounds while length(cells) > 0
        contour_arr = VT[]

        # Pick initial box
        ind, cell = first(cells)

        # Pick a starting edge
        crossing = get_first_crossing(cell)
        starting_edge = 0x01 << trailing_zeros(crossing)

        # Add the contour entry location for cell (xi_0,yi_0)
        push!(contour_arr, interpolate(x, y, z, h, ind, starting_edge, VT))

        # Start trace in forward direction
        ind_end = chase!(cells, contour_arr, x, y, z, h, ind, starting_edge, xi_range, yi_range, VT)

        if ind == ind_end
            push!(contours.lines, Curve2(contour_arr))
            continue
        end

        ind, starting_edge = advance_edge(ind, starting_edge)

        if ind[2] ∈ yi_range && ind[1] ∈ xi_range
            # Start trace in reverse direction
            chase!(cells, reverse!(contour_arr), x, y, z, h, ind, starting_edge, xi_range, yi_range, VT)
        end

        push!(contours.lines, Curve2(contour_arr))
    end

    return contours

end

end


text:
function initialize!(stochasticprogram::StochasticProgram)
    initialize!(stochasticprogram, provided_optimizer(stochasticprogram))
end

function initialize!(stochasticprogram::StochasticProgram, ::UnrecognizedOptimizerProvided)
    # Error here if an unrecognized optimizer has been provided
    error("Trying to initialize stochasticprogram with unrecognized optimizer.")
end

function initialize!(stochasticprogram::StochasticProgram, ::NoOptimizerProvided)
    @warn "Cannot initialize without optimizer. Consider [`set_optimizer`](@ref)."
    # Do not initialize anything if no optimizer is attached
    return
end

function initialize!(stochasticprogram::StochasticProgram, ::OptimizerProvided)
    # Check that all generators are there
    _check_generators(stochasticprogram)
    # If the attached optimizer is an MOI optimizer, generate
    # the deterministic equivalent problem
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache, :dep)
        # Deterministic equivalent already generated
        stochasticprogram.optimizer.optimizer = backend(cache[:dep])
        return nothing
    end
    # Generate and cache deterministic equivalent
    cache[:dep] = generate_deterministic_equivalent(stochasticprogram)
    stochasticprogram.optimizer.optimizer = backend(cache[:dep])
    return nothing
end

function initialize!(stochasticprogram::StochasticProgram, ::StructuredOptimizerProvided)
    # Check that all generators are there
    _check_generators(stochasticprogram)
    # Generate block-structured stochastic program
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    # Initialize optimizer
    stochasticprogram.optimizer.optimizer = MOI.instantiate(optimizer_constructor(stochasticprogram))
    return nothing
end

ext: jl


text:
@static if VERSION >= v"1"
    using Distributed
end

function work(numsteps, np)
  sum = 0.0
  slice = 1.0/numsteps

  sum = @distributed (+) for i = 1:numsteps
    x = (i - 0.5) * slice
    (4.0/(1.0 + x^2))
  end

  mypi = sum * slice

  return mypi
end

function picalc(numsteps)

  # Get number of processes.
  np = nworkers()

  println("Calculating PI using:")
  println("  ", numsteps, " slices")
  println("  ", np, " worker(s)")
  temp = work(1,1)
  temp = work(numsteps,np)
  start = time()

  mypi = work(numsteps,np)

  stop = time()

  elapsed = (stop - start)

  println("Obtained value of PI: ", mypi)
  println("Time taken: ", elapsed, " seconds")

end

numsteps=1000000000

if length(ARGS) > 0
  numsteps = parse(Int, ARGS[1])
end

picalc(numsteps)


text:
using StatsPlots
using Distributions
using KernelDensity

include("../scripts/examples/ou.jl")
include("../src/plotting.jl")
#### WARNING: u,v,L has to coincide with the one in ../script/examples/ou.jl inside runall
T = 10.0
ν = 1.0
μ = -5.0
L = 6
u = -1.0
v = 2.0

runall()
function empirical_sample(dt, burn, u, v, L, T)
    X = runall(true)
    tt = burn:dt:(X[end].t-0.01)
    Ξ = []
    for i in tt
         push!(Ξ, FindCoordinates(X, i))
    end
    return Ξ
end

Ξ = empirical_sample(0.1, 10.0, u, v, L, T)

function midpoint_sample(Ξ, u, v, T)
    X = []
    for j in 1:length(Ξ)
        y = u/2 + v/2 +  Λ(T/2,0,0,T)* Ξ[j].ξ[1]
        push!(X, y)
    end
    return Array{Float64,1}(X)
end

midp = midpoint_sample(Ξ, u, v, T)
#a = histogram(midp)
x = kde(midp)
aa = plot(x.x,x.density, leg = false)

#### Real density
# dX_t = ν(μ - X_t) dt + dB_t
#step
a = -5.0
b = -1.0
T2 = T/2
function β(t , a, b, v, T)
    a + 2b*exp(b*(T-t))*(v + a/b*(1 - exp(b*(T-t))))/(exp(2b*(T-t))-1)
end

function int_B(t ,a, b, T2, step, T)
    som = 0
    for i in t:step:(T2 - step)
        som += step*(b - 2b*exp(2b*(T - (i + step/2)))/(exp(2b* (T - (i + step/2))) - 1))
    end
    return som
end

function μ_(a, b, T2, step, u, v, T)
    int1 = 0.0
    for t in 0:step:T2
        int1 += step*exp(int_B(t + step/2, a, b, T2, step, T))*β(t + step/2, a, b, v, T)
    end
    int1 += u*exp(int_B(0, a, b,  T2, step, T))
end

function σ_(T2, step1, step2, a, b, T)
    int1 = 0.0
    for t in 0:step1:T2
        int1 += step1*exp(2*int_B(t + step1/2, a, b, T2, step2, T))
    end
    return int1
end

#mean
μ = μ_(a, b, T2, 0.001, u, v, T)

#variance
σ = σ_(T2, 0.001, 0.001, a, b, T)

plot!(aa, Normal(μ, sqrt(σ)), label = "true distribution")
savefig("../../output/midpoint.png")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
mutable struct DQMCStack{
        GreensElType <: Number,
        HoppingElType <: Number,
        GreensMatType <: AbstractArray{GreensElType},
        HoppingMatType <: AbstractArray{HoppingElType},
        InteractionMatType <: AbstractArray
    } <: AbstractDQMCStack

    u_stack::Vector{GreensMatType}
    d_stack::Vector{Vector{Float64}}
    t_stack::Vector{GreensMatType}

    Ul::GreensMatType
    Ur::GreensMatType
    Dl::Vector{Float64}
    Dr::Vector{Float64}
    Tl::GreensMatType
    Tr::GreensMatType
    pivot::Vector{Int64}
    tempv::Vector{GreensElType}
    tempvf::Vector{Float64}

    greens::GreensMatType
    greens_temp::GreensMatType

    tmp1::GreensMatType
    tmp2::GreensMatType

    ranges::Array{UnitRange, 1}
    n_elements::Int
    current_slice::Int # running internally over 0:mc.parameters.slices+1, where 0 and mc.parameters.slices+1 are artifcial to prepare next sweep direction.
    direction::Int

    # # -------- Global update backup
    # gb_u_stack::Array{GreensElType, 3}
    # gb_d_stack::Matrix{Float64}
    # gb_t_stack::Array{GreensElType, 3}

    # gb_greens::GreensMatType
    # gb_log_det::Float64

    # gb_conf::Array{Float64, 3}
    # # --------

    # preallocated, reused arrays
    curr_U::GreensMatType
    eV::InteractionMatType

    # hopping matrices (mu included)
    hopping_matrix::HoppingMatType
    hopping_matrix_exp::HoppingMatType
    hopping_matrix_exp_inv::HoppingMatType
    hopping_matrix_exp_squared::HoppingMatType
    hopping_matrix_exp_inv_squared::HoppingMatType

    # checkerboard hopping matrices
    checkerboard::Matrix{Int} # src, trg, bondid
    groups::Vector{UnitRange}
    n_groups::Int
    chkr_hop_half::Vector{SparseMatrixCSC{HoppingElType, Int64}}
    chkr_hop_half_inv::Vector{SparseMatrixCSC{HoppingElType, Int64}}
    chkr_hop_half_dagger::Vector{SparseMatrixCSC{HoppingElType, Int64}}
    chkr_hop::Vector{SparseMatrixCSC{HoppingElType, Int64}} # without prefactor 0.5 in matrix exponentials
    chkr_hop_inv::Vector{SparseMatrixCSC{HoppingElType, Int64}}
    chkr_hop_dagger::Vector{SparseMatrixCSC{HoppingElType, Int64}}
    chkr_mu_half::SparseMatrixCSC{HoppingElType, Int64}
    chkr_mu_half_inv::SparseMatrixCSC{HoppingElType, Int64}
    chkr_mu::SparseMatrixCSC{HoppingElType, Int64}
    chkr_mu_inv::SparseMatrixCSC{HoppingElType, Int64}

    function DQMCStack{GET, HET, GMT, HMT, IMT}() where {
            GET<:Number, HET<:Number,
            GMT<:AbstractArray{GET}, HMT<:AbstractArray{HET}, IMT<:AbstractArray
        }
        @assert isconcretetype(GET);
        @assert isconcretetype(HET);
        @assert isconcretetype(GMT);
        @assert isconcretetype(HMT);
        @assert isconcretetype(IMT);
        @assert eltype(GMT) == GET;
        @assert eltype(HMT) == HET;
        new{GET, HET, GMT, HMT, IMT}()
    end
end

# type helpers
geltype(::DQMCStack{GET, HET, GMT, HMT, IMT}) where {GET, HET, GMT, HMT, IMT} = GET
heltype(::DQMCStack{GET, HET, GMT, HMT, IMT}) where {GET, HET, GMT, HMT, IMT} = HET
gmattype(::DQMCStack{GET, HET, GMT, HMT, IMT}) where {GET, HET, GMT, HMT, IMT} = GMT
hmattype(::DQMCStack{GET, HET, GMT, HMT, IMT}) where {GET, HET, GMT, HMT, IMT} = HMT
imattype(::DQMCStack{GET, HET, GMT, HMT, IMT}) where {GET, HET, GMT, HMT, IMT} = IMT

geltype(mc::DQMC) = geltype(mc.stack)
heltype(mc::DQMC) = heltype(mc.stack)
gmattype(mc::DQMC) = gmattype(mc.stack)
hmattype(mc::DQMC) = hmattype(mc.stack)
imattype(mc::DQMC) = imattype(mc.stack)

################################################################################
### Stack Initialization
################################################################################

function initialize_stack(mc::DQMC, ::DQMCStack)
    GreensElType = geltype(mc)
    GreensMatType = gmattype(mc)
    HoppingElType = heltype(mc)
    N = length(lattice(mc))
    flv = nflavors(mc.model)

    mc.stack.n_elements = convert(Int, mc.parameters.slices / mc.parameters.safe_mult) + 1

    mc.stack.u_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:mc.stack.n_elements]
    mc.stack.d_stack = [zeros(Float64, flv*N) for _ in 1:mc.stack.n_elements]
    mc.stack.t_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:mc.stack.n_elements]

    mc.stack.greens = GreensMatType(undef, flv*N, flv*N)
    mc.stack.greens_temp = GreensMatType(undef, flv*N, flv*N)

    # used in calculate_greens
    # do not change in slice_matrices.jl or interaction_matrix_exp!
    mc.stack.Ul = GreensMatType(I, flv*N, flv*N)
    mc.stack.Ur = GreensMatType(I, flv*N, flv*N)
    mc.stack.Tl = GreensMatType(I, flv*N, flv*N)
    mc.stack.Tr = GreensMatType(I, flv*N, flv*N)
    mc.stack.Dl = ones(Float64, flv*N)
    mc.stack.Dr = ones(Float64, flv*N)
    # can be changed anywhere
    mc.stack.pivot = Vector{Int64}(undef, flv*N)
    mc.stack.tempv = Vector{GreensElType}(undef, flv*N)
    mc.stack.tempvf = Vector{Float64}(undef, flv*N)

    # can be changed anywhere
    mc.stack.tmp1 = GreensMatType(undef, flv*N, flv*N)
    mc.stack.tmp2 = GreensMatType(undef, flv*N, flv*N)

    # # Global update backup
    # mc.stack.gb_u_stack = zero(mc.stack.u_stack)
    # mc.stack.gb_d_stack = zero(mc.stack.d_stack)
    # mc.stack.gb_t_stack = zero(mc.stack.t_stack)
    # mc.stack.gb_greens = zero(mc.stack.greens)
    # mc.stack.gb_log_det = 0.
    # mc.stack.gb_conf = zero(mc.conf)

    mc.stack.ranges = UnitRange[]

    for i in 1:mc.stack.n_elements - 1
        push!(mc.stack.ranges, 1 + (i - 1) * mc.parameters.safe_mult:i * mc.parameters.safe_mult)
    end

    mc.stack.curr_U = GreensMatType(undef, flv*N, flv*N)
    mc.stack.eV = init_interaction_matrix(mc.model)

    nothing
end

# hopping
function init_hopping_matrices(mc::DQMC{M,CB}, m::Model) where {M, CB<:Checkerboard}
    init_hopping_matrix_exp(mc, m)
    CB <: CheckerboardTrue && init_checkerboard_matrices(mc, m)
    nothing
end
function init_hopping_matrix_exp(mc::DQMC, m::Model)
    N = length(lattice(m))
    flv = nflavors(m)
    dtau = mc.parameters.delta_tau

    T = hopping_matrix(mc, m)
    size(T) == (flv*N, flv*N) || error("Hopping matrix should have size "*
                                "$((flv*N, flv*N)) but has size $(size(T)) .")
    mc.stack.hopping_matrix = T
    mc.stack.hopping_matrix_exp = exp(-0.5 * dtau * T)
    mc.stack.hopping_matrix_exp_inv = exp(0.5 * dtau * T)
    mc.stack.hopping_matrix_exp_squared = mc.stack.hopping_matrix_exp * mc.stack.hopping_matrix_exp
    mc.stack.hopping_matrix_exp_inv_squared = mc.stack.hopping_matrix_exp_inv * mc.stack.hopping_matrix_exp_inv
    nothing
end

# checkerboard
rem_eff_zeros!(X::AbstractArray) = map!(e -> abs.(e)<1e-15 ? zero(e) : e,X,X)
function init_checkerboard_matrices(mc::DQMC, m::Model)
    s = mc.stack
    l = lattice(m)
    flv = nflavors(m)
    H = heltype(mc)
    N = length(l)
    dtau = mc.parameters.delta_tau
    mu = m.mu

    s.checkerboard, s.groups, s.n_groups = build_checkerboard(l)
    n_grps = s.n_groups
    cb = s.checkerboard

    T = reshape(hopping_matrix(mc, m), (N, flv, N, flv))

    s.chkr_hop_half = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
    s.chkr_hop_half_inv = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
    s.chkr_hop = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)
    s.chkr_hop_inv = Vector{SparseMatrixCSC{H, Int}}(undef, n_grps)

    for (g, gr) in enumerate(s.groups)
        Tg = zeros(H, N, flv, N, flv)
        for i in gr
            src, trg = cb[1:2,i]
            for f1 in 1:flv, f2 in 1:flv
                Tg[trg, f1, src, f2] = T[trg, f1, src, f2]
            end
        end

        Tgg = reshape(Tg, (N*flv, N*flv))
        s.chkr_hop_half[g] = sparse(rem_eff_zeros!(exp(-0.5 * dtau * Tgg)))
        s.chkr_hop_half_inv[g] = sparse(rem_eff_zeros!(exp(0.5 * dtau * Tgg)))
        s.chkr_hop[g] = sparse(rem_eff_zeros!(exp(- dtau * Tgg)))
        s.chkr_hop_inv[g] = sparse(rem_eff_zeros!(exp(dtau * Tgg)))
    end

    s.chkr_hop_half_dagger = adjoint.(s.chkr_hop_half)
    s.chkr_hop_dagger = adjoint.(s.chkr_hop)

    mus = diag(reshape(T, (N*flv, N*flv)))
    s.chkr_mu_half = spdiagm(0 => exp.(-0.5 * dtau * mus))
    s.chkr_mu_half_inv = spdiagm(0 => exp.(0.5 * dtau * mus))
    s.chkr_mu = spdiagm(0 => exp.(-dtau * mus))
    s.chkr_mu_inv = spdiagm(0 => exp.(dtau * mus))

    # hop_mat_exp_chkr = foldl(*,s.chkr_hop_half) * sqrt.(s.chkr_mu)
    # r = effreldiff(s.hopping_matrix_exp,hop_mat_exp_chkr)
    # r[find(x->x==zero(x),hop_mat_exp_chkr)] = 0.
    # println("Checkerboard - Exact ≈ ", round(maximum(absdiff(s.hopping_matrix_exp,hop_mat_exp_chkr)), 4))
    nothing
end

"""
    build_stack(mc::DQMC)

Build slice matrix stack from scratch.
"""
@bm function build_stack(mc::DQMC, ::DQMCStack)
    copyto!(mc.stack.u_stack[1], I)
    mc.stack.d_stack[1] .= one(eltype(mc.stack.d_stack[1]))
    copyto!(mc.stack.t_stack[1], I)

    @inbounds for i in 1:length(mc.stack.ranges)
        add_slice_sequence_left(mc, i)
    end

    mc.stack.current_slice = mc.parameters.slices + 1
    mc.stack.direction = -1

    nothing
end

@bm function reverse_build_stack(mc::DQMC, ::DQMCStack)
    copyto!(mc.stack.u_stack[end], I)
    mc.stack.d_stack[end] .= one(eltype(mc.stack.d_stack[end]))
    copyto!(mc.stack.t_stack[end], I)

    @inbounds for i in length(mc.stack.ranges):-1:1
        add_slice_sequence_right(mc, i)
    end

    mc.stack.current_slice = 0
    mc.stack.direction = 1

    nothing
end

################################################################################
### Slice matrix stack manipulations/updates
################################################################################

"""
    add_slice_sequence_left(mc::DQMC, idx)

Computes the next `mc.parameters.safe_mult` slice matrix products from the current `idx`
and writes them to `idx+1`. The index `idx` does not refer to the slice index,
but `mc.parameters.safe_mult` times the slice index.
"""
@bm function add_slice_sequence_left(mc::DQMC, idx::Int)
    @inbounds begin
        copyto!(mc.stack.curr_U, mc.stack.u_stack[idx])

        # println("Adding slice seq left $idx = ", mc.stack.ranges[idx])
        for slice in mc.stack.ranges[idx]
            multiply_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
        end

        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.d_stack[idx]))
        udt_AVX_pivot!(
            mc.stack.u_stack[idx + 1], mc.stack.d_stack[idx + 1], mc.stack.tmp1,
            mc.stack.pivot, mc.stack.tempv
        )
        vmul!(mc.stack.t_stack[idx + 1], mc.stack.tmp1, mc.stack.t_stack[idx])
    end
end

"""
    add_slice_sequence_right(mc::DQMC, idx)

Computes the next `mc.parameters.safe_mult` slice matrix products from the current
`idx+1` and writes them to `idx`. The index `idx` does not refer to the slice
index, but `mc.parameters.safe_mult` times the slice index.
"""
@bm function add_slice_sequence_right(mc::DQMC, idx::Int)
    @inbounds begin
        copyto!(mc.stack.curr_U, mc.stack.u_stack[idx + 1])

        for slice in reverse(mc.stack.ranges[idx])
            multiply_daggered_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
        end

        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.d_stack[idx + 1]))
        udt_AVX_pivot!(
            mc.stack.u_stack[idx], mc.stack.d_stack[idx], mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv
        )
        vmul!(mc.stack.t_stack[idx], mc.stack.tmp1, mc.stack.t_stack[idx + 1])
    end
end

################################################################################
### Green's function calculation
################################################################################

"""
    calculate_greens_AVX!(Ul, Dl, Tl, Ur, Dr, Tr, G[, pivot, temp])

Calculates the effective Greens function matrix `G` from two UDT decompositions
`Ul, Dl, Tl` and `Ur, Dr, Tr`. Additionally a `pivot` vector can be given. Note
that all inputs will be overwritten.

The UDT should follow from a set of slice_matrix multiplications, such that
`Ur, Dr, Tr = udt(B(slice)' ⋯ B(M)')` and `Ul, Dl, Tl = udt(B(slice-1) ⋯ B(1))`.
The computed Greens function is then given as `G = inv(I + Ul Dl Tl Tr Dr Ur)`
and computed here.

`Ul, Tl, Ur, Tr, G` should be square matrices, `Dl, Dr` real Vectors (from
Diagonal matrices), `pivot` an integer Vector and `temp` a Vector with the same
element type as the matrices.
"""
@bm function calculate_greens_AVX!(
        Ul, Dl, Tl, Ur, Dr, Tr, G::AbstractArray{T},
        pivot = Vector{Int64}(undef, length(Dl)),
        temp = Vector{T}(undef, length(Dl))
    ) where T
    # @bm "B1" begin
        # Used: Ul, Dl, Tl, Ur, Dr, Tr
        # TODO: [I + Ul Dl Tl Tr^† Dr Ur^†]^-1
        # Compute: Dl * ((Tl * Tr) * Dr) -> Tr * Dr * G   (UDT)
        vmul!(G, Tl, adjoint(Tr))
        vmul!(Tr, G, Diagonal(Dr))
        vmul!(G, Diagonal(Dl), Tr)
        udt_AVX_pivot!(Tr, Dr, G, pivot, temp, Val(false)) # Dl available
    # end

    # @bm "B2" begin
        # Used: Ul, Ur, G, Tr, Dr  (Ul, Ur, Tr unitary (inv = adjoint))
        # TODO: [I + Ul Tr Dr G Ur^†]^-1
        #     = [(Ul Tr) ((Ul Tr)^-1 (G Ur^†) + Dr) (G Ur)]^-1
        #     = Ur G^-1 [(Ul Tr)^† Ur G^-1 + Dr]^-1 (Ul Tr)^†
        # Compute: Ul Tr -> Tl
        #          (Ur G^-1) -> Ur
        #          ((Ul Tr)^† Ur G^-1) -> Tr
        vmul!(Tl, Ul, Tr)
        rdivp!(Ur, G, Ul, pivot) # requires unpivoted udt decompostion (Val(false))
        vmul!(Tr, adjoint(Tl), Ur)
    # end

    # @bm "B3" begin
        # Used: Tl, Ur, Tr, Dr
        # TODO: Ur [Tr + Dr]^-1 Tl^† -> Ur [Tr]^-1 Tl^†
        rvadd!(Tr, Diagonal(Dr))
    # end

    # @bm "B4" begin
        # Used: Ur, Tr, Tl
        # TODO: Ur [Tr]^-1 Tl^† -> Ur [Ul Dr Tr]^-1 Tl^†
        #    -> Ur Tr^-1 Dr^-1 Ul^† Tl^† -> Ur Tr^-1 Dr^-1 (Tl Ul)^†
        # Compute: Ur Tr^-1 -> Ur,  Tl Ul -> Tr
        udt_AVX_pivot!(Ul, Dr, Tr, pivot, temp, Val(false)) # Dl available
        rdivp!(Ur, Tr, G, pivot) # requires unpivoted udt decompostion (false)
        vmul!(Tr, Tl, Ul)
    # end

    # @bm "B5" begin
        @turbo for i in eachindex(Dr)
            Dl[i] = 1.0 / Dr[i]
        end
    # end

    # @bm "B6" begin
        # Used: Ur, Tr, Dl, Ul, Tl
        # TODO: (Ur Dl) Tr^† -> G
        vmul!(Ul, Ur, Diagonal(Dl))
        vmul!(G, Ul, adjoint(Tr))
    # end
end

"""
    calculate_greens(mc::DQMC)

Computes the effective greens function from the current state of the stack and
saves the result to `mc.stack.greens`.

This assumes the `mc.stack.Ul, mc.stack.Dl, mc.stack.Tl = udt(B(slice-1) ⋯ B(1))` and
`mc.stack.Ur, mc.stack.Dr, mc.stack.Tr = udt(B(slice)' ⋯ B(M)')`.

This should only used internally.
"""
@bm function calculate_greens(mc::DQMC, output::AbstractMatrix = mc.stack.greens)
    calculate_greens_AVX!(
        mc.stack.Ul, mc.stack.Dl, mc.stack.Tl,
        mc.stack.Ur, mc.stack.Dr, mc.stack.Tr,
        output, mc.stack.pivot, mc.stack.tempv
    )
    output
end

"""
    calculate_greens(mc::DQMC, slice[, output=mc.stack.greens, safe_mult])

Compute the effective equal-time greens function from scratch at a given `slice`.

This does not invalidate the stack, but it does overwrite `mc.stack.greens`.
"""
@bm function calculate_greens(
        mc::DQMC, slice::Int, output::AbstractMatrix = mc.stack.greens,
        conf::AbstractArray = mc.conf, safe_mult::Int = mc.parameters.safe_mult
    )
    copyto!(mc.stack.curr_U, I)
    copyto!(mc.stack.Ur, I)
    mc.stack.Dr .= one(eltype(mc.stack.Dr))
    copyto!(mc.stack.Tr, I)

    # Calculate Ur,Dr,Tr=B(slice)' ... B(M)'
    if slice+1 <= mc.parameters.slices
        start = slice+1
        stop = mc.parameters.slices
        for k in reverse(start:stop)
            if mod(k,safe_mult) == 0
                multiply_daggered_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, conf)
                vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dr))
                udt_AVX_pivot!(mc.stack.curr_U, mc.stack.Dr, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
                copyto!(mc.stack.tmp2, mc.stack.Tr)
                vmul!(mc.stack.Tr, mc.stack.tmp1, mc.stack.tmp2)
            else
                multiply_daggered_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, conf)
            end
        end
        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dr))
        udt_AVX_pivot!(mc.stack.Ur, mc.stack.Dr, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
        copyto!(mc.stack.tmp2, mc.stack.Tr)
        vmul!(mc.stack.Tr, mc.stack.tmp1, mc.stack.tmp2)
    end

    copyto!(mc.stack.curr_U, I)
    copyto!(mc.stack.Ul, I)
    mc.stack.Dl .= one(eltype(mc.stack.Dl))
    copyto!(mc.stack.Tl, I)

    # Calculate Ul,Dl,Tl=B(slice-1) ... B(1)
    if slice >= 1
        start = 1
        stop = slice
        for k in start:stop
            if mod(k,safe_mult) == 0
                multiply_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, conf)
                vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dl))
                udt_AVX_pivot!(mc.stack.curr_U, mc.stack.Dl, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
                copyto!(mc.stack.tmp2, mc.stack.Tl)
                vmul!(mc.stack.Tl, mc.stack.tmp1, mc.stack.tmp2)
            else
                multiply_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, conf)
            end
        end
        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dl))
        udt_AVX_pivot!(mc.stack.Ul, mc.stack.Dl, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
        copyto!(mc.stack.tmp2, mc.stack.Tl)
        vmul!(mc.stack.Tl, mc.stack.tmp1, mc.stack.tmp2)
    end

    return calculate_greens(mc, output)
end

################################################################################
### Stack update
################################################################################

# Green's function propagation
@inline @bm function wrap_greens!(mc::DQMC, gf, curr_slice::Int, direction::Int)
    if direction == -1
        multiply_slice_matrix_inv_left!(mc, mc.model, curr_slice - 1, gf)
        multiply_slice_matrix_right!(mc, mc.model, curr_slice - 1, gf)
    else
        multiply_slice_matrix_left!(mc, mc.model, curr_slice, gf)
        multiply_slice_matrix_inv_right!(mc, mc.model, curr_slice, gf)
    end
    nothing
end

@bm function propagate(mc::DQMC)
    @inbounds if mc.stack.direction == 1
        if mod(mc.stack.current_slice, mc.parameters.safe_mult) == 0
            mc.stack.current_slice +=1 # slice we are going to
            if mc.stack.current_slice == 1
                copyto!(mc.stack.Ur, mc.stack.u_stack[1])
                copyto!(mc.stack.Dr, mc.stack.d_stack[1])
                copyto!(mc.stack.Tr, mc.stack.t_stack[1])
                copyto!(mc.stack.u_stack[1], I)
                mc.stack.d_stack[1] .= one(eltype(mc.stack.d_stack[1]))
                copyto!(mc.stack.t_stack[1], I)
                copyto!(mc.stack.Ul, mc.stack.u_stack[1])
                copyto!(mc.stack.Dl, mc.stack.d_stack[1])
                copyto!(mc.stack.Tl, mc.stack.t_stack[1])

                calculate_greens(mc) # greens_1 ( === greens_{m+1} )

            elseif 1 < mc.stack.current_slice <= mc.parameters.slices
                idx = Int((mc.stack.current_slice - 1)/mc.parameters.safe_mult)

                copyto!(mc.stack.Ur, mc.stack.u_stack[idx+1])
                copyto!(mc.stack.Dr, mc.stack.d_stack[idx+1])
                copyto!(mc.stack.Tr, mc.stack.t_stack[idx+1])
                add_slice_sequence_left(mc, idx)
                copyto!(mc.stack.Ul, mc.stack.u_stack[idx+1])
                copyto!(mc.stack.Dl, mc.stack.d_stack[idx+1])
                copyto!(mc.stack.Tl, mc.stack.t_stack[idx+1])

                if mc.parameters.check_propagation_error
                    copyto!(mc.stack.greens_temp, mc.stack.greens)
                end

                # Should this be mc.stack.greens_temp?
                # If so, shouldn't this only run w/ mc.parameters.all_checks = true?
                wrap_greens!(mc, mc.stack.greens_temp, mc.stack.current_slice - 1, 1)

                calculate_greens(mc) # greens_{slice we are propagating to}

                if mc.parameters.check_propagation_error
                    # OPT: could probably be optimized through explicit loop
                    greensdiff = maximum(abs.(mc.stack.greens_temp - mc.stack.greens))
                    if greensdiff > 1e-7
                        push!(mc.analysis.propagation_error, greensdiff)
                        mc.parameters.silent || @printf(
                            "->%d \t+1 Propagation instability\t %.1e\n",
                            mc.stack.current_slice, greensdiff
                        )
                    end
                end

            else # we are going to mc.parameters.slices+1
                idx = mc.stack.n_elements - 1
                add_slice_sequence_left(mc, idx)
                mc.stack.direction = -1
                mc.stack.current_slice = mc.parameters.slices+1 # redundant
                propagate(mc)
            end

        else
            # Wrapping
            wrap_greens!(mc, mc.stack.greens, mc.stack.current_slice, 1)
            mc.stack.current_slice += 1
        end

    else # mc.stack.direction == -1
        if mod(mc.stack.current_slice-1, mc.parameters.safe_mult) == 0
            mc.stack.current_slice -= 1 # slice we are going to
            if mc.stack.current_slice == mc.parameters.slices
                copyto!(mc.stack.Ul, mc.stack.u_stack[end])
                copyto!(mc.stack.Dl, mc.stack.d_stack[end])
                copyto!(mc.stack.Tl, mc.stack.t_stack[end])
                copyto!(mc.stack.u_stack[end], I)
                mc.stack.d_stack[end] .= one(eltype(mc.stack.d_stack[end]))
                copyto!(mc.stack.t_stack[end], I)
                copyto!(mc.stack.Ur, mc.stack.u_stack[end])
                copyto!(mc.stack.Dr, mc.stack.d_stack[end])
                copyto!(mc.stack.Tr, mc.stack.t_stack[end])

                calculate_greens(mc) # greens_{mc.parameters.slices+1} === greens_1

                # wrap to greens_{mc.parameters.slices}
                wrap_greens!(mc, mc.stack.greens, mc.stack.current_slice + 1, -1)

            elseif 0 < mc.stack.current_slice < mc.parameters.slices
                idx = Int(mc.stack.current_slice / mc.parameters.safe_mult) + 1

                copyto!(mc.stack.Ul, mc.stack.u_stack[idx])
                copyto!(mc.stack.Dl, mc.stack.d_stack[idx])
                copyto!(mc.stack.Tl, mc.stack.t_stack[idx])
                add_slice_sequence_right(mc, idx)
                copyto!(mc.stack.Ur, mc.stack.u_stack[idx])
                copyto!(mc.stack.Dr, mc.stack.d_stack[idx])
                copyto!(mc.stack.Tr, mc.stack.t_stack[idx])

                if mc.parameters.check_propagation_error
                    copyto!(mc.stack.greens_temp, mc.stack.greens)
                end

                calculate_greens(mc)

                if mc.parameters.check_propagation_error
                    # OPT: could probably be optimized through explicit loop
                    greensdiff = maximum(abs.(mc.stack.greens_temp - mc.stack.greens))
                    if greensdiff > 1e-7
                        push!(mc.analysis.propagation_error, greensdiff)
                        mc.parameters.silent || @printf(
                            "->%d \t-1 Propagation instability\t %.1e\n",
                            mc.stack.current_slice, greensdiff
                        )
                    end
                end

                wrap_greens!(mc, mc.stack.greens, mc.stack.current_slice + 1, -1)

            else # we are going to 0
                idx = 1
                add_slice_sequence_right(mc, idx)
                mc.stack.direction = 1
                mc.stack.current_slice = 0 # redundant
                propagate(mc)
            end

        else
            # Wrapping
            wrap_greens!(mc, mc.stack.greens, mc.stack.current_slice, -1)
            mc.stack.current_slice -= 1
        end
    end
    nothing
end


text:
# YearFraction API

@inline value(yf::YearFraction) = yf.val
yearfraction(conv::DayCountConvention, day_count::Number) = YearFraction(day_count / daysperyear(conv))
yearfraction(conv::DayCountConvention, date_start::Date, date_end::Date) = yearfraction(conv, daycount(conv, date_start, date_end))
yearfraction(curve::AbstractIRCurve, maturity::Date) = yearfraction(curve_get_daycount(curve), curve_get_date(curve), maturity)

yearfractionvalue(yf::YearFraction) = value(yf)
yearfractionvalue(args...) = value(yearfraction(args...)) # old yearfraction behavior

function Base.zero(yf::YearFraction{T}) where {T}
    YearFraction(zero(T))
end

Base.:(==)(y1::YearFraction, y2::YearFraction) = value(y1) == value(y2)
Base.isless(y1::YearFraction, y2::YearFraction) = value(y1) < value(y2)

# DayCountConvention

Base.:(==)(d1::BDays252, d2::BDays252) = d1.hc == d2.hc
Base.hash(d::BDays252) = 1 + hash(d.hc)

daycount(conv::BDays252, date_start::Date, date_end::Date) = BusinessDays.bdayscount(conv.hc, date_start, date_end)
daycount(::Actual360, date_start::Date, date_end::Date) = Int(Dates.value(date_end - date_start))
daycount(::Actual365, date_start::Date, date_end::Date) = Int(Dates.value(date_end - date_start))
function daycount(::Thirty360, date_start::Date, date_end::Date)
    Y = Dates.year(date_end) - Dates.year(date_start)
    M = Dates.month(date_end) - Dates.month(date_start)

    d1 = Dates.day(date_start)
    d2 = Dates.day(date_end)
    if d1 >= 30
        d1 = 30
        if d2 >= 30
            d2 = 30
        end
    end
    D = d2-d1

    return 360Y + 30M + D
end

advancedays(conv::BDays252, date_start::Date, daycount::Int) = BusinessDays.advancebdays(conv.hc, date_start, daycount)
advancedays(::Actual360, date_start::Date, daycount::Int) = date_start + Day(daycount)
advancedays(::Actual365, date_start::Date, daycount::Int) = date_start + Day(daycount)
advancedays(::Thirty360, date_start::Date, daycount::Int) = date_start + Day(daycount)

function advancedays(conv::DayCountConvention, date_start::Date, daycount_vec::Vector{Int})
    l = length(daycount_vec)
    result = Vector{Date}(undef, l)
    for i in 1:l
        @inbounds result[i] = advancedays(conv, date_start, daycount_vec[i])
    end
    return result
end

advancedays(curve::AbstractIRCurve, daycount) = advancedays(curve_get_daycount(curve), curve_get_date(curve), daycount)

daysperyear(::BDays252) = 252
daysperyear(::Actual360) = 360
daysperyear(::Actual365) = 365
daysperyear(::Thirty360) = 360

function days_to_maturity(curve::AbstractIRCurve, maturity::Date)
    d = daycount(curve_get_daycount(curve), curve_get_date(curve), maturity)
    @assert d >= 0 "Maturity date $(maturity) should be greater than curve observation date $(curve_get_date(curve))"
    return d
end

function days_to_maturity(curve::AbstractIRCurve, maturity::YearFraction)
    return value(maturity) * daysperyear(curve_get_daycount(curve))
end


text:
@enum MoveCategory begin
    PHYSICAL
    SPECIAL
    STATUS
end

@enum NonstandardMoveType begin
    NONSTANDARD_NONE
    NONSTANDARD_PAST
    NONSTANDARD_GIGANTAMAX
    NONSTANDARD_LGPE
end
VsRecorderBase.enum_prefix(::Type{NonstandardMoveType}) = "NONSTANDARD_"

@enum MoveTarget begin
    TARGET_NORMAL
    TARGET_ADJACENT_FOE
    TARGET_ANY
    TARGET_ADJACENT_ALLY
    TARGET_FOE_SIDE
    TARGET_RANDOM_NORMAL
    TARGET_ADJACENT_ALLY_OR_SELF
    TARGET_ALLY_TEAM
    TARGET_SCRIPTED
    TARGET_ALL_ADJACENT
    TARGET_SELF
    TARGET_ALL
    TARGET_ALLIES
    TARGET_ALLY_SIDE
    TARGET_ALL_ADJACENT_FOES
end
VsRecorderBase.enum_prefix(::Type{MoveTarget}) = "TARGET_"

Base.@kwdef struct DexStats
    hp::Int
    atk::Int
    def::Int
    spa::Int
    spd::Int
    spe::Int
end

Base.@kwdef struct DexPokemon
    num::Int64
    name::String
    types::Vector{PokemonType}
    weightkg::Float64
    heightm::Float64
    abilities::Vector{AbilityID}
    base_stats::DexStats

    forme::Nullable{String} = nothing
    base_forme::Nullable{String} = nothing
    base_species::Nullable{String} = nothing
    other_formes::Nullable{Vector{String}} = nothing

    other_data::Dict{Symbol, Any} = Dict()
end

Base.@kwdef struct DexMove
    num::Int64
    accuracy::Union{Int64, Bool}
    base_power::Int64
    category::MoveCategory
    is_nonstandard::NonstandardMoveType = NONSTANDARD_NONE
    name::String
    pp::Int64
    priority::Int64
    target::MoveTarget = TARGET_NORMAL
    type::PokemonType

    other_data::Dict{Symbol, Any} = Dict()
end

const PokeDexFile = datapath("dex", "pokedex.yaml")
const MoveDexFile = datapath("dex", "movedex.yaml")

macro _dex_func(func_name, dex_name, id_type, dex_type)
    file_name = Symbol(dex_name, :File)
    T = :(OrderedDict{$id_type, $dex_type})
    quote
        const $dex_name = Ref{Nullable{$T}}(nothing)
        function $func_name(; force = false)
            @_load_data $dex_name force open($file_name) do fi
                deserialize(
                    fi,
                    $T,
                    other_key = :other_data,
                    dicttype = OrderedDict,
                    loading_dex = true
                )
            end
        end
    end |> esc
end

@_dex_func move_dex MoveDex MoveID DexMove
@_dex_func poke_dex PokeDex PokemonID DexPokemon

struct SimpleListDex{T}
    data::SimpleI18n.I18nData
end
@forward SimpleListDex.data Base.keys, Base.get, Base.haskey, Base.getindex, Base.setindex!
macro _list_dex_func(func_name, id_type, i18n_key)
    i18n_key = string(i18n_key)
    quote
        function $func_name(language = get_game_language())
            ctx = VsI18n.GlobalI18nContext[]
            lang = VsI18n.get_code(get_game_language(language))
            SimpleListDex{$id_type}(ctx.data[lang][$i18n_key])
        end
        $func_name(conf::VsConfig) = $func_name(conf.source.language)
        $func_name(ctx::PokemonContext) = $func_name(ctx.config)
        VsRecorderBase._parse(x, ::Type{$id_type}; loading_dex = false, kwargs...) = if loading_dex
            search_dex($func_name(VsI18n.EN), x) |> $id_type
        else
            $id_type(x)
        end
    end |> esc
end

@_list_dex_func nature_list Nature nature
@_list_dex_func ability_list AbilityID ability
@_list_dex_func item_list ItemID item
@_list_dex_func move_list MoveID move
@_list_dex_func poke_list PokemonID pokemon

function search_dex(dex::AbstractDict, name)
    for (id, value) in dex
        if value.name == name
            return id
        end
    end
end
function search_dex(dex::Union{SimpleListDex{T}, AbstractDict{T, String}}, name) where T
    for id in keys(dex)
        if dex[id] == name
            return T(id)
        end
    end
    T(name)
end
Base.getindex(data::SimpleI18n.I18nData, key::VsAbstractID) = data[string(key)]
Base.getindex(data::AbstractDict{T}, key::AbstractString) where T <: VsAbstractID = data[T(key)]
Base.getindex(data::SimpleI18n.I18nData, key::EnumIDTypes) = data[enum_to_string(key)]

VsRecorderBase._parse(
    dict::OrderedDict, ::Type{Vector{AbilityID}};
    _...
) = AbilityID[search_dex(ability_list(VsI18n.EN), x) for x in values(dict)]

function initialize_dex()
    move_dex()
    poke_dex()
    nothing
end


text:
# DiagnosticOptions
DiagnosticOptions() = DiagnosticOptions(create_diagnostic_opts())

"""
    create_diagnostic_opts() -> CXDiagnosticOptions
Return a pointer to a `clang::DiagnosticOptions` object.
"""
function create_diagnostic_opts()
    status = Ref{CXInit_Error}(CXInit_NoError)
    opts = clang_DiagnosticOptions_create(status)
    @assert status[] == CXInit_NoError
    return opts
end

function PrintStats(x::DiagnosticOptions)
    @check_ptrs x
    return clang_DiagnosticOptions_PrintStats(x)
end

function setShowColors(x::DiagnosticOptions, should_show::Bool)
    @check_ptrs x
    return clang_DiagnosticOptions_setShowColors(x, should_show)
end

function setShowPresumedLoc(x::DiagnosticOptions, should_show::Bool)
    @check_ptrs x
    return clang_DiagnosticOptions_setShowPresumedLoc(x, should_show)
end


text:
function get_cartpole(;
    timestep=0.1,
    gravity=[0.0; 0.0; -9.81],
    spring=0.0,
    damper=0.0,
    T=Float64)

    # Parameters
    slider_axis = [0.0; 1.0; 0.0]
    pendulum_axis = [1.0; 0.0; 0.0]
    slider_length = 1.0
    pendulum_length = 1.0
    radius = 0.075
    slider_mass = 1.0
    pendulum_mass = 0.432

    # Links
    origin = Origin{Float64}()
    slider = Capsule(1.5 * radius, slider_length, slider_mass,
        axis_offset=RotX(0.5 * π),
        color=RGBA(0.7, 0.7, 0.7, 1.0))
    pendulum = Capsule(radius, pendulum_length, pendulum_mass,
        color=RGBA(0.7, 0.7, 0.7, 1.0))
    links = [slider, pendulum]

    # Joint Constraints
    joint_origin_slider = JointConstraint(Prismatic(origin, slider, slider_axis;
        parent_vertex=szeros(Float64, 3),
        child_vertex=szeros(Float64, 3)))
    joint_slider_pendulum = JointConstraint(Revolute(slider, pendulum, pendulum_axis;
        parent_vertex=szeros(Float64, 3),
        child_vertex=[0.0; 0.0; 0.5 * pendulum_length]))
    joints = [joint_origin_slider, joint_slider_pendulum]

    # Mechanism
    mech = Mechanism(origin, links, joints,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    return mech
end

function initialize_cartpole!(mech::Mechanism{T,Nn,Ne,Nb};
        mode=:down,
        pendulum_length=1.0) where {T,Nn,Ne,Nb}

    # origin to slider
    set_maximal_configurations!(mech.origin, mech.bodies[1])
    set_maximal_velocities!(mech.bodies[1],
        v=[0.0; 0.0; 0.0],
        ω=zeros(3))

    # slider to pendulum
    if mode == :down
        set_maximal_configurations!(mech.bodies[1], mech.bodies[2],
            Δx=[0.0; 0.0; -0.5 * pendulum_length],
            Δq=RotX(π))
        set_maximal_velocities!(mech.bodies[2],
            v=zeros(3),
            ω=zeros(3))
    elseif mode == :up
        set_maximal_configurations!(mech.bodies[1], mech.bodies[2],
            Δx=[0.0; 0.0; 0.5 * pendulum_length],
            Δq=RotX(π))
        set_maximal_velocities!(mech.bodies[2],
            v=zeros(3),
            ω=zeros(3))
    end
end

function mujoco_inertia!(mech)
    mech.bodies[1].m = 1.0063
    mech.bodies[1].J = Diagonal([0.106974, 0.106974, 0.00636812])

    mech.bodies[2].m = 0.4321
    mech.bodies[2].J = Diagonal([0.0422274, 0.0422274, 0.0012155])
end


text:
#= ===================
   LOGISTIC CLASSIFIER
   =================== =#

"""
$SIGNATURES

Logistic Classifier (typically called "Logistic Regression"). This model is
a standard classifier for both binary and multiclass classification.
In the binary case it corresponds to the LogisticLoss, in the multiclass to the
Multinomial (softmax) loss. An elastic net penalty can be applied with
overall objective function

``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁``

Where `L` is either the logistic or multinomial loss and `λ` and `γ` indicate
the strength of the L2 (resp. L1) regularisation components.

## Parameters
* `penalty` (Symbol or String): the penalty to use, either `:l2`, `:l1`, `:en`
                                (elastic net) or `:none`. (Default: `:l2`)
* `lambda` (Real): strength of the regulariser if `penalty` is `:l2` or `:l1`.
                   Strength of the L2 regulariser if `penalty` is `:en`.
* `gamma` (Real): strength of the L1 regulariser if `penalty` is `:en`.
* `fit_intercept` (Bool): whether to fit an intercept (Default: `true`)
* `penalize_intercept` (Bool): whether to penalize intercept (Default: `false`)
* `solver` (Solver): type of solver to use, default if `nothing`.
"""
@with_kw_noshow mutable struct LogisticClassifier <: MMI.Probabilistic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::LogisticClassifier, nclasses::Integer) =
    LogisticRegression(m.lambda, m.gamma;
                       penalty=Symbol(m.penalty),
                       multi_class=(nclasses > 2),
                       fit_intercept=m.fit_intercept,
                       penalize_intercept=m.penalize_intercept,
                       nclasses=nclasses)

descr(::Type{LogisticClassifier}) = "Classifier corresponding to the loss function ``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the logistic loss."

#= ======================
   MULTINOMIAL CLASSIFIER
   ====================== =#

"""
$SIGNATURES

See `LogisticClassifier`, it's the same except that `multi_class` is set
to `true` by default. The other parameters are the same.
"""
@with_kw_noshow mutable struct MultinomialClassifier <: MMI.Probabilistic
    lambda::Real             = 1.0
    gamma::Real              = 0.0
    penalty::SymStr          = :l2
    fit_intercept::Bool      = true
    penalize_intercept::Bool = false
    solver::Option{Solver}   = nothing
end

glr(m::MultinomialClassifier, nclasses::Integer) =
    MultinomialRegression(m.lambda, m.gamma;
                          penalty=Symbol(m.penalty),
                          fit_intercept=m.fit_intercept,
                          penalize_intercept=m.penalize_intercept,
                          nclasses=nclasses)

descr(::Type{MultinomialClassifier}) =
    "Classifier corresponding to the loss function " *
    "``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the multinomial loss."


text:
data = map(x->split(x, " -> "),split(strip(read(joinpath(@__DIR__, "./5.txt"), String)), '\n'))
data = map(y->map(x->parse.(Int,split(x,",")),y),data)
data = map(z->reduce(vcat,z),data)

line_plot = zeros(Int,1000,1000)

function points_on_line(x1,y1,x2,y2)
  if x1 > x2
    x1, y1, x2, y2 = x2, y2, x1, y1
  end
  slope = (y2-y1)/(x2-x1)
  y_intercept = y1-slope*x1
  [(x,Int(x*slope+y_intercept)) for x in range(x1, x2) if Int(x*slope) == x*slope]
end

function plotline(x1,y1,x2,y2)
  global line_plot

  if x1 == x2
    if y1 < y2
      line_plot[x1,y1:y2] .+= 1
    else
      line_plot[x1,y2:y1] .+= 1
    end
  elseif y1 == y2
    if x1 < x2
      line_plot[x1:x2,y1] .+= 1
    else
      line_plot[x2:x1,y1] .+= 1
    end
  else
    # this is for part2 of the problem
    points = points_on_line(x2,y2,x1,y1)
    for (x,y) in points
      line_plot[x,y] += 1
    end
  end
end

function solve(data)
  global line_plot

  for d in data
    plotline(d...)
  end

  length(filter(i->i>=2,line_plot))
end

solve(data)


text:
abstract type AbstractMetabolism end
abstract type AbstractMetabolismC <: AbstractMetabolism end
abstract type AbstractMetabolismCN <: AbstractMetabolism end
abstract type AbstractMetabolismCNP <: AbstractMetabolism end
abstract type AbstractAssimilation end

@columns struct MetabolismC{KE,YEV,KM,YEM,AX,YEX,MGT} <: AbstractMetabolismC
    k_E::KE     | 1/h       | "Reserve export rate"
    y_EV::YEV   | mol/mol   | "Yield of structure on reserve"
    k_M::KM     | 1/h       | "Specific maintenance rate"
    y_EM::YEM   | mol/mol   | "Maintenance yield on reserve vs structure"
    α_X::AX     | _         | "Fraction of growth flux allocated to enzyme production"
    y_EX::YEX   | mol/mol   | "Yield on enzyme production"
    min_gt::MGT | 1/h       | "Minimum generation time"
end

for fn in fieldnames(MetabolismC)
    @eval $fn(p::MetabolismC) = p.$fn
end

growth!(r0, p::MetabolismC, E::Vector{Float64}, V::Vector{Float64}) = begin
    y_VM = y_EM(p)./y_EV(p)

    function f(r, j)
      m_E   = max(1e-8, E[j]/V[j])
      j_EC  = m_E*(k_E(p)[j] - r)
      j_EM  = k_M(p)[j]*y_EM(p)[j]
      jEM   = max(1e-8, min(j_EC, j_EM))
      jVM   = (j_EM - jEM)*y_VM[j]/y_EM(p)[j]
      j_EG  = (j_EC - jEM)/(1.0 + α_X(p)[j])
      j_G   = j_EG/y_EV(p)[j]
      res   = r - (j_G - jVM)
    end

    f1 = [r->f(r,j) for j in 1:size(r0,1)]
    r = Roots.find_zero.(f1, r0)
    #gmax = log(2)./min_gt(p)
    #r_b = min.(r, gmax)
end

growth_production!(r, p::MetabolismC, E::Vector{Float64}, V::Vector{Float64}) = begin
      y_VM  = y_EM(p)./y_EV(p)

      m_E   = E./V
      j_EC  = m_E.*(k_E(p) .- r)
      j_EM  = k_M(p).*y_EM(p)
      jEM   = min(j_EC, j_EM)
      jVM   = (j_EM - jEM).*y_VM./y_EM(p)
      j_EG  = (j_EC - jEM)./(1.0 .+ α_X(p))
      x     = α_X(p).*j_EG./y_EX(p)
      rG_CO2  = (1 .- 1 ./y_EV(p)).*j_EG.*V
      rM_CO2  = @. (jEM + jVM)*V
      rX_CO2  = (1 .- 1 ./y_EX(p)).*α_X(p).*j_EG.*V
      return x, rG_CO2, rM_CO2, rX_CO2
end

enzyme_production!(x, p::MetabolismC, V::Vector{Float64}) = begin
    J_EX = sum(x.*V, dims=1)
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using GraphMatFun;
include("setup_graphs.jl");

# Generate #include preprocessor directives for C code.
function gen_c_include(names,mkl,skip)
    blas_string=mkl ? "_MKL" : "_OpenBLAS"
    include_string=""
    for name in names
        if !(name in string.(skip))
            include_string*="#include\"$(tempdir())/$name$blas_string.c\"\n"
        end
    end
    return include_string
end

# Process code from template.
# Used for C, julia, and matlab.
function gen_main(fname,N,graphs,names,normvalues,mkl,commentsign="#",skip=[:expm_matlab,:expmpoly_matlab])

    # Common preprocessing:
    # * include files with generated code (C only);
    # * replace placeholders that appear only once; and
    # * find beginning and end of section to be repeated.
    lines=readlines(fname);
    start_repeated_line=0;
    end_repeated_line=0;
    for (i,line)=enumerate(lines)
        line=replace(line, "LOC_INC" => gen_c_include(names,mkl,skip)); # C only
        line=replace(line, "MATSIZE" => string(N));
        if (mkl)
            line=replace(line, "USING" => "using MKL;");          # Julia only
        else
            line=replace(line, "USING" => "$commentsign NO MKL"); # Julia only
        end

        lines[i]=line
        if contains(line,"START REPEATED CODE")
            start_repeated_line=i;
        end
        if contains(line,"END REPEATED CODE")
            end_repeated_line=i;
        end
    end

    repeated_code=lines[start_repeated_line+1:end_repeated_line-1];

    postprocess=lines[(end_repeated_line+1):end];
    lines=lines[1:(start_repeated_line-1)]
    for (i,g)=enumerate(graphs)
        local n;
        if (!(g in skip))
            n=names[i];
            this_sim_code=deepcopy(repeated_code);
            for (j,line)=enumerate(this_sim_code)
                if (g isa Compgraph)
                    tempdir_str=tempdir();
                    line=replace(line,"INCLUDE" =>
                      string("include(joinpath(\"$(tempdir_str)\",\"NAME.jl\"));"));
                else
                    line=replace(line,"INCLUDE" => "exp_julia=exp");
                end

                line=replace(line, "MATNORM" => string(normvalues[i]));  # Set the norm
                line=replace(line,"NAME" => n);
                line=replace(line,"FUNCTION" => n);

                this_sim_code[j]=line;
            end

            # Append
            map(code->push!(lines,code), this_sim_code)
        else
            push!(lines,"$commentsign SKIPPING $g");
            push!(lines,"");
        end
    end

    lines=[lines;postprocess]

    return lines;
end

# Asymptotic number of flops needed to evaluate the a graph.
function get_cost(graph::Compgraph)
    return count(values(graph.operations) .== :mult)+
           (4/3)*count(values(graph.operations) .== :ldiv)
end

# Generate functions for evaluating the graphs.
println("Generating graphs");
for (i,g)=enumerate(graphs)
    local n;

    n=names[i];

    print("Generating code for $n ");
    if (g isa Compgraph)
        print("cost $(get_cost(g)) ");
        gen_code(string(tempdir(),"/$n.jl"),g,funname="$n",priohelp=priohelp);
        gen_code(string(tempdir(),"/$n.m"),g,funname="$n",lang=LangMatlab(),priohelp=priohelp);

        # Distinguish between MKL and OpenBLAS in the graph generators
        # already for C. Since includes etc are different.
        fname_mkl=string(tempdir(),"/$(n)_MKL");
        fname_openblas=string(tempdir(),"/$(n)_OpenBLAS");
        gen_code("$(fname_mkl).c",g,funname="$n",lang=GraphMatFun.LangC_MKL(),priohelp=priohelp);
        gen_code("$(fname_openblas).c",g,funname="$n",lang=GraphMatFun.LangC_OpenBLAS(),priohelp=priohelp);

    end
    println();
end

n=4000
# Generate C code.
println("Generating C/MKL");
c_template="src/cputime/c_main_template.c";
lines=gen_main(c_template,n,graphs,names,normvalues,true,"//",
               [:exp_julia,:expm_matlab,:expmpoly_matlab]);
open(string(tempdir(),"/run_cputime_c_MKL.c"),"w") do io
    for line=lines
        write(io,line,"\n");
    end
end
println("Generating C/OpenBLAS");
lines=gen_main(c_template,n,graphs,names,normvalues,false,"//",
               [:exp_julia,:expm_matlab,:expmpoly_matlab]);
open(string(tempdir(),"/run_cputime_c_OpenBLAS.c"),"w") do io
    for line=lines
        write(io,line,"\n");
    end
end

# Generate Julia code.
println("Generating julia/MKL");
julia_template="src/cputime/julia_main_template.jl";
lines=gen_main(julia_template,n,graphs,names,normvalues,true);
open(string(tempdir(),"/run_cputime_julia_MKL.jl"),"w") do io
    for line=lines
        write(io,line,"\n");
    end
end
println("Generating julia/OpenBLAS");
lines=gen_main(julia_template,n,graphs,names,normvalues,false);
open(string(tempdir(),"/run_cputime_julia_OpenBLAS.jl"),"w") do io
    for line=lines
        write(io,line,"\n");
    end
end

# Generate MATLAB code.
println("Generating matlab/MKL");
lines=gen_main("src/cputime/matlab_main_template.m",n,graphs,names,normvalues,false,"%",[:exp_julia]);
open(string(tempdir(),"/run_cputime_matlab_OpenBLAS.m"),"w") do io
    for line=lines
        write(io,line,"\n");
    end
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
#using Retry
using DelimitedFiles
datadir = joinpath(@__DIR__, "..", "src", "data")
isdir(datadir) || mkdir(datadir)
@info "Downloading Bollerslev and Ghysels data..."
isfile(joinpath(datadir, "bollerslev_ghysels.txt")) || download("http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF20-1.txt", joinpath(datadir, "bollerslev_ghysels.txt"))

# @info "Downloading stock data..."
# #"DOW" is excluded because it's listed too late
# tickers = ["AAPL", "IBM", "XOM", "KO", "MSFT", "INTC", "MRK", "PG", "VZ", "WBA", "V", "JNJ", "PFE", "CSCO", "TRV", "WMT", "MMM", "UTX", "UNH", "NKE", "HD", "BA", "AXP", "MCD", "CAT", "GS", "JPM", "CVX", "DIS"]
# alldata = zeros(2786, 29)
# for (j, ticker) in enumerate(tickers)
#         @repeat 4 try
#             @info "...$ticker"
#             filename = joinpath(datadir, "$ticker.csv")
#             isfile(joinpath(datadir, "$ticker.csv")) || download("http://quotes.wsj.com/$ticker/historical-prices/download?num_rows=100000000&range_days=100000000&startDate=03/19/2008&endDate=04/11/2019", filename)
#             data = parse.(Float64, readdlm(joinpath(datadir, "$ticker.csv"), ',', String, skipstart=1)[:, 5])
#             length(data) == 2786 || error("Download failed for $ticker.")
#             alldata[:, j] .= data
#             rm(filename)
#         catch e
#             @delay_retry if 1==1 end
#         end
# end
# alldata = 100 * diff(log.(alldata), dims=1)
# open(joinpath(datadir, "dow29.csv"), "w") do io
#      writedlm(io, alldata, ',')
# end


text:
using BEAST
using CompScienceMeshes
using WiltonInts84
using Base.Test

p1 = point(0,0,0)
p2 = point(1,0,0)
p3 = point(0,1,0)

q1 = point(2,2,0)
q2 = point(2,1,0)
q3 = point(1,2,0)

τ = [p1,p2,p3]
σ = [q1,q2,q3]
m, M = BEAST.minmaxdist(τ, σ)

@show m, sqrt(2.0)
@show M, sqrt(8.0)

@test m ≈ sqrt(2.0)
@test M ≈ sqrt(8.0)

ΔR = sqrt(2)/2
r = BEAST.rings(τ,σ,ΔR)
@show r

cm = point(0.5,0.5,0.0)
cM = point(0.0,0.0,0.0)
S,V,G = WiltonInts84.wiltonints(q1,q2,q3,cm,(r[1]-2)*ΔR,(r[1]-1)*ΔR,Val{0})
@show S[3]

S,V,G = WiltonInts84.wiltonints(q1,q2,q3,cM,(r[end]-1)*ΔR,(r[end])*ΔR,Val{0})
@show S[3]

h = 0.05
Γ = meshrectangle(10.0, h, h, 3)
X = raviartthomas(Γ)

timestep = 1.0
numsteps = 10
S = timebasisc0d1(timestep, numsteps)
S = BEAST.timebasisspline2(timestep, numsteps)

c = 1.0
Ts = MWSingleLayerTDIO(c,1.0,0.0,0,0)
Th = MWSingleLayerTDIO(c,0.0,1.0,0,0)

Zs = BEAST.assemble(Ts,X,X,S)
Zh = BEAST.assemble(Th,X,X,S)

## Test the allcatestorage routine and in particular its type stability
h = 0.05
Γ = meshrectangle(10.0, h, h, 3)
X = raviartthomas(Γ)

timestep = 1.0
numsteps = 10
S = timebasisc0d1(timestep, numsteps)
S = BEAST.timebasisspline2(timestep, numsteps)

c = 1.0
Ts = MWSingleLayerTDIO(c,1.0,0.0,0,0)
Th = MWSingleLayerTDIO(c,0.0,1.0,0,0)

U = timebasiscxd0(Δt, Nt)
T = timebasisc0d1(Δt, Nt)

W = X ⊗ U
V = X ⊗ T
@code_warntype BEAST.allocatestorage(Ts,W,V)


text:
#=
Model: toy.jl

Toy model to test the implementation on certain cases.
=#
using Reachability

function compute(input_options::Pair{Symbol,<:Any}...)
    # =====================
    # Problem specification
    # =====================
    A = sparse([1, 2], [1, 2], [0.1, 0.1])

    # initial set
    X0 = BallInf([1.0, 1.0], 0.5)

    # instantiate continuous LTI system
    S = ContinuousSystem(A, X0)

    # ===============
    # Problem solving
    # ===============

    # define solver-specific options
    options = merge(Options(
        :mode => "reach",
        :property => SafeStatesProperty(HalfSpace([1., 0.], 20.)), # x1 <= 20
        :blocks => [1],
        :plot_vars => [0, 1]
        ), Options(input_options...))

    result = solve(S, options)

    # ========
    # Plotting
    # ========
    if options[:mode] == "reach"
        println("Plotting...")
        tic()
        plot(result)
        @eval(savefig(@filename_to_png))
        toc()
    en
end # function

compute(:N => 10, :T => 20.0); # warm-up
compute(:δ => 0.01, :T => 20.0); # benchmark settings (long)

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
@testset "CompoundPoissonProcess" begin

using DiffEqNoiseProcess, DiffEqBase, Test, Statistics

r = 2
rate(u,p,t) = r
W = CompoundPoissonProcess(rate,0.0,1.0)
dt = 0.1
calculate_step!(W,dt,nothing,nothing)
for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

W = CompoundPoissonProcess(rate,0.0,1.0)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

W = CompoundPoissonProcess(rate,0//1,1.0)
prob_rational = NoiseProblem(W,(0,1))
dt = 1//10
sol = solve(prob_rational;dt=dt)

dt = dt/100
for t in dt:dt:1-dt
  sol(t)
end

t = 1.0
u0 = 1.0
expected_mean = 1 + t*r
expected_variance = t*r
ensemble_prob = EnsembleProblem(prob;output_func = (sol,i)-> (sol[end],false))
sol = solve(ensemble_prob;dt=0.1,trajectories=40000)
@test abs(mean(sol) - expected_mean) < 0.4
@test abs(var(sol) - expected_variance) < 0.4 # Converges slowly

rate(du,u,p,t) = (@. du = r)
W = CompoundPoissonProcess!(rate,0.0,ones(2))
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

ensemble_prob = EnsembleProblem(prob;output_func = (sol,i)-> (sol[end],false))
sol = solve(ensemble_prob;dt=0.1,trajectories=40000)
@test abs(mean(sol) - expected_mean) < 0.4
@test abs(var(sol) - expected_variance) < 0.4 # Converges slowly
end


text:
using QuasinormalModes
using Test

function TestIfValueOfRadialWorks()
    a = -0.402546147858408 - 0.06278653515072422*im;
    Ψ = qnmfunction(s=-2,l=2,n=2,m=-2,a=0.4)
    (abs(Ψ(2.0)-a) < 10^(-9))
end

@testset "QuasinormalModes.jl" begin
    # Write your tests here.
    @test TestIfValueOfRadialWorks()
end


text:
module MyPrintf
    using Printf

    function myprintf(str, precision, value)
        fmt = str * "%." * string(precision) * "f\n"

        @eval @printf($fmt, $value)
    end
end


text:
#InspectDR: Marker control for plot widgets
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
const MARKER_COLORANT = ARGB32(0, 0, 0, 0.5)
const MARKER_LINE = LineAttributes(:solid, 2.5, MARKER_COLORANT)
const CTRLPOINT_RADIUS = Float64(4)
const CTRLPOINT_ATTR = AreaAttributes(
	line=LineAttributes(:solid, 2.5, MARKER_COLORANT), fillcolor=COLOR_TRANSPARENT
)

#==Types
===============================================================================#
struct ISMovingMarker <: InputState
	istrip::Int
	initpos::Point2D
	marker::CtrlMarker #Ref to marker
end
struct ISMovingΔInfo <: InputState
	istrip::Int
	refpos::Point2D
	initΔ::Vector2D
	marker::CtrlMarker #Ref to marker
end

#==Constructor-like functions
===============================================================================#
function CtrlMarker(x::Real, y::Real, strip::Int, ref=nothing)
	mkr = hvmarker(Point2D(x, y), MARKER_LINE, strip=strip)
	return CtrlMarker(mkr, Vector2D(0,0), BoundingBox(), ref)
end

#==Rendering functions
===============================================================================#
#TODO: move somewhere else?
function render_ctrlpoint(canvas::PCanvas2D, pt::Point2D, ixf::InputXfrm2D)
	ctx = canvas.ctx #WANTCONST
	pt = read2axis(pt, ixf)
	pt = map2dev(canvas.xf, pt)
	setlinestyle(ctx, LineStyle(CTRLPOINT_ATTR.line))
	#TODO: Deprecate workaround if behaviour is changed.
	#NOTE: Will draw line if we don't "move to" 0rad point on arc before drawing:
	#      why is this?  is this a bug? why does this not show up elsewhere?
	Cairo.move_to(ctx, pt.x+CTRLPOINT_RADIUS, pt.y) #Workaround
	cairo_circle(ctx, pt.x, pt.y, CTRLPOINT_RADIUS)
	renderfill(ctx, CTRLPOINT_ATTR.fillcolor)
	Cairo.stroke(ctx)
	return
end

function render_markercoord(canvas::PCanvas2D, pt::Point2D, font::Font,
	xfmt::NumericFormatting, yfmt::NumericFormatting, graphinfo::Graph2DInfo, strip::Int)
	xstr = formatted(pt.x, xfmt)
	ystr = formatted(pt.y, yfmt)
	str = "$xstr, $ystr"
	a = atext(str, x=pt.x, y=pt.y, xoffset=3, yoffset=-3,
		font=font, align=:tl, strip=strip)
	render(canvas, a, graphinfo)
	return
end

function render_Δinfo(canvas::PCanvas2D, p1::Point2D, p2::Point2D, Δinfo::Vector2D,
	font::Font, xfmt::NumericFormatting, yfmt::NumericFormatting, ixf::InputXfrm2D, strip::Int)
	align = ALIGN_TOP | ALIGN_LEFT #WANTCONST
	mfmt = number_fmt(ndigits=4, decfloating=true) #WANTCONST
	boxattr = AreaAttributes(
		line(style=:solid, width=1.5, color=COLOR_BLACK), COLOR_WHITE
	) #WANTCONST
	ctx = canvas.ctx #WANTCONST
	Δx = p2.x - p1.x; Δy = p2.y - p1.y
	m = Δy/Δx
	Δxstr = formatted(Δx, xfmt); Δxstr = "Δx=$Δxstr"
	Δystr = formatted(Δy, yfmt); Δystr = "Δy=$Δystr"
	mstr = formatted(m, mfmt); mstr = "Δy/Δx=$mstr"

	setfont(ctx, font)
	w, h = textextents_wh(ctx, mstr)
	Δh = h*.1 #Vertical gap
	wx, hx = textextents_wh(ctx, Δxstr)
	wy, hy = textextents_wh(ctx, Δxstr)
	w = max(w, wx, wy)
	h += hx + hy + 2*Δh

	#Figure out where to display:
	p1a = read2axis(p1, ixf); 	p2a = read2axis(p2, ixf)
	pdraw = Point2D((p1a.x + p2a.x)/2, (p1a.y + p2a.y)/2)
		pdraw = map2dev(canvas.xf, pdraw)
	x = pdraw.x + Δinfo.x; y = pdraw.y + Δinfo.y
	x -= w/2; y -= h/2
	pad = Δh
	rect = BoundingBox(x-pad, x+w+pad, y-pad, y+h+pad)
	drawrectangle(ctx, rect, boxattr)
	render(ctx, Δxstr, Point2D(x, y), font, align=align)
	y += hx + Δh
	render(ctx, Δystr, Point2D(x, y), font, align=align)
	y += hy + Δh
	render(ctx, mstr, Point2D(x, y), font, align=align)
	return rect #Want to know where we drew the box
end

function render(canvas::PCanvas2D, mg::CtrlMarkerGroup, graphinfo::Graph2DInfo, strip::Int)
	xfmt = hoverfmt(graphinfo.xfmt) #WANTCONST
	yfmt = hoverfmt(graphinfo.yfmt) #WANTCONST
	ixf = graphinfo.ixf #WANTCONST
	for elem in mg.elem
		mstrip = elem.prop.strip
		if 0 == mstrip || mstrip == strip
			render(canvas, elem.prop, graphinfo) #Render crosshairs
			render_markercoord(canvas, elem.prop.pos, mg.fntcoord, xfmt, yfmt, graphinfo, strip)
		end
	end

	#Render Δinfo on top:
	for elem in mg.elem
		mstrip = elem.prop.strip
		if 0 == mstrip || mstrip == strip
			mref = elem.ref
			if mref != nothing
				elem.Δbb = render_Δinfo(canvas, mref.prop.pos, elem.prop.pos, elem.Δinfo,
					mg.fntdelta, xfmt, yfmt, ixf, strip
				)
			end
		end
	end

	#Render control points on top:
	for elem in mg.elem
		mstrip = elem.prop.strip
		if 0 == mstrip || mstrip == strip
			render_ctrlpoint(canvas, elem.prop.pos, ixf)
		end
	end
	return
end

#==Helper functions
===============================================================================#
function hittest(marker::CtrlMarker, xf::Transform2D, ixf::InputXfrm2D, x::Float64, y::Float64)
	Δhit = Float64(CTRLPOINT_RADIUS+CTRLPOINT_ATTR.line.width/2) #WANTCONST

	pt = read2axis(marker.prop.pos, ixf)
	pt = map2dev(xf, pt)
	Δx = abs(x-pt.x); Δy = abs(y-pt.y)
	Δ = sqrt(Δx*Δx + Δy*Δy)
	return (Δ <= Δhit) #Hit successful if clicked on marker.
end
function hittest_Δinfo(marker::CtrlMarker, x::Float64, y::Float64)
	return isinside(marker.Δbb, x, y)
end

#==Action functions for bindkeys
===============================================================================#
function createmarker(pwidget::PlotWidget, ref=nothing)
	istrip = pwidget.mouseover.istrip
	pos = pwidget.mouseover.pos
	mkr = nothing
	if istrip > 0 && pos != nothing
		x, y = pos.x, pos.y
		mkr = CtrlMarker(x, y, istrip, ref)
	end
	return mkr
end

function addrefmarker(pwidget::PlotWidget)
	mkr = createmarker(pwidget)
	if mkr != nothing
		push!(pwidget.markers.elem, mkr)
		pwidget.refmarker = mkr
		refresh(pwidget, refreshdata=false)
	end
	return mkr
end

function addΔmarker(pwidget::PlotWidget, makeref::Bool=false)
	mkr = nothing
	if pwidget.refmarker != nothing
		mkr = createmarker(pwidget, pwidget.refmarker)
		if mkr != nothing
			push!(pwidget.markers.elem, mkr)
			if makeref
				pwidget.refmarker = mkr
			end
			refresh(pwidget, refreshdata=false)
		end
	end
	return mkr
end

#Add a delta marker - and make it the new reference:
addΔmarkerref(pwidget::PlotWidget) = addΔmarker(pwidget, true)

#==Other actions
===============================================================================#
function cancelmove(s::ISMovingMarker, pwidget::PlotWidget)
	s.marker.prop.pos = s.initpos
	setstate_normal(pwidget)
	refresh(pwidget, refreshdata=false)
	return
end
function deletemarker(mkr::CtrlMarker, pwidget::PlotWidget)
	markers = pwidget.markers.elem
	n = length(markers)
	for i in 1:n
		if mkr == markers[i]
			deleteat!(markers, i)
			break
		end
	end
	refresh(pwidget, refreshdata=false)
	return
end
function cancelmove_Δinfo(s::ISMovingΔInfo, pwidget::PlotWidget)
	s.marker.Δinfo = s.initΔ
	setstate_normal(pwidget)
	refresh(pwidget, refreshdata=false)
	return
end
function centerΔinfo(s::ISMovingΔInfo, pwidget::PlotWidget)
	s.marker.Δinfo = Vector2D(0,0)
	setstate_normal(pwidget)
	refresh(pwidget, refreshdata=false)
	return
end

#==Event handlers for selecting markers
===============================================================================#
#Handle mousepress if a CtrlMarker was clicked:
function handleevent_mousepress(pwidget::PlotWidget, markers::CtrlMarkerGroup,
	istrip::Int, x::Float64, y::Float64)
	xf = Transform2D(pwidget.plotinfo, istrip)
	ixf = InputXfrm2D(pwidget.plotinfo, istrip)

	#See if we clicked on a marker:
	for i in 1:length(markers.elem)
		elem = markers.elem[i]
		if hittest(elem, xf, ixf, x, y)
			initpos = elem.prop.pos
			gdk_window_set_cursor(pwidget.canvas, CURSOR_MOVE)
			pwidget.state = ISMovingMarker(istrip, initpos, elem)
			return true
		end
	end

	#Else see if we clicked on a Δ-info box:
	for i in 1:length(markers.elem)
		elem = markers.elem[i]
		if hittest_Δinfo(elem, x, y)
			refpos = Point2D(x, y) - elem.Δinfo
			gdk_window_set_cursor(pwidget.canvas, CURSOR_MOVE)
			pwidget.state = ISMovingΔInfo(istrip, refpos, elem.Δinfo, elem)
			return true
		end
	end

	return false
end

#==Event handlers for state: ISMovingMarker
===============================================================================#
function handleevent_keypress(s::ISMovingMarker, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval #Cancel move
		cancelmove(s, pwidget)
	elseif GdkKeySyms.Delete == event.keyval
		cancelmove(s, pwidget)
		deletemarker(s.marker, pwidget)
	end
end
function handleevent_mouserelease(::ISMovingMarker, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	if 1==event.button #Done moving
		setstate_normal(pwidget)
	end
end
function handleevent_mousemove(s::ISMovingMarker, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	handleevent_plothover(pwidget, event, s.istrip)

	xf = Transform2D(pwidget.plotinfo, s.istrip)
	ixf = InputXfrm2D(pwidget.plotinfo, s.istrip)
	pt = map2axis(xf, Point2D(event.x, event.y))
	pt = axis2read(pt, ixf)
	s.marker.prop.pos = pt
	refresh(pwidget, refreshdata=false)
end

#==Event handlers for state: ISMovingΔInfo
===============================================================================#
function handleevent_keypress(s::ISMovingΔInfo, pwidget::PlotWidget, event::Gtk.GdkEventKey)
	if GdkKeySyms.Escape == event.keyval #Cancel move
		cancelmove_Δinfo(s, pwidget)
	elseif Int('0') == event.keyval
		centerΔinfo(s, pwidget)
	end
end
function handleevent_mouserelease(::ISMovingΔInfo, pwidget::PlotWidget, event::Gtk.GdkEventButton)
	if 1==event.button #Done moving
		setstate_normal(pwidget)
	end
end
function handleevent_mousemove(s::ISMovingΔInfo, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	handleevent_plothover(pwidget, event, s.istrip)

	curpos = Point2D(event.x, event.y)
	s.marker.Δinfo = Vector2D(curpos - s.refpos)
	refresh(pwidget, refreshdata=false)
end

#Last line


text:
using DelftBlue
DelftBlue.main(4,2,1,0.001,0.002)
DelftBlue.main(4,2,1,0.001,0.002,test="Cylinder",mesh_file="mesh_coarse.msh")


text:
using LoopVectorization

#--------------------------------------------------------------------#

function Generate2DRectMesh(nx,ny,xmin,xmax,ymin,ymax)
    # Generates a 2D rectangular mesh of nx*ny cells
    ncell      = nx*ny
    # Twice resolution mesh
    nx2 = 2*nx + 1
    ny2 = 2*ny + 1
    # 1D axis
    x2  = LinRange(xmin,xmax,nx2)
    y2  = LinRange(ymin,ymax,ny2)
    # 2D mesh: fake ngrid for P-T space ;)
    x2d   = repeat(x2, 1, length(y2))
    y2d   = repeat(y2, 1, length(x2))'
    nodes = zeros(Int64, size(x2d))
    # Number active dofs (midfaces) - Face node numbering
    inum = 0
    for i=1:nx2
        for j=1:ny2
            if (mod(j,2)==1 && mod(i,2)==0) || (mod(j,2)==0 && mod(i,2)==1)
                inum += 1
                nodes[i,j] = inum
            end
        end
    end
    # Face nodes
    xp = zeros(inum)
    yp = zeros(inum)
    tp = zeros(Int64,inum) # node type - 1 boundary
    for i=1:nx2
        for j=1:ny2
            if (mod(j,2)==1 && mod(i,2)==0) || (mod(j,2)==0 && mod(i,2)==1)
                xp[nodes[i,j]] = x2d[i,j]
                yp[nodes[i,j]] = y2d[i,j]
                if (i==1 || i==nx2 || j==1 || j==ny2)
                    tp[nodes[i,j]] = 1
                end
            end
        end
    end
    # # South/North nodes
    # @avx for i=2:2:nx2-1
    #     for j=1:2:ny2
    #         xp[nodes[i,j]] = x2d[i,j]
    #         yp[nodes[i,j]] = y2d[i,j]
    #     end
    # end
    # # West/East nodes
    # @avx for i=1:2:nx2
    #     for j=2:2:ny2-1
    #         xp[nodes[i,j]] = x2d[i,j]
    #         yp[nodes[i,j]] = y2d[i,j]
    #     end
    # end
    # Cell 2 face numbering
    cell = zeros(4, ncell)
    @avx for i=1:nx
        for j=1:ny
            k  = i + (j-1)*nx
            ic = i+1 + (i-1)*1
            jc = j+1 + (j-1)*1
            cell[1,k] =  nodes[ic,jc-1] # South
            cell[2,k] =  nodes[ic+1,jc] # East
            cell[3,k] =  nodes[ic,jc+1] # North
            cell[4,k] =  nodes[ic-1,jc] # West
        end
    end
    return cell,xp,yp,tp
end

#--------------------------------------------------------------------#

# Create sides of mesh
xmin, xmax = 0, 1
ymin, ymax = 0, 1
nx, ny     = 5, 5  # number of elements

@time cell,xp,yp,tp = Generate2DRectMesh(nx,ny,xmin,xmax,ymin,ymax)
@time cell,xp,yp,tp = Generate2DRectMesh(nx,ny,xmin,xmax,ymin,ymax)
@time cell,xp,yp,tp = Generate2DRectMesh(nx,ny,xmin,xmax,ymin,ymax)
@time cell,xp,yp,tp = Generate2DRectMesh(nx,ny,xmin,xmax,ymin,ymax)


text:
using .MLJ

export deterministic, probabilistic

# deterministic model api
function deterministic(model, Xs::MLJ.Source, ys::MLJ.Source)
    model = scitype(model) <: MLJ.UnsupervisedScitype ?
            transform(machine(model, Xs), Xs) :
            predict(machine(model, Xs, ys), Xs)
    mach = machine(Deterministic(), Xs, ys, predict=model)
    return mach
end

# deterministic network api
function deterministic(model, Xs::MLJ.Source, ys::MLJ.Source, scores...)
    model = transform(machine(model), scores...)
    mach = machine(Deterministic(), Xs, ys, predict=model)
    return mach
end

# probabilistic model api
function probabilistic(model, Xs::MLJ.Source, ys::MLJ.Source)
    model = scitype(model) <: MLJ.UnsupervisedScitype ?
            transform(machine(model, Xs), Xs) :
            predict(machine(model, Xs, ys), Xs)
    mach = machine(Probabilistic(), Xs, ys, predict=model)
    return mach
end

# probabilistic network api
function probabilistic(model, Xs::MLJ.Source, ys::MLJ.Source, scores...)
    model = transform(machine(model), scores...)
    mach = machine(Probabilistic(), Xs, ys, predict=model)
    return mach
end