text:
include("simulate_states.jl")
include("forecast_plots.jl")

# parameters to pass to the main function
const file_date = ARGS[1]

states = [
    "ACT",
    # "NSW",
    # "NT",
    "QLD",
    "SA",
    "TAS",
    # "VIC",
    "WA",
]

# set seed for consistent plots (NOTE: this is not useful when multithreading
# enabled as we use separate seeds but the simulation pool should handle that)
rng = Random.Xoshiro(2022)

jurisdiction_assumptions = JurisdictionAssumptions()

# get the latest onset date
latest_start_date = Dates.Date(
    maximum(v for v in values(jurisdiction_assumptions.simulation_start_dates))
)

(dates, local_case_dict, import_case_dict) = read_in_cases(file_date, rng)
last_date_in_data = dates[end]
forecast_end_date = last_date_in_data + Dates.Day(35)

# create vector for dates
onset_dates = latest_start_date:Dates.Day(1):forecast_end_date

# merge all the simulation and TP files for states into single CSV
merge_simulation_files(file_date)
plot_all_forecast_intervals(file_date, states, dates, local_case_dict)


text:
function get_model(model_type)
    s = Symbol(model_type)
    return getfield(_PM, s)
end

function get_solver(pm)

    optimizer = pm["pm_solver"]
    nl = pm["pm_nl_solver"]
    mip = pm["pm_mip_solver"]
    log_level = pm["pm_log_level"]
    time_limit = pm["pm_time_limit"]
    nl_time_limit = pm["pm_nl_time_limit"]
    mip_time_limit = pm["pm_mip_time_limit"]
    tol = pm["pm_tol"]

    if optimizer == "gurobi"
        solver = JuMP.optimizer_with_attributes(
            Gurobi.Optimizer,
            "TimeLimit" => time_limit,
            "OutputFlag" => log_level,
            "FeasibilityTol" => tol,
            "OptimalityTol" => tol,
        )
    end

    if optimizer == "ipopt"
        solver = JuMP.optimizer_with_attributes(
            Ipopt.Optimizer,
            "print_level" => log_level,
            "max_cpu_time" => time_limit,
            "tol" => tol,
        )
    end

    if optimizer == "juniper" && nl == "ipopt" && mip == "cbc"
        mip_solver = JuMP.optimizer_with_attributes(
            Cbc.Optimizer,
            "logLevel" => log_level,
            "seconds" => mip_time_limit,
        )
        nl_solver = JuMP.optimizer_with_attributes(
            Ipopt.Optimizer,
            "print_level" => log_level,
            "max_cpu_time" => nl_time_limit,
            "tol" => 1e-4,
        )
        solver = JuMP.optimizer_with_attributes(
            Juniper.Optimizer,
            "nl_solver" => nl_solver,
            "mip_solver" => mip_solver,
            "log_levels" => [],
            "time_limit" => time_limit,
        )
    end

    if optimizer == "juniper" && nl == "gurobi" && mip == "cbc"
        mip_solver = JuMP.optimizer_with_attributes(
            Cbc.Optimizer,
            "logLevel" => log_level,
            "seconds" => mip_time_limit,
        )
        nl_solver = JuMP.optimizer_with_attributes(
            Gurobi.Optimizer,
            "TimeLimit" => nl_time_limit,
            "FeasibilityTol" => tol,
            "OptimalityTol" => tol,
        )
        solver = JuMP.optimizer_with_attributes(
            Juniper.Optimizer,
            "nl_solver" => nl_solver,
            "mip_solver" => mip_solver,
            "log_levels" => [],
            "time_limit" => time_limit,
        )
    end

    if optimizer == "juniper" && nl == "gurobi" && mip == "gurobi"
        mip_solver = JuMP.optimizer_with_attributes(
            Gurobi.Optimizer,
            "TimeLimit" => mip_time_limit,
            "FeasibilityTol" => tol,
            "OptimalityTol" => tol,
        )
        nl_solver = JuMP.optimizer_with_attributes(
            Gurobi.Optimizer,
            "TimeLimit" => nl_time_limit,
            "FeasibilityTol" => tol,
            "OptimalityTol" => tol,
        )
        solver = JuMP.optimizer_with_attributes(
            Juniper.Optimizer,
            "nl_solver" => nl_solver,
            "mip_solver" => mip_solver,
            "log_levels" => [],
            "time_limit" => time_limit,
        )
    end

    if optimizer == "knitro"
        solver = JuMP.optimizer_with_attributes(KNITRO.Optimizer, "tol" => tol)
    end

    if optimizer == "cbc"
        solver = JuMP.optimizer_with_attributes(
            Cbc.Optimizer,
            "seconds" => time_limit,
            "tol" => tol,
        )
    end

    if optimizer == "scip"
        solver = JuMP.optimizer_with_attributes(SCIP.Optimizer, "tol" => tol)
    end

    return solver

end

function load_pm_from_json(json_path)
    pm = Dict()
    open(json_path, "r") do f
        pm = JSON.parse(f)
    end
    for (idx, gen) in pm["gen"]
        if gen["model"] == 1
            pm["gen"][idx]["cost"] = convert(Array{Float64,1}, gen["cost"])
        end
    end
    return pm
end

ext: jl


text:
struct AdaptiveMetropolis{d, FT <: AbstractFloat, T <:
    Union{MVector{d,FT}, Vector{FT}}, SST <: StepSize,
    CT <: Cholesky} <: AdaptState
    m::T
    L::CT
    sc::FT
    dx::T
    step::SST
end

"""
    r = AdaptiveMetropolis(x0, [sc, [step]])

Constructor for AdaptiveMetropolis state.

# Arguments
- `x0`: The initial state vector
- `sc`: Scaling parameter; default `2.38/sqrt(d)` where `d` is dimension.
- `step`: Step size object; default `PolynomialStepSize(1.0)`

If `s` is `RWMState`, then proposal samples may be drawn calling
 `draw!(s, r)` and adaptation is performed with `adapt!(r, s)` or
 `adapt_rb!(r, s, α)`.
"""
function AdaptiveMetropolis(x0::T,
                            sc::FT = FT(2.38/sqrt(length(x0))),
                            step = PolynomialStepSize(one(FT))) where {
                            FT <: AbstractFloat, T<:AbstractVector{FT}}
    d = length(x0)
    L = Cholesky(eye(FT,d), :L, 0)
    AdaptiveMetropolis{d,FT,Vector{FT},typeof(step),typeof(L)}(
    Vector{FT}(x0), L, sc, Vector{FT}(x0), step)
end

# Specialisation to StaticArrays
function AdaptiveMetropolis(x0::T,
                            sc::FT = FT(2.38/sqrt(length(x0))),
                            step = PolynomialStepSize(one(FT))) where {d,
                            FT <: AbstractFloat, T<:StaticArray{Tuple{d}, FT}}
    L = Cholesky(MMatrix{d,d}(eye(FT,d)), :L, 0)
    AdaptiveMetropolis{d,FT,MVector{d,FT},typeof(step),typeof(L)}(
      MVector{d,FT}(x0), L, sc, MVector{d,FT}(x0), step)
end

@inline adapt!(sa, sr, α, k) = adapt!(sa, sr, k)
# Update function -- not yet implemented in the best possible way...
@inbounds function adapt!(sa::AdaptiveMetropolis{d},
      sr::RWMState{d}, k::Real) :: Nothing where {d}
    gamma = get(sa.step, k)
    sa.dx .= sr.x .- sa.m
    sa.m .+= gamma.*sa.dx
    rmul!(sa.dx, sqrt(gamma))
    rmul!(sa.L.factors, sqrt(1.0-gamma))
    lowrankupdate!(sa.L, sa.dx)
    nothing
end

@inbounds function draw!(sr::RWMState{d},
      s::AdaptiveMetropolis{d}) :: Nothing where {d}
    #draw!(sr, s.sc[], s.L.L)
    draw!(sr, s.L, s.sc)
    nothing
end

