text:
"""
    mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)

returns an efficient implementation of:
` sum(A[i,j] * x[i] * x[j] * op(p[i],p[j]) * op_asym(x[i],x[j],A_asym[i,j])) for i = 1:n , j = 1:n)`
where `op(p[i],p[j]) == op(p[j],p[i])` , op_asym doesn't follow this symmetry.

"""
function mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)
    #@show length(x)
    N = length(x)
    checkbounds(A, N, N)
    checkbounds(A_asym, N, N)
    @boundscheck checkbounds(p, N)
    @inbounds begin
        res1 = zero(eltype(x))

        for i = 1:N
            x[i] != 0 && begin
                res1 += p[i] * x[i]^2
                for j = 1:i - 1
                    res1 +=
                        2 *
                        x[i] *
                        x[j] *
                        op(p[i], p[j]) *
                        A[i, j] *
                        op_asym(x[i], x[j], A_asym[i, j])
                end
            end
        end
    end
    return res1
end

_gerg_asymetric_mix_rule(xi, xj, b) = b * (xi + xj) / (xi * b^2 + xj)

function T_scale(model::GERG2008,x=SA[1.])
    return mixing_rule_asymetric(
        (a,b)->sqrt(a * b),
        _gerg_asymetric_mix_rule,
        x,
        model.Tc,
        model.gamma_T,
        model.beta_T,
    )
end

# function T_scale(model::GERG2008)
#     return only(model.Tc)
# end

function T_scales(model::GERG2008,x=SA[1.])
    return model.Tc
end

# function p_scale(model::GERG2008)
#     return only(model.Pc)
# end

function p_scale(model::GERG2008,x=SA[1.])
    return dot(x,model.Pc)
end

# function p_scales(model::GERG2008,x=SA[1.])
#     return only(model.Pc)
# end

function _v_scale(model::GERG2008,x=SA[1.])
    return mixing_rule_asymetric(
        (a,b) -> ((cbrt(a) + cbrt(b))*0.5)^3,
        _gerg_asymetric_mix_rule,
        x,
        model.vc,
        model.gamma_v,
        model.beta_v,
    )
end

# function _v_scale(model::GERG2008)
#     return only(model.vc)
# end

# function lb_volume(model::GERG2008)
#     return only(model.lb_v)
# end

function lb_volume(model::GERG2008,x=SA[1.];phase=:unknown)
    return dot(x,model.lb_v)
end

function _delta(model::GERG2008, rho, T, x=SA[1.])
    vcmix = _v_scale(model,x)
    return rho * vcmix
end

function _tau(model::GERG2008, rho, T, x=SA[1.])
    Tcmix  = T_scale(model,x)
    return Tcmix / T
end

# function _delta(model::GERG2008, rho, T)
#     vcmix = _v_scale(model)
#     return rho * vcmix
# end

# function _tau(model::GERG2008, rho, T)
#     Tcmix  = T_scale(model)
#     return Tcmix / T
# end

function _f0(model::GERG2008, ρ, T, x)
    RR = 8.314472 / 8.314510
    common_type = promote_type(typeof(ρ), typeof(T), eltype(x))
    _0 = zero(common_type)
    res = _0
    x0 = _0  #for comparison
    ao2 = _0
    ao_zero = _0

    ρc = model.rhoc
    Tc = model.Tc
    nr = model.nr
    zeta = model.zeta
    for i in model.ideal_iters[1]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
            nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
            nr[5, i] * log(cosh(zeta[2, i] * τ)) +
            nr[6, i] * log(abs(sinh(zeta[3, i] * τ))) -
            nr[7, i] * log(cosh(zeta[4, i] * τ))

            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[2]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
                nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
                nr[5, i] * log(cosh(zeta[2, i] * τ)) +
                nr[6, i] * log(abs(sinh(zeta[3, i] * τ)))
            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[3]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
                nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
                nr[5, i] * log(cosh(zeta[2, i] * τ))
            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[4]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao3 = log(δ)
            a0 = RR * (ao1) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    return res
end

function _f0(model::GERG2008, rho, T)

    RR = 8.314472 / 8.314510
    common_type = promote_type(typeof(rho), typeof(T))

    res = zero(common_type)
    ao1 = zero(common_type)
    ao2 = zero(common_type)
    ao_zero = zero(common_type)

    delta = rho / model.rhoc[1]
    tau = model.Tc[1] / T

    if !iszero(length(model.ideal_iters[1]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau)) +
            model.nr[6, 1] * log(abs(sinh(model.zeta[3, 1] * tau))) -
            model.nr[7, 1] * log(cosh(model.zeta[4, 1] * tau))

        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[2]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau)) +
            model.nr[6, 1] * log(abs(sinh(model.zeta[3, 1] * tau)))
        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[3]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau))
        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[4]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao3 = log(delta)
        a0 = RR * (ao1) + ao3
        res +=  a0
    end

    return res
end

function _fr1(model::GERG2008, delta, tau)
    common_type = promote_type(typeof(delta), typeof(tau))
    res = zero(common_type)
    res0 = zero(common_type)
    res1 = zero(common_type)

    res1 = res0
    for k = 1:model.k_pol_ik[1]
        res1 += model.n0ik[1][k] * (delta^model.d0ik[1][k]) * (tau^model.t0ik[1][k])
    end

    for k = (model.k_pol_ik[1]+1):(model.k_exp_ik[1]+model.k_pol_ik[1])
        res1 +=
            model.n0ik[1][k] *
            (delta^model.d0ik[1][k]) *
            (tau^model.t0ik[1][k]) *
            exp(-delta^model.c0ik[1][k-model.k_pol_ik[1]])
    end
    res += res1

    return res
end

function _fr1(model::GERG2008, delta, tau,x)
    common_type = promote_type(typeof(delta), typeof(tau), eltype(x))
    res = zero(common_type)
    res0 = zero(common_type)
    res1 = zero(common_type)
    x0 = zero(common_type)

    for i = 1:length(model.components)
        x[i] != x0 && begin
            res1 = res0
            for k = 1:model.k_pol_ik[i]
                res1 += model.n0ik[i][k] * (delta^model.d0ik[i][k]) * (tau^model.t0ik[i][k])
            end

            for k = (model.k_pol_ik[i]+1):(model.k_exp_ik[i]+model.k_pol_ik[i])
                res1 +=
                    model.n0ik[i][k] *
                    (delta^model.d0ik[i][k]) *
                    (tau^model.t0ik[i][k]) *
                    exp(-delta^model.c0ik[i][k-model.k_pol_ik[i]])
            end
            res += x[i] * res1

        end
    end

    return res
end

function _fr2(model::GERG2008, delta, tau, x)
    common_type = promote_type(typeof(delta), typeof(tau), eltype(x))
    res = zero(common_type)
    res0 = zero(common_type)
    x0 = zero(common_type)

    for kk = 1:length(model.Aij_indices)
        i1, i2, i0 = model.Aij_indices[kk]

        x[i1] != x0 && x[i2] != x0 && begin
            res1 = res0
            for j = 1:model.k_pol_ijk[i0]
                res1 +=
                    model.nijk[i0][j] * (delta^model.dijk[i0][j]) * (tau^model.tijk[i0][j])
            end

            for j = (model.k_pol_ijk[i0]+1):(model.k_pol_ijk[i0]+model.k_exp_ijk[i0])

                idx = j - model.k_pol_ijk[i0]

                res1 +=
                    model.nijk[i0][j] *
                    (delta^model.dijk[i0][j]) *
                    (tau^model.tijk[i0][j]) *
                    exp(
                        -model.etaijk[i0][idx] * (delta - model.epsijk[i0][idx])^2 -
                        model.betaijk[i0][idx] * (delta - model.gammaijk[i0][idx]),
                    )
            end
            res += res1 * x[i1] * x[i2] * model.fij[i0]

        end
    end
    return res
end

function a_ideal(model::GERG2008, v, T, z=SA[1.0])
    N = sum(z)
    len = length(z)
    x = z/N
    rho = 1.0e-3 / (v/N)
    if len == 1
        return _f0(model, rho, T)
    else
        return _f0(model, rho, T, x)
    end
end

function a_res(model::GERG2008, v, T, z=SA[1.0])
    N = sum(z)
    len = length(z)
    rho = 1.0e-3 / (v/N)
    x = z/N
    delta = _delta(model, rho, T, x)
    tau = _tau(model, rho, T, x)
    if len == 1
        return _fr1(model, delta, tau)
    else
        return (_fr1(model, delta, tau, x) + _fr2(model, delta, tau, x))
    end
end

function eos(model::GERG2008, v, T, z=SA[1.0])
    negative_vt(v,T) && return nan_num(v,T,z)
    N = sum(z)
    len = length(z)
    rho = 1.0e-3 / (v/N)
    x = z/N
    R = R̄
    delta = _delta(model, rho, T, x)
    tau = _tau(model, rho, T, x)
    if len == 1
        return N * R * T * (_f0(model, rho, T) + _fr1(model, delta, tau))
    else
        return N * R * T * (_f0(model, rho, T, x)+ _fr1(model, delta, tau, x) + _fr2(model, delta, tau, x))
    end
end

# function ρ_TP(species::EoSModel, T, P)
#     v = volume(species,P,T)
#     return 1/v
# end
#=
julia> @time [ρ_TP(thermo_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
#4.271689 seconds (2.22 M allocations: 130.980 MiB, 1.10% gc time, 33.34% compilation time

julia> @benchmark [ρ_TP(thermo_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];

BenchmarkTools.Trial:
  memory estimate:  14.73 MiB
  allocs estimate:  320005
  --------------
  minimum time:     2.903 s (0.00% GC)
  median time:      2.994 s (0.00% GC)
  mean time:        2.994 s (0.00% GC)
  maximum time:     3.085 s (0.00% GC)
  --------------
  samples:          2
  evals/sample:     1

julia> @time [ρ_TP(CP_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
  0.538447 seconds (115.41 k allocations: 4.366 MiB, 6.34% compilation time)

julia> @benchmark [ρ_TP(CP_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
BenchmarkTools.Trial:
  memory estimate:  859.66 KiB
  allocs estimate:  50005
  --------------
  minimum time:     690.530 ms (0.00% GC)  median time:      708.036 ms (0.00% GC)  mean time:        832.741 ms (0.00% GC)  maximum time:     1.171 s (0.00% GC)
  --------------
  samples:          6
  evals/sample:     1
  =#

#julia> @time [ρ_TP(w, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
#  2.781332 seconds (889.17 k allocations: 61.773 MiB, 0.73% gc time, 2.29% compilation time)

#julia> @benchmark [ρ_TP(Clapeyron_H2O, T, P) for T in LinRange(280, 330, 100), P in LinRange(8e4, 12e5, 100)];
#=
BenchmarkTools.Trial:
  memory estimate:  57.23 MiB
  allocs estimate:  805979
  --------------
  minimum time:     2.841 s (0.31% GC)
  median time:      2.843 s (0.30% GC)
  mean time:        2.843 s (0.30% GC)
  maximum time:     2.844 s (0.30% GC)
  --------------
  samples:          2
  evals/sample:     1
  =#


text:
include("decrement.jl")

const ixVector = Vector{Float64}

mutable struct InterestRate <: Decrement
    ix::Vector{Float64}
    ifx
    # constructor with a predefined vector
    function InterestRate(v::ixVector)
        new(v,(x -> v[x]))
    end

    # constructor with a real number argument converts to a function that will produce an interest rate
    function InterestRate(i::Real)
        new(Float64[],(x -> i))
    end

    # constructor with (any) argument assumes a function that will produce an interest rate
    function InterestRate(f)
        new(Float64[],f)
    end
    # constructor with (any) argument assumes a function that will produce an interest rate
    # and provided rates for the first n... time periods
    function InterestRate(f, x...)
        new(collect(x),f)
    end
end

# the interst during time x
function i(iv::InterestRate,x)
    if x < 0 # looking for nth item from end of vector
        return iv.ix[end + Int(x) + 1]
    elseif length(iv.ix) < x
        append!(iv.ix, iv.ifx(length(iv.ix)+1))
        return i(iv,x)
    else
        return iv.ix[Int(x)]
    end
end

# the discount rate from time x to x+t
function tvx(iv::InterestRate,t,x, v=1.0)
    if t > 0
        return tvx(iv,t-1,x+1,v / (1.0 + i(iv,x)))
    else
        return v
    end
end

# ω (omega), the ultimate end of the given table
function ω(i::InterestRate)
    if length(i.ix) == 0
        return Inf
    else
        return length(i.ix)
    end
end
w(iv::InterestRate) = ω(iv)

### Convienence functions

# the discount rate at time x
vx(iv::InterestRate,x) = tvx(iv,1.0,x,1.0)

# the discount rate at time 1
v(iv::InterestRate) = vx(iv,1.0)

ext: jl


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: TKEBasedVerticalDiffusivity

#####
##### Parameters
#####

"""
    ocean_simulation(; kwargs...)

Return an `Oceananigans.Simulation` of a column model initialized with
mutable array surface boundary conditions and a linear density stratification.

Arguments
=========

    * Nz: Number of vertical grid points
    * Lz: [m] Vertical extent of domain,
    * f:  [s-1] Coriolis parameter,
    * g:  [m s-2] Gravitational acceleration,
    * T₀: [ᵒC], sea surface temperature,
    * S₀: [psu], sea surface salinity,
    * α: Thermal expansion coefficient
    * β: Haline contraction coefficient
"""
function ocean_simulation(;
    Nz = 64,  # Number of vertical grid points
    Lz = 512, # Vertical extent of domain
    f = 1e-4, # Coriolis parameter
    g = 9.81, # Gravitational acceleration
    α = 2e-4, # Thermal expansion coefficient
    β = 8e-5, # Haline contraction coefficient
)

    grid = RegularRectilinearGrid(
        size = (1, 1, Nz),
        x = (0, 1),
        y = (0, 1),
        z = (-Lz, 0),
        topology = (Periodic, Periodic, Bounded),
    )

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition([0.0]))
    v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition([0.0]))
    T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition([0.0]))
    S_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition([0.0]))

    eos = LinearEquationOfState(α = α, β = β)

    model = HydrostaticFreeSurfaceModel(
        grid = grid,
        tracers = (:T, :S, :e),
        free_surface = ImplicitFreeSurface(gravitational_acceleration = g),
        buoyancy = SeawaterBuoyancy(gravitational_acceleration = g, equation_of_state = eos),
        coriolis = FPlane(f = f),
        boundary_conditions = (T = T_bcs, S = S_bcs, u = u_bcs, v = v_bcs),
        closure = TKEBasedVerticalDiffusivity(),
    )

    simulation = Oceananigans.Simulation(model, Δt = 0.02, stop_iteration = 1)

    # Initialize the ocean state with a linear temperature and salinity stratification
    α = simulation.model.buoyancy.model.equation_of_state.α
    β = simulation.model.buoyancy.model.equation_of_state.β
    Tᵢ(x, y, z) = 16 + α * g * 5e-5 * z
    Sᵢ(x, y, z) = 35 - β * g * 5e-5 * z
    set!(simulation.model, T = Tᵢ, S = Sᵢ)

    # collect data (needs optimising)
    ocean_data = []
    data = (T = deepcopy(simulation.model.tracers.T), time = simulation.model.clock.time)
    push!(ocean_data, data)

    return simulation, ocean_data
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using CMIM
using Base.Test

# Testing data
X = Bool[true true true false false false true; true false true false false false true; false false false true false false true; false true true false true false true; true false false true false false true; false true true true true false true; true false false false false false true; true true true true true false false; false true false true false false false; true true true true false false false; true false false false true false true; false false false false false true true; true true true false false true true; true true false true false true false; false false true true false true true; false true false false true false true; true false false false false true true]
y = Bool[true, true, true, false, true, false, true, true, false, false, false, true, true, true, false, true, true]

# Assign a score to all features
features = cmim(size(X, 2), X, y)

# Compare the results against R implementation of CMIM algorithm
# https://www.rdocumentation.org/packages/praznik/versions/5.0.0/topics/CMIM
@test first.(features) == [1, 3, 5, 4, 2, 6, 7]
@test ≈(last.(features), [0.0736965, 0.0426288, 0.041282, 0.0395759, 0.00908506, 0.00341897, 0.0], atol=0.00001)


text:
"""
```
compute_meansbands(m, input_type, cond_type, output_vars; forecast_string = "",
    verbose = :low, kwargs...)

compute_meansbands(m, input_type, cond_type, output_var, df; forecast_string = "",
    population_data = DataFrame(), population_forecast = DataFrame(),
    verbose = :none, kwargs...)

compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
    forecast_string = "", population_data = DataFrame(),
    population_forecast = DataFrame(), verbose = :low,
    kwargs...)
```

Compute means and bands for pseudo-observables, observables, and shocks, and write
the results to a file. Other methods are for one `output_var` and one `var_name` respectively.

### Keyword Arguments

- `forecast_string::String = ""`: forecast identifier string (the value
  \"fcid=value\" in the forecast output filename). Required when
  `input_type == :subset`

- `density_bands::Vector{Float64} = [.5, .6, .7, .8 .9]`: a vector of percent values (between 0 and 1) for
  which to compute density bands.

- `minimize::Bool = false`: if `true`, choose shortest interval, otherwise just chop off
  lowest and highest (percent/2)

- `verbose::Symbol = :low`: level of error messages to be printed to screen. One of `:none`,
  `:low`, `:high`

- `bdd_fcast::Bool = true`: if true, calculate bounded forecasts

- `skipnan::Bool = false`: if true, remove any NaNs found in the raw forecast output series

- `df::DataFrame = DataFrame()`: if an empty DataFrame, then the function will attempt to
    load the data using `load_data`.

- `check_empty_columns::Bool = true`: if true, throw an error if
    calling load_data yields an empty column.

- `pseudo2data::AbstractDict{Symbol, Symbol} = Dict()`: Maps the names of pseudo-observables
    which require historical data during transformation to the name of the required historical series.
"""
function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol,
                            cond_type::Symbol, output_vars::Vector{Symbol};
                            forecast_string::String = "",
                            verbose::Symbol = :low, df::DataFrame = DataFrame(),
                            check_empty_columns::Bool = true,
                            bdd_fcast::Bool = true, skipnan::Bool = false,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}(),
                            kwargs...)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        output_dir = workpath(m, "forecast")
        println()
        @Base.info "Computing means and bands for input_type = $input_type, cond_type = $cond_type..."
        println("Start time: $(now())")
        println("Means and bands will be saved in $output_dir")
    end
    elapsed_time = @elapsed let
        # Determine full set of output_vars necessary for plotting desired result
        output_vars = add_requisite_output_vars(output_vars; bdd_fcast = bdd_fcast)
        if input_type == :prior
            output_vars = setdiff(output_vars, [:bddforecastobs])
        end
        # Load population data and main dataset (required for some transformations)
        if all(var -> get_product(var) == :irf, output_vars)
            population_data, population_forecast = DataFrame(), DataFrame()
        else
            population_data, population_forecast = load_population_growth(m, verbose = verbose)
            isempty(df) && (df = load_data(m, cond_type = cond_type,
                                           check_empty_columns = check_empty_columns, verbose = :none))
        end
        for output_var in output_vars
            prod = get_product(output_var)
            if VERBOSITY[verbose] >= VERBOSITY[:high]
                if prod in [:shockdec, :irf]
                    println("Computing " * string(output_var) * " for shocks:")
                else
                    print("Computing " * string(output_var) * "... ")
                end
            end

            # Compute means and bands
            mb = compute_meansbands(m, input_type, cond_type, output_var, df;
                                    forecast_string = forecast_string,
                                    population_data = population_data,
                                    population_forecast = population_forecast,
                                    skipnan = skipnan, pseudo2data = pseudo2data,
                                    verbose = verbose,
                                    kwargs...)
            GC.gc()
        end
    end
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        total_mb_time     = elapsed_time
        total_mb_time_min = total_mb_time/60

        println("\nTotal time to compute means and bands: " * string(total_mb_time_min) * " minutes")
        println("Computation of means and bands complete: " * string(now()))
    end
end

function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                            output_var::Symbol, df::DataFrame;
                            forecast_string::String = "",
                            population_data::DataFrame = DataFrame(),
                            population_forecast::DataFrame = DataFrame(),
                            skipnan::Bool = false,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}(),
                            verbose::Symbol = :none,
                            kwargs...)

    # Determine class and product
    class   = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast metadata
    metadata = get_mb_metadata(m, input_type, cond_type, output_var; forecast_string = forecast_string)

    date_list      = product == :irf ? Date[] : collect(keys(metadata[:date_inds]))
    variable_names = collect(keys(metadata[:indices]))
    pop_growth     = get_mb_population_series(product, population_data, population_forecast, date_list)

    # Compute means and bands
    if product in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                   :bddforecast, :bddforecastut, :bddforecast4q, :dettrend, :trend]
                   # :forecastlvl, :histlvl, :bddforecastlvl]
        # Get to work!
        # pmap produces an error for trendobs sometimes, so just doing this iteratively
        if product == :trend
            mb_vec = Vector{Any}(undef,length(variable_names))
            for i in 1:length(mb_vec)
                mb_vec[i] = compute_meansbands(m, input_type, cond_type, output_var,
                                               variable_names[i], df; pop_growth = pop_growth,
                                               forecast_string = forecast_string,
                                               skipnan = skipnan, pseudo2data = pseudo2data,
                                               kwargs...)
            end
        else
            mb_vec = pmap(var_name -> compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
                                                         pop_growth = pop_growth, forecast_string = forecast_string,
                                                         pseudo2data = pseudo2data,
                                                         skipnan = skipnan, kwargs...),
                          variable_names)
        end

        # Re-assemble pmap outputs
        means = DataFrame(date = date_list)
        bands = Dict{Symbol,DataFrame}()

        for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
            means[!,var_name] = var_means
            bands[var_name] = var_bands
            bands[var_name][!,:date] = date_list
        end

    elseif product in [:shockdec, :irf]
        means = product == :irf ? DataFrame() : DataFrame(date = date_list)
        bands = Dict{Symbol, DataFrame}()

        # Get to work!
        for shock_name in keys(metadata[:shock_indices])
            println(verbose, :high, "  * " * string(shock_name))

            mb_vec = pmap(var_name -> compute_meansbands(m, input_type, cond_type, output_var, var_name, df;
                                          pop_growth = pop_growth, shock_name = Nullables.Nullable(shock_name),
                                                         forecast_string = forecast_string,
                                                         skipnan = skipnan, pseudo2data = pseudo2data, kwargs...),
                          variable_names)

            # Re-assemble pmap outputs
            for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
                means[!, Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_means
                bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_bands
                if product != :irf
                    bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)][!, :date] = date_list
                end
            end
        end

    else
        error("Invalid product: $product")
    end

    mb = MeansBands(metadata, means, bands)

    # Write to file
    filepath = get_meansbands_output_file(m, input_type, cond_type, output_var,
                                          forecast_string = forecast_string)
    dirpath = dirname(filepath)
    isdir(dirpath) || mkpath(dirpath)
    JLD2.jldopen(filepath, true, true, true, IOStream) do file
        write(file, "mb", mb)
    end

    sep = prod in [:shockdec, :irf] ? "  " : ""
    println(verbose, :high, sep * "wrote " * basename(filepath))

    return mb
end

function compute_meansbands(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                            output_var::Symbol, var_name::Symbol, df::DataFrame;
                            forecast_string::String = "",
                            pop_growth::AbstractVector{Float64} = Float64[],
                            skipnan::Bool = false,
                            shock_name::Nullable{Symbol} = Nullables.Nullable{Symbol}(),
                            density_bands::Vector{Float64} = [0.5,0.6,0.7,0.8,0.9],
                            minimize::Bool = false,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}(),
                            compute_shockdec_bands::Bool = false)

    # Return only one set of bands if we read in only one draw
    if input_type in [:init, :mode, :mean]
        density_bands = [.5]
    end

    # Determine class and product
    class = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast draws
    fcast_series, transform = read_forecast_output(m, input_type, cond_type,
                                                   output_var, var_name, shock_name,
                                                   forecast_string = forecast_string)

    # Reverse transform
    y0_index = get_y0_index(m, product) # this should be index_forecast_start(m) - 4, so 4 quarters before forecast
    data = if haskey(pseudo2data, var_name) && product != :irf
        Float64.(collect(Missings.replace(Vector{Union{Missing, Float64}}(df[!,pseudo2data[var_name]]), NaN)))
    elseif class == :obs && product != :irf
        Float64.(collect(Missings.replace(Vector{Union{Missing, Float64}}(df[!,var_name]), NaN)))
    else
        fill(NaN, size(df, 1))
    end

    # data = class == :obs && product != :irf ? Float64.(collect(Missings.replace(df[:,var_name], NaN))) : fill(NaN, size(df, 1))
    transformed_series = mb_reverse_transform(fcast_series, transform, product, class,
                                              y0_index = y0_index, data = data,
                                              pop_growth = pop_growth, use_data = haskey(pseudo2data, var_name))

    # Handle NaNs
    if skipnan && output_var != :histobs && any(isnan.(transformed_series))
        # Remove rows with NaNs
        nanrows = vec(mapslices(x -> any(isnan.(x)), transformed_series, dims = Int[2]))
        transformed_series = transformed_series[.!nanrows, :]
    end

    # Compute means and bands
    means = vec(mean(transformed_series, dims = 1))
    bands = if product in [:shockdec, :dettrend, :trend] && !compute_shockdec_bands
        Dict{Symbol,DataFrame}()
    else
        find_density_bands(transformed_series, density_bands, minimize = minimize)
    end
    return means, bands
end

"""
```
compute_meansbands(models, input_types, cond_types, output_vars; forecast_strings = [],
    verbose = :low, kwargs...)

compute_meansbands(models, input_types, cond_types, output_vars, df; forecast_strings = [],
    population_data = DataFrame(), population_forecast = DataFrame(),
    verbose = :none, kwargs...)

compute_meansbands(models, input_types, cond_types, output_var, var_name, df;
    forecast_strings = [], population_data = DataFrame(),
    population_forecast = DataFrame(), verbose = :low,
    kwargs...)
```

Compute means and bands for pseudo-observables, observables, and shocks from multiple forecasts, and write
the results to a file. Other methods are for one `output_var` and one `var_name` respectively.

### Keyword Arguments

- `forecast_strings::Vector{String} = []`: forecast identifier strings for each forecast (the value
  \"fcid=value\" in the forecast output filename for individual forecasts). Required when
  `input_type == :subset`

- `density_bands::Vector{Float64} = [.5, .6, .7, .8 .9]`: a vector of percent values (between 0 and 1) for
  which to compute density bands.

- `minimize::Bool = false`: if `true`, choose shortest interval, otherwise just chop off
  lowest and highest (percent/2)

- `verbose::Symbol = :low`: level of error messages to be printed to screen. One of `:none`,
  `:low`, `:high`

- `bdd_fcast::Bool = true`: if true, calculate bounded forecasts

- `skipnan::Bool = false`: if true, remove any NaNs found in the raw forecast output series

- `ndraws::Int = 20000`: number of draws from the underlying scenarios

- `df::DataFrame = DataFrame()`: if an empty DataFrame, then the function will attempt to
    load the data using `load_data`.

- `check_empty_columns::Bool = true`: if true, throw an error if
    calling load_data yields an empty column.

- `pseudo2data::AbstractDict{Symbol, Symbol} = Dict()`: Maps the names of pseudo-observables
    which require historical data during transformation to the name of the required historical series.
"""
function compute_meansbands(models::Vector{<: AbstractDSGEModel},
                            input_types::Vector{Symbol},
                            cond_types::Vector{Symbol},
                            output_vars::Vector{Symbol};
                            forecast_strings::Vector{String} = Vector{String}("", 0),
                            weights::Vector{Float64} = Vector{Float64}(undef, 0),
                            combo_forecast_string::String = join(forecast_strings) * join(map(x->string(x), weights), "_"),
                            df::DataFrame = DataFrame(),
                            check_empty_columns::Bool = true,
                            variable_names::Vector{Symbol} = Vector{Symbol}(undef, 0),
                            shock_names::Vector{Symbol} = Vector{Symbol}(undef, 0),
                            skipnan::Bool = false, ndraws::Int = 20000,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}(),
                            verbose::Symbol = :low,
                            kwargs...)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        output_dir = workpath(models[1], "forecast")
        println()
        @Base.info "Computing means and bands for input_types = $input_types, cond_types = $cond_types..."
        println("Start time: $(now())")
        println("Means and bands will be saved in $output_dir")
    end
    elapsed_time = @elapsed let
        # Determine full set of output_vars necessary for plotting desired result
        output_vars = add_requisite_output_vars(output_vars)
        if input_types[1] == :prior
            output_vars = setdiff(output_vars, [:bddforecastobs])
        end
        # Load population data and main dataset (required for some transformations)
        if all(var -> get_product(var) == :irf, output_vars)
            population_data, population_forecast = DataFrame(), DataFrame()
        else
            population_data, population_forecast = load_population_growth(models[1], verbose = verbose)
            isempty(df) && (df = load_data(models[1], check_empty_columns = check_empty_columns, verbose = :none))
        end
        for output_var in output_vars
            prod = get_product(output_var)
            if VERBOSITY[verbose] >= VERBOSITY[:high]
                if prod in [:shockdec, :irf]
                    println("Computing " * string(output_var) * " for shocks:")
                else
                    print("Computing " * string(output_var) * "... ")
                end
            end

            # Compute means and bands
            mb = compute_meansbands(models, input_types,
                                    cond_types, output_var, df;
                                    forecast_strings = forecast_strings,
                                    combo_forecast_string = combo_forecast_string,
                                    weights = weights,
                                    population_data = population_data,
                                    population_forecast = population_forecast,
                                    variable_names = variable_names,
                                    shock_names = shock_names,
                                    skipnan = skipnan, ndraws = ndraws,
                                    pseudo2data = pseudo2data,
                                    verbose = verbose,
                                    kwargs...)
            GC.gc()
        end
    end
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        total_mb_time     = elapsed_time
        total_mb_time_min = total_mb_time/60

        println("\nTotal time to compute means and bands: " * string(total_mb_time_min) * " minutes")
        println("Computation of means and bands complete: " * string(now()))
    end
end

function compute_meansbands(models::Vector{<: AbstractDSGEModel},
                            input_types::Vector{Symbol},
                            cond_types::Vector{Symbol},
                            output_var::Symbol, df::DataFrame;
                            forecast_strings::Vector{String} = Vector{String}("", 0),
                            combo_forecast_string::String = "",
                            weights::Vector{Float64} = Vector{Float64}(undef, 0),
                            population_data::DataFrame = DataFrame(),
                            population_forecast::DataFrame = DataFrame(),
                            variable_names::Vector{Symbol} = Vector{Symbol}(undef, 0),
                            shock_names::Vector{Symbol} = Vector{Symbol}(undef, 0),
                            skipnan::Bool = false, ndraws::Int = 20000,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}(),
                            verbose::Symbol = :none,
                            kwargs...)

    # Determine class and product
    class   = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast metadata
    metadata = get_mb_metadata(models[1], input_types[1], cond_types[1], output_var; forecast_string = forecast_strings[1])
    metadata[:forecast_string] = combo_forecast_string

    date_list      = product == :irf ? Date[] : collect(keys(metadata[:date_inds]))
    if isempty(variable_names)
        variable_names = collect(keys(metadata[:indices]))
    end
    pop_growth     = get_mb_population_series(product, population_data, population_forecast, date_list)

    # Compute means and bands
    if product in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                   :bddforecast, :bddforecastut, :bddforecast4q, :dettrend, :trend]
        # Get to work!
        # pmap produces an error for trendobs sometimes, so just doing this iteratively
        if product == :trend
            mb_vec = Vector{Any}(undef,length(variable_names))
            for i in 1:length(mb_vec)
                mb_vec[i] = compute_meansbands(models, input_types,
                                               cond_types, output_var,
                                               variable_names[i], df; pop_growth = pop_growth,
                                               forecast_strings = forecast_strings,
                                               weights = weights, skipnan = skipnan,
                                               pseudo2data = pseudo2data,
                                               ndraws = ndraws, kwargs...)
            end
        else
            mb_vec = pmap(var_name -> compute_meansbands(models, input_types, cond_types, output_var, var_name, df;
                                                         pop_growth = pop_growth,
                                                         forecast_strings = forecast_strings, weights = weights,
                                                         skipnan = skipnan, ndraws = ndraws,
                                                         pseudo2data = pseudo2data, kwargs...),
                          variable_names)
        end

        # Re-assemble pmap outputs
        means = DataFrame(date = date_list)
        bands = Dict{Symbol,DataFrame}()

        for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
            means[!,var_name] = vec(var_means)
            bands[var_name] = var_bands
            bands[var_name][!,:date] = date_list
        end

    elseif product in [:shockdec, :irf]
        means = product == :irf ? DataFrame() : DataFrame(date = date_list)
        bands = Dict{Symbol, DataFrame}()

        # Get to work!
        if isempty(shock_names)
            shock_names = keys(metadata[:shock_indices])
        end
        for shock_name in shock_names
            println(verbose, :high, "  * " * string(shock_name))

            mb_vec = pmap(var_name -> compute_meansbands(models, input_types,
                                                         cond_types, output_var,
                                                         var_name, df; pop_growth = pop_growth,
                                                         forecast_strings = forecast_strings,
                                                         weights = weights, shock_name = Nullables.Nullable(shock_name),
                                                         skipnan = skipnan, ndraws = ndraws,
                                                         pseudo2data = pseudo2data, kwargs...),
                             variable_names)

            # Re-assemble pmap outputs
            for (var_name, (var_means, var_bands)) in zip(variable_names, mb_vec)
                means[!, Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_means
                bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)] = var_bands
                if product != :irf
                    bands[Symbol(var_name, DSGE_SHOCKDEC_DELIM, shock_name)][!, :date] = date_list
                end
            end
        end

    else
        error("Invalid product: $product")
    end # of if product

    mb = MeansBands(metadata, means, bands)

    # Write to file
    filepath = get_meansbands_output_file(models[1], input_types[1], cond_types[1], output_var,
                                          forecast_string = combo_forecast_string)

    dirpath = dirname(filepath)
    isdir(dirpath) || mkpath(dirpath)
    JLD2.jldopen(filepath, true, true, true, IOStream) do file
        write(file, "mb", mb)
    end

    sep = prod in [:shockdec, :irf] ? "  " : ""
    println(verbose, :high, sep * "wrote " * basename(filepath))

    return mb
end

function compute_meansbands(models::Vector,
                            input_types::Vector{Symbol},
                            cond_types::Vector{Symbol},
                            output_var::Symbol, var_name::Symbol, df::DataFrame;
                            forecast_strings::Vector{String} = Vector{String}("", 0),
                            weights::Vector{Float64} = Vector{Float64}(undef, 0),
                            pop_growth::AbstractVector{Float64} = Float64[],
                            shock_name::Nullable{Symbol} = Nullables.Nullable{Symbol}(),
                            density_bands::Vector{Float64} = [0.5,0.6,0.7,0.8,0.9],
                            minimize::Bool = false, skipnan::Bool = false,
                            compute_shockdec_bands::Bool = false, ndraws::Int = 20000,
                            pseudo2data::AbstractDict{Symbol, Symbol} = Dict{Symbol, Symbol}())

    # Return only one set of bands if we read in only one draw
    if input_types[1] in [:init, :mode, :mean]
        density_bands = [.5]
    end

    # Determine class and product
    class = get_class(output_var)
    product = get_product(output_var)

    # Read in forecast draws
    fcast_seriess = Vector{Matrix{Float64}}(undef, length(models))
    transforms = Vector(undef, length(models))
    for i in 1:length(models)
        fcast_seriess[i], transforms[i] = read_forecast_output(models[i], input_types[i], cond_types[i],
                                                         output_var, var_name, shock_name,
                                                         forecast_string = forecast_strings[i])
    end

    # Reverse transform
    y0_index = get_y0_index(models[1], product)
    data = if haskey(pseudo2data, var_name) && product != :irf
        Float64.(collect(Missings.replace(Vector{Union{Missing, Float64}}(df[!,pseudo2data[var_name]]), NaN)))
    elseif class == :obs && product != :irf
        Float64.(collect(Missings.replace(Vector{Union{Missing, Float64}}(df[!,var_name]), NaN)))
    else
        fill(NaN, size(df, 1))
    end
    # data = class == :obs && product != :irf ? Float64.(collect(Missings.replace(df[:,var_name], NaN))) : fill(NaN, size(df, 1))

    fcast_series = Matrix{Float64}(undef, 0, size(fcast_seriess[1], 2))
    if all(map(x->size(x, 1) > 1, fcast_seriess))
        for i in 1:length(models)
            sampled = sample(1:size(fcast_seriess[i], 1), Int(ndraws * weights[i]))
            fcast_series = vcat(fcast_series, fcast_seriess[i][sampled, :])
        end
    else
        fcast_series = zeros(size(fcast_seriess[1]))
        for i in 1:length(models)
            fcast_series = fcast_series + weights[i]*fcast_seriess[i]
        end
    end

    transformed_series = mb_reverse_transform(fcast_series, transforms[1], product, class,
                                              y0_index = y0_index, data = data,
                                              pop_growth = pop_growth,  use_data = haskey(pseudo2data, var_name))

    # Handle NaNs
    if skipnan && output_var != :histobs && any(isnan.(transformed_series))
        # Remove rows with NaNs
        nanrows = vec(mapslices(x -> any(isnan.(x)), transformed_series, dims = Int[2]))
        transformed_series = transformed_series[.!nanrows, :]
    end

    # Compute means and bands
    means = vec(mean(transformed_series, dims= 1))
    bands = if product in [:shockdec, :dettrend, :trend] && !compute_shockdec_bands
        Dict{Symbol,DataFrame}()
    else
        find_density_bands(transformed_series, density_bands, minimize = minimize)
    end
    return means, bands
end

function mb_reverse_transform(fcast_series::AbstractArray, transform::Function,
                              product::Symbol, class::Symbol;
                              y0_index::Int = -1,
                              data::AbstractVector{Float64} = Float64[],
                              pop_growth::AbstractVector{Float64} = Float64[],
                              use_data::Bool = false)
                              # inflation_series::AbstractArray = Float64[])
    # No transformation
    if product in [:histut, :forecastut, :bddforecastut]
        return fcast_series
    end

    if product in [:hist4q, :forecast4q, :bddforecast4q]
        transform4q = get_transform4q(transform)
        use_data = use_data ? true : (class == :obs && product != :hist4q)

        y0s = if use_data && transform4q in [loggrowthtopct_4q_percapita, loggrowthtopct_4q]
            # Sum growth rates y_{t-3}, y_{t-2}, y_{t-1}, and y_t
            data[y0_index+1:y0_index + 3]
        elseif use_data && transform4q in [logleveltopct_4q_percapita, logleveltopct_4q]
            # Divide log levels y_t by y_{t-4}
            data[y0_index:y0_index + 3]
        else
            Float64[]
        end
        reverse_transform(fcast_series, transform4q;
                          fourquarter = true, y0s = y0s,
                          pop_growth = pop_growth)
    # elseif product in [:histlvl, :forecastlvl, :bddforecastlvl]
    #     # NEED TO ADD CASE HERE FOR HISTLVL, ETC# .
        # reverse_transform(fcast_series, transform4q;
        #                   accumulate = true, y0s = y0s,
        #                   pop_growth = pop_growth)
    else
        # Use transformation that doesn't add back population growth for
        # products which are given in deviations
        if product in [:shockdec, :dettrend, :irf]
            transform = get_nopop_transform(transform)
            y0 = 0.0
        else
            y0 = class == :obs ? data[y0_index] : NaN
        end

        reverse_transform(fcast_series, transform;
                          y0 = y0, pop_growth = pop_growth)
    end
end


text:
using TwilioSMS
using Test

@testset "TwilioSMS.jl" begin
    # Currently we don't have any tests for this package
end


text:
pkgdir = "$(Pkg.dir("toJuMP"))"

@testset "gms/4stufen.gms" begin
    name = "4stufen"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)
    @test m.objSense == :Min
    @test length(m.colVal) == 150
    @test length([i for i in m.colCat if i == :Cont]) == 102
    @test length([i for i in m.colCat if i == :Bin]) == 48

    @test getlowerbound(Variable(m, 13)) == 5
    @test getvalue(Variable(m, 33)) == 176.503

end

@testset "gms/abel.gms" begin
    name = "abel"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)
    @test m.objSense == :Min
    @test length(m.colVal) == 31
    @test length([i for i in m.colCat if i == :Cont]) == 31

    @test getlowerbound(Variable(m, 2)) == 387.9
    @test getupperbound(Variable(m, 2)) == 387.9
    @test getlowerbound(Variable(m, 10)) == 85.3
    @test getupperbound(Variable(m, 10)) == 85.3
    @test getvalue(Variable(m, 3)) == 387.9

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, 225.1945831861349;atol=1e-4)
end

@testset "gms/alkyl.gms" begin
    name = "alkyl"
    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)
    @test m.objSense == :Min
    @test length(m.colVal) == 15
    @test length([i for i in m.colCat if i == :Cont]) == 15
    @test getvalue(Variable(m,1)) == -0.9
    @test getvalue(Variable(m,15)) == 1
    @test getlowerbound(Variable(m, 2)) == 0.0
    @test getupperbound(Variable(m, 2)) == 2
    @test getlowerbound(Variable(m, 3)) == 0.0
    @test getupperbound(Variable(m, 3)) == 1.6
    @test getlowerbound(Variable(m, 4)) == 0.0
    @test getupperbound(Variable(m, 4)) == 1.2
    @test getlowerbound(Variable(m, 7)) == 0.85
    @test getupperbound(Variable(m, 9)) == 12
    @test getlowerbound(Variable(m, 14)) == 0.9
    @test getupperbound(Variable(m, 15)) == 1.01010101010101

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, -1.7650001749159285;atol=1e-4)
end

@testset "gms/autocorr_bern20-03.gms" begin
    name = "autocorr_bern20-03"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)
    @test m.objSense == :Min
    @test length(m.colVal) == 21
    @test length([i for i in m.colCat if i == :Cont]) == 1
    @test length([i for i in m.colCat if i == :Bin]) == 20
end

@testset "gms/bearing.gms" begin
    name = "bearing"
    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)
    @test m.objSense == :Min
    @test length(m. colVal) == 14
    @test length([i for i in m.colCat if i == :Cont]) == 14
    @test getvalue(Variable(m, 2)) == 6
    @test getvalue(Variable(m, 3)) == 5
    @test getvalue(Variable(m, 10)) == 50
    @test getlowerbound(Variable(m, 2)) == 1
    @test getupperbound(Variable(m, 2)) == 16
    @test getlowerbound(Variable(m, 6)) == 1
    @test getupperbound(Variable(m, 6)) == 1000
    @test getupperbound(Variable(m, 2)) == 16

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, 1.9517331852884277;atol=1e-4)
end

@testset "gms/camcge.gms" begin
    name = "camcge"
    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 280
    @test length([i for i in m.colCat if i == :Cont]) == 280
    @test length([i for i in m.colCat if i == :Bin]) == 0

    @test getlowerbound(Variable(m, 101)) == 0.0
    @test getlowerbound(Variable(m, 172)) == 0.0
    @test getlowerbound(Variable(m, 200)) == 0.0

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, -191.73462423686072;atol=1e-4)
end

@testset "gms/ex4_1_3.gms" begin
    name = "ex4_1_3"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 2
    @test length([i for i in m.colCat if i == :Cont]) == 2
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, -443.67170474112413;atol=1e-4)
end

@testset "gms/ex7_3_5.gms" begin
    name = "ex7_3_5"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 14
    @test length([i for i in m.colCat if i == :Cont]) == 14
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    status = solve(m)
    @test status == :Error
end

@testset "gms/ex14_2_4.gms" begin
    name = "ex14_2_4"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 6
    @test length([i for i in m.colCat if i == :Cont]) == 6
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, 7.562126836459202e-9;atol=1e-4)
end

@testset "gms/ex14_2_8.gms" begin
    name = "ex14_2_8"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 5
    @test length([i for i in m.colCat if i == :Cont]) == 5
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, 2.624689270273677e-9;atol=1e-4)
end

@testset "gms/st_rv1.gms" begin
    name = "st_rv1"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 11
    @test length([i for i in m.colCat if i == :Cont]) == 11
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    solve(m)
    @test isapprox(m.objVal, -59.904345906908;atol=1e-4)
end

@testset "gms/water.gms" begin
    name = "water"

    gmspath = joinpath(pkgdir, "test", "examples", "gms", "$(name).gms")
    jlpath = joinpath(pkgdir, ".jls", "$(name).jl")

    gms2jump(gmspath)
    m = include(jlpath)
    @test isa(m, JuMP.Model)

    @test m.objSense == :Min
    @test length(m.colVal) == 42
    @test length([i for i in m.colCat if i == :Cont]) == 42
    @test length([i for i in m.colCat if i == :Bin]) == 0

    setsolver(m, IpoptSolver(print_level=0))
    status = solve(m)
    @test status == :Error
end


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Arbitrary-physical-units
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(
    ncols=3, width="15cm", height="2.5in",
    wspace=("10pt", "20pt"), right="10mm",
)
panel = axs[3].panel_axes("r", width="2em")

# format each ax instead of the python context change using `with`
[ax.format(small="12px", large="15px") for ax in axs]
panel.format(small="12px", large="15px")

[ax.format(
    suptitle="Arguments with arbitrary units",
    xlabel="x axis", ylabel="y axis"
) for ax in axs]
f.savefig("figs/Arbitrary-physical-units.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
module SUniward
import Base: setindex!,getindex,copy!
using Images
using OffsetArrays
import ..STC: matchlsb!,readlsb,readlsbq,LSBcostfun, AbstractImageDistortion, trypmone, ind2sub

include("incrementalsuniward.jl")
include("suniwardadd.jl")
include("convolutions.jl")

const hpdf = [-0.0544158422, 0.3128715909, -0.6756307363, 0.5853546837, 0.0158291053, -0.2840155430, -0.0004724846, 0.1287474266, 0.0173693010, -0.0440882539, -0.0139810279, 0.0087460940, 0.0048703530, -0.0003917404, -0.0006754494, -0.0001174768];
const lpdf = (-1).^(0:length(hpdf)-1).*hpdf[end:-1:1];
const flts=[lpdf*hpdf',hpdf*lpdf',hpdf*hpdf'];
const rotflts=rot180.(flts)
const wetCost = 10^8;

function suniwardcosts(cover::Array{T},sigma::T=T(1)) where {T<:AbstractFloat}
    coverPadded = padarray(cover,Pad(:symmetric,[16,16],[16,16]))
    xi = Array{Array{Float64,2},1}(undef, 3)
    R = Array{Array{Float64,2},1}(undef, 3)
    for f = 1:3
        R[f] = sameconv2(coverPadded.parent, flts[f])
        xi[f] = sameconv2(1 ./ (abs.(R[f]) .+ sigma), abs.(rotflts[f]))[16:end-17,16:end-17]
    end
    costs = xi[1] .+ xi[2] .+ xi[3]
    return(costs)
end

function suniwardcosts(cover::Array{T},stego::Array{T},sigma::T=T(1)) where {T<:AbstractFloat}
    df=cover-stego
    coverPadded= padarray(cover,Pad(:symmetric,[16,16],[16,16]))
    xi=Array{Array{Float64,2},1}(undef, 3)
    invr = padarray(zeros(eltype(cover),size(cover)),Pad(:symmetric,[16,16],[16,16]))
    for f = 1:3
        invr.parent .= 1 ./ (abs.(sameconv2(coverPadded.parent, flts[f])) .+ sigma)
        xi[f] = sameconv2(df,rotflts[f]) .* invr.parent[17:end-16,17:end-16]
    end
    return(ceiledsumabs(xi[1],1e8) + ceiledsumabs(xi[2],1e8) + ceiledsumabs(xi[3],1e8))
end

suniwardcosts(cover,stego,sigma = 1) = suniwardcosts(Float64.(cover), Float64.(stego), Float64(sigma))

end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to IGRF model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/geomagnetic_field_models
# ==========================================

# Load the test files.
#
# The test values were created using the original `igrf12syn` function. The
# source code of the program that generated those files can be seen on:
#
#   https://github.com/SatelliteToolbox/IGRF_Test
#

igrf12_geocentric_test = readdlm("./IGRF12_test_geocentric.txt")
igrf12_geodetic_test   = readdlm("./IGRF12_test_geodetic.txt")

# Function igrf12syn
# ------------------

@testset "Function igrf12syn" begin
    # Testing the geocentric part of the algorithm.
    for i = 1:size(igrf12_geocentric_test, 1)
        date  = igrf12_geocentric_test[i,1]
        r     = igrf12_geocentric_test[i,2]
        colat = igrf12_geocentric_test[i,3]
        elong = igrf12_geocentric_test[i,4]
        xt    = igrf12_geocentric_test[i,5]
        yt    = igrf12_geocentric_test[i,6]
        zt    = igrf12_geocentric_test[i,7]
        ft    = igrf12_geocentric_test[i,8]

        # Call IGRF with the same inputs as those in the test.

        if date < 2020
            (x, y, z, f) = igrf12syn(0, date, 2, r, colat, elong)
        else
            (x, y, z, f) = @test_logs (:warn, "The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.") #=
                =# igrf12syn(0, date, 2, r, colat, elong)
        end

        # Test the values.
        @test x ≈ xt atol=3e-1
        @test y ≈ yt atol=3e-1
        @test z ≈ zt atol=3e-1
        @test f ≈ ft atol=3e-1
    end

    # Testing the geodetic part of the algorithm.
    for i = 1:size(igrf12_geodetic_test, 1)
        date  = igrf12_geodetic_test[i,1]
        h     = igrf12_geodetic_test[i,2]
        colat = igrf12_geodetic_test[i,3]
        elong = igrf12_geodetic_test[i,4]
        xt    = igrf12_geodetic_test[i,5]
        yt    = igrf12_geodetic_test[i,6]
        zt    = igrf12_geodetic_test[i,7]
        ft    = igrf12_geodetic_test[i,8]

        # Call IGRF with the same inputs as those in the test.
        if date < 2020
            (x, y, z, f) = igrf12syn(0, date, 1, h, colat, elong)
        else
            (x, y, z, f) = @test_logs (:warn, "The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.") #=
                =# igrf12syn(0, date, 1, h, colat, elong)
        end

        # Test the values.
        @test x ≈ xt atol=3e-1
        @test y ≈ yt atol=3e-1
        @test z ≈ zt atol=3e-1
        @test f ≈ ft atol=3e-1
    end
end

# Function igrf12
# ---------------

@testset "Function igrf12" begin
    # Testing the geocentric part of the algorithm.
    for i = 1:size(igrf12_geocentric_test, 1)
        date  = igrf12_geocentric_test[i,1]
        r     = igrf12_geocentric_test[i,2]
        colat = igrf12_geocentric_test[i,3]
        elong = igrf12_geocentric_test[i,4]
        xt    = igrf12_geocentric_test[i,5]
        yt    = igrf12_geocentric_test[i,6]
        zt    = igrf12_geocentric_test[i,7]
        ft    = igrf12_geocentric_test[i,8]

        # Call IGRF with the same inputs as those in the test.
        (elong > 180) && (elong = elong-360)

        if date < 2020
            B = igrf12(date, r*1000, (90-colat)*pi/180, elong*pi/180)
        else
            B = @test_logs (:warn, "The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.") #=
                =# igrf12(date, r*1000, (90-colat)*pi/180, elong*pi/180)
        end

        x, y, z = B[:]
        f       = norm(B)

        # Test the values.
        @test x ≈ xt atol=3e-1
        @test y ≈ yt atol=3e-1
        @test z ≈ zt atol=3e-1
        @test f ≈ ft atol=3e-1
    end

    # Testing the geodetic part of the algorithm.
    for i = 1:size(igrf12_geodetic_test, 1)
        date  = igrf12_geodetic_test[i,1]
        h     = igrf12_geodetic_test[i,2]
        colat = igrf12_geodetic_test[i,3]
        elong = igrf12_geodetic_test[i,4]
        xt    = igrf12_geodetic_test[i,5]
        yt    = igrf12_geodetic_test[i,6]
        zt    = igrf12_geodetic_test[i,7]
        ft    = igrf12_geodetic_test[i,8]

        # Call IGRF with the same inputs as those in the test.
        (elong > 180) && (elong = elong-360)

        if date < 2020
            B = igrf12(date, h*1000, (90-colat)*pi/180, elong*pi/180, Val{:geodetic})
        else
            B = @test_logs (:warn, "The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.") #=
                =# igrf12(date, h*1000, (90-colat)*pi/180, elong*pi/180, Val{:geodetic})
        end

        x, y, z = B[:]
        f       = norm(B)

        # Test the values.
        @test x ≈ xt atol=3e-1
        @test y ≈ yt atol=3e-1
        @test z ≈ zt atol=3e-1
        @test f ≈ ft atol=3e-1
    end
end


text:
include("../Preprocessing/Encoder.jl")

"""
    SVC(; alpha=0.01, ni=1000)
Support Vector Machine Classifier.

This model learns classifiers using One-Vs-Rest.
One-Vs-Rest generates two-class classifiers divided into one class and the other classes using Logistic Regression, adopting the most likely one among all classifiers.

# Parameters:
- `K`: number of class
- `max`: maximum number of repitition
- `th`: converge threshold

# Example
```jldoctest classification
julia> model = SVC()
SVC(0.01, 1000, Logistic[])

julia> fit!(model, ct)
3-element Vector{Logistic}:
 Logistic(0.01, 1000, [0.8116709490679518 1.188329050932049; 1.7228257190036231 0.2771742809963788; -0.1519960725403138 2.1519960725403116])
 Logistic(0.01, 1000, [0.9863693439936144 1.0136306560063886; 0.8838433946106077 1.11615660538939; 1.4431044559203794 0.5568955440796174])
 Logistic(0.01, 1000, [1.262510641510418 0.7374893584895849; 0.5242383002319192 1.4757616997680822; 1.864635796779504 0.135364203220495])

julia> println(predict(model, x))
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
```
"""
mutable struct SVC
    α::Float64
    n_iter::Int64
    classifiers::Array{Logistic, 1}
    SVC(; alpha=0.01, ni=1000) = new(alpha, ni, Array{Logistic}(undef, 0))
end

function fit!(model::SVC, x, t)
    alpha, n_iter = model.α, model.n_iter
    check_size(x, t)
    c = size(t, 2)
    w = Array{Float64}(undef, 0, 0)
    classifiers = Array{Logistic}(undef, 0)
    Threads.@threads for i in 1 : c
        classifier = Logistic(alpha = alpha, ni = n_iter)
        OHE = OneHotEncoder()
        fit!(classifier, x, OHE(t[:, i]))
        push!(classifiers, classifier)
    end
    model.classifiers = classifiers
end

function (model::SVC)(x)
    p = Array{Float64}(undef, size(x, 1), length(model.classifiers))
    for i in 1 : length(model.classifiers)
        p[:, i] = forecast(model.classifiers[i], x)[:, 2]
    end
    return [findfirst(p[i, :] .== maximum(p[i, :])) for i in 1:size(p, 1)]
end


text:
#=
This shows some data wrangling using the Card data
on returns to education
=#

##
using DelimitedFiles, CSV, DataFrames, DataFramesMeta
card = CSV.read("../Data/card.csv")
display(head(card))

##
# add some variables
card[:lnwage] = round.(log.(card[:wage]), digits=4) # round to save storage space
card[:expsq] = (card[:exper].^2)/100
card[:agesq] = card[:age].^2

##
# select the ones we want
cooked = @select(card, :lnwage, :educ, :exper, :expsq, :black, :south, :smsa, :age, :agesq, :nearc4 )
display(head(cooked))

##
# write as CSV
CSV.write("cooked.csv", cooked)

##
# write as plain ASCII
data = convert(Array{Float64}, cooked)
writedlm("cooked.txt", data)

# for binary write, see JLD.jl and Feather.jl


text:
module ModelStatisticsLH

using ArgCheck, DocStringExtensions, Random
using CommonLH, StructLH

export VarInfo, check_var_info, check_var
export get_option, get_options, bounds, isscalar, var_name, value_type
export init_values
export make_test_var_info, make_test_var_infos, make_test_values

export ModelStats
export n_vars, var_names, var_meta, get_meta
export make_test_model_stats, validate_stats
export has_variable, add_variable!, delete_variable!, get_values, set_values!

export StatsCollection
export add_mstats!, delete_mstats!, replace_mstats!
export get_groups, has_mstats, get_mstats, get_stats, set_stats!
export groups_from_name, statname_from_name
export make_test_scollection

export reduce_one_field, reduce_stats

include("var_info.jl");
include("group_stats.jl");
include("stats_collection.jl");
include("reductions.jl");

format_values(v :: Vector{I}) where I <: Integer =
    string.(v);

format_values(v :: Vector{R}) where R <: Real =
    string.(round.(v; digits = 2));

format_values(v) = string.(v);

end #module


text:
#
# @description Defines the stream-cipher abstract type and the helper methods that act on this abstract type.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export StreamCipher

# @description Defines the stream-cipher abstract type. A stream-cipher is a symmetric key cipher where plaintext
# digits are combined with a pseudorandom cipher digit stream.
abstract type StreamCipher end

# @description Returns the encryption function linked to the given stream-cipher.
#
# @param {StreamCipher} self The cipher data struture.
#
# @return {Function} A function that only takes one {Vector{UInt8}} argument.
function encryptor(self::StreamCipher)
  return plaintext -> encrypt!(self, plaintext)
end

# @description Returns the decryption function linked to the given stream-cipher.
#
# @param {StreamCipher} self The cipher data struture.
#
# @return {Function} A function that only takes one {Vector{UInt8}} argument.
function decryptor(self::StreamCipher)
  return ciphertext -> decrypt!(self, ciphertext)
end


text:
using ParallelIncrements
using Test

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


text:
const logkappa_spl = opal_test_spl()

function drdm(m, r, l, P, T, rho, mu)
    #=
    Derivative of r with respect to m, from mass conservation.
    =#
    return 1 / (4 * pi * r .^ 2 .* rho)
end # drdm

function dPdm(m, r, l, P, T)
    #=
    Derivative of P with respect to m, from hydrostatic equilibrium.
    =#
    return -G * m / (4 * pi * r .^ 4)
end # dPdm

function dldm(m, r, l, P, T, rho, X, Y)
    #=
    Derivative of l with respect to m, from energy conservation.
    =#
    return epsilon_pp(rho, T, X) + epsilon_cno(rho, T, X, Y)
end # dldm

function dTdm(m, r, l, P, T, mu, kappa)
    #=
    Derivative of T with respect to m, from... thermodynamics.
    =#
    grad_rad = ∇_rad(m, l, P, T, kappa)
    grad_ad = ∇_ad()
    stable =  grad_rad < grad_ad
    if stable
        # @debug("dTdm: stable to convection")
        grad = grad_rad
    else
        # @debug("dTdm: unstable to convection")
        grad = grad_ad
    end # if
    factor = - G .* m .* T ./ (4 * pi .* r .^ 4 .* P)
     return factor .* grad
end # dTdm


text:
## Exercise 4-11
## The letters of the alphabet can be constructed from a moderate number of basic elements, like vertical and horizontal lines and a few curves. Design an alphabet that can be drawn with a minimal number of basic elements and then write functions that draw the letters.

## You should write one function for each letter, with names draw_a, draw_b, etc., and put your functions in a file named letters.jl.
println("Ans: ")

println("End.")


text:
##########################################
# Master file for running all test cases #
##########################################
using Turing; turnprogress(false)
using Test

@debug("[runtests.jl] runtests.jl loaded")

path = dirname(@__FILE__)
cd(path); include("utility.jl")
@debug("[runtests.jl] utility.jl loaded")
@debug("[runtests.jl] testing starts")

# run all tests
runtests()

@debug("[runtests.jl] all tests finished")


text:
module CoordinateSplittingPTrees

using LinearAlgebra, SparseArrays
using AbstractTrees  # for display of tree structures
import BlossomV      # we use matching to select dimension-pairs (see cs2.jl)
using DataStructures # for PriorityQueue
import Parametron    # for solving the lower-bound model. `import` prevents a conflict with Parametron.value
using Parametron: Model, Variable, Parameter, Minimize, @expression, @constraint, @objective, solve!

export Box, World
export position, boxbounds, meta, value, addpoint!
export getroot, getleaf, find_leaf_at, degree, isroot, isleaf, chaintop, iscomplete
export leaves, splits, chain, splitprint, splitprint_colored, print_tree, print_metabox

include("types.jl")
include("tree.jl")
include("polynomials.jl")
include("cs2.jl")
include("posdef.jl")

"""
CoordinateSplittingPTrees implements coordinate-splitting trees of
degree `p` (CSp-trees). These trees are reminiscent of [k-d
trees](https://en.wikipedia.org/wiki/K-d_tree), and are useful for
multidimensional full-degree polynomial interpolation.

The world bounds are specified by creating a [`World`](@ref) object.

The tree is initialized with `root = Box{p}(world)`, and new leaves are
added with `Box(parent, dims, xs, metas)` (see [`Box`](@ref)).

You can interact with boxes using the following API:

- [`position`](@ref): retrieve the position of a `Box`
- [`meta`](@ref): retrieve the metadata associated with a `Box`
- [`value`](@ref): retrieve the function value associated with a `Box`
- [`boxbounds`](@ref): retrieve the edges of a `Box`
- [`addpoint!`](@ref): split boxes to generate new positions
- [`isroot`](@ref): determine whether a given `Box` is the root-node
- [`isleaf`](@ref): determine whether a given `Box` is a leaf-node
- [`getroot`](@ref): get the root node associated with a tree
- [`getleaf`](@ref): get the leaf node associated with a (possibly-parent) box
- [`find_leaf_at`](@ref): find the leaf node containing a particular position

- [`leaves`](@ref): returns an iterator for visiting all leaf-nodes

- [`print_tree`](@ref): display the tree structure (from AbstractTrees.jl)
- [`splitprint`](@ref): a more compact display of the tree structure
- [`splitprint_colored`](@ref): similar to `splitprint` but with highlights
"""
CoordinateSplittingPTrees

"""
    z = value(box)

Retreive the function value associated with `box`. `value` calls
[`meta`](@ref), so whatever type `M` you use for metadata needs to
have a `value(::M)` function defined. `value(z::Real) = z` is
pre-defined.
"""
value(box::Box) = value(meta(box))
value(z::Real) = z

Base.isless(b1::Box, b2::Box) = isless(value(b1), value(b2))

Base.minimum(box::Box) = minimum(leaves(box))
Base.maximum(box::Box) = maximum(leaves(box))
Base.extrema(box::Box) = extrema(leaves(box))

"""
    addpoint!(box::Box, x, dimlists, metagen::Function)

Create new evaluation points by splitting existing boxes. `x`, the
final evaluation position, must be inside `box`. Dimensions are split
in the order specified by `dimlists`, which should be a
list-of-dimension-lists. For example, for a CS2 tree one might use

    dimlists = [(1,3), (2,4), (5,9), ...]

to indicate that the first split is for dimensions 1 and 3, the second
for dimensions 2 and 4, etc. All `n` dimensions should appear once in
`dimlists`; if `n` is not an integer multiple of `p`, then the final
entry should contain only the "leftover" dimensions. (Fictive
dimension(s) will be added automatically to make all splits of length
`p`.  If fictive dimensions are being used, some leaf-nodes will
appear to be duplicates of others, but in reality they are "separated"
(with size 0) along fictive dimension(s), which are not displayed, and
share the same metadata as other boxes with the same position in the
non-fictive dimensions.)

`metagen(y)` must return the metadata to be associated with position `y`,
where `y` is a vector.
"""
function addpoint!(box::Box, x, dimlists::Union{AbstractVector,Tuple}, metagen::Function, metaboxgen::Function=x->nothing)
    # Validate dimlists
    n = ndims(box)
    covered = falses(n)
    for dimlist in dimlists
        for d in dimlist
            covered[d] && error("dimension $d was duplicated")
            covered[d] = true
        end
    end
    all(covered) || error("all dimensions must be covered, got $dimlists")

    if !isleaf(box)
        box = find_leaf_at(box, x)
    end
    b = position(box)
    y = copy(b)
    L = maxchildren(box)-1
    metas = Vector{typeof(meta(box))}(undef, L)
    l = falses(degree(box))
    for dimlist in dimlists
        box0 = box
        dimlistv = [dimlist...]
        empty!(metas)
        length(dimlist) < degree(box) && resize!(l, length(dimlist))
        for i = 1:2^length(dimlist)-1
            l.chunks[1] = i  # convert the counter into a logical bit pattern
            dl = dimlistv[l]
            y0 = y[dl]
            y[dl] = x[dl]
            push!(metas, metagen(y))
            y[dl] = y0
        end
        others = Box(box, dimlist, (d->x[d]).(dimlist), (metas...,))
        box0.split.self.metabox = metaboxgen(box0.split.self)
        for bx in others
            bx.metabox = metaboxgen(bx)
        end
        # Ensure the returned box (and entire chain) is not displaced
        # along fictive dimensions
        cindex = 0x00
        for d in dimlist
            cindex = (cindex << 0x01) | (d <= n)
        end
        box = others[cindex]
        y[dimlistv] = x[dimlistv]
    end
    return box
end

"""
    addpoint!(box::Box, x, metagen::Function)

Create new evaluation points by splitting existing boxes. `x`, the
final evaluation position, must be inside `box`, and typically `box`
is the root node. Choice of splitting dimensions is automatic.

`metagen(y)` must return the metadata to be associated with position `y`,
where `y` is a vector.
"""
function addpoint!(root::Box, x, metagen::Function, metaboxgen::Function=x->nothing)
    leaf = find_leaf_at(root, x)
    dimlists = choose_dimensions(leaf)
    addpoint!(root, x, dimlists, metagen, metaboxgen)
end

"""
    addpoint_distinct!(box::Box, x, metagen::Function)

Create new evaluation points by splitting existing boxes. `x` is the
desired evaluation point; however, for any `x[i]` that are identical
to the existing target `leaf = find_leaf_at(root, x)`, a new value is
chosen using `partition_interval`.

See [`addpoint!`](@ref) for additional information.
"""
function addpoint_distinct!(root, x, metagen::Function, metaboxgen::Function=x->nothing)
    xcopy = copy(x)
    leaf = find_leaf_at(root, x)
    xleaf = position(leaf)
    for i = 1:ndims(root)
        if x[i] == xleaf[i]
            xcopy[i] = partition_interval(x[i], leaf, i)
        end
    end
    addpoint!(root, xcopy, metagen, metaboxgen)
end

"""
    y = partition_interval(x::Real, box, splitdim)

Along dimension `splitdim`, split the largest interval between the
evaluation point `x` and the edges, so that the gap between points and
edges is cut to one-third of the original (thus making all "new"
intervals equal in size):

```
old: |                 x  |
new: |     y     |     x  |
```

When the largest interval is infinite, there are two cases:

- if both edges are infinitely far away, an arbitrary point is
  generated (or select a distinct value of `world.position[splitdim]`,
  if it was initialized with a 2-tuple)
- if only one interval has infinite extent, the gap between the points
  and the new dividing edge doubles that of the interval on the
  opposite side:

```
old: |     x
new: |     x          |          y
```
"""
function partition_interval(x::Real, box::Box, splitdim)
    bb = boxbounds(box, splitdim)
    if isfinite(bb[1]) && isfinite(bb[2])
        if bb[2] - x > x - bb[1]
            return (2*bb[2] + x)/3
        else
            return (x + 2*bb[1])/3
        end
    elseif isfinite(bb[2])
        return x - 4*(bb[2]-x)
    elseif isfinite(bb[1])
        return x + 4*(x-bb[1])
    end
    w = box.world
    return world_newposition(x, w.position[splitdim], w.lower[splitdim], w.upper[splitdim])
end

function nanzero!(A::AbstractArray)
    z = zero(eltype(A))
    for I in eachindex(A)
        if isnan(A[I])
            A[I] = z
        end
    end
    return A
end

nanzero!(A::SymmetricArray) = nanzero!(A.data)

function nanzero!(A::SparseMatrixCSC)
    nanzero!(A.nzvals)
    return A
end

function nanzero!(A::SymTridiagonal)
    nanzero!(A.ev)
    nanzero!(A.dv)
    return A
end

function randleaf(root)
    L = maxchildren(root)
    box = root
    while !isleaf(box)
        split = box.split
        children = (split.self, split.others.children...)
        box = children[rand(1:L)]
    end
    return box
end

function init_plotting()
    push!(LOAD_PATH, @__DIR__)
    @eval Main using CSpPlots
    pop!(LOAD_PATH)
    nothing
end

end # module


text:
export allelecount, alleles, uniquealleles
export locidataframe, locimatrix, phasedmatrix

"""
    allelecount(locus::T) where T<:GenoArray
Return the number of unique alleles present at a locus.
"""
@inline function allelecount(locus::T) where T<:GenoArray
    unique(locus) |> skipmissing |> Base.Iterators.flatten |> unique |> length
end

"""
    alleles(locus::T) where T<:GenoArray
Return an array of all the non-missing alleles of a locus.
"""
@inline function alleles(locus::T) where T<:GenoArray
    if all(ismissing.(locus))
        int_type = eltype(typeof(locus)) |> nonmissingtype |> eltype
        return Vector{Union{Missing, int_type}}(undef, length(locus))
    end
    return Base.Iterators.flatten(skipmissing(locus)) |> collect
end

"""
    alleles(locus::T, miss::Bool = false) where T<:GenoArray
Return an array of all the non-missing alleles of a locus. Use the second positional
argument as `true` to include missing values.
"""
@inline function alleles(locus::T, miss::Bool) where T<:GenoArray
    int_type = eltype(typeof(locus)) |> nonmissingtype |> eltype
    if all(ismissing.(locus))
        return Vector{Union{Missing, int_type}}(undef, length(locus))
    end
    alle_out = Vector{Union{Missing, int_type}}(Base.Iterators.flatten(skipmissing(locus)) |> collect)
    if miss == true
        append!(alle_out, locus[locus .=== missing])
    end
    return alle_out
end

"""
    uniquealleles(locus::T) where T<:GenoArray
Return an array of all the unique non-missing alleles of a locus.
"""
@inline function uniquealleles(locus::GenoArray)
    unique(alleles(locus))
end

"""
    locidataframe(data::PopData)
Return a wide `DataFrame` of samples as columns, ommitting population information.

**Example**
```
julia> locidataframe(@nancycats)
9×237 DataFrame. Omitted printing of 232 columns
│ Row │ N215       │ N216       │ N217       │ N218       │ N219       │
│     │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │
├─────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│ 1   │ missing    │ missing    │ (135, 143) │ (133, 135) │ (133, 135) │
│ 2   │ (136, 146) │ (146, 146) │ (136, 146) │ (138, 138) │ (140, 146) │
│ 3   │ (139, 139) │ (139, 145) │ (141, 141) │ (139, 141) │ (141, 145) │
│ 4   │ (116, 120) │ (120, 126) │ (116, 116) │ (116, 126) │ (126, 126) │
│ 5   │ (156, 156) │ (156, 156) │ (152, 156) │ (150, 150) │ (152, 152) │
│ 6   │ (142, 148) │ (142, 148) │ (142, 142) │ (142, 148) │ (142, 148) │
│ 7   │ (199, 199) │ (185, 199) │ (197, 197) │ (199, 199) │ (193, 199) │
│ 8   │ (113, 113) │ (113, 113) │ (113, 113) │ (91, 105)  │ (113, 113) │
│ 9   │ (208, 208) │ (208, 208) │ (210, 210) │ (208, 208) │ (208, 208) │
```
"""
function locidataframe(data::PopData)
    unstack(select(data.genodata, Not(:population)), :name, :genotype)[:, Not(:locus)]
end

#TODO make a SMatrix instead?
"""
    locimatrix(data::PopData)
Return a matrix of genotypes with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> locimatrix(@nancycats)
237×9 Array{Union{Missing, Tuple{Int16,Int16}},2}:
 missing     (136, 146)  (139, 139)  …  (199, 199)  (113, 113)  (208, 208)
 missing     (146, 146)  (139, 145)     (185, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (141, 141)     (197, 197)  (113, 113)  (210, 210)
 (133, 135)  (138, 138)  (139, 141)     (199, 199)  (91, 105)   (208, 208)
 (133, 135)  (140, 146)  (141, 145)     (193, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (145, 149)  …  (193, 195)  (91, 113)   (208, 208)
 (135, 135)  (136, 146)  (139, 145)     (199, 199)  (105, 113)  (208, 208)
 (135, 143)  (136, 146)  (135, 149)     (193, 197)  (91, 91)    (208, 212)
 (137, 143)  (136, 146)  (139, 139)     (197, 197)  (105, 113)  (208, 212)
 (135, 135)  (132, 132)  (141, 145)     (197, 197)  (91, 105)   (208, 208)
 (137, 141)  (130, 136)  (137, 145)  …  (193, 199)  (91, 91)    (182, 182)
 (129, 133)  (130, 136)  (135, 145)     (193, 199)  (91, 113)   (182, 208)
 ⋮                                   ⋱
 (133, 135)  (136, 136)  (135, 139)  …  (199, 199)  (113, 113)  (182, 182)
 (133, 141)  (136, 136)  (135, 139)     (197, 197)  (113, 113)  (182, 208)
 (133, 141)  (130, 146)  (141, 141)     (191, 199)  missing     (208, 208)
 (123, 133)  (138, 138)  (141, 145)     (191, 197)  missing     (208, 208)
 (123, 133)  (138, 138)  (139, 139)     (197, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 139)  …  (197, 197)  missing     (208, 208)
 (133, 141)  (130, 136)  (139, 145)     (191, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 145)     (199, 199)  missing     (208, 220)
 (133, 143)  (130, 130)  (135, 145)     (197, 197)  missing     (208, 208)
 (135, 141)  (136, 144)  (143, 143)     (191, 197)  (113, 117)  (208, 208)
 (137, 143)  (130, 136)  (135, 145)  …  (193, 199)  (113, 117)  (208, 208)
 (135, 141)  (130, 146)  (135, 139)     (197, 197)  missing     (208, 208)
 ```
"""
function locimatrix(data::PopData)
    dims = size(data)
    sort_df = issorted(data.genodata, [:name, :locus]) ? sort(data.genodata, [:name, :locus]) : data.genodata
    reshape(sort_df.genotype, (dims.samples, dims.loci)) |> collect
end

"""
    phasedmatrix(data::PopData)
Return a `Vector` of length `ploidy` composed of allele matrices with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> mtx = phasedmatrix(@nancycats)
2-element Array{Array{Union{Missing, Int16},2},1}:
 [missing 136 … 113 208; missing 146 … 113 208; … ; 137 130 … 113 208; 135 130 … missing 208]
 [missing 146 … 113 208; missing 146 … 113 208; … ; 143 136 … 117 208; 141 146 … missing 208]

julia> mtx[1]
237×9 Array{Union{Missing, Int16},2}:
    missing  136  139  116         156  142  199  113         208
    missing  146  139  120         156  142  185  113         208
 135         136  141  116         152  142  197  113         210
 133         138  139  116         150  142  199   91         208
 133         140  141  126         152  142  193  113         208
 135         136  145  120         150  148  193   91         208
 135         136  139  116         152  142  199  105         208
 135         136  135  120         154  142  193   91         208
 137         136  139  116         150  142  197  105         208
 135         132  141  120         150  148  197   91         208
 137         130  137  128         152  142  193   91         182
 129         130  135  126         144  140  193   91         182
   ⋮                                      ⋮
 133         136  135     missing  146  142  199  113         182
 133         136  135     missing  150  142  197  113         182
 133         130  141     missing  148  142  191     missing  208
 123         138  141     missing  148  142  191     missing  208
 123         138  139     missing  150  142  197     missing  208
 133         136  139     missing  150  142  197     missing  208
 133         130  139     missing  152  142  191     missing  208
 133         136  139     missing  150  142  199     missing  208
 133         130  135     missing  148  142  197     missing  208
 135         136  143     missing  144  142  191  113         208
 137         130  135     missing  150  142  193  113         208
 135         130  135     missing  150  142  197     missing  208
```
"""
function phasedmatrix(data::PopData)
    dims = size(data)
    ploidy = unique(data.sampleinfo.ploidy)
    ploidy = length(ploidy) != 1 ? error("Phasing will not work on mixed-ploidy samples") : ploidy[1]
    sort_df = issorted(data.genodata, [:name, :locus]) ? sort(data.genodata, [:name, :locus]) : data.genodata
    matrices = map(j -> map(i -> ismissing(i) ? missing : i[j] , sort_df.genotype), 1:ploidy)
    map(i -> collect(reshape(i, (dims.samples, dims.loci))), matrices)
end


text:
using JuMP
using CPLEX
include("parser.jl")

function plans_coupants(MyFileName::String)
  n, s, t, S, d1, d2, p, ph, d, D = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  # probleme maitre
  ## Variables
  @variable(m, x[1:n,1:n], Bin) # certains arcs n'existent pas, ils ne sont pas utilises dans les contraintes
  @variable(m, y[1:n], Bin)
  @variable(m, z >= 0)

  ## Constraints
  @constraint(m, sum(x[s,j] for j in 1:n if d[s,j] != 0) == 1)
  @constraint(m, sum(x[i,t] for i in 1:n if d[i,t] != 0) == 1)
  @constraint(m, [v in 1:n; v != s && v != t], sum(x[i,v] for i in 1:n if d[i,v] != 0) == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m, [v in 1:n; v != t], y[v] == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m, [i in 1:n; i != t], y[t] == sum(x[i,t] for i in 1:n if d[i,t] != 0))
  @constraint(m, sum(p[i]*y[i] for i in 1:n) <= S)
  @constraint(m, z >= sum(d[i,j]*x[i,j] for i in 1:n, j in 1:n if d[i,j] != 0)) # objectif robuste reformule

  ## Objective
  @objective(m, Min, z)

  #resolution
  optimize!(m)

  z_star = JuMP.objective_value.(m)
  x_star = JuMP.getvalue.( m[:x] )
  y_star = JuMP.getvalue.( m[:y] )

  z1 = z_star + 1
  z2 = S + 1

  while z2 >= S + 1e-3 || z1 > z_star + 1e-3 || z1 < z_star - 1e-3

    #### sous-probleme 1 ####
    # Create the model
    m1 = Model(CPLEX.Optimizer)

    ## Variables
    @variable(m1, delta1[1:n,1:n] >= 0)

    ## Constraints
    @constraint(m1,[i in 1:n, j in 1:n], delta1[i,j] <= D[i,j])
    @constraint(m1, sum(delta1[i,j] for i in 1:n, j in 1:n) <= d1)

    ## Objective
    @objective(m1, Max, sum(d[i,j]*(1 + delta1[i,j])*x_star[i,j] for i in 1:n, j in 1:n if d[i,j] != 0))

    #resolution
    optimize!(m1)
    z1 = JuMP.objective_value.(m1)
    delta_1 = JuMP.getvalue.( m1[:delta1] )

    #### sous-probleme 2 ####
    # Create the model
    m2 = Model(CPLEX.Optimizer)

    ## Variables
    @variable(m2, delta2[1:n] >= 0)

    ## Constraints
    @constraint(m2,[v in 1:n], delta2[v] <= 2)
    @constraint(m2, sum(delta2[v] for v in 1:n) <= d2)

    ## Objective
    @objective(m2, Max, sum((p[i] + delta2[i]*ph[i])*y_star[i] for i in 1:n))

    #resolution
    optimize!(m2)
    z2 = JuMP.objective_value.(m2)
    delta_2 = JuMP.getvalue.( m2[:delta2] )

    if z2 >= S + 1e-3 || z1 > z_star + 1e-3 || z1 < z_star - 1e-3 # pour gerer la premiere iteration
      #### adding Constraints
      @constraint(m, z >= sum(x[i,j]*d[i,j]*(1+ delta_1[i,j]) for i in 1:n, j in 1:n if d[i,j] != 0))
      @constraint(m, sum(y[v]*(p[v] + ph[v]*delta_2[v]) for v in 1:n) <= S)

      optimize!(m)
      z_star = JuMP.objective_value.(m)
      x_star = JuMP.getvalue.( m[:x] )
      y_star = JuMP.getvalue.( m[:y] )
    end
  end
end

max_stars_repo_name: boumpteryx/ECMA

max_issues_repo_name: boumpteryx/ECMA

max_forks_repo_name: boumpteryx/ECMA


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "regex" begin
    function collect_eachmatch(re, str; overlap=false)
        [m.match for m in collect(eachmatch(re, str, overlap = overlap))]
    end

    @test collect_eachmatch(r"a?b?", "asbd") == ["a","","b","",""] ==
        collect_eachmatch(r"""a?b?""", "asbd")
    @test collect_eachmatch(r"a?b?", "asbd", overlap=true) == ["a","","b","",""]
    @test collect_eachmatch(r"\w+", "hello", overlap=true) == ["hello","ello","llo","lo","o"]
    @test collect_eachmatch(r".\s", "x \u2200 x \u2203 y") == ["x ", "∀ ", "x ", "∃ "]
    @test collect_eachmatch(r"(\w+)(\s*)", "The dark side of the moon") ==
        ["The ", "dark ", "side ", "of ", "the ", "moon"]
    @test collect_eachmatch(r"", "") == [""]
    @test collect_eachmatch(r"", "", overlap=true) == [""]
    @test collect_eachmatch(r"aa", "aaaa") == ["aa", "aa"]
    @test collect_eachmatch(r"aa", "aaaa", overlap=true) == ["aa", "aa", "aa"]
    @test collect_eachmatch(r"", "aaa") == ["", "", "", ""]
    @test collect_eachmatch(r"", "aaa", overlap=true) == ["", "", "", ""]
    @test collect_eachmatch(r"GCG","GCGCG") == ["GCG"]
    @test collect_eachmatch(r"GCG","GCGCG",overlap=true) == ["GCG","GCG"]

    # Issue 8278
    target = """71.163.72.113 - - [30/Jul/2014:16:40:55 -0700] "GET emptymind.org/thevacantwall/wp-content/uploads/2013/02/DSC_006421.jpg HTTP/1.1" 200 492513 "http://images.search.yahoo.com/images/view;_ylt=AwrB8py9gdlTGEwADcSjzbkF;_ylu=X3oDMTI2cGZrZTA5BHNlYwNmcC1leHAEc2xrA2V4cARvaWQDNTA3NTRiMzYzY2E5OTEwNjBiMjc2YWJhMjkxMTEzY2MEZ3BvcwM0BGl0A2Jpbmc-?back=http%3A%2F%2Fus.yhs4.search.yahoo.com%2Fyhs%2Fsearch%3Fei%3DUTF-8%26p%3Dapartheid%2Bwall%2Bin%2Bpalestine%26type%3Dgrvydef%26param1%3D1%26param2%3Dsid%253Db01676f9c26355f014f8a9db87545d61%2526b%253DChrome%2526ip%253D71.163.72.113%2526p%253Dgroovorio%2526x%253DAC811262A746D3CD%2526dt%253DS940%2526f%253D7%2526a%253Dgrv_tuto1_14_30%26hsimp%3Dyhs-fullyhosted_003%26hspart%3Dironsource&w=588&h=387&imgurl=occupiedpalestine.files.wordpress.com%2F2012%2F08%2F5-peeking-through-the-wall.jpg%3Fw%3D588%26h%3D387&rurl=http%3A%2F%2Fwww.stopdebezetting.com%2Fwereldpers%2Fcompare-the-berlin-wall-vs-israel-s-apartheid-wall-in-palestine.html&size=49.0KB&name=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&p=apartheid+wall+in+palestine&oid=50754b363ca991060b276aba291113cc&fr2=&fr=&tt=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&b=0&ni=21&no=4&ts=&tab=organic&sigr=13evdtqdq&sigb=19k7nsjvb&sigi=12o2la1db&sigt=12lia2m0j&sign=12lia2m0j&.crumb=.yUtKgFI6DE&hsimp=yhs-fullyhosted_003&hspart=ironsource" "Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/36.0.1985.125 Safari/537.36"""
    pat = r"""([\d\.]+) ([\w.-]+) ([\w.-]+) (\[.+\]) "([^"\r\n]*|[^"\r\n\[]*\[.+\][^"]+|[^"\r\n]+.[^"]+)" (\d{3}) (\d+|-) ("(?:[^"]|\")+)"? ("(?:[^"]|\")+)"?"""
    match(pat, target)

    # issue #26829
    @test map(m -> m.match, eachmatch(r"^$|\S", "ö")) == ["ö"]

    # issue #26199
    @test map(m -> m.match, eachmatch(r"(\p{L}+)", "Tú")) == ["Tú"]
    @test map(m -> m.match, eachmatch(r"(\p{L}+)", "Tú lees.")) == ["Tú", "lees"]
    @test map(m -> m.match, eachmatch(r"(\p{L}+)", "¿Cuál es tu pregunta?")) == ["Cuál", "es", "tu", "pregunta"]

    # Issue 9545 (32 bit)
    buf = PipeBuffer()
    show(buf, r"")
    @test read(buf, String) == "r\"\""

    # see #10994, #11447: PCRE2 allows NUL chars in the pattern
    @test occursin(Regex("^a\0b\$"), "a\0b")

    # regex match / search string must be a String
    @test_throws ArgumentError match(r"test", GenericString("this is a test"))
    @test_throws ArgumentError findfirst(r"test", GenericString("this is a test"))

    # Issue 27125
    msg = "#Hello# from Julia"
    re = r"#(.+)# from (?<name>\w+)"
    subst = s"FROM: \g<name>\n MESSAGE: \1"
    @test replace(msg, re => subst) == "FROM: Julia\n MESSAGE: Hello"

    # findall
    @test findall(r"\w+", "foo bar") == [1:3, 5:7]
    @test findall(r"\w+", "foo bar", overlap=true) == [1:3, 2:3, 3:3, 5:7, 6:7, 7:7]
    @test all(findall(r"\w*", "foo bar") .=== [1:3, 4:3, 5:7, 8:7]) # use === to compare empty ranges
    @test all(findall(r"\b", "foo bar") .=== [1:0, 4:3, 5:4, 8:7])  # use === to compare empty ranges

    # count
    @test count(r"\w+", "foo bar") == 2
    @test count(r"\w+", "foo bar", overlap=true) == 6
    @test count(r"\w*", "foo bar") == 4
    @test count(r"\b", "foo bar") == 4

    # Named subpatterns
    let m = match(r"(?<a>.)(.)(?<b>.)", "xyz")
        @test (m[:a], m[2], m["b"]) == ("x", "y", "z")
        @test sprint(show, m) == "RegexMatch(\"xyz\", a=\"x\", 2=\"y\", b=\"z\")"
    end

    # Backcapture reference in substitution string
    @test replace("abcde", r"(..)(?P<byname>d)" => s"\g<byname>xy\\\1") == "adxy\\bce"
    @test_throws ErrorException replace("a", r"(?P<x>)" => s"\g<y>")

    # Proper unicode handling
    @test  match(r"∀∀", "∀x∀∀∀").match == "∀∀"

    # 'a' flag to disable UCP
    @test match(r"\w+", "Düsseldorf").match == "Düsseldorf"
    @test match(r"\w+"a, "Düsseldorf").match == "D"

    # Regex behaves like a scalar in broadcasting
    @test occursin.(r"Hello", ["Hello", "World"]) == [true, false]

    @test startswith("abc", r"a")
    @test endswith("abc", r"c")
    @test !startswith("abc", r"b")
    @test !startswith("abc", r"c")
    @test !endswith("abc", r"a")
    @test !endswith("abc", r"b")

    @test !startswith("abc", r"A")
    @test startswith("abc", r"A"i)
    @test !endswith("abc", r"C")
    @test endswith("abc", r"C"i)

    @testset "multiplication & exponentiation" begin
        @test *(r"a") == r"a"

        @test r"a" * r"b" == r"(?:a)(?:b)"
        @test r"a" * "b"  == r"(?:a)\Qb\E"
        @test r"a" * 'b'  == r"(?:a)\Qb\E"
        @test "a"  * r"b" == r"\Qa\E(?:b)"
        @test 'a'  * r"b" == r"\Qa\E(?:b)"
        for a = (r"a", "a", 'a'),
            b = (r"b", "b", 'b'),
            c = (r"c", "c", 'c')
            a isa Regex || b isa Regex || c isa Regex || continue
            @test match(a * b * c, "abc") !== nothing
        end
        for s = ["thiscat", "thishat", "thatcat", "thathat"]
            @test match(r"this|that" * r"cat|hat", s) !== nothing
        end

        @test r"a"i * r"b"i == r"(?:a)(?:b)"i
        @test r"a"i * "b"   == r"(?:a)\Qb\E"i
        @test r"a"i * 'b'   == r"(?:a)\Qb\E"i
        @test "a"   * r"b"i == r"\Qa\E(?:b)"i
        @test 'a'   * r"b"i == r"\Qa\E(?:b)"i

        @test r"a"i  * r"b"m  == r"(?i:a)(?m:b)"
        @test r"a"im * r"b"m  == r"(?i:a)(?:b)"m
        @test r"a"im * r"b"im == r"(?:a)(?:b)"im
        @test r"a"im * r"b"i  == r"(?m:a)(?:b)"i

        r = r"" * raw"a\Eb|c"
        @test match(r, raw"a\Eb|c").match == raw"a\Eb|c"
        @test match(r, raw"c") == nothing

        # error for really incompatible options
        @test_throws ArgumentError r"a" * Regex("b", Base.DEFAULT_COMPILER_OPTS & ~Base.PCRE.UCP, Base.DEFAULT_MATCH_OPTS)
        @test_throws ArgumentError r"a" * Regex("b", Base.DEFAULT_COMPILER_OPTS, Base.DEFAULT_MATCH_OPTS & ~Base.PCRE.NO_UTF_CHECK)

        @test r"this|that"^2 == r"(?:this|that){2}"
    end

    # Test that PCRE throws the correct kind of error
    # TODO: Uncomment this once the corresponding change has propagated to CI
    #@test_throws ErrorException Base.PCRE.info(C_NULL, Base.PCRE.INFO_NAMECOUNT, UInt32)

    # test that we can get the error message of negative error codes
    @test Base.PCRE.err_message(Base.PCRE.ERROR_NOMEMORY) isa String
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
@testset "Special cases" begin
    T = Float64
    f(l) = n -> gamma(l + 1 / 2) / gamma(2l) * gamma(n + 2l) / gamma(n + l + 1 / 2)

    for (P, Q, fn) in (
        (Gegenbauer{1 / 2,T}, Legendre{T}, n -> 1.0),
        (Gegenbauer{1 / 4,T}, Jacobi{1 / 4 - 1 / 2,1 / 4 - 1 / 2,T}, f(1 / 4)),
        (Gegenbauer{3 / 4,T}, Jacobi{3 / 4 - 1 / 2,3 / 4 - 1 / 2,T}, f(3 / 4)),
        (Gegenbauer{1,T}, Jacobi{1 - 1 / 2,1 - 1 / 2,T}, f(1)),
    )
        for i in 1:5
            @test convert(Polynomial, Polynomials.basis(P, i)) ≈
                  fn(i) * convert(Polynomial, Polynomials.basis(Q, i))
        end
    end
end