text:
"""
    HelperFuncs

Generic utils.
"""
module HelperFuncs

# We can work on qualifying these later
export vertical_interpolation,
    nc_fetch_interpolate,
    fetch_interpolate_transform,
    get_height,
    normalize_profile,
    nc_fetch,
    is_face_variable,
    get_stats_path,
    compute_mse,
    penalize_nan,
    serialize_struct,
    deserialize_struct,
    jld2_path,
    scm_init_path,
    scm_output_path,
    scm_val_init_path,
    scm_val_output_path,
    ekobj_path,
    write_versions,
    expand_dict_entry,
    get_entry,
    change_entry!

using NCDatasets
using Statistics
using Interpolations
using LinearAlgebra
using Glob
using JSON
using Random

"""
    vertical_interpolation(
        var_name::String,
        filename::String,
        z_scm::Vector{FT};
    ) where {FT <: AbstractFloat}

Returns the netcdf variable var_name interpolated to heights z_scm.

Inputs:
 - `var_name` :: Name of variable in the netcdf dataset.
 - `filename` :: nc filename
 - `z_scm` :: Vertical coordinate vector onto which var_name is interpolated.
Output:
 - The interpolated vector.
"""
function vertical_interpolation(var_name::String, filename::String, z_scm::Vector{FT};) where {FT <: AbstractFloat}
    z_ref = get_height(filename, get_faces = is_face_variable(filename, var_name))
    var_ = nc_fetch(filename, var_name)
    if length(size(var_)) == 2
        # Create interpolant
        nodes = (z_ref, 1:size(var_, 2))
        var_itp = extrapolate(interpolate(nodes, var_, (Gridded(Linear()), NoInterp())), Line())
        # Return interpolated vector
        return var_itp(z_scm, 1:size(var_, 2))
    elseif length(size(var_)) == 1
        # Create interpolant
        nodes = (z_ref,)
        var_itp = LinearInterpolation(nodes, var_; extrapolation_bc = Line())
        # Return interpolated vector
        return var_itp(z_scm)
    end
end

"""
    nc_fetch_interpolate(
        var_name::String,
        filename::String,
        z_scm::Union{Vector{<:Real}, Nothing};
    )

Returns the netcdf variable var_name, possibly interpolated to heights z_scm.

Inputs:
 - `var_name` :: Name of variable in the netcdf dataset.
 - `filename` :: nc filename
 - `z_scm` :: Vertical coordinate vector onto which var_name is interpolated.
Output:
 - The interpolated vector.
"""
function nc_fetch_interpolate(var_name::String, filename::String, z_scm::Union{Vector{<:Real}, Nothing};)
    if !isnothing(z_scm)
        return vertical_interpolation(var_name, filename, z_scm)
    else
        return nc_fetch(filename, var_name)
    end
end

"""
    fetch_interpolate_transform(
        var_name::String,
        filename::String,
        z_scm::Union{Vector{<:Real}, Nothing};
    )

Returns the netcdf variable var_name, possibly interpolated to heights z_scm. If the
variable needs to be transformed to be equivalent to an SCM variable, applies the
transformation as well.

Inputs:
 - `var_name` :: Name of variable in the netcdf dataset.
 - `filename` :: nc filename
 - `z_scm` :: Vertical coordinate vector onto which var_name is interpolated.
Output:
 - The interpolated and transformed vector.
"""
function fetch_interpolate_transform(var_name::String, filename::String, z_scm::Union{Vector{<:Real}, Nothing};)
    # PyCLES vertical fluxes are per volume, not mass
    if occursin("resolved_z_flux", var_name)
        var_ = nc_fetch_interpolate(var_name, filename, z_scm)
        rho_half = nc_fetch_interpolate("rho0_half", filename, z_scm)
        var_ = var_ .* rho_half
    elseif var_name == "horizontal_vel"
        u_ = nc_fetch_interpolate("u_mean", filename, z_scm)
        v_ = nc_fetch_interpolate("v_mean", filename, z_scm)
        var_ = sqrt.(u_ .^ 2 + v_ .^ 2)
    else
        var_ = nc_fetch_interpolate(var_name, filename, z_scm)
    end
    return var_
end

"""
    get_height(filename::String; get_faces::Bool = false)

Returns the vertical cell centers or faces of the given configuration.

Inputs:
 - filename :: nc filename.
 - get_faces :: If true, returns the coordinates of cell faces. Otherwise,
    returns the coordinates of cell centers.
Output:
 - z: Vertical level coordinates.
"""
function get_height(filename::String; get_faces::Bool = false)
    if get_faces
        return nc_fetch(filename, ("zf", "z"))
    else
        return nc_fetch(filename, ("zc", "z_half"))
    end
end

"""
    normalize_profile(profile_vec, n_vars, var_vec)

Perform normalization of n_vars profiles contained in profile_vec
using the variance associated with each variable, contained
in var_vec.
"""
function normalize_profile(profile_vec, n_vars, var_vec)
    y = deepcopy(profile_vec)
    var_dof = Integer(length(profile_vec) / n_vars)
    for i in 1:n_vars
        y[(var_dof * (i - 1) + 1):(var_dof * i)] = y[(var_dof * (i - 1) + 1):(var_dof * i)] ./ sqrt(var_vec[i])
    end
    return y
end

"""
    nc_fetch(filename::String, var_names::NTuple{N, Tuple}) where {N}
    nc_fetch(filename::String, var_name::String)

Returns the data for a variable `var_name` (or
tuple of strings, `varnames`), looping through
all dataset groups.
"""
function nc_fetch(filename::String, var_names::Tuple)
    NCDataset(filename) do ds
        for var_name in var_names
            if haskey(ds, var_name)
                return Array(ds[var_name])
            else
                for group_option in ["profiles", "reference", "timeseries"]
                    haskey(ds.group, group_option) || continue
                    if haskey(ds.group[group_option], var_name)
                        return Array(ds.group[group_option][var_name])
                    end
                end
            end
        end
        error("Variables $var_names not found in the output netCDF file $filename.")
    end
end
nc_fetch(filename::String, var_name::String) = nc_fetch(filename, (var_name,))

"""
    is_face_variable(filename::String, var_name::String)

A `Bool` indicating whether the given variable is defined in a face,
or not (cell center).

TurbulenceConvection data are consistent, meaning that variables at
cell faces (resp. centers) have as dim `zf` (resp., `zc`).

PyCLES variables are inconsistent. All variables have as a dim `z`,
the cell face locations, but most of them (except for the statistical
moments of w) are actually defined at cell centers (`z_half`).
"""
function is_face_variable(filename::String, var_name::String)
    # PyCLES cell face variables
    pycles_face_vars = ["w_mean", "w_mean2", "w_mean3"]
    NCDataset(filename) do ds
        for group_option in ["profiles", "reference", "timeseries"]
            haskey(ds.group, group_option) || continue
            if haskey(ds.group[group_option], var_name)
                var_dims = dimnames(ds.group[group_option][var_name])
                if ("zc" in var_dims) | ("z_half" in var_dims)
                    return false
                elseif ("zf" in var_dims) | (var_name in pycles_face_vars)
                    return true
                elseif ("z" in var_dims) # "Inconsistent" PyCLES variables
                    return false
                else
                    error("Variable $var_name does not contain a vertical coordinate.")
                end
            end
        end
    end
end

"""
    get_stats_path(dir)

Given directory to standard LES or SCM output, fetch path to stats file.
"""
function get_stats_path(dir)
    stats = joinpath(dir, "stats")
    if !ispath(stats)
        nc_path = joinpath(dir, "*.nc")
        stat_files = glob(relpath(abspath(nc_path)))
        @assert length(stat_files) == 1 "$(length(stat_files)) stats files found with paths $nc_path"
        return stat_files[1]
    end
    try
        nc_path = joinpath(stats, "*.nc")
        stat_files = glob(relpath(abspath(nc_path)))
        @assert length(stat_files) == 1 "$(length(stat_files)) stats files found with paths $nc_path"
        return stat_files[1]
    catch e
        if isa(e, AssertionError)
            @warn "No unique stats netCDF file found in $stats. Extending search to other files."
            try
                stat_files = readdir(stats, join = true) # WindowsOS/julia relpath bug
                if length(stat_files) == 1
                    return stat_files[1]
                else
                    @error "No unique stats file found at $dir. The search returned $(length(stat_files)) results."
                end
            catch f
                if isa(f, Base.IOError)
                    @warn "Extended search errored with: $f"
                    return ""
                else
                    throw(f)
                end
            end
        else
            @warn "An error occurred retrieving the stats path at $dir. Throwing..."
            throw(e)
        end
    end
end

"""
    compute_mse(g_arr::Vector{Vector{FT}}, y::Vector{FT})::Vector{FT}
    compute_mse(g_mat::Matrix{FT}, y::Vector{FT})::Vector{FT}

Computes the L2-norm error of each vector, column or row of an array
with respect to a vector y.

Output:
 - The mse for each ensemble member.
"""
function compute_mse(g_arr::Vector{Vector{FT}}, y::Vector{FT})::Vector{FT} where {FT <: Real}
    diffs = [g - y for g in g_arr]
    errors = map(x -> dot(x, x) / length(x), diffs)
    return errors
end
function compute_mse(g_mat::Matrix{FT}, y::Vector{FT})::Vector{FT} where {FT <: Real}
    # Data are columns
    if size(g_mat, 1) == length(y)
        diffs = [g - y for g in eachcol(g_mat)]
        return map(x -> dot(x, x) / length(x), diffs)
        # Data are rows
    elseif size(g_mat, 2) == length(y)
        diffs = [g - y for g in eachrow(g_mat)]
        return map(x -> dot(x, x) / length(x), diffs)
    else
        throw(BoundsError("Dimension mismatch between inputs to `compute_error`."))
    end
end

"""
    penalize_nan(arr::Vector{FT}; penalization::FT = 1.0e5) where {FT <: AbstractFloat}

Substitutes all NaN entries in `arr` by a penalization factor.
"""
function penalize_nan(arr::Vector{FT}; penalization::FT = 1.0e5) where {FT <: AbstractFloat}
    return map(elem -> isnan(elem) ? penalization : elem, arr)
end

"""
    serialize_struct(s::T) where {T}

Serializes the given structure as a dictionary to
allow storage in JLD2 format.
"""
function serialize_struct(s::T) where {T}
    keys = propertynames(s)
    vals = getproperty.(Ref(s), keys)
    return Dict(zip(string.(keys), vals))
end

"""
    deserialize_struct(dict::Dict{String}, ::Type{T})

Deserializes the given dictionary and constructs a struct
of the given type with the dictionary values.
"""
deserialize_struct(dict::Dict{String}, ::Type{T}) where {T} = T(map(fn -> dict["$fn"], fieldnames(T))...)

"""
    jld2_path(root::String, identifier::Union{String, Int}, prefix::String)

Generates a JLD2 path, given a root path, an identifier and a prefix.
"""
function jld2_path(root::String, identifier::Union{String, Int}, prefix::String)
    return joinpath(root, "$(prefix)$(identifier).jld2")
end

scm_init_path(root, version; prefix = "scm_initializer_") = jld2_path(root, version, prefix)
scm_output_path(root, version; prefix = "scm_output_") = jld2_path(root, version, prefix)
scm_val_init_path(root, version; prefix = "scm_val_initializer_") = jld2_path(root, version, prefix)
scm_val_output_path(root, version; prefix = "scm_val_output_") = jld2_path(root, version, prefix)
ekobj_path(root, iter; prefix = "ekobj_iter_") = jld2_path(root, iter, prefix)

"""
    write_versions(versions::Vector{Int}, iteration::Int; outdir_path::String = pwd())

Writes versions associated with an EnsembleKalmanProcess iteration to a text file.
"""
function write_versions(versions::Vector{Int}, iteration::Int; outdir_path::String = pwd())
    open(joinpath(outdir_path, "versions_$(iteration).txt"), "w") do io
        for version in versions
            write(io, "$(version)\n")
        end
    end
end

"Returns the N-vector stored in `dict[key]`, or an N-vector of `nothing`"
function expand_dict_entry(dict, key, N)
    val = get(dict, key, nothing)
    r = isnothing(val) ? repeat([nothing], N) : val
    @assert length(r) == N
    r
end

"""
    get_entry(dict, key, default)

Calls `get` but logs whether the default is used.
"""
function get_entry(dict, key, default)
    if haskey(dict, key)
        return dict[key]
    else
        @info "Key $key not found in dictionary. Returning default value $default."
        return default
    end
end

"""
    change_entry!(dict, keys_and_value)

Changes the entry of a nested dictionary, giving a tuple of all its keys and the new value

Inputs:
 - `dict`           :: Parent dictionary with an arbitrary number of nested dictionaries.
 - `keys_and_value` :: Tuple of keys from the parent dictionary to the entry to be modified,
                     and the value to use to modify it.
"""
function change_entry!(dict, keys_and_value)
    function get_last_nested_dict(dict, keys)
        if length(keys) > 1
            return get_last_nested_dict(dict[keys[1]], keys[2:end])
        else
            return dict
        end
    end
    # Unpack
    value = keys_and_value[end]
    keys = keys_and_value[1:(end - 1)]
    # Modify entry
    last_dict = get_last_nested_dict(dict, keys)
    last_dict[keys[end]] = value
end

end # module


text:
using Test
@test true

ext: jl


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

struct EPOut{T<:Real} <: AbstractMetState
    μ::Vector{T}
    σ::Vector{T}
    av::Vector{T}
    va::Vector{T}
    # stoierr::T
    sol
    status::Symbol
    iter::Integer
end

status(out::EPOut) = out.status
curriter(out::EPOut) = out.iter


text:
export Activity,
    AT_GAME,
    AT_STREAMING,
    AT_LISTENING

"""
An [`Activity`](@ref)'s type. Available values are `AT_GAME`, `AT_STREAMING`, and
`AT_LISTENING`.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-types).
"""
@enum ActivityType AT_GAME AT_STREAMING AT_LISTENING AT_UNKNOWN  # Supposed to only go to 2.
@boilerplate ActivityType :lower

"""
Flags which indicate what an [`Activity`](@ref) payload contains.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-flags).
"""
@enum ActivityFlags begin
    AF_INSTANCE=1<<0
    AF_JOIN=1<<1
    AF_SPECTATE=1<<2
    AF_JOIN_REQUEST=1<<3
    AF_SYNC=1<<4
    AF_PLAY=1<<5
end
@boilerplate ActivityFlags :lower

"""
The start and stop times of an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-timestamps).
"""
struct ActivityTimestamps
    start::Union{DateTime, Missing}
    stop::Union{DateTime, Missing}
end
@boilerplate ActivityTimestamps :dict :docs :lower :merge

"""
The current party of an [`Activity`](@ref)'s player.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-party).
"""
struct ActivityParty
    id::Union{String, Missing}
    size::Union{Vector{Int}, Missing}
end
@boilerplate ActivityParty :dict :docs :lower :merge

"""
Images and hover text for an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-assets).
"""
struct ActivityAssets
    large_image::Union{String, Missing}
    large_text::Union{String, Missing}
    small_image::Union{String, Missing}
    small_text::Union{String, Missing}
end
@boilerplate ActivityAssets :dict :docs :lower :merge

"""
Secrets for Rich Presence joining and spectating of an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-secrets).
"""
struct ActivitySecrets
    join::Union{String, Missing}
    spectate::Union{String, Missing}
    match::Union{String, Missing}
end
@boilerplate ActivitySecrets :dict :docs :lower :merge

"""
A [`User`](@ref) activity.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object).
"""
struct Activity
    name::String
    type::ActivityType
    url::Union{String, Nothing, Missing}
    timestamps::Union{ActivityTimestamps, Missing}
    application_id::Union{Snowflake, Missing}
    details::Union{String, Nothing, Missing}
    state::Union{String, Nothing, Missing}
    party::Union{ActivityParty, Missing}
    assets::Union{ActivityAssets, Missing}
    secrets::Union{ActivitySecrets, Missing}
    instance::Union{Bool, Missing}
    flags::Union{Int, Missing}
end
@boilerplate Activity :dict :docs :lower :merge


text:
module NMR

"Integration range in the frequency domain, values given in ppm."
const Intrng = Tuple{Float64, Float64}

"""
Processed (frequency domain) NMR spectrum. Parameters follow the Bruker conventions:
- **re_ft**: real part of Fourier transform
- **im_ft**: imaginary part of Fourier transform
- **params**: dictionary of processing parameters, mostly from `proc` file for Bruker data
- **intrng**: list of integration ranges
- **title**: Spectrum title
"""
mutable struct ProcessedSpectrum
    re_ft :: Vector{Float64}
    im_ft :: Vector{Float64}
    params :: Dict{String, Any}
    intrng :: Array{Intrng,1}
    procno :: Int
    title :: String
end

Base.getindex(p::ProcessedSpectrum, n::Int) = p.intrng[n]
Base.getindex(p::ProcessedSpectrum, param::AbstractString) = p.params[param]
Base.getindex(p::ProcessedSpectrum, ::Colon) = p.re_ft
Base.getindex(p::ProcessedSpectrum, a::AbstractArray) = p.re_ft[a]
Base.view(p::ProcessedSpectrum, v) = @view p.re_ft[v]

"""
Unprocessed NMR experiment: FID and any associated processed spectra.

Following Bruker convention:
- **fid**: Acquired free induction decay (FID) signal
- **acqu**: Acquisition parameters, mostly from `acqu` file for Bruker data
- **procs**: Dictionary of associated processed spectra
- **default_proc**: Default processed spectrum to be used in operations that require a processed spectrum
- **name**: Experiment name
- **expno**: Experiment number if, e.g., part of a Bruker dataset
"""
mutable struct Spectrum
    fid :: Vector{Float64}
    acqu :: Dict{String, Any}
    procs :: Dict{Int,ProcessedSpectrum}
    default_proc :: Int
    name :: String
    expno :: Int
end

Base.getindex(s::Spectrum, i::Int) = s.procs[i]
Base.getindex(s::Spectrum, ::Colon) = s.procs[s.default_proc].re_ft
Base.getindex(s::Spectrum, a::AbstractArray) = s.procs[s.default_proc].re_ft[a]
Base.getindex(s::Spectrum, rng::Tuple{Float64,Float64}) = s[ppmtoindex(s,rng)]
Base.getindex(s::Spectrum, δ::Float64) = s[s.default_proc].re_ft[ppmtoindex(s, δ)]
Base.view(s::Spectrum, v) = @view s.procs[s.default_proc][v]
function Base.getindex(s::Spectrum, param::AbstractString)
    try
        s.acqu[param]
    catch err
        s[s.default_proc][param]
    end
end
Base.setindex!(s::Spectrum, d::AbstractArray, ::Colon) = (s[s.default_proc].re_ft .= d)
Base.setindex!(s::Spectrum, d::AbstractArray, r::UnitRange) = (s[s.default_proc].re_ft[r] .= d)
Base.setindex!(s::Spectrum, d, rng::Tuple{Float64, Float64}) = (s[s.default_proc].re_ft[ppmtoindex(s,rng)]=d)

Spectrum(fid :: Vector{Float64}, acqu :: Dict{Any, Any}, proc :: ProcessedSpectrum) = Spectrum(fid, acqu, Dict(1=>proc), 1, "", "")

include("utils.jl")
include("baseline.jl")
include("bruker/read.jl")
include("bruker/write.jl")
include("composition.jl")
include("decomposition.jl")
include("integration.jl")
include("interpolation.jl")
include("plotting.jl")
include("show.jl")

export Spectrum, ProcessedSpectrum, dump, # constructors, I/O
       plot, plot!, # plotting
       lsq_analyze, candidates, decompose, # decomposition
       ppmtoindex, hztoindex, ppmtohz, # unit conversion
       baseline_correct!, # processing
       limits, chemical_shifts, extract, copy!, # utility functions
       interpolate, resample, # interpolation
       integrate, # integration
       intrng, intrng_data, intrng_indices, intrng_shifts # integration regions
end


text:
# # 第4章
# 情報量基準の問題。殆ど証明でコーディングの問題は少ないです。
# 説明変数の組み合わせを求めるにあたって、Combinatrics.jlパッケージ（pythonでいうitertools）がとても役立ちます。
# また、最小値とそのインデックスを返すfindmin関数も有用です。
# ## 例題45
# ボストン不動産データセットを用いた線形回帰分析でのAIC・BICの例。
# 次のような関数をchap4.jl内で定義しています。
# ```julia
# function AIC(x,y)
#     N, k = size(x)
#     @assert N == length(y)
#     return N*log(RSS(x,y)/N) + 2k
# end
# function AIC_min(x,y)
#     _,p = size(x)
#     aic_min = Array{Float64}(undef,p)
#     aic_min_com = Array{Any}(undef,p)
#     for i in 1:p
#         combi = combinations(1:p,i) |> collect
#         aic_min[i],aic_index = [AIC(x[:,c],y) for c in combi] |> findmin
#         aic_min_com[i] = combi[aic_index]
#     end
#     return aic_min, aic_min_com
# end
# ```
using Joe
using ScikitLearn, Plots
using Plots.PlotMeasures # hide
Plots.reset_defaults() # hide
default( # hide
    titlefont  = font("JuliaMono", default(:titlefontsize),  ),  # hide
    guidefont  = font("JuliaMono",  default(:guidefontsize),  ), # hide
    tickfont   = font("JuliaMono", default(:tickfontsize),   ), # hide
    legendfont = font("JuliaMono",  default(:legendfontsize), ), # hide
    left_margin = 30px, # hide
    bottom_margin = 30px # hide
) # hide
@sk_import datasets: load_boston
df = load_boston()
X = df["data"]
X2 =X[:,vcat([1,3,5:end]...)]
y = df["target"] #目的変数

aicmin, aicindex = Joe.AIC_min(X,y)
bicmin, bicindex = Joe.BIC_min(X,y)
aicmin2, aicindex2 = Joe.AIC_min(X2,y)
bicmin2, bicindex2 = Joe.BIC_min(X2,y)

p=plot(xlabel="number of variables",ylabel="AIC or BIC values")
scatter!(p,aicmin,label="AIC")
scatter!(p,bicmin,label="BIC")
savefig(p,joinpath(@OUTPUT,"p45.svg")) # hide
# \fig{p45}
# ## 問46
# ボストン不動産データでのAIC値を最小にする組み合わせを記述したもの。
min_aic_boston,aicmin_index  = findmin(aicmin2)
println("Minimum AIC value:", min_aic_boston, ", set: ", aicindex2[aicmin_index])
# ## 問47
# BIC関数とAR2関数を定義して利用します。変数の数毎の最小値を計算する関数は
# AIC_minと同じ方法で定義しています。ただしAR2の場合は最大化しなければならないので注意が必要です。
# ```julia
# function BIC(x,y)
#     N,k = size(x)
#     @assert N == length(y)
#     return N*log(RSS(x,y)/N) + k*log(N)
# end
#
# function AR2(x,y)
#     N,k = size(x)
#     @assert N == length(y)
#     return 1 - RSS(x,y)/(N-k-1)/TSS(y)*(N-1)
# end
# ```
# どの基準で計算しても、今回の場合同じ変数の組み合わせが最適という結論になります。

#BIC
min_bic_boston,bicmin_index  = findmin(bicmin2)
println("Minimum BIC value: ", min_bic_boston, ", set: ", bicindex2[aicmin_index])
#AR2
ar2max2, ar2index2 = Joe.AR2_max(X2,y)
max_ar2_boston,ar2max_index  = findmax(ar2max2)
println("Maximum AR2 value: ", max_ar2_boston, ", set: ", ar2index2[aicmin_index])
# ## 問48
# 例題45でやったので、割愛します。


text:
# Command: Get Current URL
"""
	current_url(session::Session)::String

Get Current URL.
"""
function current_url(session::Session)::String
    @unpack addr, id = session
    response = HTTP.get("$addr/session/$id/url", [("Content-Type" => "application/json")])
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
using DataStructures
using Sockets
using Printf

CLOSED_MSG = "ECONNREFUSED"
RESET_MSG = "ECONNRESET"

Host = Union{IPAddr,String}

Banner = Union{String,Nothing}

struct PortResult
    status::String
    banner::Banner
end

struct Prefix
    symbol::String
    color
end

PREFIXES = Dict(
    "OPEN" => Prefix("[+]", :green),
    "CLOSED" => Prefix("[x]", :red),
    "FILTERED" => Prefix("[-]", :yellow)
)

function get_banner(socket)
    try
        @async (sleep(0.5); close(socket))
        banner = String(readuntil(socket, 0x0d))
        if banner === ""
            banner = "None"
        end
        banner = replace(banner, r"\n"  => s"")
        return banner
    catch e
        return "None"
    end
end

function connect_and_grab_banner(target, port)
    try
        socket = Sockets.connect(target, port)
        banner = get_banner(socket)
        close(socket)
        return PortResult("OPEN", banner)
    catch e
        if isa(e, Base.IOError)
            # Closed port
            if occursin(CLOSED_MSG, e.msg)
                return PortResult("CLOSED", nothing)
            end
            # Filtered port
            if occursin(RESET_MSG, e.msg)
                return PortResult("FILTERED", nothing)
            end
        end
    end
end

function timeout()
    sleep(1.5)
    return PortResult("FILTERED", nothing)
end

function scan_ports(target::Host, start::Int, finish::Int, only_open::Bool)
    openPorts = Dict{Int64,PortResult}()
    open_count = 0
    closed_count = 0
    filtered_count = 0
    start_time = time()
    @printf("PORT   STATUS          BANNER\n")
    for port = start:finish
        c = Channel(2)
        @async put!(c, connect_and_grab_banner(target, port))
        @async put!(c, timeout())
        result = take!(c)
        if result.status == "OPEN"
            open_count += 1
            push!(openPorts, port => result)
        end
        if result.status == "CLOSED"
            closed_count += 1
        end
        if result.status == "FILTERED"
            filtered_count += 1
        end
        if only_open && result.status != "OPEN"
            continue
        end
        prefix_buffer = IOBuffer();
        prefix = PREFIXES[result.status]
        printstyled(IOContext(prefix_buffer, :color => true), prefix.symbol, color=prefix.color)
        @printf("%-5d  %s %-8s    %s\n", port, String(take!(prefix_buffer)), result.status, result.banner)
    end
    elapsed_time = time() - start_time
    if (open_count + closed_count + filtered_count) == 0
        print("Looks like the host is down :(\n\n")
        exit(0)
    end
    print("\nTCP port scanning complete!\n")
    print("$open_count open ports, $closed_count closed and $filtered_count filtered\n")
    println("\nFinished in $elapsed_time seconds")
    return openPorts
end

max_stars_repo_name: thecrypters/nis

max_issues_repo_name: thecrypters/nis

max_forks_repo_name: thecrypters/julia-host-scanner


text:
using MeshCat, MeshCatMechanisms, RigidBodyDynamics

include(joinpath(pwd(),"models/kuka/kuka_utils.jl"))
urdf_path = joinpath(pwd(), "models/kuka/temp/kuka.urdf")
kuka = MeshCatMechanisms.parse_urdf(urdf_path, remove_fixed_tree_joints = true)
kuka_visuals = MeshCatMechanisms.URDFVisuals(urdf_path)

state = MechanismState(kuka)
state_cache = StateCache(kuka)
result = DynamicsResult(kuka)
result_cache = DynamicsResultCache(kuka)

vis = Visualizer()
mvis = MechanismVisualizer(kuka, kuka_visuals, vis[:base])
render(vis)

vis_el = visual_elements(kuka, kuka_visuals)
# set_alpha!(vis_el, 0.5)

mvis2 = MechanismVisualizer(kuka, vis[Symbol("shadow")])
MeshCatMechanisms._set_mechanism!(mvis2, vis_el)
MeshCatMechanisms._render_state!(mvis2)

setvisible!(vis[:base], false)
setvisible!(vis[:shadow], true)

function set_alpha!(visuals::Vector{VisualElement}, α)
    for el in visuals
        c = el.color
        c_new = RGBA(red(c),green(c),blue(c),α)
        el.color = c_new
    end
end

# Kuka iiwa arm parsed from URDF using RigidBodyDynamics.jl
struct Kuka{I, T} <: Model{I, T}
	n::Int
	m::Int
	d::Int

    state_cache1
    state_cache2
    state_cache3

    result_cache1
    result_cache2
    result_cache3
end

n = 14
m = 7
d = 0

results_cache1 = DynamicsResultCache(kuka)
results_cache2 = DynamicsResultCache(kuka)
results_cache3 = DynamicsResultCache(kuka)

state_cache1 = StateCache(kuka)
state_cache2 = StateCache(kuka)
state_cache3 = StateCache(kuka)

function f(model::Kuka, x, u, w)
	@error "not implemented"
end

model = Kuka{Midpoint, FixedTime}(
	n, m, d,
    state_cache1, state_cache2, state_cache3,
    results_cache1, results_cache2, results_cache3)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
module AxUtil

using Requires

include("arr.jl")
include("dpmeans.jl")
include("logging.jl")
include("math.jl")
include("flux.jl")  # must come after math
include("mcdiag.jl")
include("mcmc.jl")
include("misc.jl")
include("random.jl")

@init @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("cuda.jl")

end # module


text:
# various functions for analyzing the performance characteristics of RP trees

function local_cov_dim(eigs, eps)
    d = 1
    while !has_local_cov_dim(eigs, d, eps)
        d += 1
    end
    return d
end

has_local_cov_dim(eigs, d, eps) = sum(eigs[end-(d-1):end].^2) >= (1 - eps)*sum(eigs.^2)


text:
# function testNHBubbleX()
using CGDycore

OrdPoly = 4
OrdPolyZ=1
nx=20;
ny=2;
nz=40;
NF = nx *ny
# Caching
cache=CGDycore.Cache(OrdPoly, OrdPolyZ, nz, NF)
# Physical parameters
Param=CGDycore.PhysParameters(cache);

# Grid
lx=20000;
ly=2000;
x0=-10000;
y0=0;
Boundary = (;WE="Period", BT="Period")
Param.hS="";
Param.Grid=CGDycore.CartGrid(nx,ny,lx,ly,x0,y0,CGDycore.OrientFaceCart,Boundary,Param);
Param.TopoS="";

Param.H=10000;
Param.Grid.nz = nz
zP=zeros(nz,1);
z=zeros(nz+1,1);
dz=Param.H/nz;
zP[1]=dz/2;
for i=2:nz
  zP[i]=zP[i-1]+dz;
end
for i=2:nz+1
  z[i]=z[i-1]+dz;
end

Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.zP=zP;
Param.Grid.dz=dz;

# Discretization
OrdPoly=4;
OrdPolyZ=1;
(CG,Param)=CGDycore.Discretization(OrdPoly,OrdPolyZ,CGDycore.JacobiDG3,Param);
Param.HyperVisc=true;
Param.HyperDCurl=1.e4;
Param.HyperDGrad=1.e4;
Param.HyperDDiv=1.e4;
Param.Upwind=false

# Model
Param.ModelType="Curl";
#Param.ModelType="Div";
Param.Coriolis=false;
Param.Thermo="";
Param.Source=false;
Param.Damping=false;
Param.Flat=false;
Param.Coriolis=false;
Param.Buoyancy=true;
Param.xC0=0;
Param.zC0=2000;
Param.rC0=2000;
Param.Th0=300;
Param.DeltaTh=2;
Param.uMax=20;
Param.vMax=0;
Param.NBr=1.e-2;
Param.Equation="Compressible";
Param.RefProfile=false

# Initial conditions
Param.NumV=5;
U=zeros(CG.NumG,nz,Param.NumV);
Param.ProfRho="WarmBubble2D";
Param.ProfTheta="WarmBubble2D";
Param.ProfVel="Const";
Param.RhoPos=1;
Param.uPos=2;
Param.vPos=3;
Param.wPos=4;
Param.ThPos=5;
U[:,:,Param.RhoPos]=CGDycore.Project(CGDycore.fRho,CG,Param);
(U[:,:,Param.uPos],U[:,:,Param.vPos])=CGDycore.ProjectVec(CGDycore.fVel,CG,Param);
U[:,:,Param.ThPos]=CGDycore.Project(CGDycore.fTheta,CG,Param).*U[:,:,Param.RhoPos];

# Output
Param.vtkFileName="Bubble";
Param.vtk=0;
vtkGrid=CGDycore.vtkCGGrid(CG,CGDycore.TransCart,CGDycore.Topo,Param);
Param.cNames = [
  "Rho",
  "u",
  "v",
  "w",
  "Th"
]

time=[0];
#IntMethod="Rosenbrock";
IntMethod="RungeKutta";
if strcmp(IntMethod,"Rosenbrock")
  dtau=500;
else
  dtau=.4;
end
nIter=10000;
Param.RK=CGDycore.RungeKuttaMethod("RK4");
Param.ROS=CGDycore.RosenbrockMethod("ROSRK3");
CFL=0.125;
time=0;
SimTime=1000;
PrintTime=100;
nIter=SimTime/dtau;
PrintInt=PrintTime/dtau;
# Print initial conditions
Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
#

OP=CG.OrdPoly+1;
Param.CacheF1=zeros(OP,OP,NF,nz+1);
Param.CacheF2=zeros(OP,OP,NF,nz+1);
Param.CacheF3=zeros(OP,OP,NF,nz+1);
Param.CacheF4=zeros(OP,OP,NF,nz+1);
Param.CacheF5=zeros(OP,OP,NF,nz+1);
Param.CacheF6=zeros(OP,OP,NF,nz+1);
Param.CacheC1 = view(Param.CacheF1,:,:,:,1:nz)
Param.CacheC2 = view(Param.CacheF2,:,:,:,1:nz)
Param.CacheC3 = view(Param.CacheF3,:,:,:,1:nz)
Param.CacheC4 = view(Param.CacheF4,:,:,:,1:nz)
Param.CacheC5 = view(Param.CacheF5,:,:,:,1:nz)
Param.CacheC6 = view(Param.CacheF6,:,:,:,1:nz)
Param.Cache1=zeros(CG.NumG,nz)
Param.Cache2=zeros(CG.NumG,nz)
Param.Cache3=zeros(CG.NumG,nz)
Param.Cache4=zeros(CG.NumG,nz)
Param.Pres=zeros(OP,OP,NF,nz)
Param.KE=zeros(OP,OP,NF,nz)
Param.FCG=zeros(OP,OP,NF,nz,size(U,3))
Param.fV=zeros(size(U)..., Param.RK.nStage);
Param.Vn=zeros(size(U));
Param.RhoCG=zeros(OP,OP,NF,nz)
Param.v1CG=zeros(OP,OP,NF,nz)
Param.v2CG=zeros(OP,OP,NF,nz)
Param.wCG=zeros(OP,OP,NF,nz+1)
Param.wCCG=zeros(OP,OP,NF,nz+1)
Param.ThCG=zeros(OP,OP,NF,nz)

str = IntMethod
if str == "RungeKutta"
    for i=1:nIter
      @info "Iteration: $i"
      CGDycore.RungeKuttaExplicit!(U,dtau,CGDycore.FcnNHCurlVec!,CG,Param);
      if mod(i,PrintInt)==0
        Param.vtk=CGDycore.vtkOutput(U,vtkGrid,CG,Param);
      end
    end
end


text:
# wrapper functions for convienence
nnls(A,b) = nonneg_lsq(A,b;alg=:fnnls)
cdsc(A,b) = nonneg_lsq(A,b;alg=:cd,max_iter = 2000)

#This test fails for this method...
# Solve A*x = b for x, subject to x >=0
# A = [ 0.53879488  0.65816267
#       0.12873446  0.98669198
#       0.24555042  0.00598804
#       0.80491791  0.32793762 ]

# b = [0.888,  0.562,  0.255,  0.077]

# # Test that nnls produces the same solution as scipy
# x = [0.15512102, 0.69328985] # approx solution from scipy
# @test norm(cdsc(A,b)-x) < 1e-5

## A second test case
A2 = [ -0.24  -0.82   1.35   0.36   0.35
       -0.53  -0.20  -0.76   0.98  -0.54
        0.22   1.25  -1.60  -1.37  -1.94
       -0.51  -0.56  -0.08   0.96   0.46
        0.48  -2.25   0.38   0.06  -1.29 ]
b2 = [-1.6,  0.19,  0.17,  0.31, -1.27]
x2 = [2.2010416, 1.19009924, 0.0, 1.55001345, 0.0]
@test norm(cdsc(A2,b2)-x2) < 1e-5

#X - Test a bunch of random cases against python
# Methods internal to NNLS.jl have already been validated

# Test a bunch of random cases against nnls
for i = 1:10
  m,n = rand(1:10),rand(1:10)
  A4 = randn(m,n)
  b4 = randn(m)
  x4 = nnls(A4,b4)
  @test norm(cdsc(A4,b4)-x4) < 1e-5
end


text:
@testset "pwmeva.jl" begin

    θ = [0.0, 0.0, 0.1]

    pd = GeneralizedExtremeValue(θ[1], exp(θ[2]), θ[3])

    y = [0]

    fm = Extremes.pwmEVA{BlockMaxima, GeneralizedExtremeValue}(Extremes.BlockMaxima(Variable("y", y)), θ)

    @testset "getdistribution(fittedmodel)" begin
        @test Extremes.getdistribution(fm)[] == pd
    end

    @testset "quantile(fm, p)" begin
        # p outside of [0, 1] throws
        @test_throws AssertionError Extremes.quantile(fm, -1)

        # Test with known values
        @test quantile(fm, .99)[] ≈ quantile(pd,.99)

    end

    @testset "returnlevel(fm, returnPeriod)" begin

        # returnPeriod < 0 throws
        @test_throws AssertionError returnlevel(fm, -1)

        # Test with known values
        @test returnlevel(fm, 100).value[] ≈ quantile(pd, 1-1/100)
    end

    @testset "parametervar(fm, nboot)" begin
        # nboot < 0 throws
        @test_throws AssertionError Extremes.parametervar(fm, -1)

        # Tested in cint(fm)
    end

    @testset "quantilevar(fm, level, nboot)" begin
        # Tested in cint(fm)
    end

    @testset "cint(fm)" begin
        # returnPeriod < 0 throws
        @test_throws AssertionError cint(ReturnLevel(Extremes.BlockMaximaModel(fm), -1, [1.0]), 0.95)

        # confidencelevel not in [0, 1]
        @test_throws AssertionError cint(ReturnLevel(Extremes.BlockMaximaModel(fm), 1, [1.0]), -1)

        niter = 1000

        cover = falses(niter,3)

        Threads.@threads  for i in 1:niter
            y = rand(pd, 100)
            fm = gevfitpwm(y)
            ci = cint(fm, .95 ,1000)
            for j in eachindex(θ)
                cover[i,j] = ci[j][1] < θ[j] < ci[j][2]
            end
        end

        @test all(count.(eachcol(cover))/niter .> .85)
    end

    @testset "cint(fm, returnPeriod, confidencelevel)" begin
        # confidencelevel not in [0, 1]
        @test_throws AssertionError cint(ReturnLevel(Extremes.BlockMaximaModel(fm), 1, [1.0]), -1)

        # Test with known values

        θ = [100; log(5); .1]

        pd = GeneralizedExtremeValue(θ[1],exp(θ[2]),θ[3])

        q = quantile(pd, 1-1/100)

        niter = 1200

        cover = falses(niter)

        Threads.@threads  for i in 1:niter
            y = rand(pd, 200)
            fm = gevfitpwm(y)
            rl = returnlevel(fm, 100)
            ci = cint(rl, .95 , 500)
            cover[i] = ci[][1] < q < ci[][2]
        end

        @test count(cover)/niter > .85
    end

    @testset "showfittedEVA(io, obj, prefix)" begin
        # print does not throw
        buffer = IOBuffer()
        @test_logs Extremes.showfittedEVA(buffer, fm, prefix = "\t")

    end

end


text:
# This is a convenience function for libraries like MaBLAS and PaddedMatrices to use explicitly, and for others to force more type inference in precompilation.

const GEMMLOOPSET = loopset(
    :(for n ∈ 1:N, m ∈ 1:M
          Cₘₙ = zero(eltype(C))
          for k ∈ 1:K
              Cₘₙ += A[m,k] * B[k,n]
          end
          C[m,n] += Cₘₙ
      end)
);

function matmul_params(rs::Int, rc::Int, cls::Int)
    set_hw!(GEMMLOOPSET, rs, rc, cls)
    order = choose_order(GEMMLOOPSET)
    order[5], last(order)
end
@generated function matmul_params(::StaticInt{RS}, ::StaticInt{RC}, ::StaticInt{CLS}) where {RS,RC,CLS}
    mᵣ, nᵣ = matmul_params(RS, RC, CLS)
    Expr(:tuple, Expr(:call, Expr(:curly, :StaticInt, mᵣ)), Expr(:call, Expr(:curly, :StaticInt, nᵣ)))
end
matmul_params() = matmul_params(register_size(), register_count(), cache_linesize())


text:
using Test
using QuantumPropagators
using QuantumControlBase
using QuantumControlBase.Shapes: flattop
using UnicodePlots

@testset "propagate TLS" begin

    # from the first Krotov.jl example

    SHOWPLOT = false

    #! format: off
    σ̂_z = ComplexF64[1 0; 0 -1];
    σ̂_x = ComplexF64[0 1; 1  0];
    #! format: on

    ϵ(t) = 0.2 * flattop(t, T=5, t_rise=0.3, func=:blackman)

    function hamiltonian(Ω=1.0, ϵ=ϵ)
        Ĥ₀ = -0.5 * Ω * σ̂_z
        Ĥ₁ = σ̂_x
        return (Ĥ₀, (Ĥ₁, ϵ))
    end

    function ket(label)
        #! format: off
        result = Dict(
            "0" => Vector{ComplexF64}([1, 0]),
            "1" => Vector{ComplexF64}([0, 1]),
        )
        #! format: on
        return result[string(label)]
    end

    H = hamiltonian()

    obj = Objective(initial_state=ket(0), generator=H, target_state=ket(1))

    tlist = collect(range(0, 5, length=500))

    states = propagate_objective(obj, tlist, storage=true)

    pops = abs.(states) .^ 2
    pop0 = pops[1, :]

    SHOWPLOT && println(lineplot(tlist, pop0, ylim=[0, 1], title="0 population"))

    @test length(pop0) == length(tlist)
    @test pop0[1] ≈ 1.0
    @test abs(pop0[end] - 0.95146) < 1e-4
    @test minimum(pop0) < 0.9

end


text:
const _RotationEncodingTables = Dict()

"""
    rotation_encoding_table(T, nbits)
    rotation_encoding_table(T, X::AbstractUnitRange)

Build the rotation-invariant encoding table. If global runtime cache exists, then reuse it.
"""
rotation_encoding_table(::Type{T}, nbits::Int) where T = rotation_encoding_table(T, 0:2^nbits-1)
function rotation_encoding_table(::Type{T}, X::AbstractUnitRange) where T
    d = _RotationEncodingTables
    k = (T, minimum(X), maximum(X))
    haskey(d, k) && return d[k]
    lookup = d[k] = RotationEncodingTable{T}(X)
    return lookup
end

"""
    RotationEncodingTable{T<:Unsigned}(X)

The lookup table for the quotient space constructed using `bitrotate` equivalance relation
on bit representation of datatype `T`.

By definition, ``a`` is equivalant to ``b``(``a ~ b``) if there exists ``n`` such that
`bitrotate(a, n) == b`.

This equivalance class gives a unique encoding map from ``x∈X`` to ``a``, where ``a`` is the
minimal value of the equivalance class ``[x]``. For instance, the equivalance class of
`0b10101001` (UInt8) is

```julia
0b10101001 # 169
0b01010011 # 83
0b10100110 # 166
0b01001101 # 77
0b10011010 # 154
0b00110101 # 53
0b01101010 # 106
0b11010100 # 212
```

thus `f(169) == 0x35 == 53`.

```jldoctest
julia> using LocalBinaryPatterns: RotationEncodingTable

julia> X = RotationEncodingTable{UInt8}(0:255);

julia> X[169] # 53
0x35

julia> X = RotationEncodingTable{UInt8}(128:255);

julia> X[169] # 154
0x9a
```

The values in the lookup table is implemented in lazy cache manner so the real computation
only happens once at the first retrieval.

This lookup table is used to build an efficient implementation of rotation-invariant local
binary pattern [1].

# References

- [1] T. Ojala, M. Pietikäinen, and T. Mäenpää, “A Generalized Local Binary Pattern Operator for Multiresolution Gray Scale and Rotation Invariant Texture Classification,” in _Advances in Pattern Recognition — ICAPR 2001, vol. 2013, S. Singh, N. Murshed, and W. Kropatsch, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg_, 2001, pp. 399–408. doi: 10.1007/3-540-44732-6_41.
"""
struct RotationEncodingTable{T<:Unsigned,AT<:AbstractUnitRange,LAT} <: AbstractVector{T}
    X::AT
    Y::LAT
end
function RotationEncodingTable{T}(X::AbstractUnitRange) where T
    if !((maximum(X) <= typemax(T)) && (minimum(X) >= typemin(T)))
        throw(ArgumentError("range $(first(X)):$(last(X)) can't be uniquely represented by type $T."))
    end

    if log2(length(X)) >= 26
        datasize = round(length(X)/1024/1024/1024, digits=3)
        @warn "$(datasize)GB data will be created."
    end

    Y = Vector{Union{T,Missing}}(missing, length(X))
    RotationEncodingTable{T,typeof(X),typeof(Y)}(X, Y)
end

@inline Base.axes(X::RotationEncodingTable) = (first(X.X):last(X.X), )
@inline Base.size(X::RotationEncodingTable) = (length(X.X)-1)
@inline function Base.getindex(X::RotationEncodingTable{T}, x::Int) where T
    i = x-first(axes(X, 1))+1
    v = X.Y[i]
    if ismissing(v)
        v = minimum_bitrotate(T(x), X.X)
        X.Y[i] = v
    end
    return T(v)
end

function minimum_bitrotate(x::T, X::AbstractUnitRange) where T<:Unsigned
    function bit_rotate_project_to(x, n, X)
        v = bitrotate(x, n)
        return v in X ? v : x
    end
    minimum(bit_rotate_project_to(x, n, X) for n in 0:8sizeof(T)-1)
end


text:
###############################################################################
#
# Calculate the electron transport rate
#
###############################################################################
"""
    leaf_ETR!(photo_set::AbstractPhotoModelParaSet{FT},
              leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}

Update the electron transport variables in the leaf struct, given
- `photo_set` [`C3ParaSet`](@ref) or [`C4ParaSet`](@ref) type struct
- `leaf` [`Leaf`](@ref) type struct
"""
function leaf_ETR!(
            photo_set::C3Cytochrome{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    @unpack APAR, K_P1, Vqmax, n_C, n_L, p_i, α_1, Γ_star, φ_P1_max = leaf;
    @unpack Eff_1, Eff_2 = photo_set;

    # TODO note to run this function before rubisco_limited_rate!
    # TODO Is Vqmax necessary here, or we use Jmax as well
    # TODO Johnson and Berry used PAR here
    leaf.J_P700_j = Vqmax * APAR / (Vqmax / α_1 / φ_P1_max + APAR);
    leaf.η = 1 - n_L / n_C + (3*p_i + 7*Γ_star) / (Eff_1*p_i + Eff_2*Γ_star) / n_C;
    leaf.J_P680_j = leaf.J_P700_j / leaf.η;

    leaf.J_pot = leaf.J_P680_j;
    leaf.J = leaf.J_P680_j;

    return nothing
end

function leaf_ETR!(
            photo_set::C3ParaSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    @unpack APAR, maxPSII, Jmax, PSII_frac = leaf;

    _Jp = PSII_frac * maxPSII * APAR;
    _J  = min(_Jp, Jmax);

    leaf.J_pot = _Jp;
    leaf.J     = _J;

    return nothing
end

function leaf_ETR!(
            photo_set::C4ParaSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    @unpack APAR, maxPSII, PSII_frac = leaf

    _Jp = PSII_frac * maxPSII * APAR;

    leaf.J_pot = _Jp;
    leaf.J     = _Jp;

    return nothing
end


text:
using Statistics, SpecialFunctions

"""
    tnvar(a, b)

Variance of the standard normal distribution truncated to the interval [a,b].
"""
function tnvar(a::Real, b::Real)
    if a == b
        return zero(middle(a, b))
    elseif a < b
        m1 = tnmean(a, b)
        m2 = √tnmom2(a, b)
        return (m2 - m1) * (m2 + m1)

        # m1 = tnmean(a, b)
        # @assert a ≤ m1 ≤ b
        # return tnmom2c(m1, a, b)
    else
        return oftype(middle(a, b), NaN)
    end
end

"""
    tnvar(a, b, μ, σ)

Variance of the truncated normal distribution, where μ, σ are the mean and
standard deviation of the untruncated distribution.
"""
function tnvar(a, b, μ, σ)
    α = (a - μ) / σ
    β = (b - μ) / σ
    return tnvar(α, β) * σ^2
end


text:
# +
import Revise
import GLRenderer
import PoseComposition
import Rotations
import Images
I = Images

R = Rotations
P = PoseComposition
GL = GLRenderer
# -

obj_path = joinpath(@__DIR__, "035_power_drill/textured_simple.obj")
texture_path = joinpath(@__DIR__, "035_power_drill/texture_map.png")

# +
camera_intrinsics = GL.CameraIntrinsics(
    640, 480,
    1000.0, 1000.0,
    320.0, 240.0,
    0.01, 5.0
)

renderer = GL.setup_renderer(camera_intrinsics, GL.TextureMode())

mesh_data = GL.get_mesh_data_from_obj_file(obj_path; tex_path=texture_path)
GL.load_object!(renderer, mesh_data)

@time rgb_image, depth_image = GL.gl_render(
    renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))],
    P.IDENTITY_POSE;
#     x_min=100, x_max=camera_intrinsics.width, y_min=100, y_max=camera_intrinsics.height
)
GL.view_depth_image(depth_image)
GL.view_rgb_image(rgb_image)

# +
camera_intrinsics = GL.CameraIntrinsics(
    200, 200,
    1000.0, 1000.0,
    100.0, 100.0,
    0.01, 5.0
)
renderer = GL.setup_renderer(camera_intrinsics, GL.TextureMode())

mesh_data = GL.get_mesh_data_from_obj_file(obj_path; tex_path=texture_path)
GL.load_object!(renderer, mesh_data)

@time rgb_image, depth_image = GL.gl_render(
    renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))],
    P.IDENTITY_POSE;
#     x_min=100, x_max=camera_intrinsics.width, y_min=100, y_max=camera_intrinsics.height
)
GL.view_depth_image(depth_image)
GL.view_rgb_image(rgb_image)

# +
camera_intrinsics = GL.CameraIntrinsics(
    640, 480,
    1000.0, 1000.0,
    320.0, 240.0,
    0.01, 5.0
)

renderer = GL.setup_renderer(camera_intrinsics, GL.TextureMode())

mesh_data = GL.get_mesh_data_from_obj_file(obj_path; tex_path=texture_path)
GL.load_object!(renderer, mesh_data)

@time rgb_image, depth_image = GL.gl_render(
    renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))],
    P.IDENTITY_POSE;
#     x_min=100, x_max=camera_intrinsics.width, y_min=100, y_max=camera_intrinsics.height
)
GL.view_depth_image(depth_image)
GL.view_rgb_image(rgb_image)

# +
camera_intrinsics = GL.CameraIntrinsics(
    640, 480,
    1000.0, 1000.0,
    320.0, 240.0,
    0.01, 5.0
)

GL.set_intrinsics!(renderer,camera_intrinsics)
@time rgb_image, depth_image = GL.gl_render(
    renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))],
    P.IDENTITY_POSE;
#     x_min=100, x_max=camera_intrinsics.width, y_min=100, y_max=camera_intrinsics.height
)
GL.view_depth_image(depth_image)
GL.view_rgb_image(rgb_image)
# -

GL.view_rgb_image(rgb_image)[200:400,200:400]

# +
camera_intrinsics = GL.CameraIntrinsics(
    100, 100,
    1000.0, 1000.0,
    320.0 - 200.0, 240.0 - 200.0,
    0.01, 5.0
)

@time GL.set_intrinsics!(renderer,camera_intrinsics)
@time rgb_image2, depth_image2 = GL.gl_render(
    renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))],
    P.IDENTITY_POSE;
#     x_min=100, x_max=camera_intrinsics.width, y_min=100, y_max=camera_intrinsics.height
)
GL.view_rgb_image(rgb_image2)
# -


text:
using DataStructures
function move(cmds)
    x = y = z = 0
    for cmd in cmds
        if cmd == "e"
            x += 1
            y -= 1
        elseif cmd == "se"
            z += 1
            y -= 1
        elseif cmd == "sw"
            x -= 1
            z += 1
        elseif cmd == "w"
            x -= 1
            y += 1
        elseif cmd == "nw"
            y += 1
            z -= 1
        elseif cmd == "ne"
            x += 1
            z -= 1
        end
    end
    return (x,y,z)
end

function parseline(str)
    cmd :: Array{String, 1} = []
    i = 1
    while i <= length(str)
        if str[i] == 'n' || str[i] == 's'
            push!(cmd, str[i:i+1])
            i += 2
        else
            push!(cmd, str[i:i])
            i += 1
        end
    end
    cmd
end

function neighbours(x, y, z)
    return [(x+1,y-1,z), (x,y-1,z+1), (x-1,y,z+1), (x-1,y+1,z), (x,y+1,z-1), (x+1,y,z-1)]
end

# prints solution of both parts
function x()
    f = open("24.txt")
    s = Set{Tuple{Int,Int,Int}}()
    for i=1:344
        coord = move(parseline(readline(f)))
        if coord in s
            delete!(s, coord)
        else
            push!(s, coord)
        end
    end
    println(length(s))
    for i=1:100
        t = copy(s)
        for x = -100:100, y = -100:100
            z = -x-y
            N = neighbours(x, y, z)
            isblack = (x,y,z) in s
            c = 0
            for n in N
                if n in s
                    c += 1
                end
            end
            if isblack && (c == 0 || c > 2)
                delete!(t, (x,y,z))
            elseif !isblack && c == 2
                push!(t, (x,y,z))
            end
        end
        s = t
    end
    println(length(s))
end


text:
@testset "SOAP Descriptor" begin

@info("Testing SOAP Descriptor for DC Si with cutoff=6.5, n=4, l=4")
using JuLIP
using QUIP, Test

at = bulk(:Si, cubic=true)
desc = soap(at, 6.5, n_max=4, l_max=4)
# Reference is from the output of the QUIP command as follows
# quip atoms_filename=test.xyz descriptor_str="soap cutoff=6.5 l_max=4 n_max=4 atom_sigma=0.5 normalise=T n_Z=1 Z={14} n_species=1 species_Z={14}"
soap_ref = [ 0.23864044377742688E+000,  0.64340429976653843E-038,  0.92117979653018819E-038,
             0.48626301215687845E-005,  0.18876337954462125E-005,  0.29409765463423732E+000,
            -0.11518886345720650E-036,  0.11050218320476954E-035,  0.11763958794274931E-002,
             0.45589948982910106E-003,  0.18122123201805027E+000,  0.11109715229581757E-033,
             0.20740200612211439E-033,  0.14230028097299613E+000,  0.55054202072415959E-001,
             0.51605689877931338E+000, -0.13873220359922910E-035,  0.53089215594650063E-036,
            -0.50839942118860665E-003, -0.50804299476442307E-003,  0.44970747195745586E+000,
             0.56664079661428650E-034, -0.62628772367525336E-034, -0.86970568361713990E-001,
            -0.86763453783360533E-001,  0.55798321223812897E+000,  0.17093065575762176E-033,
             0.32536812908176952E-031,  0.26577177885526940E-001,  0.68368050294466742E-001,
             0.37916397221481701E-002,  0.29040731747973391E-037, -0.11050630810221818E-037,
             0.11135545482166505E-004,  0.46298939173893863E-005,  0.33041486666568242E-002,
            -0.88957883103660356E-036, -0.32224333298567884E-036,  0.19049288398824812E-002,
             0.79069210886284307E-003,  0.57978317308968498E-002, -0.67817971363747986E-035,
             0.24028731016100991E-034, -0.82324728991149839E-003, -0.88112752017744061E-003,
             0.30121742054713143E-004,  0.23815358022366506E-036,  0.18522830192896921E-037,
             0.12750339845444976E-004,  0.56779863069817660E-005]
soap_now = vcat(desc[1,:]...)
println("QUIP.jl SOAP:",soap_now)
println("Reference:",soap_ref)
println(@test soap_now  ≈  soap_ref)

end


text:
# Entry decision

# When capacities are not constrained, set to this value
const CapacityInf = 1e8;

"""
	$(SIGNATURES)

Abstract type for entry decision protocols.
"""
abstract type AbstractEntryDecision{F1 <: Real} <: ModelObject end

"""
	$(SIGNATURES)

Abstract type for switches governing entry decision protocols.
"""
abstract type AbstractEntrySwitches{F1 <: Real} end

## ------------  Sequential assignment

"""
	$(SIGNATURES)

Switches governing default entry protocol. This has multiple locations, capacity constraints, and sequential entry decisions (simultaneous with work/study decisions).

`valueLocal` must be fixed at zero if there is only one location.
"""
Base.@kwdef mutable struct EntryDecisionSwitches{F1} <: AbstractEntrySwitches{F1}
    nTypes :: Int
    nColleges :: Int
    nLocations :: Int
    # Min entry prob for each (feasible) college
    minEntryProb :: F1 = 0.01
    # Max entry prob for all colleges jointly
    maxEntryProb :: F1 = 0.99
    "Fix entry rates to match data?"
    fixEntryProbs :: Bool = false
    "Preference shock scale parameter"
    entryPrefScale :: F1 = 1.0
    "Calibrate preference shock scale parameter?"
    calEntryPrefScale :: Bool = true
    # Value of attending local college
    valueLocal :: F1 = 1.0
    calValueLocal :: Bool = true
    # Each student has this mass
    typeMass_jlM :: Matrix{F1}
    # College capacity (in units of typeMass)
    capacity_clM :: Matrix{F1}
    # Indices of colleges that can only be attended locally
    localOnlyIdxV :: Vector{Int} = Vector{Int}()
end

"""
	$(SIGNATURES)

Sequential entry protocol.
"""
mutable struct EntryDecision{F1} <: AbstractEntryDecision{F1}
    objId :: ObjectId
    pvec :: ParamVector
    entryPrefScale :: F1
    valueLocal :: F1
    switches :: EntryDecisionSwitches{F1}
end

## ----------------  Results object

"""
	$(SIGNATURES)

Abstract type for results from entry decisions.
"""
abstract type AbstractEntryResults{F1 <: Real} end

"""
	$(SIGNATURES)

Entry results for multiple locations. Records:
- `fracEnter_jlcM`: Fraction of students in (j,l) who attend college `c`
- `fracLocal_jlcM`: The same, but for the local college `c` (if any).
- `fracEnterBest_jlcM`: Fraction of students in (j,l) who attend college `c` and for who this is the best available college.
- `enroll_clM`: total enrollment at college (c, l).
- `enrollLocal_clM`: local enrollment at college (c, l).
- `enrollBest_clM`: enrollment of students for who this is the best available college.
- `eVal_jlM`: expected value of students in (j,l) at the start of the process.
"""
struct EntryResults{F1} <: AbstractEntryResults{F1}
    switches :: AbstractEntrySwitches{F1}
    fracEnter_jlcM :: Array{F1, 3}
    fracLocal_jlcM :: Array{F1, 3}
    fracEnterBest_jlcM :: Array{F1, 3}
    # Expected values of types
    eVal_jlM :: Matrix{F1}
    # Mass of enrollment in each college
    enroll_clM :: Matrix{F1}
    enrollLocal_clM :: Matrix{F1}
    enrollBest_clM :: Matrix{F1}
end

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


text:
module ModuleCamera
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using StaticArrays
using GeometryTypes
export ascertain_pose, obtain_intrinsics
include("camera.jl")
end


text:
@testset "quadratures" begin
    for T in (Float32, Float64, BigFloat)
        x, w = legendregauss(T, 1)
        @test iszero(x)
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 1, Bennu.LeftEndPoint())
        @test x ≈ [-one(T)]
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 1, Bennu.RightEndPoint())
        @test x ≈ [one(T)]
        @test w ≈ [2one(T)]

        x, w = legendregauss(T, 2, Bennu.LeftEndPoint())
        @test x ≈ [-one(T); T(1 // 3)]
        @test w ≈ [T(1 // 2); T(3 // 2)]

        x, w = legendregauss(T, 2, Bennu.RightEndPoint())
        @test x ≈ [T(-1 // 3); one(T)]
        @test w ≈ [T(3 // 2); T(1 // 2)]

        err = ErrorException("Must have at least two points for both ends.")
        @test_throws err legendregauss(T, 1, Bennu.BothEndPoint())

        a, b = Bennu.legendrecoefficients(T, 100)
        err = ErrorException("No convergence after 1 iterations " *
                             "(try increasing maxiterations)")
        @test_throws err Bennu.gaussrule(-one(T), one(T), a, b, Bennu.BothEndPoint(), 1)
    end

    x1, w1 = legendregauss(23)
    x2, w2 = FastGaussQuadrature.gausslegendre(23)

    @test x1 ≈ x2
    @test w1 ≈ w2

    x1, w1 = legendregausslobatto(33)
    x2, w2 = FastGaussQuadrature.gausslobatto(33)

    @test x1 ≈ x2
    @test w1 ≈ w2

    x1, w1 = legendregauss(13, Bennu.LeftEndPoint())
    x2, w2 = FastGaussQuadrature.gaussradau(13)

    @test x1 ≈ x2
    @test w1 ≈ w2
end


text:
using KahanSummation;
using SpecialFunctions;

# support functions aware of floating-point arithmetic
function log1mexp(x)
  if x > -log(2)
    log(-expm1(x))
  else
    log1p(-exp(x))
  end
end

function logexpm1(x)
  if x >= 37
    x
  elseif x >= 19
    x - exp(-x)
  else
    log(expm1(x))
  end
end

function log1pexp(x)
  if x <= -37
    exp(x)
  elseif x <= 18
    log1p(exp(x))
  elseif x <= 33.3
    x + exp(-x)
  else
    x
  end
end

# these functions are the log-addends in the transition probability summation
# depending on parameter values, they are designed to reduce floating-point
# rounding error as much as possible
function log_addend_equal_rates_v1(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + (i + j - 2 * h) * ω + h * log((1 - θ) / (1 + θ))
end

function log_addend_equal_rates_v2(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + (i + j - 2 * h) * ω + h * log((θ - 1) / (θ + 1))
end

function log_addend_birth_rate_greater_v1(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + i * ω +
  (i + j - 2 * h) * (log1mexp(θ) - log1mexp(ω + θ)) +
  h * log(-expm1(θ - ω) / expm1(θ + ω))
end

function log_addend_birth_rate_greater_v2(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + i * ω +
  (i + j - 2 * h) * (log1mexp(θ) - log1mexp(ω + θ)) +
  h * log(expm1(θ - ω) / expm1(θ + ω))
end

function log_addend_death_rate_greater_v1(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + j * ω +
  (i + j - 2 * h) * (log1mexp(θ) - log1mexp(ω + θ)) +
  h * log(-expm1(θ - ω) / expm1(θ + ω))
end

function log_addend_death_rate_greater_v2(h, i, j, θ, ω)
  log(i) + lgamma(i + j - h) - lgamma(1 + h) - lgamma(i + 1 - h) -
  lgamma(j + 1 - h) + j * ω +
  (i + j - 2 * h) * (log1mexp(θ) - log1mexp(ω + θ)) +
  h * log(expm1(θ - ω) / expm1(θ + ω))
end

# naive summation algorithm with compensated summation
function log_trans_prob(
  i::F,
  j::F,
  t::F,
  λ::F,
  μ::F
)::F where {
  F <: AbstractFloat
}
  # anything lower than ϵ (subnormal numbers) are treated as zero
  ϵ = floatmin(F)

  # a = max(i, j)
  # b = min(i, j)
  a, b = (j <= i) ? (i, j) : (j, i)

  if (t < ϵ) || ((λ < ϵ) && (μ < ϵ))
    # time is zero or both parameters are zero
    if i == j
      # probability is 1 => log(1) = 0
      0
    else
      # probability is 0 => by convention we set log(0) = -Inf
      -Inf
    end
  elseif λ ≈ μ
    # rates are equal for our floating-point machine
    ξ = 1 / λ
    θ = λ * t
    ω = log(θ / (1 + θ))

    if j > 0
      if t ≈ ξ
        # in this case we have a simple closed form solution

        # log(2) is automatically converted to a Float64 therefore we use
        # log(F(2)) for type stability
        lgamma(i + j) - lgamma(i) - lgamma(j + 1) - (i + j) * log(F(2))
      elseif t < ξ
        # addends are all positive and we can safely evaluate each value
        x = [log_addend_equal_rates_v1(h, i, j, θ, ω) for h = 0:b]

        # we have logarithms but we must sum exp(x)
        # in this case each and every exp(x[u]) is very close to zero (their sum
        # must be in [0, 1] by definition of probability)

        # use Kahan summation algorithm and other numerical tricks to reduce
        # rounding error
        sort!(x, rev=true)

        # log(exp(x[1]) + exp(x[2]) + ... + exp(x[n])) =
        # = x[1] + log(1 + exp(x[2] - x[1]) + ... + exp(x[n] - x[1]))
        x[1] + log1p(sum_kbn(exp.(x[2:end] .- x[1])))
      else
        # addends are alternating in sign and we have a huge problem
        # compute the positive and negative values separately
        x = [log_addend_equal_rates_v2(h, i, j, θ, ω) for h = 0:2:b]
        y = [log_addend_equal_rates_v2(h, i, j, θ, ω) for h = 1:2:b]

        # we sort them and try to subtract similar magnitude numbers
        sort!(x, rev=true)
        sort!(y, rev=true)

        # we cannot use exp(x) and exp(y) directly because they might have a
        # magnitude that is considered Inf by our machine
        # we normalize them by the absolute maximum
        M = max(x[1], y[1])

        # when b is odd, length(x) and length(y) are the same
        # when b is even, length(x) = length(y) + 1
        z = if (b % 2) != 0
          exp.(x .- M) .- exp.(y .- M)
        else
          [exp.(x[1:(end - 1)] .- M) .- exp.(y .- M); exp(x[end] - M)]
        end

        s = sum_kbn(z)

        # the sum must be positive because otherwise the probability would
        # be negative. If `s` is negative it means that numerical error is too
        # big
        if s > 0
          M + log(s)
        else
          -Inf
        end
      end
    else
      # j = 0
      if t ≈ ξ
        - i * log(F(2))
      else
        i * log(ω)
      end
    end
  elseif μ < ϵ
    # pure birth process
    if j >= i
      θ = λ * t

      lgamma(j) + (j - i) * logexpm1(θ) -
      (lgamma(i) + lgamma(j - i + 1) + j * θ)
    else
      # final population size cannot be lower than the initial population size
      -Inf
    end
  elseif λ < ϵ
    # pure death process
    if j <= i
      θ = μ * t

      lgamma(i + 1) + (i - j) * logexpm1(θ) -
      (lgamma(j + 1) + lgamma(i - j + 1) + i * θ)
    else
      # final population size cannot be higher than the initial population size
      -Inf
    end
  else
    ξ = log(λ / μ) / (λ - μ)

    if j > 0
      if t ≈ ξ
        # in this case we have a simple closed form solution
        lgamma(i + j) + i * log(μ) + j * log(λ) -
        (lgamma(i) + lgamma(j + 1) + (i + j) * log(λ + μ))
      elseif λ > μ
        θ = (μ - λ) * t
        ω = log(μ / λ)

        if t < ξ
          # addends are all positive and we can safely evaluate each value
          x = [log_addend_birth_rate_greater_v1(h, i, j, θ, ω) for h = 0:b]

          # we have logarithms but we must sum exp(x)
          # in this case each and every exp(x[u]) is very close to zero (their
          # sum must be in [0, 1] by definition of probability)

          # use Kahan summation algorithm and other numerical tricks to reduce
          # rounding error
          sort!(x, rev=true)

          # log(exp(x[1]) + exp(x[2]) + ... + exp(x[n])) =
          # = x[1] + log(1 + exp(x[2] - x[1]) + ... + exp(x[n] - x[1]))
          x[1] + log1p(sum_kbn(exp.(x[2:end] .- x[1])))
        else
          # addends are alternating in sign and we have a huge problem
          # compute the positive and negative values separately
          x = [log_addend_birth_rate_greater_v2(h, i, j, θ, ω) for h = 0:2:b]
          y = [log_addend_birth_rate_greater_v2(h, i, j, θ, ω) for h = 1:2:b]

          # we sort them and try to subtract similar magnitude numbers
          sort!(x, rev=true)
          sort!(y, rev=true)

          # we cannot use exp(x) and exp(y) directly because they might have a
          # magnitude that is considered Inf by our machine
          # we normalize them by the absolute maximum
          M = max(x[1], y[1])

          # when b is odd, length(x) and length(y) are the same
          # when b is even, length(x) = length(y) + 1
          z = if (b % 2) != 0
            exp.(x .- M) .- exp.(y .- M)
          else
            [exp.(x[1:(end - 1)] .- M) .- exp.(y .- M); exp(x[end] - M)]
          end

          s = sum_kbn(z)

          # the sum must be positive because otherwise the probability would
          # be negative. If `s` is negative it means that numerical error is too
          # big
          if s > 0
            M + log(s)
          else
            -Inf
          end
        end
      else
        θ = (λ - μ) * t
        ω = log(λ / μ)

        if t < ξ
          # addends are all positive and we can safely evaluate each value
          x = [log_addend_death_rate_greater_v1(h, i, j, θ, ω) for h = 0:b]

          # we have logarithms but we must sum exp(x)
          # in this case each and every exp(x[u]) is very close to zero (their
          # sum must be in [0, 1] by definition of probability)

          # use Kahan summation algorithm and other numerical tricks to reduce
          # rounding error
          sort!(x, rev=true)

          # log(exp(x[1]) + exp(x[2]) + ... + exp(x[n])) =
          # = x[1] + log(1 + exp(x[2] - x[1]) + ... + exp(x[n] - x[1]))
          x[1] + log1p(sum_kbn(exp.(x[2:end] .- x[1])))
        else
          # addends are alternating in sign and we have a huge problem
          # compute the positive and negative values separately
          x = [log_addend_death_rate_greater_v2(h, i, j, θ, ω) for h = 0:2:b]
          y = [log_addend_death_rate_greater_v2(h, i, j, θ, ω) for h = 1:2:b]

          # we sort them and try to subtract similar magnitude numbers
          sort!(x, rev=true)
          sort!(y, rev=true)

          # we cannot use exp(x) and exp(y) directly because they might have a
          # magnitude that is considered Inf by our machine
          # we normalize them by the absolute maximum
          M = max(x[1], y[1])

          # when b is odd, length(x) and length(y) are the same
          # when b is even, length(x) = length(y) + 1
          z = if (b % 2) != 0
            exp.(x .- M) .- exp.(y .- M)
          else
            [exp.(x[1:(end - 1)] .- M) .- exp.(y .- M); exp(x[end] - M)]
          end

          s = sum_kbn(z)

          # the sum must be positive because otherwise the probability would
          # be negative. If `s` is negative it means that numerical error is too
          # big
          if s > 0
            M + log(s)
          else
            -Inf
          end
        end
      end
    else
      # j = 0
      if t ≈ ξ
        i * log(μ / (λ + μ))
      else
        θ = (λ - μ) * t
        ω = log(λ / μ)
        i * log(expm1(θ) / expm1(ω + θ))
      end
    end
  end
end


text:
struct ConvexPolygon{D, N, T}
    pts::SVector{N, SVector{D, T}}
end

function ConvexPolygon(pts::Array{Array{T, 1}}, ::Val{D}, ::Val{N}) where {D, N, T}
    ConvexPolygon(SVector{N, SVector{D, T}}(pts...))
end

function support(poly::ConvexPolygon{D, N, T}, dir::SVector{D, T}) where {D, N, T}
    @inbounds poly.pts[argmax(Ref(dir').*poly.pts)]
end

function point(pt::SVector{D, T}) where {D, T}
    ConvexPolygon(SVector{1, SVector{D, T}}(pt))
end

macro point(ex)
    :(point(@SVector $ex))
end

function line(ptA::SVector{D, T}, ptB::SVector{D, T}) where {D, T}
    ConvexPolygon(SVector{2, SVector{D, T}}(ptA, ptB))
end

macro line(ex)
    :(line(@SVector($(ex.args[1])), @SVector($(ex.args[2]))))
end

function rect(center::SVector{2, T}, widths::SVector{2, T}) where {T}
    pt1 = center + SVector{2}( 0.5,  0.5).*widths
    pt2 = center + SVector{2}(-0.5,  0.5).*widths
    pt3 = center + SVector{2}(-0.5, -0.5).*widths
    pt4 = center + SVector{2}( 0.5, -0.5).*widths
    ConvexPolygon(SVector{4, SVector{2, T}}(pt1, pt2, pt3, pt4))
end

function rect(center::SVector{3, T}, widths::SVector{3, T}) where {T}
    pt1 = center + SVector{3}( 0.5,  0.5,  0.5).*widths
    pt2 = center + SVector{3}(-0.5,  0.5,  0.5).*widths
    pt3 = center + SVector{3}(-0.5, -0.5,  0.5).*widths
    pt4 = center + SVector{3}( 0.5, -0.5,  0.5).*widths
    pt5 = center + SVector{3}( 0.5,  0.5, -0.5).*widths
    pt6 = center + SVector{3}(-0.5,  0.5, -0.5).*widths
    pt7 = center + SVector{3}(-0.5, -0.5, -0.5).*widths
    pt8 = center + SVector{3}( 0.5, -0.5, -0.5).*widths
    ConvexPolygon(SVector{8, SVector{3, T}}(pt1, pt2, pt3, pt4,
                                            pt5, pt6, pt7, pt8))
end

macro rect(ex)
    :(rect(@SVector($(ex.args[1])), @SVector($(ex.args[2]))))
end

function square(center::SVector{D, T}, width::Real) where {D, T}
    widths = width*@SVector(ones(D))
    rect(center, widths)
end

macro square(ex)
    :(square(@SVector($(ex.args[1])), $(ex.args[2])))
end

function randpoly(center=[0., 0.], θ=0.0; n::Int=10, scale=1.0)
    # from http://cglab.ca/~sander/misc/ConvexGeneration/convex.html
    xs = sort(rand(n))
    ys = sort(rand(n))
    minx = xs[1]; maxx = xs[end]
    miny = ys[1]; maxy = ys[end]
    xv = similar(xs)
    yv = similar(ys)
    ltop = minx; lbot = minx;
    llef = miny; lrig = miny;
    for i = 2:n-1
        if rand(Bool)
            xv[i-1] = xs[i] - ltop
            ltop = xs[i]
        else
            xv[i-1] = lbot - xs[i]
            lbot = xs[i]
        end
        if rand(Bool)
            yv[i-1] = ys[i] - llef
            llef = ys[i]
        else
            yv[i-1] = lrig - ys[i]
            lrig = ys[i]
        end
    end
    xv[n-1] = maxx - ltop
    xv[n] = lbot - maxx
    yv[n-1] = maxy - llef
    yv[n] = lrig - maxy
    shuffle!(xv)
    shuffle!(yv)
    vec = [[xv[i], yv[i]] for i = 1:n]
    sort!(vec, by = x -> atan(x[2], x[1]))
    v = [0.0, 0.0]
    minv = [0.0, 0.0]
    pts = similar(vec)
    for i = 1:n
        pts[i] = v
        v += vec[i]
        minv = min.(minv, v)
    end
    shift = [minx-0.5, miny-0.5] - minv
    pts .+= Ref(shift)
    rotate = [cos(θ) -sin(θ); sin(θ) cos(θ)]
    pts = Ref(rotate).*pts
    pts *= scale
    pts .+= Ref(center)
    ConvexPolygon(pts, Val(2), Val(n))
end


text:
struct PoincareInvariant2ndTrapezoidal{DT,TT,ET,ΩT} <: AbstractPoincareInvariant2nd{DT}
    equ::ET
    ω::ΩT
    Δt::TT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{Double64,1,Vector{Double64}}
    J::OffsetArray{Double64,1,Vector{Double64}}
    ΔI::OffsetArray{Double64,1,Vector{Double64}}
    ΔJ::OffsetArray{Double64,1,Vector{Double64}}
end

function PoincareInvariant2ndTrapezoidal(f_equ::Function, f_surface::Function, ω::ΩT, Δt::TT, d::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT}

    if get_config(:verbosity) > 1
        println()
        println("Second Euler-Poincaré Integral Invariant (Trapezoidal)")
        println("======================================================")
        println()
        println(" nx    = ", nx)
        println(" ny    = ", ny)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    q₀ = reshape([f_surface(i/nx, j/ny) for i in 1:nx, j in 1:ny], nx*ny)

    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    J  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔI = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔJ = OffsetArray(zeros(Double64, nt+1), 0:nt)

    PoincareInvariant2ndTrapezoidal{DT,TT,typeof(equ),ΩT}(equ, ω, Δt, nx, ny, ntime, nsave, nt, I, J, ΔI, ΔJ)
end

function interpolate_trajectory(x, i1, j1, λ, μ, γ, nx, ny)
    @assert length(γ) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in eachindex(γ)
        γ[k] = x[k, nx*(j1-1)+i1] * (1-λ) * (1-μ) +
               x[k, nx*(j1-1)+i2] *    λ  * (1-μ) +
               x[k, nx*(j2-1)+i1] * (1-λ) *    μ  +
               x[k, nx*(j2-1)+i2] *    λ  *    μ
    end
end

function interpolate_derivative_i(x, i1, j1, λ, μ, γ̇, nx, ny)
    @assert length(γ̇) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in eachindex(γ̇)
        γ̇[k] = (x[k, nx*(j1-1)+i2] - x[k, nx*(j1-1)+i1]) * (1-μ) +
               (x[k, nx*(j2-1)+i2] - x[k, nx*(j2-1)+i1]) *    μ
    end
end

function interpolate_derivative_j(x, i1, j1, λ, μ, γ̇, nx, ny)
    @assert length(γ̇) == size(x, 1)
    @assert λ ≥ 0.
    @assert λ ≤ 1.
    @assert μ ≥ 0.
    @assert μ ≤ 1.

    @assert i1 > 0
    @assert i1 < nx
    @assert j1 > 0
    @assert j1 < ny

    i2 = i1 + 1
    j2 = j1 + 1

    for k in eachindex(γ̇)
        γ̇[k] = (x[k, nx*(j2-1)+i1] - x[k, nx*(j1-1)+i1]) * (1-λ) +
               (x[k, nx*(j2-1)+i2] - x[k, nx*(j1-1)+i2]) *    λ
    end
end

function integrate(t, γ, γ̇ᵢ, γ̇ⱼ, ω, b::Vector{TT}, c::Vector{TT}, q::Vector{DT}, vᵢ::Vector{DT}, vⱼ::Vector{DT}, B::Matrix{DT}) where {DT,TT}
    @assert length(b) == length(c)

    local result = zero(Double64)

    for i in eachindex(b)
        for j in eachindex(b)
            γ(c[i], c[j], q)
            γ̇ᵢ(c[i], c[j], vᵢ)
            γ̇ⱼ(c[i], c[j], vⱼ)
            ω(t, q, B)

            result += b[i] * b[j] * vector_matrix_vector_product(vⱼ, B, vᵢ)
        end
    end

    return result
end

function surface_integral(t, x::AbstractMatrix{DT}, ω, nx, ny) where {DT}
    local b = [0.5, 0.5]
    local c = [0.0, 1.0]

    local q  = zeros(DT, size(x,1))
    local vᵢ = zeros(DT, size(x,1))
    local vⱼ = zeros(DT, size(x,1))
    local B  = zeros(DT, size(x,1), size(x,1))
    local I  = zero(Double64)

    integrate_trapezoidal = (γ, γ̇ᵢ, γ̇ⱼ) -> integrate(t, γ, γ̇ᵢ, γ̇ⱼ, ω, b, c, q, vᵢ, vⱼ, B)

    for j in 1:ny-1
        for i in 1:nx-1
            γ  = (λ, μ, y) -> interpolate_trajectory(x, i, j, λ, μ, y, nx, ny)
            γ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(x, i, j, λ, μ, y, nx, ny)
            γ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(x, i, j, λ, μ, y, nx, ny)
            I += integrate_trapezoidal(γ, γ̇ᵢ, γ̇ⱼ)
        end
    end

    return I
end

function integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b::Vector{TT}, c::Vector{TT}, vᵢ::Vector{DT}, vⱼ::Vector{DT}) where {DT,TT}
    @assert length(b) == length(c)

    local result = zero(Double64)

    for i in eachindex(b)
        for j in eachindex(b)
            Θ̇ᵢ(c[i], c[j], vᵢ)
            γ̇ⱼ(c[i], c[j], vⱼ)
            result += b[i] * b[j] * vᵢ ⋅ vⱼ

            γ̇ᵢ(c[i], c[j], vᵢ)
            Θ̇ⱼ(c[i], c[j], vⱼ)
            result -= b[i] * b[j] * vᵢ ⋅ vⱼ
        end
    end

    return result
end

function surface_integral_canonical(q::AbstractMatrix{DT}, p::AbstractMatrix{DT}, nx, ny) where {DT}
    local b = [0.5, 0.5]
    local c = [0.0, 1.0]

    local vᵢ = zeros(DT, size(q,1))
    local vⱼ = zeros(DT, size(q,1))
    local I  = zero(Double64)

    integrate_trapezoidal = (γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ) -> integrate_canonical(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ, b, c, vᵢ, vⱼ)

    for j in 1:ny-1
        for i in 1:nx-1
            γ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(q, i, j, λ, μ, y, nx, ny)
            γ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(q, i, j, λ, μ, y, nx, ny)
            Θ̇ᵢ = (λ, μ, y) -> interpolate_derivative_i(p, i, j, λ, μ, y, nx, ny)
            Θ̇ⱼ = (λ, μ, y) -> interpolate_derivative_j(p, i, j, λ, μ, y, nx, ny)
            I += integrate_trapezoidal(γ̇ᵢ, γ̇ⱼ, Θ̇ᵢ, Θ̇ⱼ)
        end
    end

    return I
end

function evaluate_poincare_invariant(pinv::PoincareInvariant2ndTrapezoidal{DT}, sol::Solution) where {DT}
    local verbosity = get_config(:verbosity)

    for i in axes(sol.q,1)
        verbosity > 1 ? println("      it = ", i) : nothing
        @views pinv.I[i]  = surface_integral(sol.t[i], hcat(sol.q[i,:]...), pinv.ω, pinv.nx, pinv.ny)
        pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
        verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing

        if hasproperty(sol, :p)
            @views pinv.J[i] = surface_integral_canonical(hcat(sol.q[i,:]...), hcat(sol.p[i,:]...), pinv.nx, pinv.ny)
            pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
            verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
        end
    end

    return (DT.(pinv.I), DT.(pinv.J), DT.(pinv.ΔI), DT.(pinv.ΔJ))
end

function write_to_hdf5(pinv::PoincareInvariant2ndTrapezoidal, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing

    end
end


text:
import Term.Renderables: Renderable
using Term.Layout

import Term:
    RenderableText,
    Spacer,
    vLine,
    hLine,
    cleantext,
    textlen,
    chars,
    Panel,
    vstack,
    center!,
    leftalign!,
    rightalign!,
    leftalign,
    center,
    rightalign,
    lvstack,
    cvstack,
    rvstack

@testset "Layout - pad" begin
    @test pad("aaa", 20, :left) == "aaa                 "
    @test pad("aaa", 20, :right) == "                 aaa"
    @test pad("aaa", 20, :center) == "        aaa         "
    @test pad("aaa", 10, 20) == "          aaa                    "

    p = Panel(; width = 20, height = 10)
    pad!(p, 20, 20)
    @test p isa Panel
    @test size(p.measure) == (60, 10)

    p = Panel(; width = 20, height = 10)
    pad!(p; width = 30)
    @test p isa Panel
    @test size(p.measure) == (30, 10)
end

@testset "Layout - vertical pad" begin
    @test vertical_pad("ab", 5, :top) == "ab\n  \n  \n  \n  "
    @test vertical_pad("ab", 5, :bottom) == "  \n  \n  \n  \nab"
    @test vertical_pad("ab", 5, :center) == "  \n  \nab\n  \n  "

    @test vertical_pad("ab", 5, 5) == "  \n  \n  \n  \n  \nab\n  \n  \n  \n  \n  "

    p = Panel(; width = 20, height = 10)
    @test string(vertical_pad(p, 4, 4)) ==
          "                    \n                    \n                    \n                    \n\e[22m╭──────────────────╮\e[22m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[0m\e[22m│\e[22m\e[0m                  \e[0m\e[22m│\e[22m\e[0m\n\e[22m╰──────────────────╯\e[22m\e[0m\n                    \n                    \n                    \n                    "

    vertical_pad!(p, 4, 4)
    @test p isa Panel
    @test size(p.measure) == (20, 18)

    p = Panel(; width = 20, height = 10)
    vertical_pad!(p; height = 30)
    @test p isa Panel
    @test size(p.measure) == (20, 30)
end

@testset "\e[34mlayout - spacer" begin
    sizes = [(22, 1), (44, 123), (21, 1), (4334, 232)]
    for (w, h) in sizes
        spacer = Spacer(w, h)
        @test size(spacer.measure) == (w, h)
    end
end

@testset "\e[34mlayout - justification" begin
    make_panels() = (Panel(; width = 5), Panel(; width = 10), Panel(; width = 15))

    # right justify
    p1, p2, p3 = make_panels()

    rightalign!(p1, p2, p3)
    @test p1 isa Panel
    @test p1.measure.w == p2.measure.w == p3.measure.w

    # right justify
    p1, p2, p3 = make_panels()

    center!(p1, p2, p3)
    @test p1 isa Panel
    @test p1.measure.w == p2.measure.w == p3.measure.w

    # right justify
    p1, p2, p3 = make_panels()
    leftalign!(p1, p2, p3)
    @test p1 isa Panel
    @test p1.measure.w == p2.measure.w == p3.measure.w

    # convenience functions
    p1, p2, p3 = make_panels()
    pp = lvstack(p1, p2, p3)
    @test pp isa Renderable
    @test pp.measure.w == 15
    @test string(pp) ==
          "\e[22m╭───╮\e[22m          \n\e[22m╰───╯\e[22m\e[0m          \n\e[22m╭────────╮\e[22m     \n\e[22m╰────────╯\e[22m\e[0m     \n\e[22m╭─────────────╮\e[22m\n\e[22m╰─────────────╯\e[22m\e[0m"
    @test p1 isa Panel
    @test p1.measure.w == 5

    p1, p2, p3 = make_panels()
    pp = cvstack(p1, p2, p3)
    @test pp isa Renderable
    @test pp.measure.w == 15
    @test string(pp) ==
          "     \e[22m╭───╮\e[22m     \n     \e[22m╰───╯\e[22m\e[0m     \n  \e[22m╭────────╮\e[22m   \n  \e[22m╰────────╯\e[22m\e[0m   \n\e[22m╭─────────────╮\e[22m\n\e[22m╰─────────────╯\e[22m\e[0m"
    @test p1 isa Panel
    @test p1.measure.w == 5

    p1, p2, p3 = make_panels()
    pp = rvstack(p1, p2, p3)
    @test pp isa Renderable
    @test pp.measure.w == 15
    @test string(pp) ==
          "          \e[22m╭───╮\e[22m\n          \e[22m╰───╯\e[22m\e[0m\n     \e[22m╭────────╮\e[22m\n     \e[22m╰────────╯\e[22m\e[0m\n\e[22m╭─────────────╮\e[22m\n\e[22m╰─────────────╯\e[22m\e[0m"
    @test p1 isa Panel
    @test p1.measure.w == 5
end

@testset "\e[34mlayout - vLine " begin
    for h in [1, 22, 55, 11]
        line = vLine(h)
        @test length(line.segments) == h
        @test line.measure.h == h
    end

    lines = [(22, "bold"), (55, "red on_green")]
    for (h, style) in lines
        line = vLine(h; style = style)
        @test length(line.segments) == h
        @test line.measure.h == h
    end
    line = vLine(5; style = "red")
    @test line.segments[1].text == "\e[31m│\e[39m\e[0m"

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test vLine(22; box = box).measure.h == 22
    end

    panel = Panel(; width = 20, height = 5)
    @test length(vLine(panel).segments) == 5
    @test vLine().measure.h == displaysize(stdout)[1]
end

@testset "\e[34mlayout - hLine " begin
    for w in [1, 342, 433, 11, 22]
        line = hLine(w)
        @test length(line.segments) == 1
        @test textlen(line.segments[1].text) == w
        @test line.measure.w == w
    end

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test hLine(22; box = box).measure.w == 22
        @test hLine(22, "title"; box = box).measure.w == 22
    end

    for style in ("bold", "red on_green", "blue")
        @test textlen(hLine(11; style = style).segments[1].text) == 11
        @test textlen(hLine(11, "ttl"; style = style).segments[1].text) == 11
    end

    panel = Panel(; width = 20, height = 5)
    @test hLine().measure.w == displaysize(stdout)[2]
    @test textlen(hLine(panel).segments[1].text) == 20
end

@testset "\e[34mlayout - stack strings" begin
    s1 = "."^50
    s2 = ".\n"^5 * "."
    @test s1 / s2 isa String
end

function testlayout(p, w, h)
    _p = string(p)
    widths = textwidth.(cleantext.(split(_p, "\n")))
    @test length(unique(widths)) == 1

    @test p.measure.w == w
    @test textlen(cleantext(p.segments[1].text)) == w
    @test length(chars(cleantext(p.segments[1].text))) == w

    @test p.measure.h == h
    @test length(p.segments) == h
end

@testset "\e[34mlayout - renderable" begin
    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 / r2
    @test size(r.measure) == (50, 6)

    h1 = hLine(22)
    h2 = hLine(33)
    @test size(Measure(h1 / h2)) == (33, 2)

    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 * r2
    @test size(r.measure) == (75, 4)

    # stack other renderables
    h1 = vLine(22)
    h2 = vLine(33)
    @test size(Measure(h1 * h2)) == (2, 33)
end

@testset "\e[34mlayout - panels" begin
    p1 = Panel()
    p2 = Panel(; width = 24, height = 3)
    p3 = Panel("this {red}panel{/red}"^5; width = 30, fit = false)

    testlayout(p1 * p2, 112, 3)
    @test string(p1 * p2) ==
          "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\e[22m╭──────────────────────╮\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m\e[0m\e[0m\e[22m│\e[22m\e[0m                      \e[0m\e[22m│\e[22m\e[0m\n                                                                                \e[22m╰──────────────────────╯\e[22m\e[0m"

    testlayout(p1 / p2, 88, 5)
    @test string(p1 / p2) ==
          "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m\e[0m\n\e[22m╭──────────────────────╮\e[22m                                                                \n\e[0m\e[22m│\e[22m\e[0m                      \e[0m\e[22m│\e[22m\e[0m                                                                \n\e[22m╰──────────────────────╯\e[22m\e[0m                                                                "

    testlayout(p2 * p1, 112, 3)
    testlayout(p2 / p1, 88, 5)

    testlayout(p1 * p2 * p3, 142, 5)
    testlayout(p1 / p2 / p3, 88, 10)
    testlayout(p3 * p1 * p2, 142, 5)
    testlayout(p3 / p1 / p2, 88, 10)
end

@testset "\e[34mlayout - placeholder" begin
    ph = PlaceHolder(4, 2)
    @test length(ph.segments) == 2
    @test size(ph.measure) == (4, 2)
    @test string(ph) == "\e[2m╲ ╲ \e[22m\n\e[2m ╲ ╲\e[22m"

    @test string(PlaceHolder(25, 12)) ==
          "\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ ╲ \e[22m"

    p = Panel(; width = 8, height = 4)
    ph = PlaceHolder(p)
    @test size(ph.measure) == size(p.measure)
end

@testset "Layout - h/vstack with pad" begin
    p1, p2, p3 = PlaceHolder(5, 5), PlaceHolder(5, 5), PlaceHolder(3, 6)

    @test string(
        hstack(PlaceHolder(5, 5), PlaceHolder(5, 5), PlaceHolder(3, 6); pad = 10),
    ) ==
          "\e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m          \e[2m ╲ ╲ \e[22m          \e[2m ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m          \e[2m ╲ ╲ \e[22m          \e[2m ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲ ╲\e[22m          \e[2m╲ ╲\e[22m\n                              \e[2m ╲ \e[22m"
    @test string(
        hstack(PlaceHolder(5, 5), PlaceHolder(5, 5), PlaceHolder(3, 6); pad = 0),
    ) ==
          "\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\e[2m ╲ ╲ \e[22m\e[2m ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\e[2m ╲ ╲ \e[22m\e[2m ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲ ╲\e[22m\e[2m╲ ╲\e[22m\n          \e[2m ╲ \e[22m"

    @test string(
        vstack(PlaceHolder(5, 5), PlaceHolder(5, 5), PlaceHolder(3, 6); pad = 0),
    ) ==
          "\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  \n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  \n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  "
    @test string(
        vstack(PlaceHolder(5, 5), PlaceHolder(5, 5), PlaceHolder(3, 6); pad = 3),
    ) ==
          "\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n     \n     \n     \n     \n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n\e[2m ╲ ╲ \e[22m\n\e[2m╲ ╲ ╲\e[22m\n     \n     \n     \n     \n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  \n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  \n\e[2m╲ ╲\e[22m  \n\e[2m ╲ \e[22m  "
end

@testset "Layout - GRID" begin
    @test size(grid(layout = (3, 3), pad = (0, 0)).measure) == (48, 27)

    @test size(grid(layout = (3, 3)).measure) == (48, 27)

    @test size(grid(layout = (3, 3), pad = 2).measure) == (52, 31)

    @test size(grid(layout = (3, 3), pad = (5, 1)).measure) == (58, 29)

    @test size(grid(layout = (3, 3), pad = (5, 3)).measure) == (58, 33)

    # test passing renderables
    rens = repeat([PlaceHolder(10, 5)], 9)

    @test size(grid(rens; aspect = 1).measure) == (30, 15)

    @test size(grid(rens).measure) == (50, 10)

    @test size(grid(rens; pad = (2, 1)).measure) == (58, 11)

    @test size(grid(rens; pad = (2, 1), aspect = 0.5).measure) == (22, 29)

    @test size(grid(rens; pad = (2, 1), aspect = 1.5).measure) == (46, 17)

    g = grid(
        rens;
        pad = (2, 1),
        aspect = (12, 12),
        placeholder = PlaceHolder(10, 5; style = "red"),
    )
    @test size(g.measure) == (34, 17)

    w, h = 20, 10
    panels = collect(
        Panel("{on_$c} {/on_$c}", width = w, height = h) for c in (
            :bright_red,
            :bright_green,
            :bright_blue,
            :bright_yellow,
            :bright_magenta,
            :bright_cyan,
            :bright_black,
            :bright_white,
        )
    )

    # auto layout (default placeholder)
    for i in 2:length(panels)
        g = grid(panels[1:i])
        nc, nr = if i ≤ 3
            (i, 1)
        elseif i ≤ 5
            (3, 2)
        else
            (4, 2)
        end
        @test size(g.measure) == (w * nc, h * nr)
    end

    # matrix, explicit
    grid(reshape(panels[1:4], 2, 2))

    # vector, half explicit
    grid(panels, layout = (nothing, 4))
    grid(panels, layout = (2, nothing))

    # vector, explicit
    grid(panels; layout = (2, 4))
end


text:
module PolynomialsWrapper

using JavaCall

const Class = @jimport java.lang.Class
const DerivativeStructure = @jimport org.hipparchus.analysis.differentiation.DerivativeStructure
const Field = @jimport org.hipparchus.Field
const FieldPolynomialFunction = @jimport org.hipparchus.analysis.polynomials.FieldPolynomialFunction
const FieldPolynomialSplineFunction = @jimport org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction
const FieldUnivariateFunction = @jimport org.hipparchus.analysis.FieldUnivariateFunction
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const PolynomialFunction = @jimport org.hipparchus.analysis.polynomials.PolynomialFunction
const PolynomialFunctionLagrangeForm = @jimport org.hipparchus.analysis.polynomials.PolynomialFunctionLagrangeForm
const PolynomialFunctionNewtonForm = @jimport org.hipparchus.analysis.polynomials.PolynomialFunctionNewtonForm
const PolynomialSplineFunction = @jimport org.hipparchus.analysis.polynomials.PolynomialSplineFunction
const PolynomialsUtils = @jimport org.hipparchus.analysis.polynomials.PolynomialsUtils
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RealFieldUnivariateFunction = @jimport org.hipparchus.analysis.RealFieldUnivariateFunction

export FieldPolynomialFunction
export FieldPolynomialSplineFunction
export PolynomialFunction
export PolynomialFunctionLagrangeForm
export PolynomialFunctionNewtonForm
export PolynomialSplineFunction
export PolynomialsUtils
export add
export anti_derivative
export create_chebyshev_polynomial
export create_hermite_polynomial
export create_jacobi_polynomial
export create_laguerre_polynomial
export create_legendre_polynomial
export degree
export equals
export evaluate
export get_centers
export get_class
export get_coefficients
export get_field
export get_interpolating_points
export get_interpolating_values
export get_knots
export get_n
export get_newton_coefficients
export get_polynomials
export hash_code
export integrate
export is_valid_point
export multiply
export negate
export notify
export notify_all
export polynomial_derivative
export polynomial_spline_derivative
export shift
export subtract
export to_real_field_univariate_function
export to_string
export value
export verify_interpolation_array
export wait

include("field_polynomial_function.jl")
include("field_polynomial_spline_function.jl")
include("polynomial_function.jl")
include("polynomial_function_lagrange_form.jl")
include("polynomial_function_newton_form.jl")
include("polynomial_spline_function.jl")
include("polynomials_utils.jl")

end


text:
## Add number of workers of one's choosing:
#addprocs(20)

using ModelConstructors, DSGE, SMC
@everywhere using ModelConstructors, DSGE, SMC
@everywhere import ModelConstructors: ParameterVector
@everywhere import SMC: smc

## Initialize model object
m = AnSchorfheide()
m <= Setting(:λ, 2.0)
m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:resampler_smc, :systematic)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:use_fixed_schedule, true)

data = df_to_matrix(m, load_data(m))

## Specify settings
verbose = :low
use_chand_recursion = true
filestring_addl     = Vector{String}()

## Specify output paths
loadpath            = rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl)
savepath            = rawpath(m, "estimate", "smc_cloud.jld2", filestring_addl)
particle_store_path = rawpath(m, "estimate", "smcsave.h5",     filestring_addl)

## Define likelihood function in correct form, on all workers
DSGE.sendto(workers(), m = m)
function my_loglikelihood(parameters::ParameterVector, data::Matrix{Float64})::Float64
    DSGE.update!(m, parameters)
    DSGE.likelihood(m, data; sampler = false, catch_errors = true,
                    use_chand_recursion = true, verbose = :low)
end
@everywhere function my_loglikelihood(parameters::ParameterVector, data::Matrix{Float64})::Float64
    DSGE.update!(m, parameters)
    DSGE.likelihood(m, data; sampler = false, catch_errors = true,
                    use_chand_recursion = true, verbose = :low)
end

## Run SMC
smc(my_loglikelihood, m.parameters, data; verbose = verbose, data_vintage = data_vintage(m),
    n_parts = 400, n_Φ = 100, λ = 2.0, parallel = true, save_intermediate = false,
    particle_store_path = particle_store_path, loadpath = loadpath,
    savepath = savepath)


text:
# ---------------------------------------------------------------------------
# This is experimental, and currently unused
# ---------------------------------------------------------------------------

using DataFlow

export
	edges

function nodename(v)
    val = v.value
    isconst = isa(val, DataFlow.Constant)
    val = isconst ? val.value : val
    val, isconst
end

function edges(v::DVertex, source = Int[], destiny = Int[], names = Symbol[], cache = DataFlow.ODict())
    name, isconst = nodename(v)
    haskey(cache, isconst ? name : v) && return source,destiny,names
    push!(names, name)
    didx = length(cache)+1
    cache[isconst ? name : v] = didx
	for v′ in DataFlow.inputs(v)
		edges(v′, source, destiny, names, cache)
        sname, sisconst = nodename(v′)
		push!(source, cache[sisconst ? sname : v′])
		push!(destiny, didx)
	end
	source, destiny, names
end
edges(g::DataFlow.SyntaxGraph, args...) = edges(g.output, args...)

using PlotRecipes
@recipe function f(g::Union{DataFlow.SyntaxGraph,DVertex})
	source, destiny, names = edges(g)
    names --> map(string, names)
    method --> :tree
    root --> :left
    PlotRecipes.GraphPlot((source, destiny))
end


text:
foo() = rand()             # repeated calls change the result with each call
repeat([foo()], outer=5)   # but this only calls foo() once, clones that first value

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export isapproxzero

Base.iszero(v::AbstractVariable) = false
Base.iszero(m::AbstractMonomial) = false
Base.iszero(t::AbstractTerm) = iszero(coefficient(t))
Base.iszero(t::AbstractPolynomial) = iszero(nterms(t))

# See https://github.com/blegat/MultivariatePolynomials.jl/issues/22
# avoids the call to be transfered to eqconstant
Base.:(==)(α::Nothing, x::APL) = false
Base.:(==)(x::APL, α::Nothing) = false
Base.:(==)(α::Dict, x::APL) = false
Base.:(==)(x::APL, α::Dict) = false
Base.:(==)(α::Nothing, x::RationalPoly) = false
Base.:(==)(x::RationalPoly, α::Nothing) = false
Base.:(==)(α::Dict, x::RationalPoly) = false
Base.:(==)(x::RationalPoly, α::Dict) = false

function polyeqterm(p::AbstractPolynomial, t)
    if iszero(p)
        iszero(t)
    else
        # terms/nterms ignore zero terms
        nterms(p) == 1 && leadingterm(p) == t
    end
end
polyeqterm(p::APL, t) = polyeqterm(polynomial(p), t)

eqconstant(α, v::AbstractVariable) = false
eqconstant(v::AbstractVariable, α) = false
function _termeqconstant(t::AbstractTermLike, α)
    if iszero(t)
        iszero(α)
    else
        α == coefficient(t) && isconstant(t)
    end
end
eqconstant(α, t::AbstractTermLike) = _termeqconstant(t, α)
eqconstant(t::AbstractTermLike, α) = _termeqconstant(t, α)
eqconstant(α, p::APL) = polyeqterm(p, α)
eqconstant(p::APL, α) = polyeqterm(p, α)

function Base.:(==)(mono::AbstractMonomial, v::AbstractVariable)
    return isone(degree(mono)) && variable(mono) == v
end
function Base.:(==)(v::AbstractVariable, mono::AbstractMonomial)
    return isone(degree(mono)) && v == variable(mono)
end
function Base.:(==)(t::AbstractTerm, mono::AbstractMonomialLike)
    return isone(coefficient(t)) && monomial(t) == mono
end
function Base.:(==)(mono::AbstractMonomialLike, t::AbstractTerm)
    return isone(coefficient(t)) && mono == monomial(t)
end

function Base.:(==)(t1::AbstractTerm, t2::AbstractTerm)
    c1 = coefficient(t1)
    c2 = coefficient(t2)
    if iszero(c1)
        iszero(c2)
    else
        c1 == c2 && monomial(t1) == monomial(t2)
    end
end
Base.:(==)(p::AbstractPolynomial, t::AbstractTerm) = polyeqterm(p, t)
Base.:(==)(t::AbstractTerm, p::AbstractPolynomial) = polyeqterm(p, t)

function compare_terms(p1::AbstractPolynomial, p2::AbstractPolynomial, isz, op)
    i1 = 1
    i2 = 1
    t1 = terms(p1)
    t2 = terms(p2)
    while true
        while i1 <= length(t1) && isz(coefficient(t1[i1]))
            i1 += 1
        end
        while i2 <= length(t2) && isz(coefficient(t2[i2]))
            i2 += 1
        end
        if i1 > length(t1) && i2 > length(t2)
            return true
        end
        if i1 > length(t1) || i2 > length(t2)
            return false
        end
        if !op(t1[i1], t2[i2])
            return false
        end
        i1 += 1
        i2 += 1
    end
end

# Can there be zero term in TypedPolynomials ?
#function (==)(p1::AbstractPolynomial, p2::AbstractPolynomial)
#    nterms(p1) != nterms(p2) && return false
#    for (t1, t2) in zip(terms(p1), terms(p2))
#        @assert !iszero(t1) && !iszero(t2) # There should be no zero term
#        if t1 != t2
#            return false
#        end
#    end
#    return true
#end
Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) = compare_terms(p1, p2, iszero, ==)

Base.:(==)(p::RationalPoly, q::RationalPoly) = p.num*q.den == q.num*p.den
# Solve ambiguity with (::PolyType, ::Any)
Base.:(==)(p::APL, q::RationalPoly) = p*q.den == q.num
Base.:(==)(q::RationalPoly, p::APL) = p == q
Base.:(==)(α, q::RationalPoly) = α*q.den == q.num
Base.:(==)(q::RationalPoly, α) = α == q

# α could be a JuMP affine expression
isapproxzero(α; ztol::Real=0.) = false
function isapproxzero(α::Number; ztol::Real=Base.rtoldefault(α, α, 0))
    abs(α) <= ztol
end

isapproxzero(m::AbstractMonomialLike; kwargs...) = false
isapproxzero(t::AbstractTermLike; kwargs...) = isapproxzero(coefficient(t); kwargs...)
isapproxzero(p::APL; kwargs...) = all(isapproxzero.(terms(p); kwargs...))
isapproxzero(p::RationalPoly; kwargs...) = isapproxzero(p.num; kwargs...)

Base.isapprox(t1::AbstractTerm, t2::AbstractTerm; kwargs...) = isapprox(coefficient(t1), coefficient(t2); kwargs...) && monomial(t1) == monomial(t2)
function Base.isapprox(p1::AbstractPolynomial{S}, p2::AbstractPolynomial{T};
                       atol=0, ztol::Real=iszero(atol) ? Base.rtoldefault(S, T, 0) : atol, kwargs...) where {S, T}
    return compare_terms(p1, p2, t -> isapproxzero(t; ztol=ztol),
                         (x, y) -> isapprox(x, y; atol=atol, kwargs...))
end

Base.isapprox(p::RationalPoly, q::RationalPoly; kwargs...) = isapprox(p.num*q.den, q.num*p.den; kwargs...)
Base.isapprox(p::RationalPoly, q::APL; kwargs...) = isapprox(p.num, q*p.den; kwargs...)
Base.isapprox(p::APL, q::RationalPoly; kwargs...) = isapprox(p*q.den, q.num; kwargs...)
Base.isapprox(q::RationalPoly{C}, α; kwargs...) where {C} = isapprox(q, constantterm(α, q.den); kwargs...)
Base.isapprox(α, q::RationalPoly{C}; kwargs...) where {C} = isapprox(constantterm(α, q.den), q; kwargs...)