"""
adapt_rb!(s, r, α)

Rao-Blackwellised adaptation step of Adaptive Metropolis as suggested by
Andrieu & Thoms (Statist. Comput. 2008).

# Arguments:
- `s`: AdaptiveMetropolis object
- `r`: RWMState object
- `α`: Acceptance rate

NB: This function should be called *before* calling accept!(r).
"""
function adapt_rb!(sa::AdaptiveMetropolis{d},
      sr::RWMState{d}, α::FT, k::Real) :: Nothing where {d, FT<:AbstractFloat}
    gamma = get(sa.step, k)
    wx = (one(FT)-α)*gamma
    wy = α*gamma
    one_gamma = one(FT) - gamma
    # Update covariance Cholesky factor:
    rmul!(sa.L.factors, sqrt(one_gamma))
    sa.dx .= sr.x .- sa.m
    rmul!(sa.dx, sqrt(wx))
    lowrankupdate!(sa.L, sa.dx)
    sa.dx .= sr.y .- sa.m
    rmul!(sa.dx, sqrt(wy))
    lowrankupdate!(sa.L, sa.dx)
    # Update mean:
    sa.m .= one_gamma.*sa.m .+ wx.*sr.x .+ wy.*sr.y
    nothing
end


text:
# The Dual Arrhenius and Michealis-Menten (DAMM) model, Davidson et al. 2012
function DAMM(x, p)
     # Constants
	R = 8.314472e-3 # Universal gas constant, kJ K⁻¹ mol⁻¹
	O₂ₐ = 0.209 # Volume of O₂ in the air, L L⁻¹
	Dₗᵢ = 3.17 # Diffusion coefficient of substrate in liquid phase, dimensionless
	Dₒₐ = 1.67 # Diffusion coefficient of oxygen in air, dimensionless
	pₛₓ = 0.024 # Fraction of soil carbon that is considered soluble
     # Independent variables
	Tₛ = x[:, 1] # Soil temperature, °C
	θ = x[:, 2] # Soil moisture, m³ m⁻³
     # Parameters
	αₛₓ = p[1] # Pre-exponential factor, mgC cm⁻³ h⁻¹
     	Eaₛₓ = p[2] # Activation energy, kJ mol⁻¹
	kMₛₓ = p[3] # Michaelis constant, gC cm⁻³
    	kMₒ₂ = p[4] # Michaelis constant for O₂, L L⁻¹
	porosity = p[5] # 1 - soil buld density / soil particle density
	Sxₜₒₜ = p[6] # Total soil carbon, gC cm⁻³
     # DAMM model
	Vmax = @. (αₛₓ * exp(-Eaₛₓ/(R * (273.15 + Tₛ)))) # Maximum potential rate of respiration
     	Sₓ = @. pₛₓ * Sxₜₒₜ * Dₗᵢ * θ^3 # All soluble substrate, gC cm⁻³
	MMₛₓ = @. Sₓ / (kMₛₓ + Sₓ) # Availability of substrate factor, 0-1
	O₂ = @. Dₒₐ * O₂ₐ * ((porosity - θ)^(4/3)) # Oxygen concentration
	MMₒ₂ = @. O₂ / (kMₒ₂ + O₂) # Oxygen limitation factor, 0-1
	Resp = @. Vmax * MMₛₓ * MMₒ₂ # Respiration, mg C cm⁻³ hr⁻¹
	Respc = Resp .* 2314.8148 # Respiration, μmol CO₂ m⁻² s⁻¹
end


text:
using LOWDER

# -----------------------------------------------
# Problem definition
# -----------------------------------------------

function f(x)

    return - 0.5 * x[1] + 10 * x[2]

end

function g(x)

    return x[1] ^ 2.0 + x[2] ^ 2.0

end

fmin_list = [f, g]
x = [1.0, 1.0]
a = [0.0, 0.0]
b = [5.0, 5.0]

# -----------------------------------------------
# LOWDER parameters
# -----------------------------------------------

δinit = 1.0
Δinit = 1.2
n_points = 3
verbosity = 3
path_to_file = "./data_files/simple_runtest_01.dat"

# -----------------------------------------------
# LOWDER call
# -----------------------------------------------

sol = LOWDER.lowder(fmin_list, x, a, b; δ = δinit, Δ = Δinit, m = n_points, verbose = verbosity, filename = path_to_file)


text:
mutable struct ProposeContext{T <: Trace, P <: Parameters} <: ExecutionContext
    tr::T
    score::Float64
    visited::Visitor
    params::P
    ProposeContext(tr::T) where T <: Trace = new{T, EmptyParameters}(tr, 0.0, Visitor(), Parameters())
    ProposeContext(tr::T, params::P) where {T <: Trace, P} = new{T, P}(tr, 0.0, Visitor(), params)
end
Propose() = ProposeContext(Trace())
Propose(params) = ProposeContext(Trace(), params)

# ------------ Convenience ------------ #

