text:
# Scale-Location Student-t Stein Distribution
#
# Represents a mu + sigma * X where X is a univariate t-distribution with nu
# degrees of freedom
using Distributions: TDist
import Distributions

type SteinScaleLocationStudentT <: SteinDistribution
    nu::Float64
    mu::Float64
    sigma::Float64
    c1::Float64
    c2::Float64
    c3::Float64
end

# Constructor with canonical setting of (c1,c2,c3) = (1,1,1)
SteinScaleLocationStudentT(nu::Float64, mu::Float64, sigma::Float64) =
SteinScaleLocationStudentT(nu, mu, sigma, 1.0, 1.0, 1.0);

# Constructor for standard t-distribution with canonical constants c1, c2, c3
SteinScaleLocationStudentT(nu::Float64) =
SteinScaleLocationStudentT(nu, 0.0, 1.0, 1.0, 1.0, 1.0);

# Draw n independent samples from distribution
function rand(d::SteinScaleLocationStudentT, n::Int64)
    d.mu + d.sigma * Distributions.rand(TDist(d.nu), n)'';
end

function supportlowerbound(d::SteinScaleLocationStudentT, j::Int64)
    -Inf
end

function supportupperbound(d::SteinScaleLocationStudentT, j::Int64)
    Inf
end

# Gradient of log density of distribution evaluated at each
# entry of n x 1 matrix X
function gradlogdensity(d::SteinScaleLocationStudentT, X::Array{Float64,2})
    -(d.nu + 1) * (X - d.mu) ./ (d.nu * d.sigma^2 + (X - d.mu).^2);
end

# Cumulative distribution function
function cdf(d::SteinScaleLocationStudentT, t)
    Distributions.cdf(TDist(d.nu),(t-d.mu)/d.sigma);
end

function numdimensions(d::SteinScaleLocationStudentT)
    1
end


text:
function get_derivatives(u::Vector{Float64},p::Vector{Float64})
    # derivatives: dF/d[bifurcation_param]
    dFdp::Vector{Float64} = zeros(V.NUM)

    dFdp[V.M] = p[C.kM]*p[C.KS]/(p[C.KS] + p[C.S])^2
    dFdp[V.CD] = p[C.kCDS]*p[C.KS]/(p[C.KS] + p[C.S])^2

    return dFdp
end

ext: jl

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


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

mutable struct MonetaryCredit <: SwaggerModel
    Name::Any # spec type: Union{ Nothing, String } # spec name: Name
    EffectiveDate::Any # spec type: Union{ Nothing, DateTime } # spec name: EffectiveDate

    function MonetaryCredit(;Name=nothing, EffectiveDate=nothing)
        o = new()
        validate_property(MonetaryCredit, Symbol("Name"), Name)
        setfield!(o, Symbol("Name"), Name)
        validate_property(MonetaryCredit, Symbol("EffectiveDate"), EffectiveDate)
        setfield!(o, Symbol("EffectiveDate"), EffectiveDate)
        o
    end
end # type MonetaryCredit

const _property_map_MonetaryCredit = Dict{Symbol,Symbol}(Symbol("Name")=>Symbol("Name"), Symbol("EffectiveDate")=>Symbol("EffectiveDate"))
const _property_types_MonetaryCredit = Dict{Symbol,String}(Symbol("Name")=>"String", Symbol("EffectiveDate")=>"DateTime")
Base.propertynames(::Type{ MonetaryCredit }) = collect(keys(_property_map_MonetaryCredit))
Swagger.property_type(::Type{ MonetaryCredit }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_MonetaryCredit[name]))}
Swagger.field_name(::Type{ MonetaryCredit }, property_name::Symbol) =  _property_map_MonetaryCredit[property_name]

const _allowed_MonetaryCredit_Name = ["Recurring Charge", "Monetary Commitment", "Monetary Credit"]

function check_required(o::MonetaryCredit)
    true
end

function validate_property(::Type{ MonetaryCredit }, name::Symbol, val)
    if name === Symbol("Name")
        Swagger.validate_param(name, "MonetaryCredit", :enum, val, _allowed_MonetaryCredit_Name)
    end
end


text:
function encode(input::S) where S <: AbstractString
    flipalpha.(codeunits(input)) |> skipmissing |> collect |> s->StringPadIterator(5,s) |> collect |> String |> rstrip
    #flipalpha.(codeunits(input)) |> skipmissing |> s->padtruth(5,s)
end

decode(input::T) where T <:AbstractString = flipalpha(input) |> skipmissing |> collect |> T

struct StringPadIterator{T<:Union{AbstractString,AbstractArray}}; n::Int; s::T; end

# Indexing directly into strings like this is not a good idea, e.g. "⊗j[2]"
# Luckily we drop all difficult characters in flipalpha so we needn't worry
#
# Ideally this would accept an iterator in place of s so we didn't need to collect it needlessly
@inline Base.iterate(it::StringPadIterator) = (it.s[1],2)
@inline Base.iterate(it::StringPadIterator,pos) = begin
    if (length(it.s) < pos - pos÷(it.n+1))
        nothing
    # todo: if this is the last character, don't return it
    elseif rem(pos,it.n+1) == 0
        (' ',pos+1)
    else
        (it.s[pos - pos÷(it.n+1)], pos+1)
    end
end
Base.length(it::StringPadIterator) = length(it.s) + length(it.s)÷it.n
Base.eltype(::Type{StringPadIterator{T}}) where T = Char

# This is slower than the iterator but it was much easier to write;
# I used it to check the correctness of the function above
function padtruth(n,s)
    a = ""
    for (i,c) in enumerate(s)
        a *= rem(i,n)==0 ? "$(Char(c)) " : Char(c)
    end
    a
end

const LOWERA = codepoint('a')
const LOWERZ = codepoint('z')
const UPPERA = codepoint('A')
const UPPERZ = codepoint('Z')

flipalpha(s::AbstractString) = flipalpha.(codeunits(s))
function flipalpha(c::T)::Union{UInt8,Missing} where T <: Unsigned
    if LOWERA <= c <= LOWERZ
        T(LOWERA+LOWERZ) - c
    elseif UPPERA <= c <= UPPERZ
        T(LOWERA+LOWERZ) - c - T(LOWERA-UPPERA) # convert to lowercase
    elseif codepoint('0') <= c <= codepoint('9')
        c
    else
        missing
    end
end
# For testing:
# using BenchmarkTools
# @benchmark encode(s) setup=(s=rand('a':'z',100_000)|>join)


text:
# Autogenerated wrapper script for CImGui_jll for x86_64-linux-gnu
export compile_commands, libcimgui, libimgui

JLLWrappers.@generate_wrapper_header("CImGui")
JLLWrappers.@declare_file_product(compile_commands)
JLLWrappers.@declare_library_product(libcimgui, "libcimgui.so")
JLLWrappers.@declare_library_product(libimgui, "libimgui-cpp.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        compile_commands,
        "share/compile_commands.json",
    )

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

    JLLWrappers.@init_library_product(
        libimgui,
        "lib/libimgui-cpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Production-Destruction
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/ProductionDestruction.ipynb)
#
#md # !!! note "Overview"
#md #     System type: rational function\
#md #     State dimension: 3\
#md #     Application domain: Chemical kinetics
#
# ## Model description
#
# A production-destruction system consists of an ordinary differential equation and two constraints:
# positivity and conservativity. It means that the system states are quantities which are always positive
# and the sum of these quantities is constant. This particular family of systems is often used to test
# the stability of integration schemes.
#
# As proposed by **Michaelis-Menten** theory [^KM18], a model for three quantities can be defined as follows:
#
# ```math
#    \left\{ \begin{array}{lcl}
#    \dot{x} &=& \dfrac{-xy}{1+x}\\ \\
#    \dot{y} &=& \dfrac{xy}{1+x} - a y\\ \\
#    \dot{z} &=& ay
#    \end{array} \right.
# ```
# with nominal values ``a = 0.3`` and an initial condition ```x(0) = 9.98``, ``y(0) = 0.01`` and
# ``z(0) = 0.01``. In this model, ``x`` is the nutrients, ``y`` the phytoplankton
# and ``z`` the detritus. In this notebook we analyze the behavior of the system
# subject to uncertainties in the initial condition and/or the model's parameter `a`
# using forward reachability analysis.
#
# The solution should satisfy the following two constraints:
#
# - **Positivity:** ``x(t)``, ``y(t)``, ``z(t)`` are positive, and
# - **Conservativity:** ``x(t) + y(t) + z(t) = 10`` for all ``t \ge 0``.

# In this notebook we are interested in computing the reachable tube until the time horizon ``T=100``,
# and to verify (or embed) the  constraints. Three setups are considered, depending on the source
# of bounded uncertainties which is taken into account:
#
# - **Case I**: ``x(0) \in [9.5, 10.0]``, i.e., uncertainty on the initial condition;
# - **Case P**: ``a \in [0.296, 0.304]``, i.e., uncertainty on the parameter;
# - **Case I & P**: ``x(0) \in [9.5, 10.0]`` and ``a \in [0.296, 0.304]``, i.e., both uncertainties are mixed.

# In terms of objectives, at ``T = 100``, the constraints ``10 \in x+y+z`` and
# ``x, y, z \geq 0`` have to be verified.

# The reachability settings considered above are a variation of those in [^ARCHCOMP20_NLN].
# Please note that variables $x$ and $y$ converge towards zero, so the final volume can be used
# as a quality measure of the overapproximation.

# In the next section we explain how to check whether the constraints are satisfied
# using `LazySets.jl`.

# ## Checking satisfiability with sets

# Two measures of quality of the approximation are considered: the volume of the box
# (``x \times y \times z``) enclosing the final state (at ``T = 100``) and the total
# time of computation for evolution and verification.
# All of these results are obtained for the three setups.

# Let us discuss the implementation of each constraint satisfaction problem.
# We begin by defining the positive orthant in three-dimensional space,
# `positive_orthant`, as an unbounded polyhedron in constraint representation.

using ReachabilityAnalysis, Symbolics, Plots

@variables x y z
const positive_orthant = HPolyhedron([x >= 0, y >= 0, z >= 0], [x, y, z])

# Given a set $X \subseteq \mathbb{R}^n$, to check whether the positivity constraint
# holds correspnods to checking wheter ``X`` is included in the positive orthant.
# This computation can be done efficiently using support functions, and it is available
# in `LazySets.jl`. Multiple dispatch takes care based on the types of the arguments
# in the call `X ⊆ positive_orthant` depending on the type of ``X``.

# On the other hand, given a set ``X``, a quick way to check the conservativity constraint
# is to first overapproximate ``X`` with a box, then represent such box as a
# product-of-intervals (`IntervalBox`) `B`, and finally take the Minkowski
# sum of each line segment, using `sum(B)`. This check is only sufficient; splitting
# ``X`` can be used to refine the check.

# The verification problem is summarized in the function `prod_dest_verif`.
# It receives the solution of a reachability problem represented with
# a Taylor model flowpipe and two optional arguments that specify the time horizon
# `T` and the conservativity condition given a target state `target`,
# which defaults to `10`. The function returns the tuple `(flag, vol)` where `flag`
# is true if and only if the positivity and conservativity conditions are satisfied at the time horizon
# `T`, which defaults to `100`. The second output, `vol`, corresponds to the volume
# of the box overapproximation of the final reach-set.

function prod_dest_verif(sol; T=100.0, target=10.0)

    ## convert to a zonotopic flowpipe
    solz = overapproximate(sol, Zonotope)

    ## project the final reach-set onto the space variables x, y, z
    X = project(solz(T), vars=(1, 2, 3))

    ## check that all variables are nonnegative
    nonnegative = X ⊆ positive_orthant

    ## compute the volume of the last reach-set
    H = overapproximate(X, Hyperrectangle)
    vol = volume(H)

    ## check that that target belongs to the minkowski sum of the reach-sets projected in each coordinate
    B = convert(IntervalBox, H) # get the product-of-intervals representation
    contains_target = target ∈ sum(B)

    return nonnegative && contains_target, vol
end

# ## Case I

# Case I corresponds to initial states uncertainty. We begin by writing the system
# of ODEs in the function `prod_dest_I!`.

@taylorize function prod_dest_I!(du, u, params, t)
    local a = 0.3
    x, y, z = u[1], u[2], u[3]

    du[1] = - (x * y) / (1 + x)
    du[2] = (x * y) / (1 + x) - a * y
    du[3] = a * y
    return du
end

# We define the initial states as the cartesian product of intervals.

X0 = (9.5 .. 10.0) × (0.01 .. 0.01) × (0.01 .. 0.01)
prob = @ivp(x'= prod_dest_I!(x), dim:3, x(0) ∈ X0)

solI = solve(prob, T=100.0, alg=TMJets(abstol=1e-11, orderT=7, orderQ=1));

# Verifying that the specification holds:
property, vol = prod_dest_verif(solI)

# Now we plot ``z`` (in the ``[0, 11]`` range) w.r.t. time (in the ``[0, 100]`` range).

plot(solI, vars=(0, 3), linecolor=:orange, color=:orange, alpha=0.3, lab="I")

# ## Case P

# For the case of an uncertain parameter, we add a new state variable that corresponds
# to `a`, with constant (zero) dynamics.

@taylorize function prod_dest_IP!(du, u, params, t)
    x, y, z, a = u[1], u[2], u[3], u[4]

    du[1] = - (x * y) / (1 + x)
    du[2] = (x * y) / (1 + x) - a * y
    du[3] = a * y
    du[4] = zero(x)
    return du
end

X0 = (9.98 .. 9.98) × (0.01 .. 0.01) × (0.01 .. 0.01) × (0.296 .. 0.304)
prob = @ivp(x'= prod_dest_IP!(x), dim:4, x(0) ∈ X0)

solP = solve(prob, T=100.0, alg=TMJets(abstol=1e-12, orderT=7, orderQ=1));

# Verifying that the specification holds:
property, vol = prod_dest_verif(solP)

# Now we plot ``z`` (in the ``[0, 11]`` range) w.r.t. time (in the ``[0, 100]`` range).

plot(solP, vars=(0, 3), linecolor=:blue, color=:blue, alpha=0.3, lab="P")

# ## Case I & P

# When uncertainty in both the intial states and the paramters are present, we can
# reuse the function `prod_dest_IP!`, but setting an uncertain initial condition
# and an uncertain parameter.
# Recall that we are interested in ``x(0) \in [9.5, 10.0]`` and ``a \in [0.296, 0.304]``.

X0 = (9.5 .. 10.0) × (0.01 .. 0.01) × (0.01 .. 0.01) × (0.296 .. 0.304)
prob = @ivp(x'= prod_dest_IP!(x), dim:4, x(0) ∈ X0)

solIP = solve(prob, T=100.0, alg=TMJets(abstol=1e-11, orderT=7, orderQ=1));

# Verifying that the specification holds:
property, vol = prod_dest_verif(solIP)

# Now we plot ``z`` (in the ``[0, 11]`` range) w.r.t. time (in the ``[0, 100]`` range).

plot(solIP, vars=(0, 3), linecolor=:red, color=:red, alpha=0.3, lab="I & P")

# ## Performance optimization

# We have observed that the specifications are satisfied in all cases. The question remains
# if we can tune the problem implementation such that we can verify these properties in less time.
# In this section we put into practice the ideas explained in [Some common gotchas](@ref).
# First, consider the function `prod_dest_I!` defined above. The term `du[1] = -x*y / (1+x)`
# has the numerator and denominator in common with respect to `du[2]`, so we can store
# those values in an auxiliary variable. Similar comments apply to the function
# with varying parameter `a`. The final functions are presented below.

@taylorize function prod_dest_I_optimized!(du, u, params, t)
    local a = 0.3
    x, y, z = u[1], u[2], u[3]

    num = x * y
    den = 1 + x
    aux = num/den
    aux2 = a * y
    du[1] = -aux
    du[2] = aux - aux2
    du[3] = aux2
    return du
end

@taylorize function prod_dest_IP_optimized!(du, u, params, t)
    x, y, z, a = u[1], u[2], u[3], u[4]

    num = x * y
    den = 1 + x
    aux = num/den
    aux2 = a * y
    du[1] = -aux
    du[2] = aux - aux2
    du[3] = aux2
    du[4] = zero(x)
    return du
end

# With respect to parameter tuning, we have observed that for all three cases
# a good parameter choice is ``n_T = 7``, ``n_Q = 1``, and an adaptive absolute
# tolerance with initial value ``10^{-11}`` (resp ``10^{-12}``) for ``I``
# and ``I`` and ``P`` (resp. ``P``). This choice offers a good balance between runtime
# and precision.

# ## References

# [^ARCHCOMP20_NLN]: *ARCH-COMP20 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics.* Luca Geretti, Julien Alexandre dit Sandretto, Matthias Althoff, Luis Benet, Alexandre Chapoutot, Xin Chen, Pieter Collins, Marcelo Forets, Daniel Freire, Fabian Immler, Niklas Kochdumper, David P. Sanders and Christian Schilling (2020) ARCH20. To appear in 7th International Workshop on Applied Verification of Continuous and Hybrid Systems. 7th International Workshop on Applied Verification of Continuous and Hybrid Systems (ARCH20), vol 74, pages 49--75. 10.29007/zkf6.

# [^KM18]: Kopecz, Stefan, and Andreas Meister. *On order conditions for modified Patankar–Runge–Kutta schemes.* Applied Numerical Mathematics 123 (2018): 159-179.


text:
mutable struct Block <: AbstractBlock
  piece::Piece{Vector{Block}}
  # no children

  x::Integer
  y::Integer
end


text:
module Solution
  export opposite
  opposite(number) = -number
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# Local IDW estimator using local anisotropies


text:
function plotmarginaline(fluxobs, posterior, σ2)

  linemeansigma = marginaline(posterior, σ2)

  # project flux obs on line to find minimum, maximum t for plotting purposes

  D = round(Int, length(posterior)/2)

  ā, b̄ = mean(posterior)[1:D], mean(posterior)[D+1:2*D]

  proj(x) = dot(ā, x-b̄) / dot(ā, ā)

  tproj = [proj(fluxobs[:,i]) for i=1:size(fluxobs, 2)]

  tmin, tmax = minimum(tproj), maximum(tproj)

  tmin, tmax = sign(tmin) * abs(tmin * 1.2), sign(tmax) * abs(tmax * 1.2)

  trange = collect(LinRange(tmin, tmax, 100))

  predictions     = map(linemeansigma, trange)
  meanpredictions = reduce(hcat, [p[1] for p in predictions])
  covpredictions  = reduce(hcat, [diag(p[2]) for p in predictions])

  figure()

  # verification
  # plot3D(ā[1].*trange .+ b̄[1], ā[2].*trange .+ b̄[2], ā[3].*trange .+ b̄[3], "c-", label="mean fit")
  # map(x -> plot3D(x...,"k."), meanpredictions)

  plot3D(fluxobs[1,:], fluxobs[2,:], fluxobs[3,:], "bo", label="observed total flux")

  plot3D(vec(meanpredictions[1,:]), vec(meanpredictions[2,:]), vec(meanpredictions[3,:]), "k-", label="mean fit")

  plot3D(vec(meanpredictions[1,:]) .+ 2*vec(sqrt.(covpredictions[1,:])),
    vec(meanpredictions[2,:]),
    vec(meanpredictions[3,:]), "k--", label="std2")

  plot3D(vec(meanpredictions[1,:]),
    vec(meanpredictions[2,:] .+ 2*vec(sqrt.(covpredictions[2,:]))),
    vec(meanpredictions[3,:]), "k--", label="std2")

  plot3D(vec(meanpredictions[1,:]),
    vec(meanpredictions[2,:]),
    vec(meanpredictions[3,:]) .+ 2*vec(sqrt.(covpredictions[3,:])), "k--", label="std2")

end


text:
import UnicodePlots
import RecipesBase

function UnicodePlots.heatmap(
    field::Fields.SpectralElementField2D;
    width = 80,
    height = 40,
    kwargs...,
)
    if !(eltype(field) <: Real)
        error("Can only plot heatmaps of scalar fields")
    end

    space = axes(field)
    mesh = space.topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2

    Nu = max(div(width, n1), div(height, n2))
    M = Operators.matrix_interpolate(field, Nu)
    m1, m2 = size(M)
    domain = Meshes.domain(mesh)
    x1min = Geometry.component(domain.x1x2min, 1)
    x2min = Geometry.component(domain.x1x2min, 2)
    x1max = Geometry.component(domain.x1x2max, 1)
    x2max = Geometry.component(domain.x1x2max, 2)

    CT = Domains.coordinate_type(domain)
    X1CT = Geometry.coordinate_type(CT, 1)
    X2CT = Geometry.coordinate_type(CT, 2)
    X1CTname = Base.typename(X1CT).name
    X2CTname = Base.typename(X1CT).name

    UnicodePlots.heatmap(
        M',
        xlabel = "$(X1CTname)",
        ylabel = "$(X2CTname)",
        xoffset = x1min,
        xscale = (x1max - x1min) / (m1 - 1),
        yoffset = x2min,
        yscale = (x2max - x2min) / (m2 - 1),
        width = width,
        height = height,
        kwargs...,
    )
end

function UnicodePlots.lineplot(
    field::Fields.FiniteDifferenceField;
    name::Union{Nothing, Symbol} = nothing,
    width = 80,
    height = 40,
    kwargs...,
)
    if name !== nothing
        field = getproperty(field, name)
    else
        name = :x
    end
    space = axes(field)

    xlabel = repr(name) * " value"
    xdata = Array(parent(field))[:, 1]

    ydata = Array(parent(Spaces.coordinates_data(space)))[:, 1]
    ylabel = if field isa Spaces.FaceFiniteDifferenceSpace
        ":y faces"
    else
        ":y centers"
    end
    # fix the ylim to the column space (domain)
    ylim = extrema(ydata)
    UnicodePlots.lineplot(
        xdata,
        ydata;
        title = "Column",
        xlabel = xlabel,
        ylabel = ylabel,
        ylim = [ylim[1], ylim[2]],
        width = width,
        height = height,
        kwargs...,
    )
end

RecipesBase.@recipe function f(field::Fields.FiniteDifferenceField)
    # unwrap the data to plot
    space = axes(field)
    CT = Meshes.coordinate_type(space.mesh)
    CTName = Base.typename(CT).name

    xdata = parent(field)[:, 1]
    ydata = parent(Spaces.coordinates_data(space))[:, 1]

    # set the plot attributes
    title --> "Column"
    xguide --> ":x value"
    yguide --> (
        if field isa Spaces.FaceFiniteDifferenceSpace
            "$(CTname) faces"
        else
            "$(CTName) centers"
        end
    )

    # fix the ylim to the column space (domain)
    ylims --> extrema(ydata)
    linewidth --> 2
    legend --> false

    (xdata, ydata)
end

RecipesBase.@recipe function f(
    field::Fields.SpectralElementField2D;
    interpolate = 10,
)
    @assert interpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    # compute the interpolated data to plot
    space = axes(field)
    mesh = space.topology.mesh
    n1 = mesh.n1
    n2 = mesh.n2

    Nu = interpolate
    M = Operators.matrix_interpolate(field, Nu)

    domain = Meshes.domain(mesh)
    x1min = Geometry.component(domain.x1x2min, 1)
    x2min = Geometry.component(domain.x1x2min, 2)
    x1max = Geometry.component(domain.x1x2max, 1)
    x2max = Geometry.component(domain.x1x2max, 2)

    r1 = range(x1min, x1max, length = n1 * Nu + 1)
    r1 = r1[1:(end - 1)] .+ step(r1) ./ 2
    r2 = range(x2min, x2max, length = n2 * Nu + 1)
    r2 = r2[1:(end - 1)] .+ step(r2) ./ 2

    coord_symbols = propertynames(Fields.coordinate_field(space))

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (r1, r2, M')
end

RecipesBase.@recipe function f(
    field::Fields.ExtrudedFiniteDifferenceField;
    hinterpolate = 0,
)
    data = Fields.field_values(field)
    Ni, Nj, _, Nv, Nh = size(data)

    space = axes(field)
    hmesh = Spaces.topology(space).mesh
    vmesh = space.vertical_mesh

    #TODO: assumes VIFH layout
    @assert Nj == 1 "plotting only defined for 1D extruded fields"

    coord_symbols = propertynames(Fields.coordinate_field(space))
    hcoord_field = getproperty(Fields.coordinate_field(space), 1)
    vcoord_field = getproperty(Fields.coordinate_field(space), 2)

    if hinterpolate ≥ 1
        Nu = hinterpolate
        M_hcoord = Operators.matrix_interpolate(hcoord_field, Nu)
        M_data = Operators.matrix_interpolate(field, Nu)

        hcoord = vec(M_hcoord[1, :])
        vcoord = vec(parent(vcoord_field)[:, 1, 1, 1])
        data = M_data
    else
        hcoord = vec(parent(hcoord_field)[1, :, 1, :])
        vcoord = vec(parent(vcoord_field)[:, 1, 1, 1])
        data = reshape(parent(data), (Nv, Ni * Nh))
    end

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (hcoord, vcoord, data)
end

function play(
    timesteps::Vector;
    name::Union{Nothing, Symbol} = nothing,
    fps::Real = 5,
    width = 80,
    height = 40,
    paused = true,
)
    #TODO: compute extrema over timetep fields
    io_w, io_h = width, height
    nframes = length(timesteps)

    setraw!(io, raw) =
        ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid}, Int32), io.handle, raw)
    ansi_moveup(n::Int) = string("\e[", n, "A")
    ansi_movecol1 = "\e[1G"
    ansi_enablecursor = "\e[?25h"
    ansi_disablecursor = "\e[?25l"

    frame = lastframe = 1
    finished = false
    first_print = true
    actual_fps = 0

    # column y limits are consistent based on the space
    # compute the extrema of over field values for all timesteps

    #= TODO: need to implement field iteration
    field_limits = [extrema(timesteps[1])...]
    if length(timesteps) > 1
        for t in timesteps[2:end]
            tlim = extrema(t)
     if tlim[1] < xlim[1]
    field_limits[1] = tlim[1]
     end
     if tlim[2] > xlim[2]
    field_limits[2] = tlim[2]
     end
    end
    end
    =#

    keytask = @async begin
        try
            setraw!(stdin, true)
            while !finished
                keyin = read(stdin, Char)
                if UInt8(keyin) == 27
                    keyin = read(stdin, Char)
                    if UInt8(keyin) == 91
                        keyin = read(stdin, Char)
                        # left & up arrows
                        if UInt8(keyin) in [68, 65]
                            frame = frame <= 1 ? 1 : frame - 1
                        end
                        # right & down arrows
                        if UInt8(keyin) in [67, 66]
                            frame = frame >= nframes ? nframes : frame + 1
                        end
                    end
                end
                # seek to some pct of simulation time (1 start, 0 end)
                if keyin in '0':'9'
                    k = parse(Int, keyin)
                    pct = (k == 0 ? 10 : k) / 10.0
                    frame = round(Int, pct * nframes)
                end
                # p or spacebar to pause
                if keyin in ['p', ' ']
                    paused = !paused
                end
                # esc or q to quite
                if keyin in ['\x03', 'q']
                    finished = true
                end
            end
        catch
        finally
            setraw!(stdin, false)
        end
    end
    try
        print(ansi_disablecursor)
        setraw!(stdin, true)
        while !finished
            tim = Timer(1 / fps)
            t = @elapsed begin
                if first_print || lastframe != frame
                    field = timesteps[frame]
                    field_plot = if field isa Fields.FiniteDifferenceField
                        UnicodePlots.lineplot(
                            field,
                            title = "Column t=$(frame-1)",
                            width = io_w,
                            height = io_h,
                            #xlim = field_limits,
                        )
                    elseif field isa Fields.SpectralElementField2D
                        UnicodePlots.heatmap(field, width = io_w, height = io_h)
                    else
                        error("unknown field type: $(summary(field))")
                    end
                    # TODO: dig through structure to just compute the number of lines directly
                    nlines = countlines(IOBuffer(string(field_plot)))
                    str = sprint(; context = (:color => true)) do ios
                        println(ios, field_plot)
                        if paused
                            println(ios, "Frame: $frame/$nframes", " "^15)
                        else
                            println(
                                ios,
                                "Frame: $frame/$nframes FPS: $(round(actual_fps, digits=1))",
                                " "^5,
                            )
                        end
                    end
                    first_print ? print(str) :
                    print(ansi_moveup(nlines + 1), ansi_movecol1, str)
                end
                first_print = false
                lastframe = frame
                if !paused
                    if frame == nframes
                        frame = 1
                    else
                        frame += 1
                    end
                end
                wait(tim)
            end
            actual_fps = 1 / t
        end
    catch e
        isa(e, InterruptException) || rethrow()
    finally
        print(ansi_enablecursor)
        finished = true
        @async Base.throwto(keytask, InterruptException())
        wait(keytask)
    end
    return
end


text:
using SafeTestsets
const LONGER_TESTS = false

@time @safetestset "Stan_String" begin include("test_stan_string.jl") end
@time @safetestset "Stan" begin include("test_one_parameter_stan.jl") end
@time @safetestset "Stan" begin include("test_four_parameter_stan.jl") end


text:
import FileIO

include("model.jl")

function load_params!(gen_func, filename)
    data = FileIO.load(filename)
    for (name, param) in data
        gen_func.params[Symbol(name)] = param
    end
end

function do_inference(input_image)

    # construct trace containing observed image
    observations = DynamicAssignment()
    observations["image"] = input_image

    # fill in latent variables using proposal
    latents = get_assignment(simulate(proposal, (observations,)))
    x = latents["x"]
    y = latents["y"]
    s = latents["size"]
    letter = latents["letter"]
    angle = latents["angle"]
    println("x: $x, y: $y, size: $s, letter: $letter, angle: $angle")

    # predicted image given latents
    (trace, _) = generate(model, (), latents)
    predicted = get_call_record(trace).retval
    predicted =  min.(fill(1, size(predicted)), max.(zero(predicted), predicted))
end

load_params!(proposal, "params.jld2")

println("do inference..")
for i=1:100
    input_image = convert(Matrix{Float64}, load("prior/img-029.png"))
    predicted = do_inference(input_image)
    output_filename = @sprintf("proposed/%03d.png", i)
    FileIO.save(output_filename, colorview(Gray, predicted))
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
"""
    SrcDir(; file=".ion/templates/package/src/module.jl")

Creates a module entrypoint.

## Keyword Arguments
- `file::AbstractString`: Template file for `src/<module>.jl`.
"""
@plugin mutable struct SrcDir <: FilePlugin
    file::String = default_file("src", "module.jl")
    destination::String = ""
end

source(p::SrcDir) = p.file
destination(p::SrcDir) = p.destination
view(::SrcDir, ::PackagePlan, pkg::AbstractString) = Dict("PKG" => pkg)

# Update the destination now that we know the package name.
# Kind of hacky, but oh well.
function prehook(p::SrcDir, ::PackagePlan, pkg_dir::AbstractString)
    p.destination = joinpath("src", basename(pkg_dir) * ".jl")
end


text:
struct Atoms

    total_N::Array{Int64,1}
    local_N::Array{Int64,1}

    rho0::Array{Float64,1}

    num_atomtype::Array{Int64,1}

    types::Array{Int64,1}
    types_all::Array{Int64,1}

    masses::Array{Float64,1}
    masses_all::Array{Float64,1}

    coords::Vector{Array{Float64,1}}
    coords_all::Vector{Array{Float64,1}}

    vels::Vector{Array{Float64,1}}
    vels_all::Vector{Array{Float64,1}}

    forces::Vector{Array{Float64,1}}
    forces_all::Vector{Array{Float64,1}}

    energies::Array{Float64,1}
    energies_all::Array{Float64,1}

    bonds::Vector{Array{Int64,1}}
    bonds_all::Vector{Array{Int64,1}}

    function Atoms()
        total_N=[0]
        local_N=[0]
        rho0=[0.0]
        num_atomtype=[0]

        types=Vector{Int64}(undef,0)
        types_all=Vector{Int64}(undef,0)

        masses=Vector{Float64}(undef,0)
        masses_all=Vector{Float64}(undef,0)

        coords=Array{Array{Float64,1},1}(undef,0)
        coords_all=Array{Array{Float64,1},1}(undef,0)

        vels=Array{Array{Float64,1},1}(undef,0)
        vels_all=Array{Array{Float64,1},1}(undef,0)

        forces=Array{Array{Float64,1},1}(undef,0)
        forces_all=Array{Array{Float64,1},1}(undef,0)

        energies=Vector{Float64}(undef,0)
        energies_all=Vector{Float64}(undef,0)

        bonds=Vector{Vector{Int64}}(undef,0)
        bonds_all=Vector{Vector{Int64}}(undef,0)

        new(total_N::Array{Int64,1},
            local_N::Array{Int64,1},
            rho0::Array{Float64,1},
            num_atomtype::Array{Int64,1},
            types::Array{Int64,1},
            types_all::Array{Int64,1},

            masses::Array{Float64,1},
            masses_all::Array{Float64,1},

            coords::Vector{Array{Float64,1}},
            coords_all::Vector{Array{Float64,1}},

            vels::Vector{Array{Float64,1}},
            vels_all::Vector{Array{Float64,1}},

            forces::Vector{Array{Float64,1}},
            forces_all::Vector{Array{Float64,1}},

            energies::Array{Float64,1},
            energies_all::Array{Float64,1},

            bonds::Vector{Array{Int64,1}},
            bonds_all::Vector{Array{Int64,1}})
    end
end

function init_atoms!(atoms::Atoms,
                    types_all::Array{Int64,1},
                    masses_all::Array{Float64,1},
                    coords_all::Array{Float64,2},
                    forces_all::Array{Float64,2},
                    vels_all::Array{Float64,2},
                    energies_all::Array{Float64,1})
    for i in 1:atoms.total_N[1]
        push!(atoms.types_all,types_all[i])
        push!(atoms.masses_all,masses_all[i])
        push!(atoms.coords_all,coords_all[i,1:3])
        push!(atoms.forces_all,forces_all[i,1:3])
        push!(atoms.vels_all,vels_all[i,1:3])
        push!(atoms.energies_all,energies_all[i])
    end
end

function init_bonds!(atoms::Atoms,
                    bonds_all::Array{Int64,2})
    for bondi in 1:size(bonds_all)[1]
        push!(atoms.bonds_all,bonds_all[bondi,1:3])
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.332774675925122 <= q <= 6.4153373403467775)

                     add_NL_constraint(m, :(sigmoid(0.1742606373223392 + -0.05222015296022331*$(x[1]) + 0.15266538164806986*$(x[2]) + 0.27567618482533085*$(x[3])) + sigmoid(0.7992746561529409 + 0.8053279528303205*$(x[1]) + -0.8864483595131469*$(x[2]) + 0.6110618799366847*$(x[3])) + sigmoid(-0.6752579159956698 + 0.27352401436271423*$(x[1]) + 0.2855995340868476*$(x[2]) + -0.33283297885270313*$(x[3])) + sigmoid(0.2430039547312175 + -0.8377319978555251*$(x[1]) + 0.39266174748446403*$(x[2]) + -0.9683058237799194*$(x[3])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
mutable struct CanonicalModel
    JuMPmodel::JuMP.AbstractModel
    variables::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    constraints::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    cost_function::JuMP.AbstractJuMPScalar
    expressions::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    parameters::Union{Nothing,Dict{Symbol, JuMP.Containers.DenseAxisArray}}
    initial_conditions::Dict{Symbol, Array{PSI.InitialCondition}}
    pm_model::Union{Nothing,PM.GenericPowerModel}
end

_variable_type(cm::CanonicalModel) = JuMP.variable_type(cm.JuMPmodel)


text:
@params struct JuniperIpoptOptions
    nt::NamedTuple
    subsolver_options
    first_order::Bool
end
function JuniperIpoptOptions(;
    first_order = true,
    linear_constraints = false,
    subsolver_options = IpoptOptions(
        first_order = first_order,
        linear_constraints = linear_constraints,
    ),
    kwargs...,
)
    first_order = !hasproperty(subsolver_options.nt, :hessian_approximation) ||
        subsolver_options.nt.hessian_approximation == "limited-memory"
    return JuniperIpoptOptions((; kwargs...), subsolver_options, first_order)
end

@params mutable struct JuniperIpoptWorkspace <: Workspace
    model::VecModel
    problem::JuMPProblem
    x0::AbstractVector
    integers::AbstractVector{<:Bool}
    options::JuniperIpoptOptions
    counter::Base.RefValue{Int}
end
function JuniperIpoptWorkspace(
    model::VecModel, x0::AbstractVector = getinit(model);
    options = JuniperIpoptOptions(), kwargs...,
)
    integers = model.integer
    @assert length(integers) == length(x0)
    nt1 = options.subsolver_options.nt
    subsolver_options = map(keys(nt1)) do k
        string(k) => nt1[k]
    end
    nl_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, Dict(subsolver_options)...)
    nt2 = options.nt
    solver_options = map(keys(nt2)) do k
        string(k) => nt2[k]
    end
    optimizer = JuMP.optimizer_with_attributes(
        Juniper.Optimizer, "nl_solver" => nl_solver, Dict(solver_options)...,
    )
    problem, counter = get_jump_problem(
        model, copy(x0); first_order = options.first_order,
        optimizer = optimizer, integers = integers,
    )
    return JuniperIpoptWorkspace(model, problem, copy(x0), integers, options, counter)
end
@params struct JuniperIpoptResult <: AbstractResult
    minimizer
    minimum
    problem
    status
    fcalls::Int
end

function optimize!(workspace::JuniperIpoptWorkspace)
    @unpack problem, options, x0, counter = workspace
    counter[] = 0
    jump_problem = workspace.problem
    jump_model = jump_problem.model
    moi_model = jump_model.moi_backend
    MOI.set.(Ref(moi_model), Ref(MOI.VariablePrimalStart()), workspace.problem.vars, x0)
    MOI.optimize!(moi_model)
    minimizer = MOI.get(moi_model, MOI.VariablePrimal(), jump_problem.vars)
    objval = MOI.get(moi_model, MOI.ObjectiveValue())
    term_status = MOI.get(moi_model, MOI.TerminationStatus())
    primal_status = MOI.get(moi_model, MOI.PrimalStatus())
    return JuniperIpoptResult(
        minimizer, objval, problem, (term_status, primal_status), counter[],
    )
end

struct JuniperIpoptAlg <: AbstractOptimizer end

function Workspace(model::VecModel, optimizer::JuniperIpoptAlg, args...; kwargs...,)
    return JuniperIpoptWorkspace(model, args...; kwargs...)
end


text:
using Test, LinearAlgebra, Distributions, AdvancedHMC, Random
using AdvancedHMC.Adaptation: WelfordVar, NaiveVar, WelfordCov, NaiveCov, add_sample!, get_var, get_cov, reset!
using DiffResults: GradientResult, value, gradient
using ForwardDiff: gradient!

# Check that the estimated variance is approximately correct.
@testset "Online v.s. naive v.s. true var/cov estimation" begin
    D = 10
    n_samples = 100_000

    var_welford = WelfordVar(D)
    var_naive = NaiveVar()
    var_estimators = [var_welford, var_naive]
    cov_welford = WelfordCov(D)
    cov_naive = NaiveCov()
    cov_estimators = [cov_welford, cov_naive]
    estimators = [var_estimators..., cov_estimators...]

    for dist in [MvNormal(zeros(D), ones(D)), Dirichlet(D, 1)]
        for _ = 1:n_samples
            s = rand(dist)
            for estimator in estimators
                add_sample!(estimator, s)
            end
        end

        @test get_var(var_welford) ≈ get_var(var_naive) atol=0.1D
        for estimator in var_estimators
            @test get_var(estimator) ≈ var(dist) atol=0.1D
        end

        @test get_cov(cov_welford) ≈ get_cov(cov_naive) atol=0.1D^2
        for estimator in cov_estimators
            @test get_cov(estimator) ≈ cov(dist) atol=0.1D^2
        end

        for estimator in estimators
            reset!(estimator)
        end
    end
end

@testset "Preconditioner constructors" begin
    θ = [0.0, 0.0, 0.0, 0.0]
    pc1 = Preconditioner(UnitEuclideanMetric) # default dim = 2
    pc2 = Preconditioner(DiagEuclideanMetric)
    pc3 = Preconditioner(DenseEuclideanMetric)

    # Var adaptor dimention should be increased to length(θ) from 2
    AdvancedHMC.adapt!(pc1, θ, 1.)
    AdvancedHMC.adapt!(pc2, θ, 1.)
    AdvancedHMC.adapt!(pc3, θ, 1.)
    @test AdvancedHMC.Adaptation.getM⁻¹(pc2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(pc3) == LinearAlgebra.diagm(0 => ones(length(θ)))
end

@testset "Stan HMC adaptors" begin
    θ = [0.0, 0.0, 0.0, 0.0]

    adaptor1 = StanHMCAdaptor(
        1000,
        Preconditioner(UnitEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor2 = StanHMCAdaptor(
        1000,
        Preconditioner(DiagEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor3 = StanHMCAdaptor(
        1000,
        Preconditioner(DenseEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )

    AdvancedHMC.adapt!(adaptor1, θ, 1.)
    AdvancedHMC.adapt!(adaptor2, θ, 1.)
    AdvancedHMC.adapt!(adaptor3, θ, 1.)
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor3) == LinearAlgebra.diagm(0 => ones(length(θ)))
end

let D=10
    function buildfuncs(target)
        ℓπ(θ) = logpdf(target, θ)

        function ∂ℓπ∂θ(θ)
            res = GradientResult(θ)
            gradient!(res, ℓπ, θ)
            return (value(res), gradient(res))
        end

        return ℓπ, ∂ℓπ∂θ
    end

    function runnuts(ℓπ, ∂ℓπ∂θ, metric; n_samples=2_000)
        n_adapts = 1_000

        θ_init = randn(D)

        h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
        prop = NUTS(Leapfrog(find_good_eps(h, θ_init)))
        adaptor = StanHMCAdaptor(
            n_adapts,
            Preconditioner(metric),
            NesterovDualAveraging(0.8, prop.integrator)
        )
        samples, stats = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; verbose=false)
        return (samples=samples, stats=stats, adaptor=adaptor)
    end

    @testset "Adapted mass v.s. true variance" begin
        Random.seed!(123)
        n_tests = 5

        @testset "DiagEuclideanMetric" begin
            for _ in 1:n_tests
                # Random variance
                σ = ones(D) + abs.(randn(D))

                # Diagonal Gaussian
                target = MvNormal(zeros(D), σ)
                ℓπ, ∂ℓπ∂θ = buildfuncs(target)

                res = runnuts(ℓπ, ∂ℓπ∂θ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ σ .^ 2 rtol=0.2

                res = runnuts(ℓπ, ∂ℓπ∂θ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.covar ≈ diagm(0 => σ.^2) rtol=0.25
            end
        end

        @testset "DenseEuclideanMetric" begin
            for _ in 1:n_tests
                # Random covariance
                m = randn(D, D)
                Σ = m' * m

                # Correlated Gaussian
                target = MvNormal(zeros(D), Σ)
                ℓπ, ∂ℓπ∂θ = buildfuncs(target)

                res = runnuts(ℓπ, ∂ℓπ∂θ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ diag(Σ) rtol=0.2

                res = runnuts(ℓπ, ∂ℓπ∂θ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.covar ≈ Σ rtol=0.25
            end
        end

    end

    @testset "Initialisation adaptor by metric" begin
        target = MvNormal(zeros(D), 1)
        ℓπ, ∂ℓπ∂θ = buildfuncs(target)

        mass_init = fill(0.5, D)
        res = runnuts(ℓπ, ∂ℓπ∂θ, DiagEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.var == mass_init

        mass_init = diagm(0 => fill(0.5, D))
        res = runnuts(ℓπ, ∂ℓπ∂θ, DenseEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.covar == mass_init
    end
end


text:
# This file was generated, do not modify it. # hide
res = KnnResult(3) # allocates memory for 10 nearest neighbors
for v in rand(db, 10)
    reuse!(res)  # reuses the res object
    @time search(G, v, res)
    @show minimum(res), maximum(res), argmin(res), argmax(res)
    @show res.id, res.dist
end

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# from a suggestion at https://discourse.julialang.org/t/install-a-packages-from-a-list/30920/3

using Pkg, Distributed

dependencies = [
    "JLD",
    "LoopVectorization",
    "Optim",
    "NLsolve",
    "StatsBase",
    "Distributions",
    "Plots",
    "Revise"
]

Pkg.add(dependencies)

# dynamically start up cores if not done by command line
if nprocs() == 1
  addprocs(length(Sys.cpu_info()) ÷ 2 - 1)
end

max_stars_repo_name: irudik/docker-julia-starter

max_issues_repo_name: irudik/docker-julia-starter

max_forks_repo_name: irudik/docker-julia-starter


text:
using OptimNewton
using Test

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


text:
using Test
include("../../utils/distancias.jl")
include("../../utils/validation.jl")
include("../../naive-classifiers/constant.jl")
include("../../learner/euclidean-1-NN.jl")

data = [-1 0; 1 1; -1 -3; 1 -2; 9 -1; 1 2; 10 20]
class = ['-', '+',  '-',   '-',  '-', '+' , '+' ]

@testset "Función 1-NN" begin
    @test OneNN([2,2], EuclideanDistance, data, class ) == '+'
    @test OneNN([-1,-3], EuclideanDistance, data, class ) == '-'

    oneNN_from_learner = LearnerEuclideanOneNN(data, class)
    @test oneNN_from_learner([2,2]) == '+'
    @test oneNN_from_learner([-1,-3]) == '-'

end

# Distribución de índices
# tamaño 7, tomado en 3 folds [3,5, 7], 123 45 67
l1 = [3,5,7]
l1_test = [CalculateIndex(3,7,i) for i in 1:3 ]
# tamaño 8, tomado en 3 folds [3,6,8] 123 456 78
l2 = [3,6,8]
l2_test = [CalculateIndex(3,8,i) for i in 1:3 ]
# tamaño 9, tomado en 4 folds [3,5, 7, 9] 123 45 67 89
l3 = [3,5,7, 9]
l3_test = [CalculateIndex(4,9,i) for i in 1:4 ]
@testset "Index distribution" begin
    @test l1 == l1_test
    @test l2 == l2_test
    @test l3 == l3_test
end

@testset "Naive Classifiers" begin
    cte = ConstantLearned([1,1,1], [1,2,3])
    @test cte(123) == 1
end

@testset "Leave one out" begin
    data = [-1 0; 1 1; -1 -3; 1 -2; 9 -1; 1 2; 10 20]
    class = ['-', '+',  '-',   '-',  '-', '+' , '+' ]
    accuracy = LeaveOneOut(
        data,
        class,
       LearnerEuclideanOneNN
    )
    @test 0<= accuracy <=100
end


text:
struct Backend{backend} end

Backend(b::Backend) = b
Backend(s::AbstractString) = Backend{Symbol(s)}()
Backend(s::Symbol) = Backend{s}()
Backend(s::Nothing) = Backend{s}()
backend(::Backend{N}) where {N} = N
Backend() = Backend{Symbol()}()

macro Backend_str(s)
  return Backend{Symbol(s)}
end

const current_backend = Ref{Union{Nothing,Backend}}(nothing)

visualize(::Backend{nothing}, args...; kwargs...) = nothing

set_backend!(::Nothing) = (current_backend[] = nothing)
function set_backend!(backend::Backend)
  original_backend = current_backend[]
  current_backend[] = backend
  return original_backend
end
set_backend!(backend::Union{Symbol,String}) = set_backend!(Backend(backend))

get_backend() = isnothing(current_backend[]) ? default_backend() : current_backend[]

function plot(::Backend{T}, args...; kwargs...) where {T}
  return error("plot not implemented for backend type $T.")
end
function draw_edge!(::Backend{T}, args...; kwargs...) where {T}
  return error("draw_edge! not implemented for backend type $T.")
end
function annotate!(::Backend{T}, args...; kwargs...) where {T}
  return error("annotate! not implemented for backend type $T.")
end

function translate_color(::Backend{T}, color) where {T}
  return error("translate_color not implemented for backend type $T and color $color")
end

point_to_line(v1, v2) = ([v1[1], v2[1]], [v1[2], v2[2]])


text:
export set_inequality_constraints!, set_terminal_inequality_constraints!, set_equality_constraints!, set_terminal_equality_constraints!

function set_inequality_constraints!( model::Model, func::Union{Nothing,Function}; dx::CallbackArg=:autodiff, du::CallbackArg=:autodiff, dp::CallbackArg=:autodiff )
    # Allow the inequality constraints to be disabled
    if isa( func, Nothing )
        model.ineqCon = false

        model.callbacks.hfct_cb     = nothing
        model.callbacks.dhdx_vec_cb = nothing
        model.callbacks.dhdu_vec_cb = nothing
        model.callbacks.dhdp_vec_cb = nothing

        return
    else
        model.ineqCon = true
        model.callbacks.hfct_cb = func
    end

    if isa( dx, Function )
        model.callbacks.dhdx_vec_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( du, Function )
        model.callbacks.dhdu_vec_cb = du
    elseif du == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(du)) for du" ) )
    end

    if isa( dp, Function )
        model.callbacks.dhdp_vec_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end
end

function set_terminal_inequality_constraints!( model::Model, func::Function; dx::CallbackArg=:autodiff, dp::CallbackArg=:autodiff, dT::CallbackArg=:autodiff )
    # Allow the terminal inequality constraints to be disabled
    if isa( func, Nothing )
        model.termIneqCon = false

        model.callbacks.hTfct_cb     = nothing
        model.callbacks.dhTdp_vec_cb = nothing
        model.callbacks.dhTdp_vec_cb = nothing
        model.callbacks.dhTdT_vec_cb = nothing

        return
    else
        model.termIneqCon = true
        model.callbacks.hTfct_cb = func
    end

    if isa( dx, Function )
        model.callbacks.dhTdx_vec_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( dp, Function )
        model.callbacks.dhTdp_vec_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end

    if isa( dT, Function )
        model.callbacks.dhTdT_vec_cb = dT
    elseif dT == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dT)) for dT" ) )
    end
end

function set_equality_constraints!( model::Model, func::Function; dx::CallbackArg=:autodiff, du::CallbackArg=:autodiff, dp::CallbackArg=:autodiff )
    # Allow the equality constraints to be disabled
    if isa( func, Nothing )
        model.eqCon = false

        model.callbacks.gfct_cb     = nothing
        model.callbacks.dgdx_vec_cb = nothing
        model.callbacks.dgdu_vec_cb = nothing
        model.callbacks.dgdp_vec_cb = nothing

        return
    else
        model.eqCon = true
        model.callbacks.gfct_cb = func
    end

    if isa( dx, Function )
        model.callbacks.dgdx_vec_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( du, Function )
        model.callbacks.dgdu_vec_cb = du
    elseif du == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(du)) for du" ) )
    end

    if isa( dp, Function )
        model.callbacks.dgdp_vec_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end
end

function set_terminal_equality_constraints!( model::Model, func::Function; dx::CallbackArg=:autodiff, dp::CallbackArg=:autodiff, dT::CallbackArg=:autodiff )
    # Allow the equality constraints to be disabled
    if isa( func, Nothing )
        model.termEqCon = false

        model.callbacks.gTfct_cb     = nothing
        model.callbacks.dgTdx_vec_cb = nothing
        model.callbacks.dgTdp_vec_cb = nothing
        model.callbacks.dgTdT_vec_cb = nothing

        return
    else
        model.termEqCon = true
        model.callbacks.gTfct_cb = func
    end

    if isa( dx, Function )
        model.callbacks.dgTdx_vec_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( dp, Function )
        model.callbacks.dgTdp_vec_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end

    if isa( dT, Function )
        model.callbacks.dgTdT_vec_cb = dT
    elseif dT == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dT)) for dT" ) )
    end
end


text:
const vertex_codes = [
    "\x15\x45\x51\x54\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA",
    "\x15\x45\x51\x55\x56\x59\x5A\x65\x66\x6A\x95\x96\x9A\xA6\xAA",
    "\x01\x05\x11\x15\x41\x45\x51\x55\x56\x5A\x66\x6A\x96\x9A\xA6\xAA",
    "\x01\x15\x16\x45\x46\x51\x52\x55\x56\x5A\x66\x6A\x96\x9A\xA6\xAA\xAB",
    "\x15\x45\x54\x55\x56\x59\x5A\x65\x69\x6A\x95\x99\x9A\xA9\xAA",
    "\x05\x15\x45\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xAA",
    "\x05\x15\x45\x55\x56\x59\x5A\x66\x6A\x96\x9A\xAA",
    "\x05\x15\x16\x45\x46\x55\x56\x59\x5A\x66\x6A\x96\x9A\xAA\xAB",
    "\x04\x05\x14\x15\x44\x45\x54\x55\x59\x5A\x69\x6A\x99\x9A\xA9\xAA",
    "\x05\x15\x45\x55\x56\x59\x5A\x69\x6A\x99\x9A\xAA",
    "\x05\x15\x45\x55\x56\x59\x5A\x6A\x9A\xAA",
    "\x05\x15\x16\x45\x46\x55\x56\x59\x5A\x5B\x6A\x9A\xAA\xAB",
    "\x04\x15\x19\x45\x49\x54\x55\x58\x59\x5A\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x05\x15\x19\x45\x49\x55\x56\x59\x5A\x69\x6A\x99\x9A\xAA\xAE",
    "\x05\x15\x19\x45\x49\x55\x56\x59\x5A\x5E\x6A\x9A\xAA\xAE",
    "\x05\x15\x1A\x45\x4A\x55\x56\x59\x5A\x5B\x5E\x6A\x9A\xAA\xAB\xAE\xAF",
    "\x15\x51\x54\x55\x56\x59\x65\x66\x69\x6A\x95\xA5\xA6\xA9\xAA",
    "\x11\x15\x51\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\xA5\xA6\xAA",
    "\x11\x15\x51\x55\x56\x5A\x65\x66\x6A\x96\xA6\xAA",
    "\x11\x15\x16\x51\x52\x55\x56\x5A\x65\x66\x6A\x96\xA6\xAA\xAB",
    "\x14\x15\x54\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x99\xA5\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x9A\xA6\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x96\x9A\xA6\xAA\xAB",
    "\x15\x16\x55\x56\x5A\x66\x6A\x6B\x96\x9A\xA6\xAA\xAB",
    "\x14\x15\x54\x55\x59\x5A\x65\x69\x6A\x99\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x9A\xAA",
    "\x15\x16\x55\x56\x59\x5A\x66\x6A\x6B\x9A\xAA\xAB",
    "\x14\x15\x19\x54\x55\x58\x59\x5A\x65\x69\x6A\x99\xA9\xAA\xAE",
    "\x15\x19\x55\x59\x5A\x69\x6A\x6E\x99\x9A\xA9\xAA\xAE",
    "\x15\x19\x55\x56\x59\x5A\x69\x6A\x6E\x9A\xAA\xAE",
    "\x15\x1A\x55\x56\x59\x5A\x6A\x6B\x6E\x9A\xAA\xAB\xAE\xAF",
    "\x10\x11\x14\x15\x50\x51\x54\x55\x65\x66\x69\x6A\xA5\xA6\xA9\xAA",
    "\x11\x15\x51\x55\x56\x65\x66\x69\x6A\xA5\xA6\xAA",
    "\x11\x15\x51\x55\x56\x65\x66\x6A\xA6\xAA",
    "\x11\x15\x16\x51\x52\x55\x56\x65\x66\x67\x6A\xA6\xAA\xAB",
    "\x14\x15\x54\x55\x59\x65\x66\x69\x6A\xA5\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\xA6\xAA",
    "\x15\x16\x55\x56\x5A\x65\x66\x6A\x6B\xA6\xAA\xAB",
    "\x14\x15\x54\x55\x59\x65\x69\x6A\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\xAA",
    "\x15\x16\x55\x56\x59\x5A\x65\x66\x69\x6A\x6B\xAA\xAB",
    "\x14\x15\x19\x54\x55\x58\x59\x65\x69\x6A\x6D\xA9\xAA\xAE",
    "\x15\x19\x55\x59\x5A\x65\x69\x6A\x6E\xA9\xAA\xAE",
    "\x15\x19\x55\x56\x59\x5A\x65\x66\x69\x6A\x6E\xAA\xAE",
    "\x15\x55\x56\x59\x5A\x66\x69\x6A\x6B\x6E\x9A\xAA\xAB\xAE\xAF",
    "\x10\x15\x25\x51\x54\x55\x61\x64\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x11\x15\x25\x51\x55\x56\x61\x65\x66\x69\x6A\xA5\xA6\xAA\xBA",
    "\x11\x15\x25\x51\x55\x56\x61\x65\x66\x6A\x76\xA6\xAA\xBA",
    "\x11\x15\x26\x51\x55\x56\x62\x65\x66\x67\x6A\x76\xA6\xAA\xAB\xBA\xBB",
    "\x14\x15\x25\x54\x55\x59\x64\x65\x66\x69\x6A\xA5\xA9\xAA\xBA",
    "\x15\x25\x55\x65\x66\x69\x6A\x7A\xA5\xA6\xA9\xAA\xBA",
    "\x15\x25\x55\x56\x65\x66\x69\x6A\x7A\xA6\xAA\xBA",
    "\x15\x26\x55\x56\x65\x66\x6A\x6B\x7A\xA6\xAA\xAB\xBA\xBB",
    "\x14\x15\x25\x54\x55\x59\x64\x65\x69\x6A\x79\xA9\xAA\xBA",
    "\x15\x25\x55\x59\x65\x66\x69\x6A\x7A\xA9\xAA\xBA",
    "\x15\x25\x55\x56\x59\x5A\x65\x66\x69\x6A\x7A\xAA\xBA",
    "\x15\x55\x56\x5A\x65\x66\x69\x6A\x6B\x7A\xA6\xAA\xAB\xBA\xBB",
    "\x14\x15\x29\x54\x55\x59\x65\x68\x69\x6A\x6D\x79\xA9\xAA\xAE\xBA\xBE",
    "\x15\x29\x55\x59\x65\x69\x6A\x6E\x7A\xA9\xAA\xAE\xBA\xBE",
    "\x15\x55\x59\x5A\x65\x66\x69\x6A\x6E\x7A\xA9\xAA\xAE\xBA\xBE",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x6B\x6E\x7A\xAA\xAB\xAE\xBA\xBF",
    "\x45\x51\x54\x55\x56\x59\x65\x95\x96\x99\x9A\xA5\xA6\xA9\xAA",
    "\x41\x45\x51\x55\x56\x59\x5A\x65\x66\x95\x96\x99\x9A\xA5\xA6\xAA",
    "\x41\x45\x51\x55\x56\x5A\x66\x95\x96\x9A\xA6\xAA",
    "\x41\x45\x46\x51\x52\x55\x56\x5A\x66\x95\x96\x9A\xA6\xAA\xAB",
    "\x44\x45\x54\x55\x56\x59\x5A\x65\x69\x95\x96\x99\x9A\xA5\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x65\x6A\x95\x96\x99\x9A\xA6\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x66\x6A\x95\x96\x99\x9A\xA6\xAA\xAB",
    "\x45\x46\x55\x56\x5A\x66\x6A\x96\x9A\x9B\xA6\xAA\xAB",
    "\x44\x45\x54\x55\x59\x5A\x69\x95\x99\x9A\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x69\x6A\x95\x96\x99\x9A\xA9\xAA\xAE",
    "\x45\x55\x56\x59\x5A\x6A\x95\x96\x99\x9A\xAA",
    "\x45\x46\x55\x56\x59\x5A\x6A\x96\x9A\x9B\xAA\xAB",
    "\x44\x45\x49\x54\x55\x58\x59\x5A\x69\x95\x99\x9A\xA9\xAA\xAE",
    "\x45\x49\x55\x59\x5A\x69\x6A\x99\x9A\x9E\xA9\xAA\xAE",
    "\x45\x49\x55\x56\x59\x5A\x6A\x99\x9A\x9E\xAA\xAE",
    "\x45\x4A\x55\x56\x59\x5A\x6A\x9A\x9B\x9E\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x56\x59\x65\x66\x69\x95\x96\x99\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x59\x65\x66\x6A\x95\x96\x9A\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x5A\x65\x66\x6A\x95\x96\x9A\xA5\xA6\xAA\xAB",
    "\x51\x52\x55\x56\x5A\x66\x6A\x96\x9A\xA6\xA7\xAA\xAB",
    "\x54\x55\x56\x59\x65\x69\x6A\x95\x99\x9A\xA5\xA6\xA9\xAA",
    "\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA",
    "\x15\x45\x51\x55\x56\x59\x5A\x65\x66\x6A\x95\x96\x9A\xA6\xAA\xAB",
    "\x55\x56\x5A\x66\x6A\x96\x9A\xA6\xAA\xAB",
    "\x54\x55\x59\x5A\x65\x69\x6A\x95\x99\x9A\xA5\xA9\xAA\xAE",
    "\x15\x45\x54\x55\x56\x59\x5A\x65\x69\x6A\x95\x99\x9A\xA9\xAA\xAE",
    "\x15\x45\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xA6\xA9\xAA\xAB\xAE",
    "\x55\x56\x59\x5A\x66\x6A\x96\x9A\xA6\xAA\xAB",
    "\x54\x55\x58\x59\x5A\x69\x6A\x99\x9A\xA9\xAA\xAD\xAE",
    "\x55\x59\x5A\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x55\x56\x59\x5A\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x55\x56\x59\x5A\x6A\x9A\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x65\x66\x69\x95\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x65\x66\x69\x6A\x95\x96\xA5\xA6\xA9\xAA\xBA",
    "\x51\x55\x56\x65\x66\x6A\x95\x96\xA5\xA6\xAA",
    "\x51\x52\x55\x56\x65\x66\x6A\x96\xA6\xA7\xAA\xAB",
    "\x54\x55\x59\x65\x66\x69\x6A\x95\x99\xA5\xA6\xA9\xAA\xBA",
    "\x15\x51\x54\x55\x56\x59\x65\x66\x69\x6A\x95\xA5\xA6\xA9\xAA\xBA",
    "\x15\x51\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x9A\xA5\xA6\xA9\xAA\xAB\xBA",
    "\x55\x56\x5A\x65\x66\x6A\x96\x9A\xA6\xAA\xAB",
    "\x54\x55\x59\x65\x69\x6A\x95\x99\xA5\xA9\xAA",
    "\x15\x54\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x99\x9A\xA5\xA6\xA9\xAA\xAE\xBA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x9A\xA6\xA9\xAA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x96\x9A\xA6\xAA\xAB",
    "\x54\x55\x58\x59\x65\x69\x6A\x99\xA9\xAA\xAD\xAE",
    "\x55\x59\x5A\x65\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\x99\x9A\xA9\xAA\xAE",
    "\x15\x55\x56\x59\x5A\x66\x69\x6A\x9A\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x61\x64\x65\x66\x69\x95\xA5\xA6\xA9\xAA\xBA",
    "\x51\x55\x61\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xB6\xBA",
    "\x51\x55\x56\x61\x65\x66\x6A\xA5\xA6\xAA\xB6\xBA",
    "\x51\x55\x56\x62\x65\x66\x6A\xA6\xA7\xAA\xAB\xB6\xBA\xBB",
    "\x54\x55\x64\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xB9\xBA",
    "\x55\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x55\x56\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x55\x56\x65\x66\x6A\xA6\xAA\xAB\xBA\xBB",
    "\x54\x55\x59\x64\x65\x69\x6A\xA5\xA9\xAA\xB9\xBA",
    "\x55\x59\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x15\x55\x56\x59\x5A\x65\x66\x69\x6A\xA5\xA6\xA9\xAA\xBA",
    "\x15\x55\x56\x5A\x65\x66\x69\x6A\xA6\xAA\xAB\xBA\xBB",
    "\x54\x55\x59\x65\x68\x69\x6A\xA9\xAA\xAD\xAE\xB9\xBA\xBE",
    "\x55\x59\x65\x69\x6A\xA9\xAA\xAE\xBA\xBE",
    "\x15\x55\x59\x5A\x65\x66\x69\x6A\xA9\xAA\xAE\xBA\xBE",
    "\x55\x56\x59\x5A\x65\x66\x69\x6A\xAA\xAB\xAE\xBA\xBF",
    "\x40\x41\x44\x45\x50\x51\x54\x55\x95\x96\x99\x9A\xA5\xA6\xA9\xAA",
    "\x41\x45\x51\x55\x56\x95\x96\x99\x9A\xA5\xA6\xAA",
    "\x41\x45\x51\x55\x56\x95\x96\x9A\xA6\xAA",
    "\x41\x45\x46\x51\x52\x55\x56\x95\x96\x97\x9A\xA6\xAA\xAB",
    "\x44\x45\x54\x55\x59\x95\x96\x99\x9A\xA5\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\xA6\xAA",
    "\x45\x46\x55\x56\x5A\x95\x96\x9A\x9B\xA6\xAA\xAB",
    "\x44\x45\x54\x55\x59\x95\x99\x9A\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\xAA",
    "\x45\x46\x55\x56\x59\x5A\x95\x96\x99\x9A\x9B\xAA\xAB",
    "\x44\x45\x49\x54\x55\x58\x59\x95\x99\x9A\x9D\xA9\xAA\xAE",
    "\x45\x49\x55\x59\x5A\x95\x99\x9A\x9E\xA9\xAA\xAE",
    "\x45\x49\x55\x56\x59\x5A\x95\x96\x99\x9A\x9E\xAA\xAE",
    "\x45\x55\x56\x59\x5A\x6A\x96\x99\x9A\x9B\x9E\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x65\x95\x96\x99\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x65\x66\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x51\x55\x56\x65\x66\x95\x96\x9A\xA5\xA6\xAA",
    "\x51\x52\x55\x56\x66\x95\x96\x9A\xA6\xA7\xAA\xAB",
    "\x54\x55\x59\x65\x69\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x45\x51\x54\x55\x56\x59\x65\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x45\x51\x55\x56\x59\x5A\x65\x66\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xAB\xEA",
    "\x55\x56\x5A\x66\x6A\x95\x96\x9A\xA6\xAA\xAB",
    "\x54\x55\x59\x65\x69\x95\x99\x9A\xA5\xA9\xAA",
    "\x45\x54\x55\x56\x59\x5A\x65\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xAE\xEA",
    "\x45\x55\x56\x59\x5A\x6A\x95\x96\x99\x9A\xA6\xA9\xAA",
    "\x45\x55\x56\x59\x5A\x66\x6A\x95\x96\x99\x9A\xA6\xAA\xAB",
    "\x54\x55\x58\x59\x69\x95\x99\x9A\xA9\xAA\xAD\xAE",
    "\x55\x59\x5A\x69\x6A\x95\x99\x9A\xA9\xAA\xAE",
    "\x45\x55\x56\x59\x5A\x69\x6A\x95\x96\x99\x9A\xA9\xAA\xAE",
    "\x45\x55\x56\x59\x5A\x6A\x96\x99\x9A\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x65\x95\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x65\x66\x95\x96\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x65\x66\x95\x96\xA5\xA6\xAA",
    "\x51\x52\x55\x56\x65\x66\x95\x96\xA5\xA6\xA7\xAA\xAB",
    "\x54\x55\x59\x65\x69\x95\x99\xA5\xA6\xA9\xAA",
    "\x51\x54\x55\x56\x59\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xBA\xEA",
    "\x51\x55\x56\x65\x66\x6A\x95\x96\x9A\xA5\xA6\xA9\xAA",
    "\x51\x55\x56\x5A\x65\x66\x6A\x95\x96\x9A\xA5\xA6\xAA\xAB",
    "\x54\x55\x59\x65\x69\x95\x99\xA5\xA9\xAA",
    "\x54\x55\x59\x65\x69\x6A\x95\x99\x9A\xA5\xA6\xA9\xAA",
    "\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA",
    "\x55\x56\x59\x5A\x65\x66\x6A\x95\x96\x9A\xA6\xA9\xAA\xAB",
    "\x54\x55\x58\x59\x65\x69\x95\x99\xA5\xA9\xAA\xAD\xAE",
    "\x54\x55\x59\x5A\x65\x69\x6A\x95\x99\x9A\xA5\xA9\xAA\xAE",
    "\x55\x56\x59\x5A\x65\x69\x6A\x95\x99\x9A\xA6\xA9\xAA\xAE",
    "\x55\x56\x59\x5A\x66\x69\x6A\x96\x99\x9A\xA6\xA9\xAA\xAB\xAE\xAF",
    "\x50\x51\x54\x55\x61\x64\x65\x95\xA5\xA6\xA9\xAA\xB5\xBA",
    "\x51\x55\x61\x65\x66\x95\xA5\xA6\xA9\xAA\xB6\xBA",
    "\x51\x55\x56\x61\x65\x66\x95\x96\xA5\xA6\xAA\xB6\xBA",
    "\x51\x55\x56\x65\x66\x6A\x96\xA5\xA6\xA7\xAA\xAB\xB6\xBA\xBB",
    "\x54\x55\x64\x65\x69\x95\xA5\xA6\xA9\xAA\xB9\xBA",
    "\x55\x65\x66\x69\x6A\x95\xA5\xA6\xA9\xAA\xBA",
    "\x51\x55\x56\x65\x66\x69\x6A\x95\x96\xA5\xA6\xA9\xAA\xBA",
    "\x51\x55\x56\x65\x66\x6A\x96\xA5\xA6\xAA\xAB\xBA\xBB",
    "\x54\x55\x59\x64\x65\x69\x95\x99\xA5\xA9\xAA\xB9\xBA",
    "\x54\x55\x59\x65\x66\x69\x6A\x95\x99\xA5\xA6\xA9\xAA\xBA",
    "\x55\x56\x59\x65\x66\x69\x6A\x95\x9A\xA5\xA6\xA9\xAA\xBA",
    "\x55\x56\x5A\x65\x66\x69\x6A\x96\x9A\xA5\xA6\xA9\xAA\xAB\xBA\xBB",
    "\x54\x55\x59\x65\x69\x6A\x99\xA5\xA9\xAA\xAD\xAE\xB9\xBA\xBE",
    "\x54\x55\x59\x65\x69\x6A\x99\xA5\xA9\xAA\xAE\xBA\xBE",
    "\x55\x59\x5A\x65\x66\x69\x6A\x99\x9A\xA5\xA6\xA9\xAA\xAE\xBA\xBE",
    "\x55\x56\x59\x5A\x65\x66\x69\x6A\x9A\xA6\xA9\xAA\xAB\xAE\xBA",
    "\x40\x45\x51\x54\x55\x85\x91\x94\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x41\x45\x51\x55\x56\x85\x91\x95\x96\x99\x9A\xA5\xA6\xAA\xEA",
    "\x41\x45\x51\x55\x56\x85\x91\x95\x96\x9A\xA6\xAA\xD6\xEA",
    "\x41\x45\x51\x55\x56\x86\x92\x95\x96\x97\x9A\xA6\xAA\xAB\xD6\xEA\xEB",
    "\x44\x45\x54\x55\x59\x85\x94\x95\x96\x99\x9A\xA5\xA9\xAA\xEA",
    "\x45\x55\x85\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xDA\xEA",
    "\x45\x55\x56\x85\x95\x96\x99\x9A\xA6\xAA\xDA\xEA",
    "\x45\x55\x56\x86\x95\x96\x9A\x9B\xA6\xAA\xAB\xDA\xEA\xEB",
    "\x44\x45\x54\x55\x59\x85\x94\x95\x99\x9A\xA9\xAA\xD9\xEA",
    "\x45\x55\x59\x85\x95\x96\x99\x9A\xA9\xAA\xDA\xEA",
    "\x45\x55\x56\x59\x5A\x85\x95\x96\x99\x9A\xAA\xDA\xEA",
    "\x45\x55\x56\x5A\x95\x96\x99\x9A\x9B\xA6\xAA\xAB\xDA\xEA\xEB",
    "\x44\x45\x54\x55\x59\x89\x95\x98\x99\x9A\x9D\xA9\xAA\xAE\xD9\xEA\xEE",
    "\x45\x55\x59\x89\x95\x99\x9A\x9E\xA9\xAA\xAE\xDA\xEA\xEE",
    "\x45\x55\x59\x5A\x95\x96\x99\x9A\x9E\xA9\xAA\xAE\xDA\xEA\xEE",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\x9B\x9E\xAA\xAB\xAE\xDA\xEA\xEF",
    "\x50\x51\x54\x55\x65\x91\x94\x95\x96\x99\xA5\xA6\xA9\xAA\xEA",
    "\x51\x55\x91\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xE6\xEA",
    "\x51\x55\x56\x91\x95\x96\x9A\xA5\xA6\xAA\xE6\xEA",
    "\x51\x55\x56\x92\x95\x96\x9A\xA6\xA7\xAA\xAB\xE6\xEA\xEB",
    "\x54\x55\x94\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xE9\xEA",
    "\x55\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x55\x56\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x55\x56\x95\x96\x9A\xA6\xAA\xAB\xEA\xEB",
    "\x54\x55\x59\x94\x95\x99\x9A\xA5\xA9\xAA\xE9\xEA",
    "\x55\x59\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x45\x55\x56\x59\x5A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x45\x55\x56\x5A\x95\x96\x99\x9A\xA6\xAA\xAB\xEA\xEB",
    "\x54\x55\x59\x95\x98\x99\x9A\xA9\xAA\xAD\xAE\xE9\xEA\xEE",
    "\x55\x59\x95\x99\x9A\xA9\xAA\xAE\xEA\xEE",
    "\x45\x55\x59\x5A\x95\x96\x99\x9A\xA9\xAA\xAE\xEA\xEE",
    "\x55\x56\x59\x5A\x95\x96\x99\x9A\xAA\xAB\xAE\xEA\xEF",
    "\x50\x51\x54\x55\x65\x91\x94\x95\xA5\xA6\xA9\xAA\xE5\xEA",
    "\x51\x55\x65\x91\x95\x96\xA5\xA6\xA9\xAA\xE6\xEA",
    "\x51\x55\x56\x65\x66\x91\x95\x96\xA5\xA6\xAA\xE6\xEA",
    "\x51\x55\x56\x66\x95\x96\x9A\xA5\xA6\xA7\xAA\xAB\xE6\xEA\xEB",
    "\x54\x55\x65\x94\x95\x99\xA5\xA6\xA9\xAA\xE9\xEA",
    "\x55\x65\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x51\x55\x56\x65\x66\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x51\x55\x56\x66\x95\x96\x9A\xA5\xA6\xAA\xAB\xEA\xEB",
    "\x54\x55\x59\x65\x69\x94\x95\x99\xA5\xA9\xAA\xE9\xEA",
    "\x54\x55\x59\x65\x69\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x55\x56\x59\x65\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xEA",
    "\x55\x56\x5A\x66\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xAB\xEA\xEB",
    "\x54\x55\x59\x69\x95\x99\x9A\xA5\xA9\xAA\xAD\xAE\xE9\xEA\xEE",
    "\x54\x55\x59\x69\x95\x99\x9A\xA5\xA9\xAA\xAE\xEA\xEE",
    "\x55\x59\x5A\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xAE\xEA\xEE",
    "\x55\x56\x59\x5A\x6A\x95\x96\x99\x9A\xA6\xA9\xAA\xAB\xAE\xEA",
    "\x50\x51\x54\x55\x65\x95\xA1\xA4\xA5\xA6\xA9\xAA\xB5\xBA\xE5\xEA\xFA",
    "\x51\x55\x65\x95\xA1\xA5\xA6\xA9\xAA\xB6\xBA\xE6\xEA\xFA",
    "\x51\x55\x65\x66\x95\x96\xA5\xA6\xA9\xAA\xB6\xBA\xE6\xEA\xFA",
    "\x51\x55\x56\x65\x66\x95\x96\xA5\xA6\xA7\xAA\xAB\xB6\xBA\xE6\xEA\xFB",
    "\x54\x55\x65\x95\xA4\xA5\xA6\xA9\xAA\xB9\xBA\xE9\xEA\xFA",
    "\x55\x65\x95\xA5\xA6\xA9\xAA\xBA\xEA\xFA",
    "\x51\x55\x65\x66\x95\x96\xA5\xA6\xA9\xAA\xBA\xEA\xFA",
    "\x55\x56\x65\x66\x95\x96\xA5\xA6\xAA\xAB\xBA\xEA\xFB",
    "\x54\x55\x65\x69\x95\x99\xA5\xA6\xA9\xAA\xB9\xBA\xE9\xEA\xFA",
    "\x54\x55\x65\x69\x95\x99\xA5\xA6\xA9\xAA\xBA\xEA\xFA",
    "\x55\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xBA\xEA\xFA",
    "\x55\x56\x65\x66\x6A\x95\x96\x9A\xA5\xA6\xA9\xAA\xAB\xBA\xEA",
    "\x54\x55\x59\x65\x69\x95\x99\xA5\xA9\xAA\xAD\xAE\xB9\xBA\xE9\xEA\xFE",
    "\x55\x59\x65\x69\x95\x99\xA5\xA9\xAA\xAE\xBA\xEA\xFE",
    "\x55\x59\x65\x69\x6A\x95\x99\x9A\xA5\xA6\xA9\xAA\xAE\xBA\xEA",
    "\x55\x56\x59\x5A\x65\x66\x69\x6A\x95\x96\x99\x9A\xA5\xA6\xA9\xAA\xAB\xAE\xBA\xEA",
]

# Lookup Tables & Gradients

struct Vertex4D
    dx::Float32
    dy::Float32
    dz::Float32
    dw::Float32
    xsvp::Int64
    ysvp::Int64
    zsvp::Int64
    wsvp::Int64

    function Vertex4D(xsv::Int, ysv::Int, zsv::Int, wsv::Int)
        ssv = Float32((xsv + ysv + zsv + wsv) * UNSKEW_4D)
        new(-xsv - ssv, -ysv - ssv, -zsv - ssv, -wsv - ssv,
            (xsv * PRIME_X)%Int64,
            (ysv * PRIME_Y)%Int64,
            (zsv * PRIME_Z)%Int64,
            (wsv * PRIME_W)%Int64)
    end
end
Vertex4D(v) = Vertex4D(((v >> 0) & 3) - 1, ((v >> 2) & 3) - 1,
                       ((v >> 4) & 3) - 1, ((v >> 6) & 3) - 1)

const offsets_4D = Vector{UInt16}(undef, 256)

function create_lattice_map()
    num_vertices = 0
    tab = falses(256)
    for i = 1:256
        codes = vertex_codes[i]
        num_vertices += sizeof(codes)
        offsets_4D[i] = num_vertices
        for v in codeunits(codes)
            tab[v+1] = true
        end
    end
    siz = count(tab)
    inv_map = zeros(UInt8, 256)
    vertices = Vector{UInt8}(undef, siz)
    vertex_map = Vector{Vertex4D}(undef, siz)
    r = findfirst(tab)
    i = 1
    while r !== nothing
        inv_map[r] = i
        vertices[i] = r-1
        vertex_map[i] = Vertex4D(r-1)
        i += 1
        r = findnext(tab, r+1)
    end
    io = IOBuffer()
    println(io, "const vertex_ind = [")
    for codes in vertex_codes
        print(io, "    ")
        for v in codeunits(codes)
            i = inv_map[v+1]
            print(io, i<10 ? " " : "", i, ',')
        end
        println(io)
    end
    println(io, "]")
    String(take!(io)), vertex_map, vertices, inv_map
end


text:
module SeisPDF

# Read input trace

# Read response file
include("read_response.jl")

# Split trace to 1-hour then 15-minute
include("slicing.jl")

# For each 15-minute segment
# Detrend
include("detrend.jl")

# Demean
include("demean.jl")

# FFT
include("compute_fft.jl")

# Remove response
include("remove_response.jl")

# Taper
include("taper.jl")

# Calculate PSD
include("calculate_psd.jl")
# end segment

# Calculate PSD and PDF summary
include("summarize_psd.jl")
include("summarize_pdf.jl")

end


text:
# This script is intended to be a dependency for the hydropower example.
using PowerSystems
using Statistics
using Dates
using TimeSeries

# We can use some of the simple data that's been assembled for testing PowerSimulations.
PowerSystems.download(PowerSystems.TestData; branch = "master") # *note* add `force=true` to get a fresh copy
data_dir = joinpath(dirname(dirname(pathof(PowerSystems))), "data");

rawsys = PowerSystems.PowerSystemTableData(
    joinpath(data_dir, "5-bus-hydro"),
    100.0,
    joinpath(data_dir, "5-bus-hydro", "user_descriptors.yaml");
    generator_mapping_file = joinpath(data_dir, "5-bus-hydro", "generator_mapping.yaml"),
)
c_sys5_hy_uc = System(
    rawsys,
    timeseries_metadata_file = joinpath(
        data_dir,
        "forecasts",
        "5bus_ts",
        "7day",
        "timeseries_pointers_da_7day.json",
    ),
    time_series_in_memory = true,
)
transform_single_time_series!(c_sys5_hy_uc, 24, Hour(24))

c_sys5_hy_ed = System(
    rawsys,
    timeseries_metadata_file = joinpath(
        data_dir,
        "forecasts",
        "5bus_ts",
        "7day",
        "timeseries_pointers_rt_7day.json",
    ),
    time_series_in_memory = true,
)
transform_single_time_series!(c_sys5_hy_ed, 12, Hour(1))

c_sys5_hy_wk = System(rawsys, time_series_in_memory = true)

n = 3# number of days

# And an hourly system with longer time scales
MultiDay = collect(
    get_forecast_initial_times(c_sys5_hy_uc)[1]:Hour(24):(get_forecast_initial_times(
        c_sys5_hy_uc,
    )[1] + Day(n)),
);

for load in get_components(PowerLoad, c_sys5_hy_uc)
    fc_values = Vector{Float64}()
    for t in MultiDay
        fc = get_time_series_array(Deterministic, load, "max_active_power"; start_time = t)
        push!(fc_values, mean(values(fc)))
    end
    wk_fc = TimeArray(MultiDay, fc_values)
    add_time_series!(
        c_sys5_hy_wk,
        get_component(PowerLoad, c_sys5_hy_wk, get_name(load)),
        SingleTimeSeries("max_active_power", wk_fc),
    )
end

for gen in get_components(HydroGen, c_sys5_hy_uc)
    fc_values = Vector{Float64}()
    for t in MultiDay
        fc = get_time_series_array(Deterministic, gen, "max_active_power"; start_time = t)
        push!(fc_values, mean(values(fc)))
    end

    wk_fc = TimeArray(MultiDay, fc_values)

    add_time_series!(
        c_sys5_hy_wk,
        get_component(HydroGen, c_sys5_hy_wk, get_name(gen)),
        SingleTimeSeries("max_active_power", wk_fc),
    )
end

for gen in get_components(HydroEnergyReservoir, c_sys5_hy_uc)
    for label in ["storage_capacity", "inflow", "storage_target"]
        fc_values = Vector{Float64}()
        for t in MultiDay
            fc = get_time_series_array(
                Deterministic,
                gen,
                label;
                start_time = get_forecast_initial_times(c_sys5_hy_uc)[1],
            )
            push!(fc_values, mean(values(fc)))
        end

        wk_fc = TimeArray(MultiDay, fc_values)

        add_time_series!(
            c_sys5_hy_wk,
            get_component(HydroGen, c_sys5_hy_wk, get_name(gen)),
            SingleTimeSeries(label, wk_fc),
        )
    end
end
transform_single_time_series!(c_sys5_hy_wk, 2, Hour(48))


text:
#=

  Anagrams/Deranged anagrams in Julia.

  http://rosettacode.org/wiki/Anagrams/Deranged_anagrams
  """
  Two or more words are said to be anagrams if they have the same characters, but in
  a different order.

  By analogy with derangements we define a deranged anagram as two words with the same
  characters, but in which the same character does not appear in the same position
  in both words.

  Task
  Use the word list at unixdict [http://wiki.puzzlers.org/pub/wordlists/unixdict.txt]
  to find and display the longest deranged anagram.
  """

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("jl_utils.jl")
using Underscores

# are w1 and w2 deranged?
function is_deranged(w1,w2)
    for i in 1:length(w1)
        if w1[i] === w2[i]
            return false
        end
    end
    return true
end

function sorted_hash(words)
    d = Dict()
    for word in words
        s = sort(split(word,""))
        get!(d,s,[])
        push!(d[s],word)
    end
    return d
end

function deranged_angrams(words)
    d = sorted_hash(words)

    max_len = 0
    max_words = []
    for (sorted,ws) in collect(d)
        if length(ws) == 1
            continue
        end
        len = length(ws)
        wlen  = length(join(sorted,""))
        if len > 1 && wlen > max_len
            found = filter(w->is_deranged(w[1],w[2]),combinations(2,ws))
            if length(found) > 0
                if max_len == wlen
                    push!(max_words,found[1])
                else
                    max_words = found
                    max_len = wlen
                end
            end
        end
    end

    println("max_len:$max_len")
    println("max_words:$max_words")
end

# wordlist = "eng_dict.txt"
wordlist = "unixdict.txt"
# wordlist = "words_lower.txt"
words = readlines(wordlist)
deranged_angrams(words)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
import DataStructures
import StaticArrays
using BenchmarkTools

include("../src/soalight.jl")
include("../src/allocator.jl")
include("../src/list.jl")

function testbase1(n, p)
    l = DataStructures.nil()
    for j in 1:n[]
        l = DataStructures.cons((j, p) , l)
    end
    s = 0
    while !isa(l, DataStructures.Nil)
        s += l.head[1]
        l = l.tail
    end
    @assert s == n[] * (n[] + 1) / 2
end

function testbase2(n, p)
    l = nil()
    for j in 1:n[]
        l = cons((j, p) , l)
    end
    s = 0
    while l != nothing
        s += l.car[1]
        l = l.cdr
    end
    @assert s == n[] * (n[] + 1) / 2
end

function testalloc(n, p, alloc)
    l = nil(alloc)
    for j in 1:n[]
        l = cons((j, p), l, alloc)
    end
    s = 0
    while l != 0
        car, cdr = alloc[l]
        s += car[1]
        l = cdr
    end
    @assert s == n[] * (n[] + 1) / 2
    emptyend!(alloc)
end

function testfree(n, p, alloc)
    l = nil(alloc)
    for j in 1:n[]
        l = cons((j, p), l, alloc)
    end
    s = 0
    while l != 0
        car, cdr = alloc[l]
        s += car[1]
        deallocate(l, alloc)
        l = cdr
    end
    @assert s == n[] * (n[] + 1) / 2
    @assert isempty(alloc)
end

ans = nothing
for n in [100000, 10]
    println(n, " elements")
    for p in [1, 10, 100]
        println(" payload of ", p, " float(s)")

        Payload = StaticArrays.SVector{p, Float64}

        GC.gc()

        print("  DataStructures.List                    ")
        @btime ($ans = testbase1(Ref($n), zeros($Payload)))

        GC.gc()

        print("  List without allocator                 ")
        @btime ($ans = testbase2(Ref($n), zeros($Payload)))

        GC.gc()

        print("  with fixed allocator                   ")
        alloc = Allocator{ListNode{Tuple{Int, Payload}, Int}, Int}(n)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with resizable allocator               ")
        alloc = Allocator{ListNode{Tuple{Int, Payload}, Int}, Int}(nothing)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with fixed free list allocator         ")
        alloc = FreeListAllocator{ListNode{Tuple{Int, Payload}, Int}, Int}(n)
        @btime ($ans = testfree(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with resizable free list allocator     ")
        alloc = FreeListAllocator{ListNode{Tuple{Int, Payload}, Int}, Int}(nothing)
        @btime ($ans = testfree(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with fixed SOA allocator               ")
        store = TupleVector{ListNode{Tuple{Int, Payload}, Int}}(n)
        alloc = SOAllocator{ListNode{Tuple{Int, Payload}, Int}, Int, typeof(store)}(store, n)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $alloc))
        GC.gc()

        print("  with resizable SOA allocator           ")
        store = TupleVector{ListNode{Tuple{Int, Payload}, Int}}(N)
        alloc = SOAllocator{ListNode{Tuple{Int, Payload}, Int}, Int, typeof(store)}(store, nothing)
        @btime ($ans = testalloc(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with fixed free list SOA allocator     ")
        store = TupleVector{ListNode{Tuple{Int, Payload}, Int}}(n)
        alloc = FreeListSOAllocator{ListNode{Tuple{Int, Payload}, Int}, Int, typeof(store)}(store, n)
        @btime ($ans = testfree(Ref($n), zeros($Payload), $alloc))

        GC.gc()

        print("  with resizable free list SOA allocator ")
        store = TupleVector{ListNode{Tuple{Int, Payload}, Int}}(N)
        alloc = FreeListSOAllocator{ListNode{Tuple{Int, Payload}, Int}, Int, typeof(store)}(store, nothing)
        @btime ($ans = testfree(Ref($n), zeros($Payload), $alloc))
    end
end


text:
# # Decaying barotropic QG beta-plane turbulence
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_betadecay.ipynb).
#
# An example of decaying barotropic quasi-geostrophic turbulence on a beta plane.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Plots, Printf, Statistics, Random"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Plots, Printf, Random

using Statistics: mean

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution: n² grid points
stepper = "FilteredRK4"  # timestepper
     dt = 0.04           # timestep
 nsteps = 2000           # total number of time-steps
 nsubs  = 20             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
β = 10.0      # planetary PV gradient
μ = 0.0       # bottom drag
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing a viscosity
# coefficient `ν` leads to the module's default value: `ν=0`. In this example numerical instability
# due to accumulation of enstrophy at high wavenumbers is taken care with the `FilteredTimestepper`
# we picked. Thus, we choose not to do any dealiasing by providing `aliased_fraction=0`.

prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ,
                                  dt=dt, stepper=stepper, aliased_fraction=0)
nothing # hide

# and define some shortcuts
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# ## Setting initial conditions

# Our initial condition consist of a flow that has power only at wavenumbers with
# ``6 < \frac{L}{2\pi} \sqrt{k_x^2 + k_y^2} < 10`` and initial energy ``E_0``.
# `ArrayType()` function returns the array type appropriate for the device, i.e., `Array` for
# `dev = CPU()` and `CuArray` for `dev = GPU()`.

E₀ = 0.08 # energy of initial condition

K = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber

Random.seed!(1234)
q₀h = ArrayType(dev)(randn(Complex{eltype(grid)}, size(sol)))
@. q₀h = ifelse(K < 6  * 2π/L, 0, q₀h)
@. q₀h = ifelse(K > 10 * 2π/L, 0, q₀h)
@. q₀h[1, :] = 0    # remove any power from zonal wavenumber k=0
q₀h *= sqrt(E₀ / SingleLayerQG.energy(q₀h, vars, params, grid)) # normalize q₀ to have energy E₀
q₀ = irfft(q₀h, grid.nx)

SingleLayerQG.set_q!(prob, q₀)
nothing # hide

# Let's plot the initial vorticity and streamfunction. Note that when plotting, we decorate
# the variable to be plotted with `Array()` to make sure it is brought back on the CPU when
# `vars` live on the GPU.

p1 = heatmap(x, y, Array(vars.q'),
         aspectratio = 1,
              c = :balance,
           clim = (-12, 12),
          xlims = (-grid.Lx/2, grid.Lx/2),
          ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity ∂v/∂x-∂u/∂y",
     framestyle = :box)

p2 = contourf(x, y, Array(vars.ψ'),
        aspectratio = 1,
             c = :viridis,
        levels = range(-0.7, stop=0.7, length=20),
          clim = (-0.35, 0.35),
         xlims = (-grid.Lx/2, grid.Lx/2),
         ylims = (-grid.Ly/2, grid.Ly/2),
        xticks = -3:3,
        yticks = -3:3,
        xlabel = "x",
        ylabel = "y",
         title = "initial streamfunction ψ",
    framestyle = :box)

layout = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout = layout, size = (800, 360))

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps)
Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_decayingbetaturb"
plotname = "snapshots"
filename = joinpath(filepath, "decayingbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
out = Output(prob, filename, (:sol, get_sol))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction and
# their corresponding zonal mean structure.

function plot_output(prob)
  q = Array(prob.vars.q)
  ψ = Array(prob.vars.ψ)
  q̄ = Array(mean(q, dims=1)')
  ū = Array(mean(prob.vars.u, dims=1)')

  pq = heatmap(x, y, q',
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-12, 12),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ∂v/∂x-∂u/∂y",
        framestyle = :box)

  pψ = contourf(x, y, ψ',
       aspectratio = 1,
            legend = false,
                 c = :viridis,
            levels = range(-0.7, stop=0.7, length=20),
              clim = (-0.35, 0.35),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  pqm = plot(q̄, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-2.2, 2.2),
            xlabel = "zonal mean q",
            ylabel = "y")
  plot!(pqm, 0*y, y, linestyle=:dash, linecolor=:black)

  pum = plot(ū, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-0.55, 0.55),
            xlabel = "zonal mean u",
            ylabel = "y")
  plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)

  layout = @layout Plots.grid(2, 2)
  p = plot(pq, pqm, pψ, pum, layout = layout, size = (800, 720))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)

  if j % round(Int, nsteps/nsubs / 4) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = Array(vars.q)
  p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t)
  p[3][1][:z] = Array(vars.ψ)
  p[2][1][:x] = Array(mean(vars.q, dims=1)')
  p[4][1][:x] = Array(mean(vars.u, dims=1)')

  stepforward!(prob, diags, nsubs)
  SingleLayerQG.updatevars!(prob)

end

mp4(anim, "singlelayerqg_betadecay.mp4", fps=8)

# ## Save

# Finally, we can save, e.g., the last snapshot via
# ```julia
# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
# savefig(savename)
# ```


text:
# Crying baby problem described in DMU book
# State: hungry = true; not hungry = false
# Action: feed = true; do nothing = false
# Observation: crying = true; not crying = false

mutable struct BabyPOMDP <: POMDP{Bool, Bool, Bool}
    r_feed::Float64
    r_hungry::Float64
    p_become_hungry::Float64
    p_cry_when_hungry::Float64
    p_cry_when_not_hungry::Float64
    discount::Float64
end
BabyPOMDP(r_feed, r_hungry) = BabyPOMDP(r_feed, r_hungry, 0.1, 0.8, 0.1, 0.9)
BabyPOMDP() = BabyPOMDP(-5., -10.)

updater(problem::BabyPOMDP) = DiscreteUpdater(problem)

actions(::BabyPOMDP) = (true, false)
actionindex(::BabyPOMDP, a::Bool) = a + 1
states(::BabyPOMDP) = (true, false)
stateindex(::BabyPOMDP, s::Bool) = s + 1
observations(::BabyPOMDP) = (true, false)
obsindex(::BabyPOMDP, o::Bool) = o + 1

# start knowing baby is not not hungry
initialstate(::BabyPOMDP) = BoolDistribution(0.0)
initialobs(m::BabyPOMDP, s) = observation(m, s)

function transition(pomdp::BabyPOMDP, s::Bool, a::Bool)
    if a # fed
        return BoolDistribution(0.0)
    elseif s # did not feed when hungry
        return BoolDistribution(1.0)
    else # did not feed when not hungry
        return BoolDistribution(pomdp.p_become_hungry)
    end
end

function observation(pomdp::BabyPOMDP, sp::Bool)
    if sp # hungry
        return BoolDistribution(pomdp.p_cry_when_hungry)
    else
        return BoolDistribution(pomdp.p_cry_when_not_hungry)
    end
end

function reward(pomdp::BabyPOMDP, s::Bool, a::Bool)
    r = 0.0
    if s # hungry
        r += pomdp.r_hungry
    end
    if a # feed
        r += pomdp.r_feed
    end
    return r
end

discount(p::BabyPOMDP) = p.discount

# some example policies
struct Starve <: Policy end
action(::Starve, ::B) where {B} = false
updater(::Starve) = NothingUpdater()

struct AlwaysFeed <: Policy end
action(::AlwaysFeed, ::B) where {B} = true
updater(::AlwaysFeed) = NothingUpdater()

# feed when the previous observation was crying - this is nearly optimal
struct FeedWhenCrying <: Policy end
updater(::FeedWhenCrying) = PreviousObservationUpdater()
function action(::FeedWhenCrying, b::Union{Nothing, Bool})
    if b == nothing || b == false # not crying (or null)
        return false
    else # is crying
        return true
    end
end
action(::FeedWhenCrying, b::Bool) = b
action(p::FeedWhenCrying, b::Missing) = false
# assume the second argument is a distribution
action(::FeedWhenCrying, d::Any) = pdf(d, true) > 0.5

# deprecated in POMDPs 0.9
POMDPs.initialstate_distribution(::BabyPOMDP) = BoolDistribution(0.0)


text:
""" Functions for handling the modification of the assignements in the branch and bound """

# All functions are public

"""
Assign the next variable to one.
"""
function ASSIGNMENT_addVar!(assignment::Assignment, prob::BiOKP)

	assignment.lastAssigned += 1
	assignment.assign[assignment.lastAssigned] = 1

	assignment.weight += prob.weights[assignment.lastAssigned]

	assignment.profit += prob.profits[1:end, assignment.lastAssigned]
end

"""
Assign the current variable to zero.

The parameter `canAddVar` indicates if the variable has been assigned to one in the previous iteration. If `canAddVar` is true, we have to remove the weight and profit of the variable from the assignment.
"""
function ASSIGNMENT_removeVar!(assignment::Assignment, prob::BiOKP, canAddVar::Bool)

    if !canAddVar
		assignment.lastAssigned += 1
	else
		assignment.weight -= prob.weights[assignment.lastAssigned]
		assignment.profit -= prob.profits[1:end, assignment.lastAssigned]
	end

	assignment.assign[assignment.lastAssigned] = 0
end

"""
Update the `assignment` to match the version of the parent (unassigned the current variable).

NOTE : no need to take of the weight of object ... because the backtrack always happens after ASSIGNMENT_removeVar!
"""
function ASSIGNMENT_backtrack!(assignment::Assignment)

	assignment.assign[assignment.lastAssigned] = -1
	assignment.lastAssigned -= 1
end


text:
module Archimedes

using LinearAlgebra
using Statistics

using Unitful
using Luxor

const _current_luxor_drawing = Ref{Luxor.Drawing}()

export Point2D,
       PointMass,
       mass,
       getx,
       gety,
       radius,
       bbox,
       centroid,
       gravitycenter,
       draw,
       drawing,
       current_drawing,
       BoxShip,
       corners,
       labeled_point,
       buoyancycenter,
       toship,
       isocarene,
       metacenter,
       setopacity,
       puttext

struct Point2D
    x::typeof(1.0u"m")
    y::typeof(1.0u"m")
end
getx(p::Point2D) = p.x
gety(p::Point2D) = p.y
Base.:-(a::Point2D, b::Point2D) = Point2D(a.x-b.x, a.y-b.y)
Base.:+(a::Point2D, b::Point2D) = Point2D(a.x+b.x, a.y+b.y)
Base.:*(a::Point2D, b::Number) = Point2D(a.x*b, a.y*b)
Base.:*(b::Number, a::Point2D) = Point2D(a.x*b, a.y*b)
Base.:/(a::Point2D, b::Number) = Point2D(a.x/b, a.y/b)
LinearAlgebra.norm(a::Point2D) = sqrt(a.x^2 + a.y^2)

struct PointMass
    coordinates::Point2D
    mass::typeof(1.0u"kg")
    radius::typeof(1.0u"m")
end
PointMass(x,y,m,r) = PointMass(Point2D(x,y),m,r)
coords(p::PointMass) = p.coordinates
mass(p::PointMass) = p.mass
getx(p::PointMass) = getx(coords(p))
gety(p::PointMass) = gety(coords(p))
radius(p::PointMass) = p.radius

"""
Get the bounding box (xmin, xmax, ymin, ymax) of an array of pointmasses
"""
bbox(pointmasses) = (extrema(getx.(pointmasses))..., extrema(gety.(pointmasses))...)

"""
Get the centroid of an array of pointmasses
"""
centroid(pointmasses) = Point2D(mean(getx.(pointmasses)), mean(gety.(pointmasses)))

"""
Get the center of gravity
"""
function gravitycenter(pointmasses)
  m = sum(mass.(pointmasses))
  x = sum(getx.(pointmasses) .* mass.(pointmasses)) / m
  y = sum(gety.(pointmasses) .* mass.(pointmasses)) / m
  r = sqrt(sum(radius.(pointmasses) .^ 2))
  return PointMass(x, y, m, r)
end

"""
Helper to map the points with units to pixel-coordinates
"""
struct CoordMapping
  origin::Point2D
  scaling::typeof(1.0/(1.0u"m"))
end

function CoordMapping(width::Number, bbox::Tuple)
  (xmin,xmax,ymin,ymax) = bbox
  center = Point2D(mean((xmin,xmax)), mean((ymin,ymax)))
  scaling = width/(xmax-xmin)
  return CoordMapping(center, scaling)
end

scale(x, mapping::CoordMapping) = x*mapping.scaling

function remap(p, mapping)
  (xoff, yoff) = (getx(mapping.origin), gety(mapping.origin))
  return Luxor.Point(((getx(p)-xoff, -(gety(p)-yoff)).*mapping.scaling)...)
end

"""
Initialize a drawing, returing the point mapping
"""
function drawing(bbox, figwidth=300)
  (xmin,xmax,ymin,ymax) = bbox
  ar = (xmax-xmin)/(ymax-ymin)
  figheight = Int(round(figwidth/ar))
  margin = 30
  bottom_padding = 40
  _current_luxor_drawing[] = Drawing(figwidth+2*margin, figheight+2*margin+bottom_padding, :svg)
  origin()
  background("white")
  return CoordMapping(figwidth, bbox)
end

function puttext(txt, p::Point2D, mapping::CoordMapping)
  sethue("black")
  c = remap(p,mapping)
  text(txt, c, halign=:center, valign=:top)
end

"""
Draws a pointmass into the current context
"""
function draw(p::PointMass, mapping::CoordMapping; hue="black")
  sethue(hue)
  c = remap(p,mapping)
  r = scale(radius(p), mapping)
  circle(c, r, :fill)
  spacing = 5.0
  step = 10.0
  text("x = $(getx(p))", Point(c.x,c.y+r+spacing), halign=:center, valign=:top)
  text("y = $(gety(p))", Point(c.x,c.y+r+spacing+step), halign=:center, valign=:top)
  text("m = $(mass(p))", Point(c.x,c.y+r+spacing+2*step), halign=:center, valign=:top)
end

function draw(p::Point2D, mapping::CoordMapping, r=5.0; hue="black", label="")
  sethue(hue)
  c = remap(p,mapping)
  circle(c, r, :fill)
  if label != ""
    text(label, Point(c.x+1.2*r,c.y+1.2*r), halign=:center, valign=:top)
  end
end

"""
Shows the drawing and removes the temp file
"""
function current_drawing()
  finish()
  return _current_luxor_drawing[]
end

struct BoxShip
  width::typeof(1.0u"m")
  height::typeof(1.0u"m")
  draft::typeof(1.0u"m")
  KG::typeof(1.0u"m")
  heel::Float64
  vshift::typeof(1.0u"m")

  BoxShip(width, height, draft, KG, heel=0.0, vshift=0.0u"m") = new(width,height,draft,KG,heel,vshift)
end

forward_trans(p, s::BoxShip) = Point2D(cos(s.heel)*getx(p) - sin(s.heel)*gety(p),  sin(s.heel)*getx(p) + cos(s.heel)*gety(p) + s.vshift)
inverse_trans(p, s::BoxShip) = Point2D(cos(s.heel)*getx(p) + sin(s.heel)*(gety(p) - s.vshift), -sin(s.heel)*getx(p) + cos(s.heel)*(gety(p) - s.vshift))

"""
Convert position p to ship coordinates
"""
toship(p, s) = inverse_trans(p, s) - Point2D(0.0u"m", s.height/2 - s.draft)

function corners(s::BoxShip)
  x = s.width/2
  ymax = s.height-s.draft
  ymin = -s.draft
  return forward_trans.(Point2D.([-x,x,x,-x], [ymin,ymin,ymax,ymax]), Ref(s))
end

function waterline(s::BoxShip)
  x = 1.4*s.width/2
  y = 0.0u"m"
  return (Point2D(-x,y), Point2D(x,y))
end

function gravitycenter(s::BoxShip)
  return forward_trans(Point2D(0.0u"m", s.KG-s.draft), s)
end

bbox(s::BoxShip) = 1.4 .* (-s.width/2, s.width/2, -s.draft, s.height-s.draft)

function wl_intersect(p1, p2)
  x1 = getx(p1)
  y1 = gety(p1)
  x2 = getx(p2)
  y2 = gety(p2)
  r = (x2-x1)/(y2-y1)
  return Point2D(-y1*r+x1, 0.0u"m")
end

"""
Corners of the underwater part
"""
function carene(s::BoxShip)
  ship_pts = corners(s)
  N = length(ship_pts)
  pts_above = collect(Iterators.filter((p) -> gety(p[2]) > 0u"m", enumerate(ship_pts)))
  sort!(pts_above, lt = (a,b) -> getx(a[2]) < getx(b[2]))
  l1 = pts_above[1][1]
  l2 = l1 % N + 1
  wl_left = wl_intersect(ship_pts[l1], ship_pts[l2])
  r1 = pts_above[end][1]
  r2 = (r1 - 2 + N) % N + 1
  wl_right = wl_intersect(ship_pts[r1], ship_pts[r2])
  pts_below = collect(Iterators.filter((p) -> gety(p[2]) < 0u"m", enumerate(ship_pts)))
  sort!(pts_below, lt = (a,b) -> getx(a[2]) < getx(b[2]))
  return (getindex.(pts_below, 2)..., wl_right, wl_left)
end

function metacenter(s::BoxShip)
  trap = carene(s)
  w_wl = getx(trap[end-1]) - getx(trap[end])
  F = inverse_trans(Point2D((getx(trap[end-1]) + getx(trap[end]))/2, 0.0u"m"),s)
  BM = (w_wl^3 / 12)/carene_area(s)
  B = inverse_trans(buoyancycenter(s),s)
  nvec = F - B
  Mship = B + Point2D(BM*sin(s.heel), BM*cos(s.heel))
  return forward_trans(Mship, s)
end

"""
Make a 3-vector out of a Point2D
"""
vec3(p::Point2D) = [getx(p), gety(p), 0.0u"m"]

"""
Area of a triangle
"""
function area(pts)
  v = vec3.(pts)
  return abs(cross(v[3]-v[1], v[2]-v[1])[3])/2
end

function triangulate(pts::Union{AbstractArray{ET},NTuple{N,ET} where N}) where ET
  c = centroid(pts)
  triags = Array{NTuple{3,ET}}(undef,length(pts))
  N = length(pts)
  for i in 1:length(pts)
    triags[i] = (c, pts[i], pts[i%N+1])
  end
  return triags
end

function buoyancycenter(s::BoxShip)
  trap = carene(s)
  triags = triangulate(trap)
  areas = area.(triags)
  centroids = centroid.(triags)
  A = sum(areas)
  x = sum(getx.(centroids) .* areas) / A
  y = sum(gety.(centroids) .* areas) / A
  return Point2D(x,y)
end

function carene_area(s::BoxShip)
  trap = carene(s)
  triags = triangulate(trap)
  return sum(area.(triags))
end

function isocarene(s0, θ)
  θ0 = s0.heel
  T = s0.draft
  DT_l = -T
  DT_u = T
  A0 = carene_area(s0)
  for i in 1:50
    DT = (DT_u+DT_l)/2
    s1 = BoxShip(s0.width, s0.height, T, s0.KG, θ+θ0, DT)
    A1 = carene_area(s1)
    if abs(A1-A0)/A0 < 1e-8
      return s1
    end

    if A1 < A0
      DT_u = DT
    else
      DT_l = DT
    end
  end
  error("isocarene did not converge")
end

function labeled_point(p, label, hue="black", radius=5.0)
  sethue(hue)
  circle(p, radius, :fill)
  text(label, Point(p.x+1.2*radius,p.y+1.2*radius), halign=:center, valign=:top)
end

function draw(m::CoordMapping, s::BoxShip, transformation=((p,::BoxShip) -> p); showM=true, showB=true)
  ship_pts = transformation.(corners(s),Ref(s))
  p = remap.(ship_pts, Ref(m))
  sethue("black")
  line(p[1], p[2], :stroke)
  line(p[2], p[3], :stroke)
  line(p[4], p[1], :stroke)
  sethue("darkgrey")
  line(p[3], p[4], :stroke)
  labeled_point(remap(transformation(gravitycenter(s),s),m), "G")
  if showM || showB
    M = remap(transformation(metacenter(s),s),m)
    B = remap(transformation(buoyancycenter(s),s),m)
    if showM
      labeled_point(M, "M")
    end
    if showB
      labeled_point(B, "B", "red")
    end
    if showM && showB
      setdash("dot")
      line(M, B, :stroke)
      setdash("solid")
    end
  end
  wl = remap.(transformation.(waterline(s),Ref(s)), Ref(m))
  sethue("blue")
  line(wl[1], wl[2], :stroke)
  #trap = carene(s)
  #F = remap(transformation(Point2D((getx(trap[end-1]) + getx(trap[end]))/2, 0.0u"m"),s),m)
  #labeled_point(F, "F", "blue")
  return m
end

function draw(s::BoxShip, figwidth=300, transformation=((p,::BoxShip) -> p); showM=true, showB=true)
  draw(drawing(bbox(s), figwidth), s, transformation, showM=showM, showB=showB)
end

setopacity = Luxor.setopacity

end # module