text:
using Plots
using Oceananigans

using Statistics

xnodes(ϕ) = repeat(reshape(ϕ.grid.xC, ϕ.grid.Nx, 1, 1), 1, ϕ.grid.Ny, ϕ.grid.Nz)
ynodes(ϕ) = repeat(reshape(ϕ.grid.yC, 1, ϕ.grid.Ny, 1), ϕ.grid.Nx, 1, ϕ.grid.Nz)
znodes(ϕ) = repeat(reshape(ϕ.grid.zC, 1, 1, ϕ.grid.Nz), ϕ.grid.Nx, ϕ.grid.Ny, 1)

xnodes(ϕ::FaceFieldX) = repeat(reshape(ϕ.grid.xF[1:end-1], ϕ.grid.Nx, 1, 1), 1, ϕ.grid.Ny, ϕ.grid.Nz)
ynodes(ϕ::FaceFieldY) = repeat(reshape(ϕ.grid.yF[1:end-1], 1, ϕ.grid.Ny, 1), ϕ.grid.Nx, 1, ϕ.grid.Nz)
znodes(ϕ::FaceFieldZ) = repeat(reshape(ϕ.grid.zF[1:end-1], 1, 1, ϕ.grid.Nz), ϕ.grid.Nx, ϕ.grid.Ny, 1)

zerofunk(args...) = 0

function set_ic!(model; u=zerofunk, v=zerofunk, w=zerofunk, T=zerofunk, S=zerofunk)
    model.velocities.u.data .= u.(xnodes(model.velocities.u), ynodes(model.velocities.u), znodes(model.velocities.u))
    model.velocities.v.data .= v.(xnodes(model.velocities.v), ynodes(model.velocities.v), znodes(model.velocities.v))
    model.velocities.w.data .= w.(xnodes(model.velocities.w), ynodes(model.velocities.w), znodes(model.velocities.w))
    model.tracers.T.data    .= T.(xnodes(model.tracers.T),    ynodes(model.tracers.T),    znodes(model.tracers.T))
    model.tracers.S.data    .= S.(xnodes(model.tracers.S),    ynodes(model.tracers.S),    znodes(model.tracers.S))
    return nothing
end

plotxzslice(ϕ, slice=1, args...; kwargs...) = pcolormesh(
    view(xnodes(ϕ), :, slice, :), view(znodes(ϕ), :, slice, :), view(ϕ.data, :, slice, :), args...; kwargs...)

plotxyslice(ϕ, slice=1, args...; kwargs...) = pcolormesh(
    view(xnodes(ϕ), :, :, slice), view(ynodes(ϕ), :, :, slice), view(ϕ.data, :, :, slice), args...; kwargs...)

function total_kinetic_energy(model)
    return 0.5 * (
          sum(model.velocities.u.data.^2)
        + sum(model.velocities.v.data.^2)
        + sum(model.velocities.w.data.^2)
        )
end

function total_kinetic_energy(u, v, w)
    return 0.5 * (sum(u.data.^2) + sum(v.data.^2) + sum(w.data.^2))
end

function total_energy(model, N)
    b = model.tracers.T.data .- mean(model.tracers.T.data, dims=(1, 2))
    return 0.5 * (
          sum(model.velocities.u.data.^2)
        + sum(model.velocities.v.data.^2)
        + sum(model.velocities.w.data.^2)
        + sum(b.^2) / N^2
        )
end

function w_relative_error(model, w)
    w_ans = FaceFieldZ(w.(
        xnodes(model.velocities.w),
        ynodes(model.velocities.w),
        znodes(model.velocities.w),
        model.clock.time), model.grid)

    return mean(
        (model.velocities.w.data .- w_ans.data).^2) / mean(w_ans.data.^2)

end

function u_relative_error(model, u)
    u_ans = FaceFieldX(u.(
        xnodes(model.velocities.u),
        ynodes(model.velocities.u),
        znodes(model.velocities.u),
        model.clock.time), model.grid)

    return mean(
        (model.velocities.u.data .- u_ans.data).^2 ) / mean(u_ans.data.^2)
end

function T_relative_error(model, T)
    T_ans = CellField(T.(
        xnodes(model.tracers.T),
        ynodes(model.tracers.T),
        znodes(model.tracers.T),
        model.clock.time), model.grid)

    return mean(
        (model.tracers.T.data .- T_ans.data).^2 ) / mean(T_ans.data.^2)
end

"""
    make_vertical_slice_movie(model::Model, nc_writer::NetCDFOutputWriter,
                              var_name, Nt, Δt, var_offset=0, slice_idx=1)

Make a movie of a vertical slice produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀). A `slice_idx` can be specified to select
the index of the y-slice to be plotted (useful when plotting vertical slices
from a 3D model, it should be set to 1 for 2D xz-slice models).
"""
function make_vertical_slice_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0, slice_idx=1)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.contour(model.grid.xC, reverse(model.grid.zC), rotl90(var[:, slice_idx, :] .- var_offset),
                      fill=true, levels=9, linewidth=0, color=:balance,
                      clims=(-0.011, 0.011), title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
        # Plots.heatmap(model.grid.xC, model.grid.zC, rotl90(var[:, slice_idx, :]) .- var_offset,
        #               color=:balance, clims=(-0.01, 0.01), title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

"""
    make_horizontal_slice_movie(model::Model, nc_writer::NetCDFOutputWriter,
                                var_name, Nt, Δt, var_offset=0)

Make a movie of a horizontal slice produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀).
"""
function make_horizontal_slice_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.heatmap(model.grid.xC, model.grid.yC, var[:, :, 1] .- var_offset,
                      color=:balance, clims=(-0.01, 0.01),
                      title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

"""
    make_vertical_profile_movie(model::Model, nc_writer::NetCDFOutputWriter,
                                var_name, Nt, Δt, var_offset=0)

Make a movie of a vertical profile produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀).
"""
function make_vertical_profile_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.plot(var[1, 1, :] .- var_offset, model.grid.zC,
                   title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

using NetCDF

function make_avg_temperature_profile_movie()
    Nt, dt = 86400, 0.5
    freq = 3600
    N_frames = Int(Nt/freq)
    filename_prefix = "convection"
    var_offset = 273.15

    Nz, Lz = 128, 100
    dz = Lz/Nz
    zC = -dz/2:-dz:-Lz

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")

        filepath = filename_prefix * lpad(freq*n, 9, "0") * ".nc"
        field_data = ncread(filepath, "T")
        ncclose(filepath)

        T_profile = mean(field_data; dims=[1,2])

        Plots.plot(reshape(T_profile, Nz) .- var_offset, zC,
                   title="t=$(freq*n*dt) s ($(round(freq*n*dt/86400; digits=2)) days)")
    end

    mp4(animation, filename_prefix * "$(round(Int, time())).mp4", fps=30)
end


text:
module SStools

import Distributions.MvNormal

export kalman_filter, fast_state_smoother, simulate

"""
Perform Kalman filtering on the data y.
Conventions are as in Durbin and Koopman (2012).
Relevant dimensions are:
- Nt: number of time points
- Np: dimension of observation space
- Nm: dimension of state space
- Nr: dimension of state noise covariance

Parameters:
- a0: mean of prior on states: (Nm, 1) or (Nm,)
- P0: variance of prior on states: (Nm, Nm)
- Z: either a (Np, Nm, Nt) array or a (Np, Nm) array
- H: either a (Np, Np, Nt) array or a (Np, Np) array
- T: either a (Nm, Nm, Nt) array or a (Nm, Nm) array
- R: either a (Nm, Nr, Nt) array or a (Nm, Nr) array
- Q: either a (Nr, Nr, Nt) array or a (Nr, Nr) array

Returns:
- v: (Np, Nt) vector of residuals at each time
- K: (Nm, Np, Nt) Kalman gain matrix
- Finv: (Np, Np, Nt) inverse of prediction variance matrix
"""
function kalman_filter(y, a0, P0, Z, H, T, R, Q)
    Np, Nt = size(y)
    Nm = size(P0, 1)

    # preallocate arrays
    a = Array(Float64, Nm, Nt)
    P = Array(Float64, Nm, Nm, Nt)
    v = Array(Float64, Np, Nt)
    Finv = Array(Float64, Np, Np, Nt)
    K = Array(Float64, Nm, Np, Nt)

    # initialize
    a[:, 1] = a0;
    P[:, :, 1] = P0;

    # iterate
    for t in 1:Nt
        local a_t = slice(a, :, t)
        local P_t = slice(P, :, :, t)
        local y_t = slice(y, :, t)
        local Z_t = ndims(Z) < 3 ? Z : slice(Z, :, :, t)
        local H_t = ndims(H) < 3 ? H : slice(H, :, :, t)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)
        local R_t = ndims(R) < 3 ? R : slice(R, :, :, t)
        local Q_t = ndims(Q) < 3 ? Q : slice(Q, :, :, t)

        v_t = y_t - Z_t * a_t
        F_t = Z_t * P_t * Z_t' + H_t
        Finv_t = inv(F_t)
        K_t = T_t * P_t * Z_t' * Finv_t
        L_t = T_t - K_t * Z_t

        if t < Nt
            a[:, t + 1] = T_t * a_t + K_t * v_t
            P[:, :, t + 1] = T_t * P_t * L_t' + R_t * Q_t * R_t'
        end

        v[:, t] = v_t
        K[:, :, t] = K_t
        Finv[:, :, t] = Finv_t
    end

    return v, K, Finv, a, P
end

"""
Runs an efficient smoothing algorithm to get posterior mean of state
vector.

Parameters are as in kalman_filter.
"""
function fast_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)
    # infer dimensions
    Np, Nt = size(v)
    Nm = size(P0, 1)

    # preallocate
    r0 = 0.
    r = Array(Float64, Nm, Nt)
    r[:, end] = 0
    alpha = Array(Float64, Nm, Nt);

    # iterate backward
    for t in Nt:-1:1
        local Finv_t = slice(Finv, :, :, t)
        local v_t = slice(v, :, t)
        local K_t = slice(K, :, :, t)
        local r_t = slice(r, :, t)
        local Z_t = ndims(Z) < 3 ? Z : slice(Z, :, :, t)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)

        local u = Finv_t * v_t - K_t' * r_t
        local thisr = Z_t' * u + T_t' * r_t

        if t - 1 > 0
            r[:, t - 1] = thisr
        else
            r0 = thisr
        end
    end

    # run model forward
    alpha[:, 1] = a0 + P0 * r0

    for t in 1:(Nt - 1)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)
        local R_t = ndims(R) < 3 ? R : slice(R,  :, :, t)
        local Q_t = ndims(Q) < 3 ? Q : slice(Q,  :, :, t)
        local RQR = R_t * Q_t * R_t'

        alpha[:, t + 1] = T_t * alpha[:, t] + RQR * r[:, t]
    end

    return alpha
end

"""
Draw a sample of the state trajectory from the smoothed posterior.
"""
function simulate(y, a0, P0, Z, H, T, R, Q; interleaved=true)
    # get dimensions
    Np, Nt = size(y)
    Nm = size(a0, 1)
    Nr = size(Q, 1)

    # preallocate
    alpha_plus = Array(Float64, Nm, Nt)
    y_plus = Array(Float64, Np, Nt)

    # simulate data
    alpha_plus[:, 1] = rand(MvNormal(a0, P0))

    for t in 1:Nt
        local Z_t = ndims(Z) < 3 ? Z : slice(Z, :, :, t)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)
        local R_t = ndims(R) < 3 ? R : slice(R, :, :, t)
        # random number generators should take copies, not views?
        local H_t = ndims(H) < 3 ? H : H[:, :, t]
        local Q_t = ndims(Q) < 3 ? Q : Q[:, :, t]

        # draw disturbances
        local ϵ = rand(MvNormal(full(H_t)))
        y_plus[:, t] = Z_t * alpha_plus[:, t] + ϵ

        if t < Nt
            η = rand(MvNormal(full(Q_t)))
            alpha_plus[:, t + 1] = T_t * alpha_plus[:, t] + R_t * η
        end
    end

    # calculate smoothed means:
    if interleaved
        # actual data
        v, K, Finv, a, P = interleaved_kalman_filter(y, a0, P0, Z, H, T, R, Q)
        alpha_hat = interleaved_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)

        # simulated data
        v, K, Finv, a, P = interleaved_kalman_filter(y_plus, a0, P0, Z, H, T, R, Q)
        alpha_hat_plus = interleaved_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)

    else
        # actual data
        v, K, Finv, a, P = kalman_filter(y, a0, P0, Z, H, T, R, Q)
        alpha_hat = fast_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)

        # simulated data
        v, K, Finv, a, P = kalman_filter(y_plus, a0, P0, Z, H, T, R, Q)
        alpha_hat_plus = fast_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)
    end

    # combine
    alpha_draw = alpha_plus - alpha_hat_plus + alpha_hat

    return alpha_draw
end

"""
Like Kalman filter, but works by transforming a vector-valued
observation into a sequence of scalar observations.
"""
function interleaved_kalman_filter(y, a0, P0, Z, H, T, R, Q)
    # test for diagonality of H
    # H (minus its diagonals) should be the 0 matrix
    if ndims(H) == 2
        assert(H == diagm(diag(H)))
    end

    Np, Nt = size(y)
    Nm = size(P0, 1)

    # preallocate
    a = Array(Float64, Nm, Np + 1, Nt)
    P = Array(Float64, Nm, Nm, Np + 1, Nt)
    v = Array(Float64, Np, Nt)
    Finv = Array(Float64, Np, Nt)
    K = Array(Float64, Nm, Np, Nt)

    # initialize
    a[:, 1, 1] = a0
    P[:, :, 1, 1] = P0

    for t in 1:Nt, i in 1:Np
        local Z_t = ndims(Z) < 3 ? Z : slice(Z, :, :, t)
        local H_t = ndims(H) < 3 ? H : slice(H, :, :, t)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)
        local R_t = ndims(R) < 3 ? R : slice(R, :, :, t)
        local Q_t = ndims(Q) < 3 ? Q : slice(Q, :, :, t)
        z = vec(full(Z_t[i, :]))  # now a column vector

        v[i, t] = y[i, t] - dot(z, a[:, i, t])
        F = dot(z, P[:, :, i, t] * z) + H_t[i, i]
        Finv[i, t] = 1. / F
        if F != 0
            K[:, i, t] = P[:, :, i, t] * z * Finv[i, t]
        else
            K[:, i, t] = 0
        end

        a[:, i + 1, t] = a[:, i, t] + K[:, i, t] * v[i, t]
        P[:, :, i + 1, t] = P[:, :, i, t] - F * (K[:, i, t] * K[:, i, t]')

        if i == Np && t < Nt
            a[:, 1, t + 1] = T_t * a[:, Np + 1, t]
            P[:, :, 1, t + 1] = T_t * P[:, :, Np + 1, t] * T_t' + R_t * Q_t * R_t'
        end
    end

    return v, K, Finv, a, P
end

function interleaved_state_smoother(v, K, Finv, a0, P0, Z, T, R, Q)
    # infer dimensions
    Np, Nt = size(v)
    Nm = size(P0, 1)

    # preallocate
    r = Array(Float64, Nm, Np + 1, Nt)

    # initialize
    r[:, end, end] = 0

    # iterate backward
    for t in Nt:-1:1, i in (Np + 1):-1:2
        ii = i - 1  # handles offset between r and Z/K/F indices

        local Z_t = ndims(Z) < 3 ? Z : slice(Z, :, :, t)
        L = eye(Nm) - K[:, ii, t] * Z_t[ii, :]
        z = vec(full(Z_t[ii, :]))  # now a column vector

        r[:, i - 1, t] = z * (v[ii, t] * Finv[ii, t]) + L' * r[:, i, t]

        if i == 2 && t > 1
            # NOTE: t - 1 BELOW !!!
            local T_t = ndims(T) < 3 ? T : slice(T, :, :, t - 1)
            r[:, Np + 1, t - 1] = T_t' * r[:, 1, t]
        end
    end

    # initialize again
    α = Array(Float64, Nm, Nt)
    α[:, 1] = a0 + P0 * r[:, 1, 1]

    for t in 1:(Nt - 1)
        local T_t = ndims(T) < 3 ? T : slice(T, :, :, t)
        local R_t = ndims(R) < 3 ? R : slice(R, :, :, t)
        local Q_t = ndims(Q) < 3 ? Q : slice(Q, :, :, t)
        local RQR = R_t * Q_t * R_t'

        α[:, t + 1] = T_t * α[:, t] + RQR * r[:, 1, t + 1]
    end

    return α
end

end  # module

ext: jl

max_stars_repo_name: jmxpearson/labcr

max_issues_repo_name: jmxpearson/labcr

max_forks_repo_name: jmxpearson/labcr


text:
using PyCall
using StatsBase

import Base: convert, show

abstract type Layer end

convert(t::Type{T}, x::PyObject) where {T<:Layer} = error("convert(::Type{$T}, ::PyObject) not implemented.")
PyObject(l::T) where {T<:Layer} = error("PyObject(::$T) not implemented.")

function StatsBase.weights(l::T) where {T<:Layer}
    obj = PyObject(l)
    return obj[:get_weights]()
end

function weights!(l::T, W::Array) where {T<:Layer}
    obj = PyObject(l)
    obj[:set_weights](W)
end

function config(l::T) where {T<:Layer}
    obj = PyObject(l)
    return obj[:get_config]()
end

function input(l::T) where {T<:Layer}
    obj = PyObject(l)
    return Tensor(obj[:input])
end

function input(l::T, i::Int) where {T<:Layer}
    obj = PyObject(l)
    return Tensor(obj[:get_input_at](i))
end

function output(l::T) where {T<:Layer}
    obj = PyObject(l)
    return Tensor(obj[:output])
end

function output(l::T, i::Int) where {T<:Layer}
    obj = PyObject(l)
    return Tensor(obj[:get_output_at](i))
end

function input_shape(l::T) where {T<:Layer}
    obj = PyObject(l)
    return obj[:input_shape]
end

function input_shape(l::T, i::Int) where {T<:Layer}
    obj = PyObject(l)
    return obj[:get_input_shape_at](i)
end

function output_shape(l::T) where {T<:Layer}
    obj = PyObject(l)
    return obj[:output_shape]
end

function output_shape(l::T, i::Int) where {T<:Layer}
    obj = PyObject(l)
    return obj[:get_output_shape_at](i)
end

export Layer, weights!, config, input, output, input_shape, output_shape

module Layers

import PyCall: PyObject, pycall

import ..Keras
import ..Keras: PyDoc

#===========================================
Autogenerating Keras layer wrappers.
===========================================#
const keras_layers = Dict{String, Array}(
    "core" => [
        "Dense",
        "Activation",
        "Dropout",
        "Flatten",
        "Reshape",
        "Permute",
        "RepeatVector",
        "Lambda",
        "ActivityRegularization",
        "Masking",
    ],
    "convolutional" => [
        "Conv1D",
        "Conv2D",
        "SeparableConv2D",
        "Conv2DTranspose",
        "Conv3D",
        "Cropping1D",
        "Cropping2D",
        "Cropping3D",
        "UpSampling1D",
        "UpSampling2D",
        "UpSampling3D",
        "ZeroPadding1D",
        "ZeroPadding2D",
        "ZeroPadding3D",
    ],
    "pooling" => [
        "MaxPooling1D",
        "MaxPooling2D",
        "MaxPooling3D",
        "AveragePooling1D",
        "AveragePooling2D",
        "AveragePooling3D",
        "GlobalMaxPooling1D",
        "GlobalMaxPooling2D",
        "GlobalAveragePooling1D",
        "GlobalAveragePooling2D",
    ],
    "local" => [
        "LocallyConnected1D",
        "LocallyConnected2D"
    ],
    "recurrent" => [
        "Recurrent",
        "SimpleRNN",
        "GRU",
        "LSTM",
    ],
    "embeddings" => [
        "Embedding",
    ],
    "merge" => [
        "Add",
        "Multiply",
        "Average",
        "Maximum",
        "Concatenate",
        "Dot",
    ],
    "advanced_activations" => [
        "LeakyReLU",
        "PReLU",
        "ELU",
        "ThresholdedReLU",
    ],
    "normalization" => [
        "BatchNormalization",
    ],
    "noise" => [
        "GaussianNoise",
        "GaussianDropout",
    ],
    "wrappers" => [
        "TimeDistributed",
        "Bidirectional",
    ],
)

for (submod, layers) in keras_layers
    for l in layers
        layer_name = Symbol(l)

        @eval begin
            struct $layer_name <: Keras.Layer
                obj::PyObject

                @doc PyDoc(Keras._layers, Symbol($l)) function $layer_name(args...; kwargs...)
                    new(Keras._layers[Symbol($l)](args...; kwargs...))
                end
            end

            #convert(::Type{$(layer_name)}, obj::PyObject) = $layer_name(obj)
            PyObject(layer::$(layer_name)) = layer.obj
            # Base.Docs.doc(layer::$(layer_name)) = Base.Docs.doc(layer.obj)
            function (layer::$(layer_name))(args...; kwargs...)
                obj = PyObject(layer)
                return obj[:__call__](args...; kwargs...)
            end
        end
    end
end

const keras_merge_funcs = [
    "add",
    "multiply",
    "average",
    "maximum",
    "concatenate",
    "dot",
]

for l in keras_merge_funcs
    layer_name = Symbol(l)

    @eval begin
        @doc PyDoc(Keras._layers, Symbol($l)) function $layer_name(args...; kwargs...)
            Keras._layers[Symbol($l)](args...; kwargs...)
        end
    end
end

end


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export handle_request

#------------------------------------------------------------------------------
# Base Request Handler
#------------------------------------------------------------------------------

function handle_request (s::TCPSocket)
    r = nothing

    while isopen(s)
        r = form_request(s, r)

        info("[request] Formed request $r, dispatching")

        handle_request(s, r)
    end
end

#------------------------------------------------------------------------------
# Generic Request Handlers
#------------------------------------------------------------------------------

function handle_request (s::TCPSocket, r::CallRequest)
    x::Any = r.fun((r.args)...)
    info("[request] Applied $(r.fun) to $(r.args)")
    handle_response(s, x)
end

function handle_request (s::TCPSocket, r::CastRequest)
    info("[request] CastRequest received")
    noreply(s)
    x::Any = r.fun((r.args)...)
    x
end

function handle_request (s::TCPSocket, r::InfoRequest)
    info("[request] InfoRequest received")
    noreply(s)
end

function handle_request (s::TCPSocket, r::ErrorRequest)
    info("[request] ErrorRequest received")
    noreply_close(s)
    error("Error request $r signaled")
end

function handle_request (s::TCPSocket, r::StreamingCallRequest)
    info("[request] StreamingCallRequest received")
    bytes::Array{Uint8} = read_binary_stream(s)
    x::Any = r.fun((r.args)..., bytes)
    handle_response(s, x)
end

function handle_request (s::TCPSocket, r::CallbackCastRequest)
    info("[request] CallbackCastRequest received")
    noreply(s)
    x::Any = r.fun((r.args)...)
    # TODO: Initiate callback in the connected client
    msg = (:call, r.callback_mod, r.callback_fun, tuple(x))
end

#------------------------------------------------------------------------------
# Request Formation Functions
#------------------------------------------------------------------------------

function form_request (s::TCPSocket, r::Union(Request,Void))
    if r == nothing
        next_request(s)
    else
        next = next_request(s)
        modify(next, r)
    end
end

function next_request (s::TCPSocket)
    info("[request] Reading next request")
    bytes = read_erlang_term(s)
    info("[request] Read erlang term $bytes")
    msg, nbytes = decode_message(bytes)
    info("[request] Decoded erlang message $msg")
    is(msg[1], :call) ? call_request(msg[2], msg[3], tuple(msg[4]...)) :
    is(msg[1], :cast) ? cast_request(msg[2], msg[3], tuple(msg[4]...)) :
    is(msg[1], :info) ? InfoRequest(msg[2], msg[3]) :
    is(msg[1], :error) ? ErrorRequest(msg[2]) :
    error("[request] Message does not conform to protocol")
end

#------------------------------------------------------------------------------
# Request Modification Functions
#------------------------------------------------------------------------------

function modify(r1::InfoRequest, r2::InfoRequest)
    InfoRequest(r2)
end

function modify(r1::CallRequest, r2::InfoRequest)
    if is(r2.command, :stream)
        StreamingCallRequest(r2)
    else
        r1
    end
end

function modify(r1::CastRequest, r2::InfoRequest)
    if is(r2.command, :callback)
        CallbackCastRequest(r2)
    else
        r1
    end
end

max_issues_repo_name: esl/FsServer

max_forks_repo_name: esl/FsServer


text:
import Pkg, Libdl, PackageCompiler

env_to_precompile = ARGS[1]

sysimage_path = joinpath(env_to_precompile, "JuliaSysimage.$(Libdl.dlext)")

project_filename = isfile(joinpath(env_to_precompile, "JuliaProject.toml")) ? joinpath(env_to_precompile, "JuliaProject.toml") : joinpath(env_to_precompile, "Project.toml")

project = Pkg.API.read_project(project_filename)

used_packages = Symbol.(collect(keys(project.deps)))

@info "Now building a custom sysimage for the environment '$env_to_precompile'."

PackageCompiler.create_sysimage(used_packages, sysimage_path=sysimage_path, project=env_to_precompile)

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
## -------------  Logistic function

"""
	$(SIGNATURES)

Generalized logistic function object.

``
lb + (ub - lb) / (1 + lrShifter * exp(-slope * x)) ^ twister
``

Note that slope can be set to 1 if the input is `x / x0`. But it makes sense to have a slope when x0 depends on the college.

# Parameters
- `lb` and `ub` are the bounds of `f(x)`.
- `lrShifter >= 0` roughly shifts the curve left/right. It really sets the grad prob at h = 0. Higher `lrShifter` implies lower grad probs.
- `slope > 0` is the slope.
- `twister > 0` produces asymmetry (`twister = 1` is symmetric).
"""
@with_kw mutable struct GeneralizedLogistic
    lb :: Double = 0.0
    ub :: Double = 1.0
    lrShifter :: Double = 1.0
    slope :: Double = 1.0
    twister :: Double = 1.0
end

"""
	$(SIGNATURES)

Generalized logistic function.
"""
function logistic(x; lb = 0.0, ub = 1.0, lrShifter = 1.0, slope = 1.0, twister = 1.0)
    return lb .+ (ub - lb) ./ ((1 .+ lrShifter .* exp.(-slope .* x)) .^ twister);
end

function logistic(glf :: GeneralizedLogistic, x)
    return logistic(x; lb = glf.lb, ub = glf.ub, lrShifter = glf.lrShifter,
        slope = glf.slope, twister = glf.twister)
end

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


text:
"""
    arcmin_to_kpc(c::Cosmology.AbstractCosmology, θ::Unitful.AbstractQuantity, z::Real)

Convert `arcmin` to `kpc` for a given redshift `z` and cosmology `c`.
Conserves unit information.
"""
function arcmin_to_kpc(c::Cosmology.AbstractCosmology, θ::Unitful.AbstractQuantity, z::Real)

    # angular diameter distance in Mpc
    dA = angular_diameter_dist(c, z)

    return dA * 2π / 360 / 60.0u"arcminute" * θ |> u"kpc"
end

"""
    arcmin_to_kpc(c::Cosmology.AbstractCosmology, θ::Real, z::Real)

Convert `arcmin` to `kpc` for a given redshift `z` and cosmology `c`.
Deletes unit information.
"""
arcmin_to_kpc(c::Cosmology.AbstractCosmology, θ::Real, z::Real ) = arcmin_to_kpc( c, θ * 1.0u"arcminute", z) |> ustrip

"""
    arcmin_to_kpc(θ::Real, h::AbstractGadgetHeader)

Convert `arcmin` to `kpc` for a given redshift and cosmology taken from `AbstractGadgetHeader`.
"""
arcmin_to_kpc(θ::Real, h::AbstractGadgetHeader) = arcmin_to_kpc( cosmology(h), θ, h.z) |> ustrip

"""
    arcmin_to_kpc(θ::Unitful.AbstractQuantity, h::AbstractGadgetHeader)

Convert `arcmin` to `kpc` for a given redshift and cosmology taken from `AbstractGadgetHeader`.
"""
arcmin_to_kpc(θ::Unitful.AbstractQuantity, h::AbstractGadgetHeader) = arcmin_to_kpc( cosmology(h), θ, h.z)


text:
using StructJuMP

m = StructuredModel()

factories = 1:3
centers = 1:5
scenarios = ["lo", "mid", "hi"]

cap = [500, 450, 650]

demand = [150 160 170;
          100 120 135;
          250 270 300;
          300 325 350;
          600 700 800]

prob = [0.25, 0.5, 0.25]

transcost = [2.49 5.21 3.76 4.85 2.07;
             1.46 2.54 1.83 1.86 4.76;
             3.26 3.08 2.60 3.76 4.45]

prodcost = 14
price = 24
wastecost = 4

@variable(m, ship[factories,centers] >= 0)
@variable(m, product[factories] >= 0)
@variable(m, sales[centers] >= 0)
@variable(m, waste[centers] >= 0)
@variable(m, profit)

@variable(m, received[centers] >= 0)
@objective(m, Max, -sum(transcost[i,j]*ship[i,j] for i=factories, j=centers) + sum(prodcost*product[i] for i=factories))
for j in centers
    @constraint(m, received[j] == sum(ship[i,j] for i=factories))
end
for i in factories
    @constraint(m, product[i] == sum(ship[i,j] for j=centers))
end

for (s, elem) in enumerate(scenarios)
    bl = StructuredModel(parent=m, id=s)
    @variable(bl, 0 <= salesw[i=centers] <= demand[i,s])
    @variable(bl, wastew[centers] >= 0)
    @objective(bl, :Max, sum(price*prob[s]*salesw[j] for j=centers) - sum(wastecost*prob[s]*wastew[j] for j=centers))
    for j in centers
        @constraint(bl, received[j] == salesw[j]+wastew[j])
    end
end


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

mutable struct IoK8sApimachineryPkgApisMetaV1APIResource <: SwaggerModel
    categories::Any # spec type: Union{ Nothing, Vector{String} } # spec name: categories
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespaced::Any # spec type: Union{ Nothing, Bool } # spec name: namespaced
    shortNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: shortNames
    singularName::Any # spec type: Union{ Nothing, String } # spec name: singularName
    storageVersionHash::Any # spec type: Union{ Nothing, String } # spec name: storageVersionHash
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs
    version::Any # spec type: Union{ Nothing, String } # spec name: version

    function IoK8sApimachineryPkgApisMetaV1APIResource(;categories=nothing, group=nothing, kind=nothing, name=nothing, namespaced=nothing, shortNames=nothing, singularName=nothing, storageVersionHash=nothing, verbs=nothing, version=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("categories"), categories)
        setfield!(o, Symbol("categories"), categories)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("namespaced"), namespaced)
        setfield!(o, Symbol("namespaced"), namespaced)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("shortNames"), shortNames)
        setfield!(o, Symbol("shortNames"), shortNames)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("singularName"), singularName)
        setfield!(o, Symbol("singularName"), singularName)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("storageVersionHash"), storageVersionHash)
        setfield!(o, Symbol("storageVersionHash"), storageVersionHash)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResource, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIResource

const _property_map_IoK8sApimachineryPkgApisMetaV1APIResource = Dict{Symbol,Symbol}(Symbol("categories")=>Symbol("categories"), Symbol("group")=>Symbol("group"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("namespaced")=>Symbol("namespaced"), Symbol("shortNames")=>Symbol("shortNames"), Symbol("singularName")=>Symbol("singularName"), Symbol("storageVersionHash")=>Symbol("storageVersionHash"), Symbol("verbs")=>Symbol("verbs"), Symbol("version")=>Symbol("version"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIResource = Dict{Symbol,String}(Symbol("categories")=>"Vector{String}", Symbol("group")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("namespaced")=>"Bool", Symbol("shortNames")=>"Vector{String}", Symbol("singularName")=>"String", Symbol("storageVersionHash")=>"String", Symbol("verbs")=>"Vector{String}", Symbol("version")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIResource }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIResource))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIResource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIResource[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIResource }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIResource[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIResource)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("namespaced")) === nothing) && (return false)
    (getproperty(o, Symbol("singularName")) === nothing) && (return false)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1APIResource }, name::Symbol, val)
end


text:
import Zygote
struct ZygoteGradAlg <: GradAlg end
const ZygoteGrad = ZygoteGradAlg()
import ..Omega
Zygote.@nograd Omega.Space.increment

function gradmap(rv, ω::Ω)
  l = apl(rv, ω)
  params = Params(values(ω)) # We can avoid doing this every time.
  g = gradient(params) do
    rv(ω)
  end
  g
end

function lineargradient(rv, ω, ::ZygoteGradAlg)
  Zygote.gradient(ωvec -> apl(rv, unlinearize(ωvec, ω)), linearize(ω))
end

function grad(rv, ω, v, ::ZygoteGradAlg)
  grads = Zygote.gradient(Zygote.Params(v)) do
    rv(ω)
  end
end

grad(rv, ω, ::ZygoteGradAlg) = grad(rv, ω, values(ω), ZygoteGrad)

function gradarray(rv, ω, ::ZygoteGradAlg)
  vs = values(ω)
  grads_ = grad(rv, ω, vs, ZygoteGrad)
  map(v -> grads_[v], vs)
end


text:
# No info
struct Flat <: ContinuousUnivariateDistribution end

Distributions.rand(d::Flat) = rand()
Distributions.logpdf(d::Flat, x::T) where T<:Real= zero(x)
Distributions.minimum(d::Flat) = -Inf
Distributions.maximum(d::Flat) = +Inf

# For vec support
Distributions.rand(d::Flat, n::Int) = Vector([rand() for _ = 1:n])
Distributions.logpdf(d::Flat, x::AbstractVector{<:Real}) = zero(x)

# Pos
struct FlatPos{T<:Real} <: ContinuousUnivariateDistribution
    l::T
end

Distributions.rand(d::FlatPos) = rand() + d.l
Distributions.logpdf(d::FlatPos, x::Real) = x <= d.l ? -Inf : zero(x)
Distributions.minimum(d::FlatPos) = d.l
Distributions.maximum(d::FlatPos) = Inf

# For vec support
Distributions.rand(d::FlatPos, n::Int) = Vector([rand() for _ = 1:n] .+ d.l)
function Distributions.logpdf(d::FlatPos, x::AbstractVector{<:Real})
    return any(x .<= d.l) ? -Inf : zero(x)
end

# Binomial with logit
struct BinomialLogit{T<:Real, I<:Integer} <: DiscreteUnivariateDistribution
    n::I
    logitp::T
end

struct VecBinomialLogit{T<:Real, I<:Integer} <: DiscreteUnivariateDistribution
    n::Vector{I}
    logitp::Vector{T}
end

function logpdf_binomial_logit(n, logitp, k)
    logcomb = -StatsFuns.log1p(n) - SpecialFunctions.lbeta(n - k + 1, k + 1)
    return logcomb + k * logitp - n * StatsFuns.log1pexp(logitp)
end

function Distributions.logpdf(d::BinomialLogit{<:Real}, k::Int)
    return logpdf_binomial_logit(d.n, d.logitp, k)
end

function Distributions.logpdf(d::VecBinomialLogit{<:Real}, ks::Vector{<:Integer})
    return sum(logpdf_binomial_logit.(d.n, d.logitp, ks))
end


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

struct DefaultApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks whether a domain name in the cloudapp.net zone is available for use.
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: domainNameLabel::String
Return: DnsNameAvailabilityResult
"""
function checkDnsNameAvailability(_api::DefaultApi, location::String, api_version::String, subscriptionId::String; domainNameLabel=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DnsNameAvailabilityResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/locations/{location}/CheckDnsNameAvailability", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "domainNameLabel", domainNameLabel)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export checkDnsNameAvailability


text:
import RFFT
using Test, FFTW, LinearAlgebra

@testset begin
    for region in (1:2, 1, 2)
        for sz in ((5,6), (6,5))
            pair = RFFT.RCpair{Float64}(undef, sz, region)
            r = @inferred(real(pair))
            c = @inferred(complex(pair))
            b = rand(eltype(r), size(r))
            copyto!(r, b)
            RFFT.rfft!(pair)
            RFFT.irfft!(pair)
            @test r ≈ b
            pfwd = RFFT.plan_rfft!(pair)
            pinv = RFFT.plan_irfft!(pair)
            pinv(pfwd(pair))
            @test r ≈ b

            pair2 = copy(pair)
            @test real(pair2) == real(pair)
            @test complex(pair2) == complex(pair)
        end
    end
end


text:
module PETLION

using StatsBase: mean
using SciMLBase: DAEFunction, DAEProblem, step!, init
using Dierckx: Spline1D
using GeneralizedGenerated: mk_function, RuntimeFn
using LinearAlgebra: diagind, Tridiagonal, norm
using KLU: klu, klu!, KLUFactorization
using SparseArrays: sparse, findnz, SparseMatrixCSC, spzeros, spdiagm
using SparseDiffTools: matrix_colors, ForwardColorJacCache, forwarddiff_color_jacobian!
using RecursiveArrayTools: VectorOfArray
using Symbolics: @variables, Num, gradient, jacobian_sparsity, expand_derivatives, Differential, get_variables, sparsejacobian, substitute, simplify, build_function, IfElse
using RecipesBase
using SpecialFunctions: erf
using SHA: sha1

import PkgVersion
import SuiteSparse
import Sundials
import LinearAlgebra

# Must be loaded last
using BSON: @load, @save

export simulate, simulate!
export solution

export petlion
export boundary_stop_conditions, options_simulation, discretizations_per_section, options_numerical

export D_s_eff_isothermal, D_s_eff
export rxn_rate_isothermal, rxn_rate
export D_eff_linear, D_eff
export K_eff, K_eff_isothermal
export thermodynamic_factor_linear, thermodynamic_factor

export OCV_LCO,  OCV_NMC
export OCV_LiC6, OCV_LiC6_with_NMC

export rxn_BV, rxn_BV_γMod_01
export rxn_MHC

include("outputs.jl")
include("structures.jl")
include("params.jl")
include("external.jl")
include("set_variables.jl")
include("simulate.jl")
include("model_evaluation.jl")
include("checks.jl")
include("generate_functions.jl")
include("physics_equations/residuals.jl")
include("physics_equations/scalar_residual.jl")
include("physics_equations/auxiliary_states_and_coefficients.jl")
include("physics_equations/custom_functions.jl")
include("physics_equations/numerical_tools.jl")

## Backwards compatability
const Params = petlion
const run_model = simulate
const run_model! = simulate!
const options_model = options_simulation
const model_output = solution

export Params, run_model, run_model!, options_model

end # module


text:
"""
    model_infill(search_range::Vector{Tuple{Float64,Float64}},plan::AbstractArray{T,2},
    samples::AbstractArray{T,2},sm_interpolant; options::Options=Options()) where T

Infill function that takes calculates the location of new samples based
on the supplied options. The returned options are updated updated to facilitate
cycling through the infill objective functions.

...
# Arguments
- `search_range::Vector{Tuple{Float64,Float64}`:
    a vector of tuples containing the lower and upper limits
    for each dimension. `length(search_range)` is equal to number
    of dimensions.
- `plan::AbstractArray{T,2}`:
    sample locations where each column corresponds to the location of one point.
    `size(plan) = (num_dimensions,num_samples)`.
- `samples::AbstractArray{T,2}`:
    function value at each sample location. each column contains one value from
    the corresponding plan location. `size(samples) = (1,num_samples)`.
- `options=Options()`:
    all options available to customize the surrogate infill.
...
"""
function model_infill(search_range::Vector{Tuple{Float64,Float64}},plan::AbstractArray{T,2},
        samples::AbstractArray{T,2},sm_interpolant; options::Options=Options()) where T

    @unpack rbf_opt_gens, num_infill_points, trace,
            infill_funcs, infill_iterations = options

    trace && print_infill_head()

    (length(samples) != size(plan,2)) && error("plan and samples do not have the correct sizes")
    ((length(sm_interpolant(plan[:,1])) == 1) && (any(infill_funcs.==:std))) && error(":std infill criteria can not be used with one interpolant")

    # Extract the requested infill objective functions
    infill_obj_fun = infill_objective(sm_interpolant,plan,samples,infill_funcs)

    # Find the infill optimization points
    infill_plan,infill_type,infill_prediction,res_bboptim,options =
        infill_opt(search_range,infill_iterations,num_infill_points,infill_obj_fun,infill_funcs,plan,sm_interpolant,options)

    # Add additional points to fill up the desired number of infill points if required
    infill_plan,infill_type,infill_prediction =
        infill_add(sm_interpolant,samples,plan,infill_prediction,search_range,infill_plan,
        infill_type,num_infill_points,infill_obj_fun,infill_funcs,res_bboptim)

    trace && print_infill_tail(infill_type,num_infill_points,infill_plan,infill_prediction)

    return infill_plan, infill_type, infill_prediction, options
end

# function model_infill(search_range,plan,samples,sm_interpolant,vals,criteria,options)

#     @unpack rbf_opt_gens, num_interpolants, num_infill_points,
#             trace, infill_funcs, infill_iterations, categorical = options

#     if trace
#         println("Finding new infill samples ...")
#     end

#     #Infill function options
#     min_infill_fun = minimum_infill(sm_interpolant)
#     median_infill_fun = median_infill(sm_interpolant)
#     mean_infill_fun = mean_infill(sm_interpolant)
#     dist_infill_fun = distance_infill(plan,samples,sm_interpolant)
#     std_infill_fun = std_infill(sm_interpolant)

#     #Get the infill objective functions
#     call(f, x) = f(x)
#     library = Dict(
#         :min => x -> min_infill_fun(x),
#         :median => x -> median_infill_fun(x),
#         :mean => x -> mean_infill_fun(x),
#         :dist => x -> dist_infill_fun(x),
#         :std => x -> std_infill_fun(x)
#     )
#     functions_to_call = Tuple([library[s] for s in infill_funcs])
#     infill_obj_fun = function (x)
#         call.(functions_to_call, Ref(x))
#     end

#     num_infill_obj_funs = length(infill_obj_fun(plan[:,1]))

#     # Try generating pareto optimal infill points. Wrapped in try block due to
#     # the algorithm sometimes failing to generate a result.
#     infill_incomplete = true
#     j = 0
#     res_bboptim = nothing

#     while infill_incomplete && j < 51
#         try
#             res_bboptim = bboptimize(infill_obj_fun; Method=:borg_moea,
#                     FitnessScheme=ParetoFitnessScheme{num_infill_obj_funs}(is_minimizing=true),
#                     SearchRange=search_range, ϵ=0.00001,
#                     MaxFuncEvals=infill_iterations,
#                     MaxStepsWithoutProgress=20_000,TraceMode=:silent);

#             infill_incomplete = false
#         catch
#             j += 1
#         end
#     end

#     # Pick the infill points
#     infill_plan = Array{Float64,2}(undef,size(plan,1),0)
#     infill_type = Array{Symbol,1}()
#     infill_prediction = Array{Float64,1}()

#     # Cycle through the objective functions if the number of infill points is
#     # less than the number of objective functions.
#     if num_infill_points < num_infill_obj_funs

#         infill_obj_funs = Array{Int64}(undef,num_infill_points)
#         for i = 1:num_infill_points
#             infill_obj_funs[i] = first(Iterators.drop(Iterators.cycle(1:num_infill_obj_funs),criteria+num_infill_obj_funs-2+i))
#         end

#         # Find the next cyclic infill objective function to use
#         criteria = infill_obj_funs[end]+1
#     else
#         infill_obj_funs = 1:num_infill_obj_funs
#     end

#     for i in infill_obj_funs
#         pf = pareto_frontier(res_bboptim)
#         best_obj1, idx_obj1 = findmin(map(elm -> fitness(elm)[i], pf))
#         bo1_solution = BlackBoxOptim.params(pf[idx_obj1]) # get the solution candidate itself...

#         v = copy(permutedims(bo1_solution'))
#         if categorical
#             for (i,x_val) in enumerate(v)
#                 c_ind = closest_index(x_val, vals[i])
#                 v[i] = vals[i][c_ind]
#             end
#         end

#         # Add the infill point if it does not exist in the plan or infill_plan
#         if !colinmat(infill_plan,v) && !colinmat(plan,v)
#             push!(infill_prediction,median(sm_interpolant(vec(v))))
#             infill_plan = [infill_plan v]
#             push!(infill_type,infill_funcs[i])
#         end
#     end

#     # Add additional points to fill up the desired number of infill points
#     for i = 1:(num_infill_points-size(infill_plan,2))
#         infill_samples = Array{Float64,2}(undef,1,size(infill_plan,2))
#         for i = 1:size(infill_plan,2)
#             infill_samples[i] = median(sm_interpolant(infill_plan[:,i]))
#         end

#         # Find the point in the pareto front which is located furthest away from previous samples
#         dist_infill_fun = distance_infill([plan infill_plan],[samples infill_samples],sm_interpolant)
#         dist_obj = Inf
#         par_f = pareto_frontier(res_bboptim)
#         best_add_infill_solution = BlackBoxOptim.params(par_f[1])
#         for pf in par_f
#             bo1_solution = BlackBoxOptim.params(pf)
#             if (cur_val = dist_infill_fun(bo1_solution)) < dist_obj
#                 dist_obj = cur_val
#                 best_add_infill_solution = copy(bo1_solution)
#             end
#         end

#         # Add the infill point if it does not exist in the plan or infill_plan
#         v = copy(permutedims(best_add_infill_solution'))
#         if categorical
#             for (i,x_val) in enumerate(v)
#                 c_ind = closest_index(x_val, vals[i])
#                 v[i] = vals[i][c_ind]
#             end
#         end

#         if !colinmat(infill_plan,v) && !colinmat(plan,v)
#             push!(infill_prediction,median(sm_interpolant(vec(v))))
#             infill_plan = [infill_plan v]
#             push!(infill_type,:pareto)
#         end
#     end

#     # Add random points if there aren't enough in the pareto front
#     while (num_infill_points-size(infill_plan,2)) != 0

#         # Add the infill point if it does not exist in the plan or infill_plan
#         v = Array{Float64,2}(undef,size(infill_plan,1),1)
#         for i = 1:length(v)
#             v[i] = scale(rand(),search_range[i][1],search_range[i][2],old_min=0,old_max=1)
#         end

#         # Add the infill point if it does not exist in the plan or infill_plan
#         if categorical
#             for (i,x_val) in enumerate(v)
#                 c_ind = closest_index(x_val, vals[i])
#                 v[i] = vals[i][c_ind]
#             end
#         end
#         if !colinmat(infill_plan,v) && !colinmat(plan,v)
#             push!(infill_prediction,median(sm_interpolant(vec(v))))
#             infill_plan = [infill_plan v]
#             push!(infill_type,:rand)
#         end
#     end

#     if trace
#         println("Infill samples:")
#         println("---------------------------------------------------------------")
#         for i = (length(infill_type)-num_infill_points+1):length(infill_type)
#             print(@sprintf("%-15s",infill_type[i]))
#         end
#         print("\n")
#         println("---------------------------------------------------------------")
#         for j = 1:size(infill_plan,1)
#             for i = 1:size(infill_plan,2)
#                 print(@sprintf("%-15.7g",infill_plan[j,i]))
#             end
#             print("\n")
#         end
#         println("---------------------------------------------------------------")
#         _, min_loc = findmin(infill_prediction)
#         for i = 1:size(infill_plan,2)
#             if i == min_loc
#                 printstyled(@sprintf("%-15.7g",infill_prediction[i]); color=:light_green, bold=true)
#             else
#                 printstyled(@sprintf("%-15.7g",infill_prediction[i]))
#             end
#         end
#         print("\t prediction\n")
#         println("---------------------------------------------------------------")
#     end

#     return infill_plan, criteria, infill_type, infill_prediction
# end

# function model_infill_brute(search_range,plan,samples,sm_interpolant,vals,criteria,options)

#     @unpack rbf_opt_gens, num_interpolants, num_infill_points,
#             trace, infill_funcs, infill_iterations, categorical = options

#     if trace
#         println("Finding new infill samples ...")
#     end

#     #Infill function options
#     min_infill_fun = minimum_infill(sm_interpolant)
#     median_infill_fun = median_infill(sm_interpolant)
#     mean_infill_fun = mean_infill(sm_interpolant)
#     dist_infill_fun = distance_infill(plan,samples,sm_interpolant)
#     std_infill_fun = std_infill(sm_interpolant)

#     #Get the infill objective functions
#     call(f, x) = f(x)
#     library = Dict(
#         :min => x -> min_infill_fun(x),
#         :median => x -> median_infill_fun(x),
#         :mean => x -> mean_infill_fun(x),
#         :dist => x -> dist_infill_fun(x),
#         :std => x -> std_infill_fun(x)
#     )
#     functions_to_call = Tuple([library[s] for s in infill_funcs])
#     infill_obj_fun = function (x)
#         call.(functions_to_call, Ref(x))
#     end

#     num_infill_obj_funs = length(infill_obj_fun(plan[:,1]))

#     # Pick the infill points
#     infill_plan = Array{Float64,2}(undef,size(plan,1),0)
#     infill_type = Array{Symbol,1}()
#     infill_prediction = Array{Float64,1}()

#     # Cycle through the objective functions if the number of infill points is
#     # less than the number of objective functions.
#     if num_infill_points < num_infill_obj_funs

#         infill_obj_funs = Array{Int64}(undef,num_infill_points)
#         for i = 1:num_infill_points
#             infill_obj_funs[i] = first(Iterators.drop(Iterators.cycle(1:num_infill_obj_funs),criteria+num_infill_obj_funs-2+i))
#         end

#         # Find the next cyclic infill objective function to use
#         criteria = infill_obj_funs[end]+1
#     else
#         infill_obj_funs = 1:num_infill_obj_funs
#     end

#     for i in infill_obj_funs
#         fun = functions_to_call[i]
#         possible_designs = vec(collect(Base.Iterators.product(vals...)))

#         fun_out = [fun([design...]) for design in possible_designs]

#         min_locs = sortperm(fun_out)
#         found_min = false
#         j = 1
#         while found_min == false
#             v = [possible_designs[min_locs[j]]...]

#             # Add the infill point if it does not exist in the plan or infill_plan
#             if !colinmat(infill_plan,v) && !colinmat(plan,v)
#                 push!(infill_prediction,median(sm_interpolant(vec(v))))
#                 infill_plan = [infill_plan v]
#                 push!(infill_type,infill_funcs[i])
#                 found_min=true
#             end
#             j+=1
#         end
#     end

#     if trace
#         println("Infill samples:")
#         println("---------------------------------------------------------------")
#         for i = (length(infill_type)-num_infill_points+1):length(infill_type)
#             print(@sprintf("%-15s",infill_type[i]))
#         end
#         print("\n")
#         println("---------------------------------------------------------------")
#         for j = 1:size(infill_plan,1)
#             for i = 1:size(infill_plan,2)
#                 print(@sprintf("%-15.7g",infill_plan[j,i]))
#             end
#             print("\n")
#         end
#         println("---------------------------------------------------------------")
#         _, min_loc = findmin(infill_prediction)
#         for i = 1:size(infill_plan,2)
#             if i == min_loc
#                 printstyled(@sprintf("%-15.7g",infill_prediction[i]); color=:light_green, bold=true)
#             else
#                 printstyled(@sprintf("%-15.7g",infill_prediction[i]))
#             end
#         end
#         print("\t prediction\n")
#         println("---------------------------------------------------------------")
#     end

#     return infill_plan, criteria, infill_type, infill_prediction
# end


text:
StatsBase.describe(X::NullableVector) = StatsBase.describe(STDOUT, X)

function StatsBase.describe{T<:Real}(io::IO, X::NullableVector{T})
    nullcount = sum(X.isnull)
    pnull = 100nullcount/length(X)
    if pnull != 100 # describe will fail if dropnull returns an empty vector
        describe(io, dropnull(X))
    else
        println(io, "Summary Stats:")
        println(io, "Type:           $(eltype(X))")
    end
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe(io::IO, X::NullableVector)
    nullcount = sum(X.isnull)
    pnull = 100nullcount/length(X)
    println(io, "Summary Stats:")
    println(io, "Length:         $(length(X))")
    println(io, "Type:           $(eltype(X))")
    println(io, "Number Unique:  $(length(unique(X)))")
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe{T<:Real}(io::IO, X::AbstractVector{Nullable{T}})
    nullcount = sum(_isnull, X)
    pnull = 100nullcount/length(X)
    if pnull != 100 # describe will fail if dropnull returns an empty vector
        describe(io, dropnull(X))
    else
        println(io, "Summary Stats:")
        println(io, "Type:           $(eltype(X))")
    end
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe{T<:Nullable}(io::IO, X::AbstractVector{T})
    nullcount = sum(_isnull, X)
    pnull = 100nullcount/length(X)
    println(io, "Summary Stats:")
    println(io, "Length:         $(length(X))")
    println(io, "Type:           $(eltype(X))")
    println(io, "Number Unique:  $(length(unique(X)))")
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end


text:
include("../src/TrajectoryMatrices.jl")
using TrajectoryMatrices
using Base.Test

info("Beginning tests...")
data = collect(1:21)

# Make sure everything is compiled before testing
traj = TrajectoryMatrix(data, 3)

@test traj[1,3] == 3
@test traj[19, 3] == 21
@test size(traj) == (19,3)

trans_traj = traj'

@test trans_traj[3,1] == 3
@test trans_traj[3, 19] == 21
@test size(trans_traj) == (3,19)

data_mem = @allocated collect(1:21)
traj_mem = @allocated TrajectoryMatrix(data, 3)
tran_mem = @allocated traj'

@test traj_mem < (data_mem * 2)
@test tran_mem < (traj_mem * 3)

data = [1 1 1 1; 2 2 2 2; 3 3 3 3]
expected = [1 1 1 1 2 2 2 2; 2 2 2 2 3 3 3 3]
traj = TrajectoryMatrix(data, 2)

for i in eachindex(expected)
    @test expected[i] == traj[i]
end

expected = expected'
traj = traj'

for i in eachindex(expected)
    @test expected[i] == traj[i]
end

traj[1] = 8
@test traj[1] == 8

traj[1,1] = 9
@test traj[1,1] == 9

info("All tests passed.")


text:
module TestDistributions

# using Revise
using Test
using MLJBase
using CategoricalArrays
import Distributions

## UNIVARIATE NOMINAL

v = collect("asdfghjklzxc")
d = UnivariateNominal(v, [0.09, 0.02, 0.1, 0.1, 0.1, 0.1, 0.1, 0.11, 0.01, 0.1, 0.07, 0.1])
@test pdf(d, 's') ≈ 0.02
@test mode(d) == 'k'
rand(d, 5)

v = collect("abcd")
d = UnivariateNominal(v, [0.2, 0.3, 0.1, 0.4])
sample = rand(d, 10^4)
freq_given_label = Distributions.countmap(sample)
pairs  = collect(freq_given_label)
sort!(pairs, by=pair->pair[2], alg=QuickSort)
sorted_labels = first.(pairs)
# if this fails it is bug or an exceedingly rare event or a bug:
@test sorted_labels == ['c', 'a', 'b', 'd']

v=['a', 'b', 'a', 'b', 'c', 'b', 'a', 'a']
d = Distributions.fit(UnivariateNominal, v)
@test pdf(d, 'a') ≈ 0.5
@test pdf(d, 'b') ≈ 0.375
@test pdf(d, 'c') ≈ 0.125

# to check fitting to categorical returns zero prob for missing
# levels, we add and drop new level to a categorical version of v:
w = categorical(append!(v, 'f'))
w = w[1:end-1]
e = Distributions.fit(UnivariateNominal, w)
@test e.prob_given_label['f'] == 0

# arithmetic
d1 = UnivariateNominal(['a', 'b'], [0.2, 0.8])
d2 = UnivariateNominal(['b', 'c'], [0.3, 0.7])
d = average([d1, d2])
@test d.prob_given_label['a'] ≈ 0.1 && d.prob_given_label['b'] ≈ 0.55 &&  d.prob_given_label['c'] ≈ 0.35

end # module

true


text:
asparams(::Affine, ::Val{:μ}) = asℝ
asparams(::Affine, ::Val{:σ}) = asℝ₊
asparams(::Type{A}, ::Val{:μ}) where {A<:Affine} = asℝ
asparams(::Type{A}, ::Val{:σ}) where {A<:Affine} = asℝ₊

asparams(::Affine{N,M,T}, ::Val{:μ}) where {N,M,T<:AbstractArray} = as(Array, asℝ, size(d.μ))
asparams(::Affine{N,M,T}, ::Val{:σ}) where {N,M,T<:AbstractArray} = as(Array, asℝ, size(d.σ))


text:
function correlation(sample, k)
    # Stabilized correlation function χ(k) = χ(|i-j|) = ⟨xᵢxⱼ⟩ - ⟨xᵢ⟩⟨xⱼ⟩
    # Following Welford derivation from https://changyaochen.github.io/welford/

    M = length(sample)
    corr = sumx = sumy = 0.0
    for i in 1:M-k
        # ̅xᵢ₊₁ = ̅xᵢ + (xᵢ₊₁ - ̅xᵢ) / (N+1)
        # ̅yᵢ₊₁ = ̅yᵢ + (yᵢ₊₁ - ̅yᵢ) / (N+1)
        # χᵢ₊₁ = χᵢ + [(xᵢ₊₁ - ̅xᵢ)(yᵢ₊₁ - ̅yᵢ₊₁) - χᵢ] / (i + 1)

        invN = 1.0 / i
        sumy += invN * (sample[i + k] - sumy)
        sumx_delta = invN * (sample[i] - sumx)
        corr += sumx_delta * (sample[i + k] - sumy) - invN * corr
        sumx += sumx_delta
    end
    return corr
end

# QMCM eq 3.15
"""
    unbinned_tau(sample[; truncate = true, max_rel_err = 0.0, min_sample_size = 32])

Estimates the autocorrelation time `τ = ∑ₖ (1 - k/N) χₖ/χ₀` from a sample using the
correlation function `χₖ = ⟨xᵢ xᵢ₊ₖ⟩ - ⟨xᵢ⟩⟨xᵢ₊ₖ⟩` where x are values in the sample.

Note that by default this function truncated the sum over k if
`v / τ < max_rel_err / (N - k)` or if the `k > N - min_sample_size`. The former
cuts of fluctuations around 0 and can be turned of with `truncate = false`.
"""
function unbinned_tau(sample; truncate = true, max_rel_err = 0.0, min_sample_size = 32)
    tau = 0.0

    for k in 1:length(sample) - max(1, min_sample_size)
        v = (1 - k/length(sample)) * correlation(sample, k)
        tau += v

        # We assume the worst case is a constant correlation for the tail.
        # There are (length(sample) - k) values left
        # The average prefactor is 0.5 * (1 - k/length(sample))
        # If the sum of the tail is smaller than max_rel_err * tau, i.e. if tau
        # can at most increase by a factor of max_rel_error, we stop the loop
        if truncate && 0.5 * v * (length(sample) - k) < max_rel_err * tau
            @debug "Cancelled unbinned_tau summation after $k iterations."
            break
        end
    end

    return tau / correlation(sample, 0)
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestMOF

import JSON
import JSONSchema
import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOF = MOI.FileFormats.MOF

const TEST_MOF_FILE = "test.mof.json"

const SCHEMA =
    JSONSchema.Schema(JSON.parsefile(MOI.FileFormats.MOF.SCHEMA_PATH))

function _validate(filename::String)
    MOI.FileFormats.compressed_open(
        filename,
        "r",
        MOI.FileFormats.AutomaticCompression(),
    ) do io
        object = JSON.parse(io)
        ret = JSONSchema.validate(SCHEMA, object)
        if ret !== nothing
            error(
                "Unable to read file because it does not conform to the MOF " *
                "schema: ",
                ret,
            )
        end
    end
end

struct UnsupportedSet <: MOI.AbstractSet end
struct UnsupportedFunction <: MOI.AbstractFunction end

function _test_model_equality(model_string, variables, constraints; suffix = "")
    model = MOF.Model()
    MOIU.loadfromstring!(model, model_string)
    MOI.write_to_file(model, TEST_MOF_FILE * suffix)
    model_2 = MOF.Model()
    MOI.read_from_file(model_2, TEST_MOF_FILE * suffix)
    MOI.Test.util_test_models_equal(model, model_2, variables, constraints)
    return _validate(TEST_MOF_FILE * suffix)
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
struct ExprEvaluator <: MOI.AbstractNLPEvaluator
    objective::Expr
    constraints::Vector{Expr}
end
MOI.features_available(::ExprEvaluator) = [:ExprGraph]
MOI.initialize(::ExprEvaluator, features) = nothing
MOI.objective_expr(evaluator::ExprEvaluator) = evaluator.objective
MOI.constraint_expr(evaluator::ExprEvaluator, i::Int) = evaluator.constraints[i]

function HS071(x::Vector{MOI.VariableIndex})
    x1, x2, x3, x4 = x
    return MOI.NLPBlockData(
        MOI.NLPBoundsPair.([25, 40], [Inf, 40]),
        ExprEvaluator(
            :(x[$x1] * x[$x4] * (x[$x1] + x[$x2] + x[$x3]) + x[$x3]),
            [
                :(x[$x1] * x[$x2] * x[$x3] * x[$x4] >= 25),
                :(x[$x1]^2 + x[$x2]^2 + x[$x3]^2 + x[$x4]^2 == 40),
            ],
        ),
        true,
    )
end

function test_HS071()
    model = MOF.Model()
    x = MOI.add_variables(model, 4)
    for (index, variable) in enumerate(x)
        MOI.set(model, MOI.VariableName(), variable, "var_$(index)")
    end
    MOI.add_constraints(model, x, Ref(MOI.Interval(1.0, 5.0)))
    MOI.set(model, MOI.NLPBlock(), HS071(x))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.write_to_file(model, TEST_MOF_FILE)
    @test replace(read(TEST_MOF_FILE, String), '\r' => "") ==
          replace(read(joinpath(@__DIR__, "nlp.mof.json"), String), '\r' => "")
    return _validate(TEST_MOF_FILE)
end

function test_nonlinear_error_handling()
    node_list = MOF.Object[]
    string_to_variable = Dict{String,MOI.VariableIndex}()
    variable_to_string = Dict{MOI.VariableIndex,String}()
    # Test unsupported function for Expr -> MOF.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(not_supported_function(x)),
        node_list,
        variable_to_string,
    )
    # Test unsupported function for MOF -> Expr.
    @test_throws Exception MOF.convert_mof_to_expr(
        MOF.OrderedObject("type" => "not_supported_function", "value" => 1),
        node_list,
        string_to_variable,
    )
    # Test n-ary function with no arguments.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(min()),
        node_list,
        variable_to_string,
    )
    # Test unary function with two arguments.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(sin(x, y)),
        node_list,
        variable_to_string,
    )
    # Test binary function with one arguments.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(^(x)),
        node_list,
        variable_to_string,
    )
    # An expression with something other than :call as the head.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(a <= b <= c),
        node_list,
        variable_to_string,
    )
    # Hit the default fallback with an un-interpolated complex number.
    @test_throws Exception MOF.convert_expr_to_mof(
        :(1 + 2im),
        node_list,
        variable_to_string,
    )
    # Invalid number of variables.
    @test_throws Exception MOF.substitute_variables(
        :(x[1] * x[2]),
        [MOI.VariableIndex(1)],
    )
    # Function-in-Set
    @test_throws Exception MOF.extract_function_and_set(:(foo in set))
    # Not a constraint.
    @test_throws Exception MOF.extract_function_and_set(:(x^2))
    # Two-sided constraints
    @test MOF.extract_function_and_set(:(1 <= x <= 2)) ==
          MOF.extract_function_and_set(:(2 >= x >= 1)) ==
          (:x, MOI.Interval(1, 2))
    # Less-than constraint.
    @test MOF.extract_function_and_set(:(x <= 2)) == (:x, MOI.LessThan(2))
end

function test_Roundtrip_nonlinear_expressions()
    x = MOI.VariableIndex(123)
    y = MOI.VariableIndex(456)
    z = MOI.VariableIndex(789)
    string_to_var = Dict{String,MOI.VariableIndex}("x" => x, "y" => y, "z" => z)
    var_to_string = Dict{MOI.VariableIndex,String}(x => "x", y => "y", z => "z")
    for expr in [
        2,
        2.34,
        2 + 3im,
        x,
        :(1 + $x),
        :($x - 1),
        :($x + $y),
        :($x + $y - $z),
        :(2 * $x),
        :($x * $y),
        :($x / 2),
        :(2 / $x),
        :($x / $y),
        :($x / $y / $z),
        :(2^$x),
        :($x^2),
        :($x^$y),
        :($x^(2 * $y + 1)),
        :(sin($x)),
        :(sin($x + $y)),
        :(2 * $x + sin($x)^2 + $y),
        :(sin($(3im))^2 + cos($(3im))^2),
        :($(1 + 2im) * $x),
        :(ceil($x)),
        :(floor($x)),
        :($x < $y),
        :($x <= $y),
        :($x > $y),
        :($x >= $y),
        :($x == $y),
        :($x != $y),
        # :($x && $y), :($x || $y),
        :(ifelse($x > 0, 1, $y)),
    ]
        node_list = MOF.OrderedObject[]
        object = MOF.convert_expr_to_mof(expr, node_list, var_to_string)
        @test MOF.convert_mof_to_expr(object, node_list, string_to_var) == expr
    end
end

function test_nonlinear_readingwriting()
    # Write to file.
    model = MOF.Model()
    (x, y) = MOI.add_variables(model, 2)
    MOI.set(model, MOI.VariableName(), x, "var_x")
    MOI.set(model, MOI.VariableName(), y, "y")
    con = MOI.add_constraint(
        model,
        MOF.Nonlinear(:(2 * $x + sin($x)^2 - $y)),
        MOI.EqualTo(1.0),
    )
    MOI.set(model, MOI.ConstraintName(), con, "con")
    MOI.write_to_file(model, TEST_MOF_FILE)
    # Read the model back in.
    model2 = MOF.Model()
    MOI.read_from_file(model2, TEST_MOF_FILE)
    con2 = MOI.get(model2, MOI.ConstraintIndex, "con")
    foo2 = MOI.get(model2, MOI.ConstraintFunction(), con2)
    # Test that we recover the constraint.
    @test foo2.expr == :(2 * $x + sin($x)^2 - $y)
    @test MOI.get(model, MOI.ConstraintSet(), con) ==
          MOI.get(model2, MOI.ConstraintSet(), con2)
    return _validate(TEST_MOF_FILE)
end

function test_show()
    @test sprint(show, MOF.Model()) == "A MathOptFormat Model"
end

function test_nonempty_model()
    model = MOF.Model(warn = true)
    MOI.add_variable(model)
    @test !MOI.is_empty(model)
    exception = ErrorException(
        "Cannot read model from file as destination model is not empty.",
    )
    @test_throws exception MOI.read_from_file(
        model,
        joinpath(@__DIR__, "empty_model.mof.json"),
    )
    options = MOF.get_options(model)
    @test options.warn
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOI.read_from_file(model, joinpath(@__DIR__, "empty_model.mof.json"))
    options2 = MOF.get_options(model)
    @test options2.warn
end

function test_failing_models()
    @testset "$(filename)" for filename in filter(
        f -> endswith(f, ".mof.json"),
        readdir(joinpath(@__DIR__, "failing_models")),
    )
        @test_throws Exception MOI.read_from_file(
            MOF.Model(),
            joinpath(@__DIR__, "failing_models", filename),
        )
    end
end

function test_Blank_variable_name()
    model = MOF.Model()
    variable = MOI.add_variable(model)
    @test_throws Exception MOF.moi_to_object(variable, model)
    MOI.FileFormats.create_unique_names(model, warn = true)
    @test MOF.moi_to_object(variable, model) ==
          MOF.OrderedObject("name" => "x1")
end

function test_Duplicate_variable_name()
    model = MOF.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    y = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), y, "x")
    @test MOF.moi_to_object(x, model) == MOF.OrderedObject("name" => "x")
    @test MOF.moi_to_object(y, model) == MOF.OrderedObject("name" => "x")
    MOI.FileFormats.create_unique_names(model, warn = true)
    @test MOF.moi_to_object(x, model) == MOF.OrderedObject("name" => "x")
    @test MOF.moi_to_object(y, model) == MOF.OrderedObject("name" => "x_1")
end

function test_Blank_constraint_name()
    model = MOF.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    c = MOI.add_constraint(model, f, MOI.ZeroOne())
    name_map = Dict(x => "x")
    MOI.FileFormats.create_unique_names(model, warn = true)
    @test MOF.moi_to_object(c, model, name_map)["name"] == "c1"
end

function test_Duplicate_constraint_name()
    model = MOF.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    c1 = MOI.add_constraint(model, f, MOI.LessThan(1.0))
    c2 = MOI.add_constraint(model, f, MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ConstraintName(), c1, "c")
    MOI.set(model, MOI.ConstraintName(), c2, "c")
    name_map = Dict(x => "x")
    @test MOF.moi_to_object(c1, model, name_map)["name"] == "c"
    @test MOF.moi_to_object(c2, model, name_map)["name"] == "c"
    MOI.FileFormats.create_unique_names(model, warn = true)
    @test MOF.moi_to_object(c1, model, name_map)["name"] == "c_1"
    @test MOF.moi_to_object(c2, model, name_map)["name"] == "c"
end

function test_empty_model()
    model = MOF.Model()
    MOI.write_to_file(model, TEST_MOF_FILE)
    model_2 = MOF.Model()
    MOI.read_from_file(model_2, TEST_MOF_FILE)
    return MOI.Test.util_test_models_equal(model, model_2, String[], String[])
end

function test_FEASIBILITY_SENSE()
    model = MOF.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.write_to_file(model, TEST_MOF_FILE)
    model_2 = MOF.Model()
    MOI.read_from_file(model_2, TEST_MOF_FILE)
    return MOI.Test.util_test_models_equal(model, model_2, ["x"], String[])
end

function test_empty_function_term()
    model = MOF.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
        MOI.GreaterThan(1.0),
    )
    MOI.set(model, MOI.ConstraintName(), c, "c")
    MOI.write_to_file(model, TEST_MOF_FILE)
    model_2 = MOF.Model()
    MOI.read_from_file(model_2, TEST_MOF_FILE)
    return MOI.Test.util_test_models_equal(model, model_2, ["x"], ["c"])
end

function test_min_objective()
    return _test_model_equality(
        """
variables: x
minobjective: x
""",
        ["x"],
        String[],
    )
end

function test_max_objective()
    return _test_model_equality(
        """
variables: x
maxobjective: x
""",
        ["x"],
        String[],
        suffix = ".gz",
    )
end

function test_min_scalaraffine()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
""",
        ["x"],
        String[],
    )
end

function test_max_scalaraffine()
    return _test_model_equality(
        """
variables: x
maxobjective: 1.2x + 0.5
""",
        ["x"],
        String[],
        suffix = ".gz",
    )
end

function test_singlevariable_in_lower()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x >= 1.0
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_upper()
    return _test_model_equality(
        """
variables: x
maxobjective: 1.2x + 0.5
x <= 1.0
""",
        ["x"],
        String[],
        suffix = ".gz",
    )
end

function test_singlevariable_in_interval()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x in Interval(1.0, 2.0)
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_equalto()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x == 1.0
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_zeroone()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x in ZeroOne()
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_integer()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x in Integer()
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_Semicontinuous()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x in Semicontinuous(1.0, 2.0)
""",
        ["x"],
        String[],
    )
end

function test_singlevariable_in_Semiinteger()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x + 0.5
x in Semiinteger(1.0, 2.0)
""",
        ["x"],
        String[],
    )
end

function test_scalarquadratic_objective()
    return _test_model_equality(
        """
variables: x
minobjective: 1.0*x*x + -2.0x + 1.0
""",
        ["x"],
        String[],
    )
end

function test_SOS1()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in SOS1([1.0, 2.0, 3.0])
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_SOS2()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in SOS2([1.0, 2.0, 3.0])
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_Reals()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in Reals(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_Zeros()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in Zeros(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_Nonnegatives()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in Nonnegatives(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_Nonpositives()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in Nonpositives(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_PowerCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in PowerCone(2.0)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_DualPowerCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in DualPowerCone(0.5)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_GeometricMeanCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in GeometricMeanCone(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_Complements()
    return _test_model_equality(
        "variables: x, y\nc1: [x, y] in Complements(2)",
        ["x", "y"],
        ["c1"],
    )
end

function test_vectoraffine_in_zeros()
    return _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [1.0x + -3.0, 2.0y + -4.0] in Zeros(2)
""",
        ["x", "y"],
        ["c1"],
    )
end

function test_vectorquadratic_in_nonnegatives()
    return _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [1.0*x*x + -2.0x + 1.0, 2.0y + -4.0] in Nonnegatives(2)
""",
        ["x", "y"],
        ["c1"],
    )
end

function test_ExponentialCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in ExponentialCone()
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_DualExponentialCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in DualExponentialCone()
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_SecondOrderCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in SecondOrderCone(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_RotatedSecondOrderCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in RotatedSecondOrderCone(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_PositiveSemidefiniteConeTriangle()
    return _test_model_equality(
        """
variables: x1, x2, x3
minobjective: x1
c1: [x1, x2, x3] in PositiveSemidefiniteConeTriangle(2)
""",
        ["x1", "x2", "x3"],
        ["c1"],
    )
end

function test_PositiveSemidefiniteConeSquare()
    return _test_model_equality(
        """
variables: x1, x2, x3, x4
minobjective: x1
c1: [x1, x2, x3, x4] in PositiveSemidefiniteConeSquare(2)
""",
        ["x1", "x2", "x3", "x4"],
        ["c1"],
    )
end

function test_LogDetConeTriangle()
    return _test_model_equality(
        """
variables: t, u, x1, x2, x3
minobjective: x1
c1: [t, u, x1, x2, x3] in LogDetConeTriangle(2)
""",
        ["t", "u", "x1", "x2", "x3"],
        ["c1"],
    )
end

function test_LogDetConeSquare()
    return _test_model_equality(
        """
variables: t, u, x1, x2, x3, x4
minobjective: x1
c1: [t, u, x1, x2, x3, x4] in LogDetConeSquare(2)
""",
        ["t", "u", "x1", "x2", "x3", "x4"],
        ["c1"],
    )
end

function test_RootDetConeTriangle()
    return _test_model_equality(
        """
variables: t, x1, x2, x3
minobjective: x1
c1: [t, x1, x2, x3] in RootDetConeTriangle(2)
""",
        ["t", "x1", "x2", "x3"],
        ["c1"],
    )
end

function test_RootDetConeSquare()
    return _test_model_equality(
        """
variables: t, x1, x2, x3, x4
minobjective: x1
c1: [t, x1, x2, x3, x4] in RootDetConeSquare(2)
""",
        ["t", "x1", "x2", "x3", "x4"],
        ["c1"],
    )
end

function test_Indicator()
    _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [x, y] in Indicator{ACTIVATE_ON_ONE}(GreaterThan(1.0))
""",
        ["x", "y"],
        ["c1"],
    )

    return _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [x, y] in Indicator{ACTIVATE_ON_ZERO}(GreaterThan(1.0))
""",
        ["x", "y"],
        ["c1"],
    )
end

function test_NormOneCone()
    return _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [x, y] in NormOneCone(2)
""",
        ["x", "y"],
        ["c1"],
    )
end

function test_NormInfinityCone()
    return _test_model_equality(
        """
variables: x, y
minobjective: x
c1: [x, y] in NormInfinityCone(2)
""",
        ["x", "y"],
        ["c1"],
    )
end

function test_RelativeEntropyCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in RelativeEntropyCone(3)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_NormSpectralCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in NormSpectralCone(1, 2)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_NormNuclearCone()
    return _test_model_equality(
        """
variables: x, y, z
minobjective: x
c1: [x, y, z] in NormNuclearCone(1, 2)
""",
        ["x", "y", "z"],
        ["c1"],
    )
end

function test_v04()
    model = MOF.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "v0.4.mof.json"))
    model_2 = MOF.Model()
    MOI.Utilities.loadfromstring!(
        model_2,
        """
variables: x, y
minobjective: x
c: x + y >= 1.0
x in Interval(0.0, 1.0)
y in ZeroOne()
""",
    )
    MOI.Test.util_test_models_equal(model, model_2, ["x", "y"], ["c"])
    return
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # allow time for unlink to happen
    rm(TEST_MOF_FILE, force = true)
    rm(TEST_MOF_FILE * ".gz", force = true)
    return
end

end

TestMOF.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Formatting
include("ncurses.jl")

mutable struct Game
	screen_size::NamedTuple{(:x, :y),Tuple{Int64,Int64}}
	player_pos::NamedTuple{(:x, :y),Tuple{Int64,Int64}}
end

function Game_new()
	# startup
	stdscreen = initscr() # gets things started and gives us root window
	cbreak()              # normalizes the terminal
	curs_set(0)           # hide cursor
	noecho()              # dont wait until newline
	nodelay(stdscreen,1)  # this makes input non blocking and returns -1 if not key

	clear()
	max_y,max_x = getmaxyx(stdscreen)
	player =
	Game(
		(x=max_x,y=max_y),
		(x=div(max_x,2),y=div(max_y,2))
	)
end

function Game_render(game::Game)
	clear()
	mvprintw(game.player_pos.y,game.player_pos.x,"@")
	refresh()
end

function Game_run(game::Game)
	# Render first time to show something
	Game_render(game)
	while true
		i = getch()
		if i == -1
			# we have recieved nothing
			continue
		elseif i == 27
			i = getch()
			if i == -1
				# we have received an ESCAPE key
				break
			elseif i == 91
				# We might have an arrow
				i = getch()
				if i == 65
					game.player_pos = (x=game.player_pos.x,y=max(game.player_pos.y-1,0))
				elseif i == 66
					game.player_pos = (x=game.player_pos.x,y=min(game.player_pos.y+1,game.screen_size.y-1))
				elseif i == 67
					game.player_pos = (x=min(game.player_pos.x+1,game.screen_size.x-1),y=game.player_pos.y)
				elseif i == 68
					game.player_pos = (x=max(game.player_pos.x-1,0),y=game.player_pos.y)
				end
				Game_render(game)
			end
		end
	end
	endwin()
end

game = Game_new()
Game_run(game)

max_stars_repo_name: richardanaya/julia-roguelike

max_issues_repo_name: richardanaya/julia-roguelike

max_forks_repo_name: richardanaya/julia-roguelike


text:
@testset "Chernoff tests" begin
    d = Chernoff()

    cdftest = [
        0.005 0.5037916689930134;
        0.015 0.5113737159744436;
        0.025 0.5189518911337264;
        0.035 0.5265236172270374;
        0.045 0.5340863229612878;
        0.055 0.541637445363313;
        0.065 0.5491744321382676;
        0.075 0.556694744013348;
        0.085 0.5641958570651835;
        0.095 0.5716752650256268;
        0.105 0.5791304815656103;
        0.115 0.5865590425511593;
        0.125 0.5939585082711147;
        0.135 0.6013264656321359;
        0.145 0.6086605303182655;
        0.155 0.6159583489126876;
        0.165 0.6232176009794599;
        0.175 0.6304360010992436;
        0.185 0.6376113008624671;
        0.195 0.6447412908097765;
        0.205 0.6518238023240649;
        0.215 0.6588567094682105;
        0.225 0.6658379307650819;
        0.235 0.6727654309193446;
        0.245 0.6796372224828434;
        0.255 0.6864513674498119;
        0.265 0.6932059787930698;
        0.275 0.6998992219329853;
        0.285 0.7065293161345655;
        0.295 0.7130945358373242;
        0.305 0.7195932119132353;
        0.315 0.7260237328492151;
        0.325 0.7323845458573087;
        0.335 0.7386741579064688;
        0.345 0.7448911366812273;
        0.355 0.7510341114556635;
        0.365 0.7571017738953651;
        0.375 0.7630928787745712;
        0.385 0.7690062446148507;
        0.395 0.7748407542420853;
        0.405 0.7805953552614361;
        0.415 0.786269060454087;
        0.425 0.7918609480850309;
        0.435 0.7973701621351792;
        0.445 0.8027959124488315;
        0.455 0.8081374748004513;
        0.465 0.8133941908816786;
        0.475 0.8185654682017057;
        0.485 0.823650779913423;
        0.495 0.8286496645651691;
        0.505 0.8335617257530529;
        0.515 0.8383866317232886;
        0.525 0.8431241148805353;
        0.535 0.8477739712206347;
        0.545 0.8523360597092733;
        0.555 0.8568103015557271;
        0.565 0.8611966794490666;
        0.575 0.8654952367057263;
        0.585 0.8697060763547471;
        0.595 0.8738293601627617;
        0.605 0.8778653075886909;
        0.615 0.8818141946885997;
        0.625 0.88567635294178;
        0.635 0.8894521680505922;
        0.645 0.8931420786562767;
        0.655 0.8967465750281972;
        0.665 0.9002661976867428;
        0.675 0.9037015359952119;
        0.685 0.9070532266935508;
        0.695 0.9103219523993228;
        0.705 0.9135084400870754;
        0.715 0.9166134594908737;
        0.725 0.9196378215191889;
        0.735 0.9225823766240466;
        0.745 0.9254480131350427;
        0.755 0.928235655591972;
        0.765 0.9309462630294286;
        0.775 0.9335808272718444;
        0.785 0.9361403712094327;
        0.795 0.9386259470397855;
        0.805 0.9410386345416801;
        0.815 0.9433795393000659;
        0.825 0.9456497909787722;
        0.835 0.9478505415455782;
        0.845 0.9499829635303033;
        0.855 0.9520482483056357;
        0.865 0.9540476043275662;
        0.875 0.9559822554235589;
        0.885 0.9578534391151732;
        0.895 0.9596624048960603;
        0.905 0.9614104126081986;
        0.915 0.9630987307693821;
        0.925 0.9647286349803182;
        0.935 0.9663014063237074;
        0.945 0.9678183298266878;
        0.955 0.9692806929171217;
        0.965 0.9706897839442293;
        0.975 0.9720468907261572;
        0.985 0.9733532991278337;
        0.995 0.9746102916986603
    ]

    pdftest = [
        0.000 0.758345;
        0.050 0.755123;
        0.100 0.745532;
        0.150 0.729792;
        0.200 0.708260;
        0.250 0.681422;
        0.300 0.649874;
        0.350 0.614303;
        0.400 0.575464;
        0.450 0.534159;
        0.500 0.491208;
        0.550 0.447424;
        0.600 0.403594;
        0.650 0.360447;
        0.700 0.318646;
        0.750 0.278760;
        0.800 0.241264;
        0.850 0.206523;
        0.900 0.174795;
        0.950 0.146231;
        1.000 0.120880;
        1.050 0.0987031;
        1.100 0.0795821;
        1.150 0.0633362;
        1.200 0.0497372;
        1.250 0.0385246;
        1.300 0.0294208;
        1.350 0.0221441;
        1.400 0.0164201;
        1.450 0.0119902;
        1.500 0.0086186;
        1.550 0.0060955;
        1.600 0.0042401;
        1.650 0.0028995;
        1.700 0.0019485;
        1.750 0.0012861;
        1.800 0.0008334;
        1.850 0.0005300;
        1.900 0.0003306;
        1.950 0.0002022;
        2.000 0.0001212
    ]

    for i=1:size(cdftest, 1)
        @test isapprox(cdf(d, cdftest[i, 1]), cdftest[i, 2])
    end

    for i=1:size(pdftest, 1)
        @test isapprox(pdf(d, pdftest[i, 1]), pdftest[i, 2] ; atol = 1e-6)
    end
end


text:
module LRUCache

include("cyclicorderedset.jl")
export LRU

using Base.Threads
using Base: Callable

_constone(x) = 1

# Default cache size
mutable struct LRU{K,V} <: AbstractDict{K,V}
    dict::Dict{K, Tuple{V, LinkedNode{K}, Int}}
    keyset::CyclicOrderedSet{K}
    currentsize::Int
    maxsize::Int
    lock::SpinLock
    by::Any
    finalizer::Any

    function LRU{K, V}(; maxsize::Int, by = _constone, finalizer = nothing) where {K, V}
        dict = Dict{K, V}()
        keyset = CyclicOrderedSet{K}()
        new{K, V}(dict, keyset , 0, maxsize, SpinLock(), by, finalizer)
    end
end

LRU(; maxsize::Int, by = _constone, finalizer = nothing) =
    LRU{Any,Any}(maxsize = maxsize, by = by, finalizer = finalizer)

Base.show(io::IO, lru::LRU{K, V}) where {K, V} =
    print(io, "LRU{$K, $V}(; maxsize = $(lru.maxsize))")

function Base.iterate(lru::LRU, state...)
    next = iterate(lru.keyset, state...)
    if next === nothing
        return nothing
    else
        k, state = next
        v, = lru.dict[k]
        return k=>v, state
    end
end

Base.length(lru::LRU) = length(lru.keyset)
Base.isempty(lru::LRU) = isempty(lru.keyset)
function Base.sizehint!(lru::LRU, n::Integer)
    lock(lru.lock) do
        sizehint!(lru.dict, n)
    end
    return lru
end

_unsafe_haskey(lru::LRU, key) = haskey(lru.dict, key)
function Base.haskey(lru::LRU, key)
    lock(lru.lock) do
        return _unsafe_haskey(lru, key)
    end
end
function Base.get(lru::LRU, key, default)
    lock(lru.lock) do
        if _unsafe_haskey(lru, key)
            v = _unsafe_getindex(lru, key)
            return v
        else
            return default
        end
    end
end
function Base.get(default::Callable, lru::LRU, key)
    lock(lru.lock)
    if _unsafe_haskey(lru, key)
        v = _unsafe_getindex(lru, key)
        unlock(lru.lock)
        return v
    else
        unlock(lru.lock)
    end
    return default()
end
function Base.get!(lru::LRU{K, V}, key, default) where {K, V}
    evictions = Tuple{K, V}[]
    v = lock(lru.lock) do
        if _unsafe_haskey(lru, key)
            v = _unsafe_getindex(lru, key)
            return v
        end
        v = default
        _unsafe_addindex!(lru, v, key)
        _unsafe_resize!(lru, evictions)
        return v
    end
    _finalize_evictions!(lru.finalizer, evictions)
    return v
end
function Base.get!(default::Callable, lru::LRU{K, V}, key) where {K, V}
    evictions = Tuple{K, V}[]
    lock(lru.lock)
    if _unsafe_haskey(lru, key)
        v = _unsafe_getindex(lru, key)
        unlock(lru.lock)
        return v
    else
        unlock(lru.lock)
    end
    v = default()
    lock(lru.lock)
    if _unsafe_haskey(lru, key)
        # should we test that this yields the same result as default()
        v = _unsafe_getindex(lru, key)
    else
        _unsafe_addindex!(lru, v, key)
        _unsafe_resize!(lru, evictions)
    end
    unlock(lru.lock)
    _finalize_evictions!(lru.finalizer, evictions)
    return v
end

function _unsafe_getindex(lru::LRU, key)
    v, n, s = lru.dict[key]
    _move_to_front!(lru.keyset, n)
    return v
end
function Base.getindex(lru::LRU, key)
    lock(lru.lock) do
        if _unsafe_haskey(lru, key)
            v = _unsafe_getindex(lru, key)
            return v
        else
            throw(KeyError(key))
        end
    end
end
function _unsafe_addindex!(lru::LRU{K}, v, key) where K
    n = LinkedNode{K}(key)
    rotate!(_push!(lru.keyset, n))
    s = lru.by(v)::Int
    lru.currentsize += s
    lru.dict[key] = (v, n, s)
end

function Base.setindex!(lru::LRU{K, V}, v, key; dofinalize = true) where {K, V}
    evictions = Tuple{K, V}[]
    lock(lru.lock) do
        if _unsafe_haskey(lru, key)
            old_v, n, s = lru.dict[key]
            if dofinalize && lru.finalizer !== nothing
                push!(evictions, (key, old_v))
            end
            lru.currentsize -= s
            s = lru.by(v)::Int
            lru.currentsize += s
            lru.dict[key] = (v, n, s)
            _move_to_front!(lru.keyset, n)
        else
            _unsafe_addindex!(lru, v, key)
        end
        _unsafe_resize!(lru, evictions)
    end
    _finalize_evictions!(lru.finalizer, evictions)
    return lru
end

function _unsafe_resize!(lru::LRU{K, V}, evictions::Vector{Tuple{K, V}},
                         maxsize::Integer = lru.maxsize) where {K, V}
    lru.maxsize = maxsize
    while lru.currentsize > lru.maxsize
        key = pop!(lru.keyset)
        v, n, s = pop!(lru.dict, key)
        if lru.finalizer !== nothing
            push!(evictions, (key, v))
        end
        lru.currentsize -= s
    end
    return
end
function Base.resize!(lru::LRU{K, V}; maxsize::Integer = lru.maxsize) where {K, V}
    @assert 0 <= maxsize
    evictions = Tuple{K, V}[]
    lock(lru.lock) do
        _unsafe_resize!(lru, evictions, maxsize)
    end
    _finalize_evictions!(lru.finalizer, evictions)
    return lru
end

function Base.delete!(lru::LRU{K, V}, key; dofinalize = true) where {K, V}
    v = lock(lru.lock) do
        v, n, s = pop!(lru.dict, key)
        lru.currentsize -= s
        _delete!(lru.keyset, n)
        return v
    end
    if dofinalize && lru.finalizer !== nothing
        lru.finalizer(key, v)
    end
    return lru
end
function Base.pop!(lru::LRU{K, V}, key; dofinalize = true) where {K, V}
    (key, v) = lock(lru.lock) do
        v, n, s = pop!(lru.dict, key)
        lru.currentsize -= s
        _delete!(lru.keyset, n)
        return (key, v)
    end
    if dofinalize && lru.finalizer !== nothing
        lru.finalizer(key, v)
    end
    return v
end

function Base.empty!(lru::LRU{K, V}; dofinalize = true) where {K, V}
    evictions = Tuple{K, V}[]
    lock(lru.lock) do
        if !dofinalize || lru.finalizer === nothing
            lru.currentsize = 0
            empty!(lru.dict)
            empty!(lru.keyset)
        else
            sizehint!(evictions, length(lru))
            _unsafe_resize!(lru, evictions, 0)
        end
    end
    _finalize_evictions!(lru.finalizer, evictions)
    return lru
end

function _finalize_evictions!(finalizer, evictions)
    for (key, value) in evictions
        finalizer(key, value)
    end
    return
end

end # module


text:
using Test

include("../polynomial.jl")

@testset "Polynomial result" begin
  @testset "nthPower(x, power)" begin
    @test nthPower(1,1000) == 1
    @test nthPower(2,3) == 8
    @test nthPower(10,3) == 1000
    @test_throws DomainError(-3, "argument must be nonnegative") nthPower(10,-3)
    @test nthPower(10,0) == 1
  end

  @testset "calcTerm(coeffient, (value, power))" begin
    @test calcTerm(0, (2, 34)) == 0
    @test calcTerm(-1, (2, 2)) == -4
    @test calcTerm(3, (1, 200)) == 3
  end

  @testset "calcPolynomial((coeffient, value, power)::firstTermTuple, ... , (coeffient, value, power)::lastTermTuple)" begin
    # y = 4 + x^3 -2*x^2 + 3*x when x = 1
    @test calcPolynomial((4,1,0), (1,1,3), (-2,1,2), (3,1,1)) == 6
    @test calcPolynomial((8,1,0)) == 8
  end

end