function propose(fn::Function, args...)
    ctx = Propose()
    ret = ctx(fn, args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.score
end

function propose(params, fn::Function, args...)
    ctx = Propose(params)
    ret = ctx(fn, args...)
    return ret, HierarchicalCallSite(ctx.tr, ctx.score, fn, args, ret), ctx.score
end

function propose(fn::typeof(rand), d::Distribution{K}) where K
    ctx = Propose()
    addr = gensym()
    ret = ctx(fn, addr, d)
    return ret, get_top(ctx.tr, addr), ctx.score
end

function propose(fn::typeof(markov), call::Function, len::Int, args...)
    ctx = Propose()
    addr = gensym()
    ret = ctx(fn, addr, call, len, args...)
    return ret, get_sub(ctx.tr, addr), ctx.score
end

function propose(params, fn::typeof(markov), call::Function, len::Int, args...)
    ctx = Propose(params)
    addr = gensym()
    ret = ctx(fn, addr, call, len, args...)
    return ret, get_sub(ctx.tr, addr), ctx.score
end

function propose(fn::typeof(plate), call::Function, args::Vector)
    ctx = Propose(params)
    addr = gensym()
    ret = ctx(fn, addr, call, args)
    return ret, get_sub(ctx.tr, addr), ctx.score
end

function propose(params, fn::typeof(plate), call::Function, args::Vector)
    ctx = Propose(params)
    addr = gensym()
    ret = ctx(fn, addr, call, args)
    return ret, get_sub(ctx.tr, addr), ctx.score
end

function propose(fn::typeof(plate), d::Distribution{K}, len::Int) where K
    ctx = Propose()
    addr = gensym()
    ret = ctx(fn, addr, d, len)
    return ret, get_sub(ctx.tr, addr), ctx.score
end

# ------------ includes ------------ #

include("hierarchical/propose.jl")
include("plate/propose.jl")
include("markov/propose.jl")
include("factor/propose.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct ProposeContext{T <: Trace} <: ExecutionContext
    tr::T
    score::Float64
    visited::Visitor
    params::LearnableParameters
end
```

`ProposeContext` is used to propose traces for inference algorithms which use custom proposals. `ProposeContext` instances can be passed sets of `LearnableParameters` to configure the propose with parameters which have been learned by differentiable programming.

Inner constructors:

```julia
ProposeContext(tr::T) where T <: Trace = new{T}(tr, 0.0, LearnableParameters())
```

Outer constructors:

```julia
Propose() = ProposeContext(Trace())
```
""", ProposeContext)

@doc(
"""
```julia
ret, g_cl, w = propose(fn::Function, args...)
ret, cs, w = propose(fn::typeof(rand), d::Distribution{K}) where K
ret, v_cl, w = propose(fn::typeof(markov), call::Function, len::Int, args...)
ret, v_cl, w = propose(fn::typeof(plate), call::Function, args::Vector)
ret, v_cl, w = propose(fn::typeof(plate), d::Distribution{K}, len::Int) where K
```

`propose` provides an API to the `ProposeContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, a `RecordSite` instance specialized to the call, and the score `w`.
""", propose)


text:
# ---
# title: 1420. Build Array Where You Can Find The Maximum Exactly K Comparisons
# id: problem1420
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/build-array-where-you-can-find-the-maximum-exactly-k-comparisons/description/>
# hidden: true
# ---
#
# Given three integers `n`, `m` and `k`. Consider the following algorithm to
# find the maximum element of an array of positive integers:
#
# ![](https://assets.leetcode.com/uploads/2020/04/02/e.png)
#
# You should build the array arr which has the following properties:
#
#   * `arr` has exactly `n` integers.
#   * `1 <= arr[i] <= m` where `(0 <= i < n)`.
#   * After applying the mentioned algorithm to `arr`, the value `search_cost` is equal to `k`.
#
# Return _the number of ways_ to build the array `arr` under the mentioned
# conditions. As the answer may grow large, the answer  **must be**  computed
# modulo `10^9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2, m = 3, k = 1
#     Output: 6
#     Explanation: The possible arrays are [1, 1], [2, 1], [2, 2], [3, 1], [3, 2] [3, 3]
#
#
# **Example 2:**
#
#
#
#     Input: n = 5, m = 2, k = 3
#     Output: 0
#     Explanation: There are no possible arrays that satisify the mentioned conditions.
#
#
# **Example 3:**
#
#
#
#     Input: n = 9, m = 1, k = 1
#     Output: 1
#     Explanation: The only possible array is [1, 1, 1, 1, 1, 1, 1, 1, 1]
#
#
# **Example 4:**
#
#
#
#     Input: n = 50, m = 100, k = 25
#     Output: 34549172
#     Explanation: Don't forget to compute the answer modulo 1000000007
#
#
# **Example 5:**
#
#
#
#     Input: n = 37, m = 17, k = 7
#     Output: 418930126
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 50`
#   * `1 <= m <= 100`
#   * `0 <= k <= n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using MPI
using HDF5
using Test

@testset "mpio" begin

MPI.Init()
info = MPI.Info()

h = HDF5.mpihandles[sizeof(MPI.CComm)]
ccomm = reinterpret(h, MPI.CComm(MPI.COMM_WORLD))
cinfo = reinterpret(h, MPI.CInfo(info))

nprocs = MPI.Comm_size(MPI.COMM_WORLD)
myrank = MPI.Comm_rank(MPI.COMM_WORLD)

fileprop = p_create(HDF5.H5P_FILE_ACCESS)
HDF5.h5p_set_fapl_mpio(fileprop, ccomm, cinfo)
h5comm, h5info = HDF5.h5p_get_fapl_mpio(fileprop, h)

# compare the pointer in case of OpenMPI
if h == HDF5.Hmpih64
  c2 = unsafe_load(reinterpret(Ptr{Clong}, ccomm))
  c3 = unsafe_load(reinterpret(Ptr{Clong}, h5comm))
  #
  @test c2 == c3
end

# open file in parallel and write dataset
fn = String(MPI.bcast(collect(tempname()), 0, MPI.COMM_WORLD))
f = h5open(fn, "w", "fapl_mpio", (ccomm, cinfo) )
@test isopen(f)

g = g_create(f, "mygroup")
A = [ myrank + i for i = 1:10]
dset = d_create(g, "B", datatype(Int64), dataspace(10, nprocs), "chunk", (10, 1), "dxpl_mpio", HDF5.H5FD_MPIO_COLLECTIVE)
dset[:,myrank+1] = A
close(f)

MPI.Barrier(MPI.COMM_WORLD)
f = h5open(String(fn), "r", "fapl_mpio", (ccomm, cinfo))
@test isopen(f)
@test names(f) == ["mygroup"]
# read(f, name, pv...) is already taken by multi-dataset read, bool withargs is only dummy argument
B=read(f, "mygroup/B", true, "dxpl_mpio", HDF5.H5FD_MPIO_COLLECTIVE)
@test !isempty(B)
@test A == vec(B[:, myrank + 1])
B=f["mygroup/B", "dxpl_mpio", HDF5.H5FD_MPIO_COLLECTIVE]
@test !isempty(B)
@test A == vec(B[:, myrank + 1])
close(f)

MPI.Barrier(MPI.COMM_WORLD)
B = h5read(fn, "mygroup/B", "fapl_mpio", (ccomm, cinfo), "dxpl_mpio", HDF5.H5FD_MPIO_COLLECTIVE)
@test A == vec(B[:, myrank + 1])
MPI.Barrier(MPI.COMM_WORLD)
B = h5read(fn, "mygroup/B", (:, myrank + 1), "fapl_mpio", (ccomm, cinfo), "dxpl_mpio", HDF5.H5FD_MPIO_COLLECTIVE)
@test A == vec(B)

# we need to close HDF5 and finalize the info object before finalizing MPI
finalize(info)
HDF5.h5_close()
MPI.Barrier(MPI.COMM_WORLD)
MPI.Finalize()

end # testset mpio


text:
## Define the model
## parameters
γ₁, β₁, λ₁, θ₁, κ = 11.0, 1.0, 12.48, 1.6, 2.6
γ₂, β₂, λ₂, θ₂, xstar = 22.0, 1.0, 16.25, 1.0, 1.6

## state space is {11,10,01,00}
T = [
    -(γ₁ + γ₂) γ₂ γ₁ 0;
    β₂ -(γ₁ + β₂) 0 γ₁;
    β₁ 0 -(γ₂ + β₁) γ₂;
    0 β₁ β₂ -(β₂ + β₂);
    ]

C = [λ₁ - θ₁; λ₁ - θ₁; -θ₁; -θ₁]

r₁₁(x) = (λ₂ - κ) * (x .== 0) +
    (λ₂ - θ₂) * ((x .> 0) .& (x .<= xstar)) +
    λ₂ * (x .> xstar)
r₁₀(x) = -κ * (x .== 0) +
    -θ₂ * ((x .> 0) .& (x .<= xstar))
r₀₁(x) = (λ₂ - κ) * (x .== 0) +
    (λ₂ - θ₂) * ((x .> 0) .& (x .<= xstar)) +
    λ₂ * (x .> xstar)
r₀₀(x) = -κ * (x .== 0) +
    -θ₂ * ((x .> 0) .& (x .<= xstar))

R₁₁(x) = x .* (λ₂ - θ₂) .* ((x .> 0) .& (x .< xstar)) .+
    ((x .- xstar) .* λ₂ .+ xstar .* (λ₂ - θ₂)) .* (x .>= xstar)
R₁₀(x) = x .* -θ₂ * ((x .> 0) .& (x .< xstar)) + xstar .* -θ₂ * (x .>= xstar)
R₀₁(x) = x .* (λ₂ - θ₂) .* ((x .> 0) .& (x .< xstar)) .+
    ((x .- xstar) .* λ₂ .+ xstar .* (λ₂ - θ₂)) .* (x .>= xstar)
R₀₀(x) = x .* -θ₂ * ((x .> 0) .& (x .< xstar)) + xstar .* -θ₂ * (x .>= xstar)

r = (
    r = function (x)
        [r₁₁(x) r₁₀(x) r₀₁(x) r₀₀(x)]
    end,
    R = function (x)
        [R₁₁(x) R₁₀(x) R₀₁(x) R₀₀(x)]
    end,
)

## make models
simBounds = [0 Inf; -Inf Inf] # bounds for simulation only
simModel = SFFM.Model( T, C, r, Bounds = simBounds)
println("created simModel with no upper bound")
println("")

approxBounds = [0 48; -Inf Inf] # bounds for approximation only
approxModel = SFFM.Model( T, C, r, Bounds = approxBounds)
println("Created approxModel with upper bound at x=", approxModel.Bounds[1,end])
println("")

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
export minibatch, Data
import Base: iterate, eltype, length, rand, repeat, summary, show
using Base.Iterators: Cycle
using Base: @propagate_inbounds, tail
using Random: randperm
using Knet: atype

"Iterator of minibatches (x,y pairs) returned by `minibatch`."
mutable struct Data{T}; x; y; batchsize; length; partial; imax; indices; shuffle; xsize; ysize; xtype; ytype; end

"""
    minibatch(x, [y], batchsize; shuffle, partial, xtype, ytype, xsize, ysize)

Return an iterator of minibatches [(xi,yi)...] given data tensors x, y and batchsize.

The last dimension of x and y give the number of instances and should be equal. `y` is
optional, if omitted a sequence of `xi` will be generated rather than `(xi,yi)` tuples.  Use
`repeat(d,n)` for multiple epochs, `Iterators.take(d,n)` for a partial epoch, and
`Iterators.cycle(d)` to cycle through the data forever (this can be used with `converge`).
If you need the iterator to continue from its last position when stopped early (e.g. by a
break in a for loop), use `Iterators.Stateful(d)` (by default the iterator would restart
from the beginning).

Keyword arguments:

- `shuffle=false`: Shuffle the instances every epoch.
- `partial=false`: If true include the last partial minibatch < batchsize.
- `xtype=typeof(x)`: Convert xi in minibatches to this type.
- `ytype=typeof(y)`: Convert yi in minibatches to this type.
- `xsize=size(x)`: Convert xi in minibatches to this shape (with last dimension adjusted for batchsize).
- `ysize=size(y)`: Convert yi in minibatches to this shape (with last dimension adjusted for batchsize).
"""
minibatch

function minibatch(x,y,batchsize; shuffle=false,partial=false,xsize=size(x),ysize=size(y),
                   # default xtype, ytype should be robust to ndims change:
                   xtype = (eltype(x) <: AbstractFloat ? atype() : (typeof(x).name.wrapper){eltype(x)}),
                   ytype = (eltype(y) <: AbstractFloat ? atype() : (typeof(y).name.wrapper){eltype(y)}))
    nx = size(x)[end]
    if nx != size(y)[end]; throw(DimensionMismatch()); end
    x2 = reshape(x, :, nx)
    y2 = reshape(y, :, nx)
    imax = partial ? nx : nx - batchsize + 1
    Data{Tuple{xtype,ytype}}(x2,y2,batchsize,nx,partial,imax,1:nx,shuffle,xsize,ysize,xtype,ytype)
end

function minibatch(x,batchsize; shuffle=false,partial=false,xsize=size(x),
                   xtype = (eltype(x) <: AbstractFloat ? atype() : (typeof(x).name.wrapper){eltype(x)}))
    nx = size(x)[end]
    x2 = reshape(x, :, nx)
    imax = partial ? nx : nx - batchsize + 1
    Data{xtype}(x2,nothing,batchsize,nx,partial,imax,1:nx,shuffle,xsize,nothing,xtype,nothing)
end

@propagate_inbounds function iterate(d::Data, i=0)     # returns data in d.indices[i+1:i+batchsize]
    if i >= d.imax
        return nothing
    end
    if d.shuffle && i == 0
        d.indices = randperm(d.length)
    end
    nexti = min(i + d.batchsize, d.length)
    ids = d.indices[i+1:nexti]
    xbatch = try convert(d.xtype, reshape(d.x[:,ids],d.xsize[1:end-1]...,length(ids)))
    catch e; throw(DimensionMismatch("Error converting batch from $(summary(d.x)) to $(d.xtype): $e")); end
    if d.y == nothing
        return (xbatch,nexti)
    else
        ybatch = try convert(d.ytype, reshape(d.y[:,ids],d.ysize[1:end-1]...,length(ids)))
        catch; throw(DimensionMismatch("Y tensor not compatible with size=$(d.ysize) and type=$(d.ytype)")); end
        return ((xbatch,ybatch),nexti)
    end
end

eltype(::Type{Data{T}}) where T = T

function length(d::Data)
    n = d.length / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

function rand(d::Data)
    i = rand(0:(d.length-d.batchsize))
    return iterate(d, i)[1]
end

# IterTools.ncycle(data,n) for multiple epochs
# Base.Iterators.cycle(data) to go forever
# Base.Iterators.take(data,n) for partial epochs
# IterTools.takenth(itr,n) to report every n iterations

struct Repeat; data::Data; n::Int; end

function repeat(d::Data, n::Int)
    @warn "repeat(d::Data,n) is deprecated, use IterTools.ncycle instead." maxlog=1
    @assert n >= 0
    Repeat(d,n)
end

length(r::Repeat) = r.n * length(r.data)
eltype(r::Repeat) = eltype(r.data)
eltype(c::Cycle{Data}) = eltype(c.xs)
eltype(c::Cycle{Repeat}) = eltype(c.xs)

@propagate_inbounds function iterate(r::Repeat, s=(1,))
    epoch, state = s[1], tail(s)
    epoch > r.n && return nothing
    next = iterate(r.data, state...)
    next === nothing && return iterate(r, (epoch+1,))
    (next[1], (epoch, next[2]))
end

# Give length info in summary:
summary(d::Data) = "$(length(d))-element $(typeof(d))"
show(io::IO, d::Data) = print(IOContext(io,:compact=>true), summary(d))
show(io::IO, ::MIME"text/plain", d::Data) = show(io, d)


text:
struct BetterFeaturization <: SeaPearl.AbstractFeaturization end

function SeaPearl.featurize(sr::SeaPearl.DefaultStateRepresentation{BetterFeaturization,TS}) where TS
    g = sr.cplayergraph
    features = zeros(Float32, coloring_generator.n+6, nv(g))
    for i in 1:nv(g)
        cp_vertex = SeaPearl.cpVertexFromIndex(g, i)
        if isa(cp_vertex, SeaPearl.VariableVertex)
            features[1,i] = 1.
            if g.cpmodel.objective == cp_vertex.variable
                features[6, i] = 1.
            end
        end
        if isa(cp_vertex, SeaPearl.ConstraintVertex)
            features[2, i] = 1.
            constraint = cp_vertex.constraint
            if isa(constraint, SeaPearl.NotEqual)
                features[4, i] = 1.
            end
            if isa(constraint, SeaPearl.LessOrEqual)
                features[5, i] = 1.
            end
        end
        if isa(cp_vertex, SeaPearl.ValueVertex)
            features[3, i] = 1.
            value = cp_vertex.value
            features[6+value, i] = 1.
        end
    end
    features
end

function SeaPearl.feature_length(::Type{SeaPearl.DefaultStateRepresentation{BetterFeaturization, TS}}) where TS
    coloring_generator.n+6
end

function SeaPearl.global_feature_length(::Type{SeaPearl.DefaultStateRepresentation{BetterFeaturization, TS}}) where TS
    return 0
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
module gen_points

using Clustering
using Distances
using Distributions
export generate_points,
       get_centroids

function generate_points()
    #1. Define a large circle
    #2. Sample 'center' points from the large circle
    #3. For each 'center' points, sample a disk around it

    # generate 'center' points
    n_center = 20
    theta = rand(Uniform(0,2*pi), n_center)
    distance = rand(Normal(2, 0.2), n_center);

    center = zeros(n_center,2)
    for i=1:n_center
        center[i, 1] = distance[i] * cos(theta[i])
        center[i, 2] = distance[i] * sin(theta[i])
    end

    # define the covariance matrix for 2D Gaussian
    sigma = Array{Float64}(undef, 2, 2)
    sigma[1,1] = 0.05
    sigma[1, 2] = 0
    sigma[2, 1] = 0
    sigma[2, 2] = 0.05

    # sample points around centers
    n_samples = 10
    points = Matrix{Float64}(undef, 2, 0)

    for i = 1:n_center
        sampled_points = rand(MvNormal(center[i,:], sigma), n_samples)
        points = hcat(points, sampled_points)
    end
    return points

end

function get_centroids(points, n_clusters)

    # compute distance matrix
    D = pairwise(Euclidean(), points, points, dims=2);

    # clustering
    clustering = hclust(D, linkage = :single)
    clusters = cutree(clustering, k = n_clusters)

    # get centroids of each cluster
    centroids = Matrix{Float64}(undef, 2,0)
    for i =1:n_clusters
        index = findall(x -> x == i, clusters)
        c = mean(points[:, index], dims = 2)
        centroids = hcat(centroids, c)
    end

    return clusters, centroids
end

function compute_distances(points, centroids)
    # total points
    P = hcat(points, centroids)
    n_points = size(points)[2]
    n_centroids = size(centroids)[2]

    # compute distances
    D = pairwise(Euclidean(), P, P, dims=2)

    # Define submatrices
    D_R1 = D[1:n_points, 1:n_points]
    D_R2 = D[n_points+1:end, n_points+1:end]
    D_R1_R2 = D[1:n_points, n_points+1:end]
        # rows (landmarks): R1
        # columns (witness) : R2
    D_R2_R1 = D[n_points+1:end, 1:n_points];
        # rows (landmarks): R2
        # columns (witness) : R1
    return D
end

end


text:
module IntervalModel

using Common
# using NLopt
using Ipopt

function init_intervals(log_size::Float64, l0_count::Float64, level_count::Int64)
	interval = Array(Float64, 0)
	for i = 1:level_count
		# push!(interval, log_size * l0_count)
		push!(interval, log_size * l0_count * (10. ^ Float64(i - 1)))
	end
	interval
end

function calculate_wa_twolevel!(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r::Array{Float64}, wa_w::Array{Float64})
	# mem->log
	wa_r[1] = 0.
	wa_w[1] = 1.

	# log->0
	wa_r[2] = 0.
	wa_w[2] = Common.unique(X, log_size) / log_size

	# ## amortized, full destination level
	# # 0->1, 1->2, ...
	# for i in 1:(length(intervals) - 1)
	#	wa[2 + i] = Common.unique(X, intervals[i] + intervals[i + 1]) / intervals[i]
	# end
	# wa[2 + length(intervals)] = Float64(X.count) / intervals[end]

	# ## amortized, compact entire level
	# # 0->1, 1->2, ...
	# interval = 0.
	# next_interval = intervals[1]
	# for i in 1:(length(intervals) - 1)
	#	interval = interval * 0.5 + next_interval
	#	next_interval = intervals[i + 1]
	#	wa[2 + i] = unique_avg(X, interval, interval * 0.5 + next_interval) / interval
	# end
	# interval = interval * 0.5 + next_interval
	# wa[2 + length(intervals)] = Float64(X.count) / interval

	## deamortized, compact each sstable in a round-robin way
	# 0->1, 1->2, ...
	interval = 0.
	next_interval = intervals[1]
	for i in 1:(length(intervals) - 1)
		if i == 1
			# 0->1 compaction is usually a whole level
			# do not use interval_from_density() and adding extra unique() to WA that are caused by using small tables
			interval = next_interval
			next_interval = intervals[i + 1]
			wa_r[2 + i] = (Common.unique(X, log_size) * l0_count + Common.unique(X, next_interval)) / interval
			wa_w[2 + i] = Common.unique(X, interval + next_interval) / interval
		else
			interval = interval + interval_from_density(X, Common.unique(X, next_interval))
			next_interval = intervals[i + 1]
			# using additional unique(); see SizeMode.jl for details
			wa_r[2 + i] = (Common.unique(X, interval) + Common.unique(X, next_interval) + Common.unique(X, interval) * 1.) / interval
			wa_w[2 + i] = (Common.unique(X, interval + next_interval) + Common.unique(X, interval) * 1.) / interval
		end
	end
	interval = interval + interval_from_density(X, Common.unique(X, next_interval))
	# using additional unique(); see SizeMode.jl for details
	wa_r[2 + length(intervals)] = (Common.unique(X, interval) + Float64(X.count) + Common.unique(X, interval) * 1.) / interval
	wa_w[2 + length(intervals)] = (Float64(X.count) + Common.unique(X, interval) * 1.) / interval

	wa_r, wa_w
end

function calculate_wa_twolevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64})
	wa_r = Array(Float64, 2 + length(intervals))
	wa_w = Array(Float64, 2 + length(intervals))
	calculate_wa_twolevel!(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r, wa_w)
end

# function calculate_wa_twolevel_ratios(X::Distribution, log_size::Float64, l0_count::Float64, interval_ratios::Array{Float64})
#	current_interval = interval_ratios[1]
#	intervals = interval_ratios * (log_size * l0_count / current_interval)
#	return calculate_wa_twolevel(X, log_size, l0_count, intervals)
# end

function calculate_sizes_twolevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64})
	sizes = Array(Float64, length(intervals))

	for i in 1:(length(intervals) - 1)
		sizes[i] = Common.unique(X, intervals[i + 1])
	end
	sizes[length(intervals)] = Float64(X.count)

	sizes
end

function calculate_wa_multilevel!(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r::Array{Float64}, wa_w::Array{Float64})
	# TODO: wa_r

	# mem->log
	wa_r[1] = 0.
	wa_w[1] = 1.

	# log->0
	wa_r[2] = 0.
	wa_w[2] = Common.unique(X, log_size) / log_size

	# ## amortized, full destination level
	# # 0->1, 1->2, ...
	# for i in 1:(length(intervals) - 1)
	#	# level-0...1 size, level-0...2 size, ...
	#	level_size = Common.unique(X, geom_mean(intervals[(i + 1):end]))
	#	wa[2 + i] = level_size / intervals[i]
	# end
	# wa[2 + length(intervals)] = Float64(X.count) / intervals[end]

	## amortized, compact entire level (TODO: do we need to modify interval to consider "0.5" factor?)
	# 0->1, 1->2, ...
	# interval = 0.
	# next_interval = geom_mean(intervals)
	for i in 1:(length(intervals) - 1)
		wa_r[2 + i] = 0.
		wa_w[2 + i] = unique_avg(X, geom_mean(intervals[i:end]), geom_mean(intervals[i:end]) * 0.5 + geom_mean(intervals[(i + 1):end])) / intervals[i]
		# interval = interval * 0.5 + next_interval
		# next_interval = geom_mean(intervals[(i + 1):end])
		# wa[2 + i] = unique_avg(X, interval, interval * 0.5 + next_interval) / interval
	end
	wa_r[2 + length(intervals)] = 0.
	wa_w[2 + length(intervals)] = Float64(X.count) / intervals[end]

	wa_r, wa_w
end

function calculate_wa_multilevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64})
	wa_r = Array(Float64, 2 + length(intervals))
	wa_w = Array(Float64, 2 + length(intervals))
	calculate_wa_multilevel!(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r, wa_w)
end

# function calculate_wa_multilevel_ratios(X::Distribution, log_size::Float64, l0_count::Float64, interval_ratios::Array{Float64})
#	current_interval = geom_mean(interval_ratios)
#	intervals = interval_ratios * (log_size * l0_count / current_interval)
#	return calculate_wa_multilevel(X, log_size, l0_count, intervals)
# end

function calculate_sizes_multilevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64})
	sizes = Array(Float64, length(intervals))

	for i in 1:(length(intervals) - 1)
		sizes[i] = Common.unique(X, geom_mean(intervals[i + 1:end]))
	end
	sizes[length(intervals)] = Float64(X.count)

	sizes
end

function optimize_wa_twolevel(X::Distribution, log_size::Float64, l0_count::Float64, init_intervals::Array{Float64}, wa_r_factor::Float64, ftol::Float64, max_time::Float64)
	n = X.count
	level_count = length(init_intervals)

	# v2 = Array(Float64, level_count)
	# v2[1] = log_size * l0_counum

	# count = 0
	# wa_r = Array(Float64, 2 + level_count)
	# wa_w = Array(Float64, 2 + level_count)
	# f = (v, grad) -> begin
	#	count += 1
	#	v2[2:level_count] = v
	#	get_wa(wa_r_factor, calculate_wa_twolevel!(X, log_size, l0_count, v2, wa_r, wa_w))
	# end

	# v = init_intervals[2:end]

	# opt = Opt(:LN_COBYLA, level_count - 1)
	# min_objective!(opt, f)
	# # inequality_constraint!(opt, (v, grad) -> log_size * l0_count - v[1])		# <= 0
	# for i = 1:(level_count - 2)
	#	inequality_constraint!(opt, (v, grad) -> v[i] - v[i + 1])		# <= 0
	# end
	# ftol_abs!(opt, ftol)
	# maxtime!(opt, max_time)
	# @time (minf, minx, ret) = optimize(opt, v)
	# println("got $minf at $minx after $count iterations (returned $ret)")

	# cat(1, [log_size * l0_count], minx)

	#######################

	v2 = Array(Float64, level_count)
	v2[1] = log_size * l0_count

	count = 0
	wa_r = Array(Float64, 2 + level_count)
	wa_w = Array(Float64, 2 + level_count)

	eval_f = (v) -> begin
		count += 1
		v2[2:level_count] = v
		get_wa(wa_r_factor, calculate_wa_twolevel!(X, log_size, l0_count, v2, wa_r, wa_w))
	end

	eval_grad_f = (v, grad_f) -> begin
		v2[2:level_count] = v
		y = get_wa(wa_r_factor, calculate_wa_twolevel!(X, log_size, l0_count, v2, wa_r, wa_w))
		for i = 2:level_count
			diff = max(v2[i] * 0.001, 1.)
			org = v2[i]
			v2[i] += diff
			grad_f[i - 1] = (get_wa(wa_r_factor, calculate_wa_twolevel!(X, log_size, l0_count, v2, wa_r, wa_w)) - y) / diff
			v2[i] = org
		end
	end

	eval_g = (v, g) -> begin
		for i = 1:(level_count - 2)
			g[i] = v[i] - v[i + 1]
		end
	end

	# level i's interval - level i+1's interval <= 0
	eval_jac_g = (v, mode, rows, cols, values) -> begin
		if mode == :Structure
			c = 1
			for i = 1:level_count - 2
				rows[c] = i
				cols[c] = i
				c += 1
				rows[c] = i
				cols[c] = i + 1
				c += 1
			end
		else
			c = 1
			for i = 1:level_count - 2
				values[c] = 1.
				c += 1
				values[c] = -1.
				c += 1
			end
		end
	end

	v_L = [log_size * l0_count for i = 1:level_count - 1]
	v_U = [2.e19 for i = 1:level_count - 1]

	g_L = [-2.e19 for i = 1:level_count - 2]
	g_U = [0. for i = 1:level_count - 2]

	prob = createProblem(level_count - 1, v_L, v_U,
						 level_count - 2, g_L, g_U,
						 (level_count - 2) * 2, 0,
						 eval_f, eval_g, eval_grad_f, eval_jac_g)

	addOption(prob, "hessian_approximation", "limited-memory")

	addOption(prob, "tol", ftol)
	addOption(prob, "max_cpu_time", max_time)
	addOption(prob, "acceptable_iter", 1000)

	addOption(prob, "print_level", 2)

	prob.x = init_intervals[2:end]

	@time status = solveProblem(prob)

	ret = Ipopt.ApplicationReturnStatus[status]
	minf = prob.obj_val
	minx = prob.x
	println("got $minf at $minx after $count iterations (returned $ret)")

	cat(1, [log_size * l0_count], minx)
end

function optimize_wa_multilevel(X::Distribution, log_size::Float64, l0_count::Float64, init_intervals::Array{Float64}, wa_r_factor::Float64, ftol::Float64, max_time::Float64)
	n = X.count
	level_count = length(init_intervals)

	# v2 = Array(Float64, level_count)

	# count = 0
	# wa_r = Array(Float64, 2 + level_count)
	# wa_w = Array(Float64, 2 + level_count)
	# f = (v, grad) -> begin
	#	count += 1
	#	# we need to make geom_mean(cat(1, [X], v)) = log_size * l0_count
	#	# 1/X + ..  = 1 / (log_size * l0_count)
	#	# 1/X = 1 / (log_size * l0_count) - ...
	#	# X = 1 / (1 / (log_size * l0_count) - ...)
	#	#   = geom_mean(cat(1, [log_size * l0_count], -v))

	#	# v2[1] = geom_mean(cat(1, [log_size * l0_count], -v))

	#	v2[1] = -(log_size * l0_count)
	#	v2[2:level_count] = v
	#	v2[1] = -geom_mean(v2)
	#	get_wa(wa_r_factor, calculate_wa_multilevel!(X, log_size, l0_count, v2, wa_r, wa_w))
	# end

	# v = init_intervals[2:end]

	# opt = Opt(:LN_COBYLA, level_count - 1)
	# min_objective!(opt, f)
	# for i = 1:(level_count - 2)
	#	inequality_constraint!(opt, (v, grad) -> v[i] - v[i + 1])		# <= 0
	# end
	# ftol_abs!(opt, ftol)
	# maxtime!(opt, max_time)
	# @time (minf, minx, ret) = optimize(opt, v)
	# println("got $minf at $minx after $count iterations (returned $ret)")

	# x = geom_mean(cat(1, [log_size * l0_count], -minx))
	# cat(1, [x], minx)

	#######################

	v2 = Array(Float64, level_count)
	v2[1] = log_size * l0_count

	count = 0
	wa_r = Array(Float64, 2 + level_count)
	wa_w = Array(Float64, 2 + level_count)

	eval_f = (v) -> begin
		count += 1
		v2[1] = -(log_size * l0_count)
		v2[2:level_count] = v
		v2[1] = -geom_mean(v2)
		# note that v2[1] can become negative accidentally, which is not valid for unique()
		get_wa(wa_r_factor, calculate_wa_multilevel!(X, log_size, l0_count, v2, wa_r, wa_w))
	end

	eval_grad_f = (v, grad_f) -> begin
		v2[1] = -(log_size * l0_count)
		v2[2:level_count] = v
		v2[1] = -geom_mean(v2)
		y = get_wa(wa_r_factor, calculate_wa_multilevel!(X, log_size, l0_count, v2, wa_r, wa_w))
		for i = 2:level_count
			diff = max(v2[i] * 0.001, 1.)
			org = v2[i]
			v2[i] += diff
			v2[1] = -(log_size * l0_count)
			v2[1] = -geom_mean(v2)
			grad_f[i - 1] = (get_wa(wa_r_factor, calculate_wa_multilevel!(X, log_size, l0_count, v2, wa_r, wa_w)) - y) / diff
			v2[i] = org
		end
	end

	eval_g = (v, g) -> begin
		for i = 1:(level_count - 2)
			g[i] = v[i] - v[i + 1]
		end
	end

	# level i's interval - level i+1's interval <= 0
	eval_jac_g = (v, mode, rows, cols, values) -> begin
		if mode == :Structure
			c = 1
			for i = 1:level_count - 2
				rows[c] = i
				cols[c] = i
				c += 1
				rows[c] = i
				cols[c] = i + 1
				c += 1
			end
		else
			c = 1
			for i = 1:level_count - 2
				values[c] = 1
				c += 1
				values[c] = -1
				c += 1
			end
		end
	end

	v_L = [log_size * l0_count for i = 1:level_count - 1]
	v_U = [2.e19 for i = 1:level_count - 1]

	g_L = [-2.e19 for i = 1:level_count - 2]
	g_U = [0. for i = 1:level_count - 2]

	prob = createProblem(level_count - 1, v_L, v_U,
						 level_count - 2, g_L, g_U,
						 (level_count - 2) * 2, 0,
						 eval_f, eval_g, eval_grad_f, eval_jac_g)

	addOption(prob, "hessian_approximation", "limited-memory")

	addOption(prob, "tol", ftol)
	addOption(prob, "max_cpu_time", max_time)
	addOption(prob, "acceptable_iter", 1000)

	addOption(prob, "print_level", 2)

	prob.x = init_intervals[2:end]

	@time status = solveProblem(prob)

	ret = Ipopt.ApplicationReturnStatus[status]
	minf = prob.obj_val
	minx = prob.x
	println("got $minf at $minx after $count iterations (returned $ret)")

	x = geom_mean(cat(1, [log_size * l0_count], -minx))
	cat(1, [x], minx)
end

function print_twolevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r_factor::Float64)
	level_count = length(intervals)

	println("intervals = ", [iround(v) for v in intervals])
	println("exp. size = ", [iround(Common.unique(X, v)) for v in intervals])
	println("(", [round(intervals[i] / intervals[i - 1] * 100.) / 100. for i in 2:length(intervals)], " X)")
	wa = calculate_wa_twolevel(X, log_size, l0_count, intervals)
	println("WA (mem->log) = ", wa[2][1])
	println("WA (log->0) = ", wa[2][2])
	for i = 1:level_count; println("WA ($(i-1)->$i) = ", wa[2][i + 2]) end
	println("WA = ", get_wa(wa_r_factor, wa))
end

function print_multilevel(X::Distribution, log_size::Float64, l0_count::Float64, intervals::Array{Float64}, wa_r_factor::Float64)
	level_count = length(intervals)

	println("intervals = ", [iround(v) for v in intervals])
	println("(", [round(intervals[i] / intervals[i - 1] * 100.) / 100. for i in 2:length(intervals)], " X)")
	wa = calculate_wa_multilevel(X, log_size, l0_count, intervals)
	println("avg L0 intervals = ", iround(geom_mean(intervals)))
	println("WA (mem->log) = ", wa[2][1])
	println("WA (log->0) = ", wa[2][2])
	for i = 1:level_count; println("WA ($(i-1)->$i) = ", wa[2][i + 2]) end
	println("WA = ", get_wa(wa_r_factor, wa))
end

end

max_stars_repo_name: efficient/msls-eval

max_issues_repo_name: theopengroup/msls-eval

max_forks_repo_name: theopengroup/msls-eval


text:
# Implicitly solve heat equation: ∂_t T = α * ∇²(T) + β*sin(γ*z)
# Boundary conditions: ∇T(z=a) = ∇T_bottom, T(z=b) = T_top
import Plots
using LinearAlgebra
using DiffEqBase
using OrdinaryDiffEq: SplitODEProblem, solve, IMEXEuler
import SciMLBase

a,b, n = 0, 1, 10               # zmin, zmax, number of cells
n̂_min, n̂_max = -1, 1            # Outward facing unit vectors
α = 100;                        # thermal diffusivity, larger means more stiff
β, γ = 10000, π;                # source term coefficients
Δt = 1000;                      # timestep
N_t = 10;                       # number of timesteps
FT = Float64;                   # float type
Δz = FT(b-a)/FT(n)
Δz² = Δz^2;
∇²_op = [1/Δz², -2/Δz², 1/Δz²]; # row or Laplacian operator
∇T_bottom = 10;                 # Temperature gradient at the top
T_top = 1;                      # Temperature at the bottom
S(z) = β*sin(γ*z)               # source term, (sin for easy integration)
zf = range(a, b, length=n+1);
zc = map(i-> zf[i]+Δz/2, 1:n);

function T_analytic(z) # Analytic steady state solution
    c1 = ∇T_bottom-β*cos(γ*a)/(γ*α)
    c2 = T_top-(β*sin(γ*b)/(γ^2*α)+c1*b)
    return β*sin(γ*z)/(γ^2*α)+c1*z+c2
end

T = zeros(FT, n+1);

# Equations: derivation to matrix form
# ∂_t T = α * ∇²(T) + S
# (T^{n+1}-T^n) = Δt (α * ∇²(T^{n+1}) + S)
# (T^{n+1} - Δt α * ∇²(T^{n+1})) = T^n + Δt*S
# (I - Δt α * ∇²) T^{n+1} = T^n + Δt*S

# Derive Dirichlet boundary stencil & source:
# ∂_t T = α * (T[i-1]+T[b]-2*T[i])/Δz² + S
# ∂_t T = α * (T[i-1]-2*T[i])/Δz² + S + α * T[b] / Δz²

# Derive Neumann boundary stencil & source:
# ∇T_bottom*n̂= (T[g] - T[i])/(2Δz),     n̂ = [-1,1] ∈ [zmin,zmax]
# T[i] + 2*Δz*∇T_bottom*n̂ = T[g]
# ∂_t T = α * (((T[i] + 2*Δz*∇T_bottom*n̂) - T[b])/Δz - (T[b] - T[i])/Δz)/Δz + S
# ∂_t T = α * (((T[i]) - T[b])/Δz - (T[b] - T[i])/Δz)/Δz + S + α/Δz²*2*Δz*∇T_bottom*n̂
# ∂_t T = α * (2*T[i] - 2*T[b])/Δz² + S + 2α/Δz*∇T_bottom*n̂
∇² = Tridiagonal(
    ones(FT, n) .* ∇²_op[1],
    ones(FT, n+1)   .* ∇²_op[2],
    ones(FT, n) .* ∇²_op[3]
);

# Modify boundary stencil to account for BCs
∇².d[1] = -2/Δz²
∇².du[1] = +2/Δz²

# Modify boundary stencil to account for BCs
∇².du[n] = 0  # modified stencil
∇².d[n+1] = 0 # to ensure `∂_t T = 0` at `z=zmax`
∇².dl[n] = 0  # to ensure `∂_t T = 0` at `z=zmax`

# Compute boundary source: α * T[b] / Δz²
AT_b = zeros(FT, n+1);
AT_b[1] = α*2/Δz*∇T_bottom*n̂_min;
AT_b[end-1] = α*T_top/Δz²;

# Set initial condition:
T .= 1;
T[n+1] = T_top # set top BC

function rhs!(dT, T, params, t)
    n = params.n
    i = 1:n # interior domain
    dT[i] .= S.(zf[i]) .+ AT_b[i]
    return dT
end

D = α .* ∇²

params = (;n)

tspan = (FT(0), N_t*FT(Δt))

prob = SplitODEProblem(
    SciMLBase.DiffEqArrayOperator(
        D,
    ),
    rhs!,
    T,
    tspan,
    params
)

println("Solving...")

alg = IMEXEuler(linsolve=LinSolveFactorize(lu!))

sol = solve(
    prob,
    alg,
    dt = Δt,
    saveat = range(FT(0), N_t*FT(Δt), length=5),
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

T_end = sol.u[end]

p1 = Plots.plot(zf, T_analytic.(zf), label="analytic", markershape=:circle, markersize=6)
p1 = Plots.plot!(p1, zf, T_end, label="numerical", markershape=:diamond)
p1 = Plots.plot!(p1, title="T ∈ cell faces")

p2 = Plots.plot(zf, abs.(T_end .- T_analytic.(zf)), label="error", markershape=:circle, markersize=6)
p2 = Plots.plot!(p2, title="T ∈ cell faces")

Plots.plot(p1, p2)


text:
module IsConstFoldable_tests

using Test

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

end


text:
# Contains all the structures used and define the 1D and 2D types of the problem

# Define the 2 SNFE problems types for 1D and 2D cases
abstract type Problem end
abstract type OneDim <: Problem end
abstract type TwoDim <: Problem end

include("DomainConstruct.jl")  # File containing the structure of our problem domain
include("ProbStruct.jl")       # File containing the structures of probSNFE (input and output)
include("SolveStruct.jl")      # File containing the structures of solveSNFE


text:
export voidLayer

type voidLayer<:Layer
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# Port of the autocorrelation length estimation of the emcee Python package,
# under MIT License. Original authors Dan Foreman-Mackey et al.
#
# See also:
#
# https://dfm.io/posts/autocorr/
# https://github.com/dfm/emcee/
# https://github.com/dfm/emcee/issues/209

function _ac_next_pow_two(i::Integer)
    nbits = 8 * sizeof(i)
    1 << (nbits - leading_zeros(i > 0 ? i - 1 : i))
end

"""
    fft_autocov(v::AbstractVector{<:Real})
    fft_autocov(v::AbstractVectorOfSimilarVectors{<:Real})

Compute the autocovariance of of variate series `v`, separately for each
degree of freedom.

Uses FFT, in contract to `StatsBase.autocov`.
"""
function fft_autocov end

function fft_autocov(x::AbstractVector{<:Real})
    n = length(eachindex(x))
    n2 = 2 * _ac_next_pow_two(n)
    x2 = zeros(eltype(x), n2)
    idxs2 = firstindex(x2):(firstindex(x2) + n - 1)
    x2_view = view(x2, idxs2)
    x2_view .= x .- mean(x)

    x2_fft = rfft(x2)
    x2_fft .= abs2.(x2_fft)
    acf = irfft(x2_fft, size(x2, 1))[idxs2]

    acf .*= inv(n)

    acf
end

function fft_autocov(v::AbstractVectorOfSimilarVectors{<:Real})
    X = flatview(v)
    n = size(X, 2)
    n2 = 2 * _ac_next_pow_two(n)
    X2 = zeros(eltype(X), size(X, 1), n2)
    idxs = axes(X2,2)
    idxs2 = first(idxs):(first(idxs) + n - 1)
    X2_view = view(X2, :, idxs2)
    X2_view .= X .- mean(X, dims = 2)

    X2_fft = rfft(X2, 2)
    X2_fft .= abs2.(X2_fft)
    acf = irfft(X2_fft, size(X2, 2), 2)[:, idxs2]

    acf .*= inv(n)

    nestedview(acf)
end

"""
    fft_autocor(v::AbstractVector{<:Real})
    fft_autocor(v::AbstractVectorOfSimilarVectors{<:Real})

Compute the autocorrelation function (ACF) of variate series `v`, separately
for each degree of freedom.

Uses FFT, in contract to `StatsBase.autocor`.
"""
function fft_autocor end

function fft_autocor(v::AbstractVector{<:Real})
    acf = fft_autocov(v)
    acf ./= first(acf)
    acf
end

function fft_autocor(v::AbstractVectorOfSimilarVectors{<:Real})
    acf = flatview(fft_autocov(v))
    acf ./= acf[:, first(axes(acf, 2))]
    nestedview(acf)
end


text:
using  LinearAlgebra
const  IdentityMatrix = UniformScaling{Bool}

function clip(s::Real, min_thresh::Real, max_thresh::Real, min_new::Real = min_thresh, max_new::Real = max_thresh)
	s = ifelse(s < min_thresh, min_new, ifelse(s > max_thresh, max_new, s))
end

function scaled_norm(E::IdentityMatrix, v::Array, p::T = T(2)) where {T <: AbstractFloat}
	E.λ ? norm(v, p) : zero(eltype(v))
end

function scaled_norm(E::Diagonal{T}, v::Array{T}, p::T = T(2)) where {T <: AbstractFloat}
	if p == 2
		return scaled_norm2(E, v)
	elseif p == Inf
		return scaled_norm_Inf(E, v)
	elseif p == 1
		return scaled_norm1(E, v)
	else
		throw(ArgumentError("bad norm specified"))
	end
end

function scaled_norm2(E::Diagonal, v::Array)
	sum_sq  = zero(eltype(v))
	for i = 1:length(v)
		sum_sq += (E.diag[i] * v[i])^2
	end
	return sqrt(sum_sq)
end

function scaled_norm_Inf(E::Diagonal, v::Array)
	norm  = zero(eltype(v))
	for i = 1:length(v)
		norm = max(norm, abs(E.diag[i] * v[i]))
	end
	return norm
end

function scaled_norm1(E::Diagonal, v::Array)
	norm  = zero(eltype(v))
	for i = 1:length(v)
		norm += abs(E.diag[i] * v[i])
	end
	return norm
end

function col_norms!(v::Array{Tf, 1},
	A::Matrix{Tf};
	reset::Bool = true) where {Tf <: AbstractFloat}

	if reset
		fill!(v,0.)
	end

	for i = 1:size(A, 2)
		v[i] = max(v[i], norm(view(A, :, i), Inf))
	end
	return v
end

function col_norms!(v::Array{Tf, 1},
	A::SparseMatrixCSC{Tf,Ti}; reset::Bool = true) where{Tf <: AbstractFloat, Ti <: Integer}

	if reset
		fill!(v,0.)
	end

	@inbounds for i = eachindex(v)
		for j = A.colptr[i]:(A.colptr[i + 1] - 1)
			tmp = abs(A.nzval[j])
			v[i] = v[i] > tmp ? v[i] : tmp;
		end
	end
	return v
end

function row_norms!(v::Array{Tf, 1},
	A::Matrix{Tf};
	reset::Bool = true) where{Tf <: AbstractFloat}

	if reset
		fill!(v,0.)
	end

	for i = 1:size(A, 1)
		v[i] = max(v[i], norm(view(A, i, :), Inf))
	end
	return v
end

function row_norms!(v::Array{Tf, 1},
	A::SparseMatrixCSC{Tf, Ti};
	reset::Bool = true) where{Tf <: AbstractFloat, Ti <: Integer}

	if reset
		fill!(v,0.)
	end

	@inbounds for i = 1:(A.colptr[end] - 1)
		idx = A.rowval[i]
		tmp = abs(A.nzval[i])
		v[idx] = v[idx] > tmp ? v[idx] : tmp
	end
	return v
end

function scalarmul!(A::SparseMatrixCSC, c::Real)
	A.nzval .*= c
end

function scalarmul!(A::AbstractMatrix, c::Real)
	A .*= c
end

function lmul!(L::Diagonal{T}, M::SparseMatrixCSC{T}) where {T <: AbstractFloat}

	#NB : Same as:  @views M.nzval .*= D.diag[M.rowval]
	#but this way allocates no memory at all and
	#is marginally faster
	m, n = size(M)
	(m == length(L.diag)) || throw(DimensionMismatch())

	@inbounds for i = 1:(M.colptr[end] - 1)
	 		M.nzval[i] *= L.diag[M.rowval[i]]
	end
	return M
end

lmul!(L::IdentityMatrix, M::AbstractMatrix) = L.λ ? M : M .= zero(eltype(M))

function lmul!(L::Diagonal{T}, x::AbstractVector{T}) where {T <: AbstractFloat}
	(length(L.diag) == length(x)) || throw(DimensionMismatch())
	@. x = x * L.diag
	return nothing
end

lmul!(L::IdentityMatrix, x::AbstractVector{T}) where {T <: AbstractFloat} = L.λ ? x : x .= zero(eltype(x))

function rmul!(M::SparseMatrixCSC{T}, R::Diagonal{T}) where {T <: AbstractFloat}

	m, n = size(M)
	(n == length(R.diag)) || throw(DimensionMismatch())

	@inbounds for i = 1:n, j = M.colptr[i]:(M.colptr[i + 1] - 1)
		 	M.nzval[j] *= R.diag[i]
	end
	return M
end

rmul!(M::AbstractMatrix, R::IdentityMatrix) = R.λ ? R : R .= zero(eltype(R))

function lrmul!(L::Diagonal{T}, M::SparseMatrixCSC{T}, R::Diagonal{T}) where {T <: AbstractFloat}

	m, n = size(M)
	Mnzval  = M.nzval
	Mrowval = M.rowval
	Mcolptr = M.colptr
	Rd      = R.diag
	Ld      = L.diag
	(m == length(Ld) && n == length(Rd)) || throw(DimensionMismatch())

	@inbounds for i = 1:n
		for j = Mcolptr[i]:(Mcolptr[i + 1] - 1)
	 		Mnzval[j] *= Ld[Mrowval[j]] * Rd[i]
		end
	end
	return M
end

lrmul!(L::IdentityMatrix,
	M::AbstractMatrix,
	R::IdentityMatrix) = (L.λ && R.λ) ? M : M .= zero(eltype(M))

lrmul!(L::Diagonal,
	M::SparseMatrixCSC,
	R::IdentityMatrix) = R.λ ? lmul!(L, M) : M .= zero(eltype(M))

lrmul!(L::Diagonal,
	M::AbstractMatrix,
	R::Diagonal) = LinearAlgebra.lmul!(L, LinearAlgebra.rmul!(M, R))

lrmul!(L::Diagonal,
	M::AbstractMatrix,
	R::IdentityMatrix) = R.λ ? LinearAlgebra.lmul!(L, M) : M .= zero(eltype(M))

lrmul!(L::IdentityMatrix,
	M::AbstractMatrix,
	R::Diagonal) = L.λ ? LinearAlgebra.rmul!(M, R) : M .= zero(eltype(M))

"""
    symmetrize_upper!(A)

Symmetrizes the matrix A by calculating A = 0.5 * (A + A') but only performs the operation on the upper triangular part.
"""
function symmetrize_upper!(A::AbstractMatrix{T}) where {T <: AbstractFloat}
	n = size(A, 1)
	@assert(size(A, 1) == size(A, 2), "Matrix is not square.")
	@inbounds for j in 1:n, i in 1:j
 		A[i, j] = (A[i, j] + A[j, i]) / T(2)
	end
	nothing
end

"""
    symmetrize_full!(A)

Symmetrizes the matrix A by calculating A = 0.5 * (A + A') and storing the result in-place.
"""
function symmetrize_full!(A::AbstractMatrix{T}) where {T <: AbstractFloat}
	n = size(A, 1)
	@assert(size(A, 1) == size(A, 2), "Matrix is not square.")
	@inbounds for j in 1:n, i in 1:j
 		A[i, j] = (A[i, j] + A[j, i]) / T(2)
 		A[j, i] = A[i, j]
	end
	nothing
end

# this function assumes real symmetric X and only considers the upper triangular part
function is_pos_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
	# See https://math.stackexchange.com/a/13311
	@inbounds for i = 1:size(X, 1)
		X[i, i] += tol
	end
	F = cholesky!(Symmetric(X), check = false)
	return issuccess(F)
end

function is_neg_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
	@. X *= -one(T)
	return is_pos_def!(X, tol)
end

is_pos_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(copy(X), tol)
is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)

"Round x to the closest multiple of N."
function round_multiple(x::T, N::T) where {T <: Integer}
	return floor(T, x +  0.5 * N - rem(x + 0.5 * N, N))
end


text:
# https://html.spec.whatwg.org/multipage/web-sockets.html

module ExampleWebSocket

using Bukdu
const ServerHost = "localhost"
const ServerPort = 8080

struct WSController <: ApplicationController
    conn::Conn
end

function index(c::WSController)
    render(HTML, """
<!DOCTYPE html>

  <meta charset="utf-8" />
  <title>WebSocket Test</title>

WebSocket Test

julia> ws = first(Bukdu.websockets())
HTTP.WebSockets.WebSocket{Sockets.TCPSocket}(Sockets.TCPSocket(RawFD(0x0000001a) active, 0 bytes waiting), 0x01, true, UInt8[0x57, 0x65, 0x62, 0x53, 0x6f, 0x63, 0x6b, 0x65, 0x74, 0x20, 0x72, 0x6f, 0x63, 0x6b, 0x73], UInt8[], false, false)

julia> write(ws, "hello")
5

""")
end

end # module ExampleWebSocket

if PROGRAM_FILE == basename(@__FILE__)

using Bukdu
import .ExampleWebSocket: WSController, ServerHost, ServerPort, index

plug(Plug.Logger, access_log=(path=normpath(@__DIR__, "access.log"),), formatter=Plug.LoggerFormatter.datetime_message)

routes() do
    get("/", WSController, index)
    plug(Plug.WebSocket) #
end

Bukdu.start(ServerPort, host=ServerHost)

# Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu