text:
# SP Constructs #
# ========================== #
"""
    WS(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenarioaData; solver = JuMP.UnsetSolver())

Generate a **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

In other words, generate the first stage and the second stage of the `stochasticprogram` as if `scenario` is known to occur. Optionally, a capable `solver` can be supplied to `WS`. Otherwise, any previously set solver will be used.

See also: [`DEP`](@ref), [`EVP`](@ref)
"""
function WS(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Return WS model
    return _WS(generator(stochasticprogram,:stage_1), generator(stochasticprogram,:stage_2), stage_parameters(stochasticprogram, 1), stage_parameters(stochasticprogram, 2), scenario, internal_solver(supplied_solver))
end
function _WS(stage_one_generator::Function,
             stage_two_generator::Function,
             stage_one_params::Any,
             stage_two_params::Any,
             scenario::AbstractScenario,
             solver::MathProgBase.AbstractMathProgSolver)
    ws_model = Model(solver = solver)
    stage_one_generator(ws_model, stage_one_params)
    ws_obj = copy(ws_model.obj)
    stage_two_generator(ws_model, stage_two_params, scenario, ws_model)
    append!(ws_obj, ws_model.obj)
    ws_model.obj = ws_obj
    return ws_model
end
"""
    WS_decision(stochasticprogram::TwoStageStochasticProgram, scenario::AbstractScenario; solver = JuMP.UnsetSolver())

Calculate the optimizer of the **wait-and-see** (`WS`) model of the two-stage `stochasticprogram`, corresponding to `scenario`.

Optionally, supply a capable `solver` to solve the wait-and-see problem. The default behaviour is to rely on any previously set solver.

See also: [`WS`](@ref)
"""
function WS_decision(stochasticprogram::StochasticProgram{2}, scenario::AbstractScenario; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve WS model for supplied scenario
    ws_model = WS(stochasticprogram, scenario, solver = solver)
    solve(ws_model)
    # Return WS decision
    decision = ws_model.colVal[1:decision_length(stochasticprogram)]
    if any(isnan.(decision))
        warn("Optimal decision not defined. Check that the EVP model was properly solved.")
    end
    return decision
end
"""
    EWS(stochasticprogram::StochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected wait-and-see result** (`EWS`) of the `stochasticprogram`.

In other words, calculate the expectated result of all possible wait-and-see models, using the provided scenarios in `stochasticprogram`. Optionally, a capable `solver` can be supplied to solve the intermediate problems. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function EWS(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot determine EWS without a solver.")
    end
    # Solve all possible WS models and compute EWS
    return _EWS(stochasticprogram, internal_solver(supplied_solver))
end
"""
    EWS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver = JuMP.UnsetSolver(), confidence = 0.95, N::Integer = 1000)

Approximately calculate the **expected wait-and-see result** (`EWS`) of the two-stage `stochasticmodel` to the given `confidence` level, over the scenario distribution induced by `sampler`.

Supply a capable `solver` to solve the intermediate problems. `N` is the number of scenarios to sample.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function EWS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, N::Integer)
    sp = sample(stochasticmodel, sampler, N)
    𝔼WS, σ = _stat_EWS(sp, solver)
    z = quantile(Normal(0,1), confidence)
    L = 𝔼WS - z*σ/sqrt(N)
    U = 𝔼WS + z*σ/sqrt(N)
    return ConfidenceInterval(L, U, confidence)
end
function _EWS(stochasticprogram::TwoStageStochasticProgram{S,SP}, solver::MathProgBase.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    return sum([begin
                ws = _WS(stochasticprogram.generator[:stage_1],
                         stochasticprogram.generator[:stage_2],
                         stage_parameters(stochasticprogram, 1),
                         stage_parameters(stochasticprogram, 2),
                         scenario,
                         solver)
                solve(ws)
                probability(scenario)*getobjectivevalue(ws)
                end for scenario in scenarios(stochasticprogram.scenarioproblems)])
end
function _EWS(stochasticprogram::TwoStageStochasticProgram{S,SP}, solver::MathProgBase.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    partial_ews = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_ews[i] = remotecall_fetch((sp,stage_one_generator,stage_two_generator,stage_one_params,stage_two_params,solver)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return 0.0
                return sum([begin
                            ws = _WS(stage_one_generator,
                                     stage_two_generator,
                                     stage_one_params,
                                     stage_two_params,
                                     scenario,
                                     solver)
                            solve(ws)
                            probability(scenario)*getobjectivevalue(ws)
                            end for scenario in scenarioproblems.scenarios])
                end,
                w,
                stochasticprogram.scenarioproblems[w-1],
                stochasticprogram.generator[:stage_1],
                stochasticprogram.generator[:stage_2],
                stage_parameters(stochasticprogram, 1),
                stage_parameters(stochasticprogram, 2),
                solver)
        end
    end
    return sum(partial_ews)
end
function _stat_EWS(stochasticprogram::TwoStageStochasticProgram{S,SP},
                        solver::MPB.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    ws_generator = scenario -> WS(stochasticprogram, scenario; solver = solver)
    𝔼WS, σ² = welford(ws_generator, scenarios(stochasticprogram))
    return 𝔼WS, sqrt(σ²)
end
function _stat_EWS(stochasticprogram::TwoStageStochasticProgram{S,SP},
                        solver::MPB.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    partial_welfords = Vector{Tuple{Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch((sp,stage_one_generator,stage_two_generator,stage_one_params,stage_two_params,solver)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return zero(eltype(x)), zero(eltype(x))
                ws_generator = scenario -> _WS(stage_one_generator,
                                               stage_two_generator,
                                               stage_one_params,
                                               stage_two_params,
                                               scenario;
                                               solver = solver)
                return (welford(ws_generator, scenarioproblems.scenarios)..., length(scenarioproblems.scenarios))
            end,
            w,
            stochasticprogram.scenarioproblems[w-1],
            stochasticprogram.generator[:stage_1_vars],
            stochasticprogram.generator[:stage_2],
            stage_parameters(stochasticprogram, 1),
            stage_parameters(stochasticprogram, 2),
            solver)
        end
    end
    𝔼WS, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼WS, sqrt(σ²)
end
"""
    DEP(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Generate the **deterministically equivalent problem** (`DEP`) of the two-stage `stochasticprogram`.

In other words, generate the extended form the `stochasticprogram` as a single JuMP model. Optionally, a capable `solver` can be supplied to `DEP`. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`WS`](@ref)
"""
function DEP(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache,:dep)
        dep = cache[:dep]
        setsolver(dep, supplied_solver)
        return dep
    end
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Define first-stage problem
    dep_model = Model(solver = internal_solver(supplied_solver))
    generator(stochasticprogram,:stage_1)(dep_model, stage_parameters(stochasticprogram, 1))
    dep_obj = copy(dep_model.obj)
    # Define second-stage problems, renaming variables according to scenario.
    stage_two_params = stage_parameters(stochasticprogram, 2)
    visited_objs = collect(keys(dep_model.objDict))
    for (i, scenario) in enumerate(scenarios(stochasticprogram))
        generator(stochasticprogram,:stage_2)(dep_model, stage_two_params, scenario, dep_model)
        append!(dep_obj,probability(scenario)*dep_model.obj)
        for (objkey,obj) ∈ filter(kv->kv.first ∉ visited_objs, dep_model.objDict)
            newkey = if isa(obj,JuMP.Variable)
                varname = add_subscript(dep_model.colNames[obj.col], i)
                dep_model.colNames[obj.col] = varname
                dep_model.colNamesIJulia[obj.col] = varname
                newkey = Symbol(varname)
            elseif isa(obj,Array{JuMP.Variable})
                JuMP.fill_var_names(JuMP.REPLMode, dep_model.colNames, obj)
                arrayname = add_subscript(objkey, i)
                for var in obj
                    splitname = split(dep_model.colNames[var.col],"[")
                    varname = @sprintf("%s[%s", add_subscript(splitname[1],i), splitname[2])
                    dep_model.colNames[var.col] = varname
                    dep_model.colNamesIJulia[var.col] = varname
                end
                newkey = Symbol(arrayname)
            elseif isa(obj,JuMP.ConstraintRef)
                arrayname = add_subscript(objkey, i)
                newkey = Symbol(arrayname)
            elseif isa(obj,Array{ConstraintRef})
                arrayname = add_subscript(objkey, i)
                newkey = Symbol(arrayname)
            elseif isa(obj,JuMP.JuMPArray)
                newkey = if isa(obj,JuMP.JuMPArray{JuMP.ConstraintRef})
                    arrayname = add_subscript(objkey, i)
                    newkey = Symbol(arrayname)
                else
                    JuMP.fill_var_names(JuMP.REPLMode, dep_model.colNames, obj)
                    arrayname = add_subscript(dep_model.varData[obj].name, i)
                    newkey = Symbol(arrayname)
                    dep_model.varData[obj].name = newkey
                    for var in obj.innerArray
                        splitname = split(dep_model.colNames[var.col],"[")
                        varname = @sprintf("%s[%s", add_subscript(splitname[1],i), splitname[2])
                        dep_model.colNames[var.col] = varname
                        dep_model.colNamesIJulia[var.col] = varname
                    end
                    newkey
                end
            else
                continue
            end
            dep_model.objDict[newkey] = obj
            delete!(dep_model.objDict,objkey)
            push!(visited_objs,newkey)
        end
    end
    dep_model.obj = dep_obj
    # Cache DEP
    cache[:dep] = dep_model
    # Return DEP
    return dep_model
end
"""
    VRP(stochasticprogram::StochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **value of the recouse problem** (`VRP`) in `stochasticprogram`.

In other words, optimize the stochastic program and return the optimal value. Optionally, supply a capable `solver` to optimize the stochastic program. Otherwise, any previously set solver will be used.

See also: [`EVPI`](@ref), [`EWS`](@ref)
"""
function VRP(stochasticprogram::StochasticProgram; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve DEP
    optimize!(stochasticprogram, solver = solver)
    # Return optimal value
    return optimal_value(stochasticprogram)
end
"""
    VRP(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver = JuMP.UnsetSolver(), confidence = 0.95)

Return a confidence interval around the **value of the recouse problem** (`VRP`) of `stochasticmodel` to the given `confidence` level.

Optionally, supply a capable `solver` to optimize the stochastic program. Otherwise, any previously set solver will be used.

See also: [`EVPI`](@ref), [`VSS`](@ref), [`EWS`](@ref)
"""
function VRP(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95)
    ss = optimize!(stochasticmodel, sampler; solver = solver, confidence = confidence)
    return confidence_interval(ss)
end
"""
    EVPI(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected value of perfect information** (`EVPI`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `VRP` and `EWS`. Optionally, supply a capable `solver` to solve the intermediate problems. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`EWS`](@ref), [`VSS`](@ref)
"""
function EVPI(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot determine EVPI without a solver.")
    end
    # Calculate VRP
    vrp = VRP(stochasticprogram, solver=supplied_solver)
    # Solve all possible WS models and calculate EWS
    ews = _EWS(stochasticprogram, internal_solver(supplied_solver))
    # Return EVPI = EWS-VRP
    return abs(ews-vrp)
end
"""
    EVPI(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver = JuMP.UnsetSolver(), confidence = 0.95)

Approximately calculate the **expected value of perfect information** (`EVPI`) of the two-stage `stochasticmodel` to the given `confidence` level, over the scenario distribution induced by `sampler`.

In other words, calculate confidence intervals around `VRP` and `EWS`. If they do not overlap, the EVPI is statistically significant, and a confidence interval is calculated and returned. Optionally, supply a capable `solver` to solve the intermediate problems. Otherwise, any previously set solver will be used.

See also: [`VRP`](@ref), [`EWS`](@ref), [`VSS`](@ref)
"""
function EVPI(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, tol::AbstractFloat = 1e-1, kwargs...)
    # Condidence level
    α = (1-confidence)/2
    # Calculate confidence interval around VRP
    ss = optimize!(stochasticmodel, sampler; solver = solver, confidence = 1-α, tol = tol, kwargs...)
    vrp = confidence_interval(ss)
    # EWS solution of the corresponding size
    ews = EWS(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = ss.N)
    try
        evpi = ConfidenceInterval(lower(ews) - upper(vrp), upper(ews) - lower(vrp), confidence)
        lower(evpi) >= -tol || error()
        return evpi
    catch
        @warn "EVPI is not statistically significant to the chosen confidence level and tolerance"
        return ConfidenceInterval(-Inf, Inf, 1.0)
    end
end
"""
    EVP(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Generate the **expected value problem** (`EVP`) of the two-stage `stochasticprogram`.

In other words, generate a wait-and-see model corresponding to the expected scenario over all available scenarios in `stochasticprogram`. Optionally, supply a capable `solver` to `EVP`. Otherwise, any previously set solver will be used.

See also: [`EVP_decision`](@ref), [`EEV`](@ref), [`EV`](@ref), [`WS`](@ref)
"""
function EVP(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache,:evp)
        evp = cache[:evp]
        setsolver(evp, supplied_solver)
        return evp
    end
    # Create EVP as a wait-and-see model of the expected scenario
    ev_model = WS(stochasticprogram, expected(stochasticprogram), solver = supplied_solver)
    # Cache EVP
    cache[:evp] = ev_model
    # Return EVP
    return ev_model
end
"""
    EVP_decision(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the optimizer of the `EVP` of the two-stage `stochasticprogram`.

Optionally, supply a capable `solver` to solve the expected value problem. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EV`](@ref), [`EEV`](@ref)
"""
function EVP_decision(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot determine EWS without a solver.")
    end
    # Solve EVP
    evp = EVP(stochasticprogram, solver = solver)
    solve(evp)
    # Return EVP decision
    decision = evp.colVal[1:decision_length(stochasticprogram)]
    if any(isnan.(decision))
        warn("Optimal decision not defined. Check that the EVP model was properly solved.")
    end
    return decision
end
"""
    EV(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the optimal value of the `EVP` of the two-stage `stochasticprogram`.

Optionally, supply a capable `solver` to solve the expected value problem. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EVP_decision`](@ref), [`EEV`](@ref)
"""
function EV(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP model
    evp = EVP(stochasticprogram; solver = solver)
    solve(evp)
    # Return optimal value
    return getobjectivevalue(evp)
end
"""
    EEV(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **expected value of the expected value solution** (`EEV`) of the two-stage `stochasticprogram`.

In other words, evaluate the `EVP` decision. Optionally, supply a capable `solver` to solve the intermediate problems. The default behaviour is to rely on any previously set solver.

See also: [`EVP`](@ref), [`EV`](@ref)
"""
function EEV(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP model
    evp_decision = EVP_decision(stochasticprogram; solver = solver)
    # Calculate EEV by evaluating the EVP decision
    eev = evaluate_decision(stochasticprogram, evp_decision; solver = solver)
    # Return EEV
    return eev
end
"""
    EEV(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver = JuMP.UnsetSolver(), confidence = 0.95, N::Integer = 100, Ñ::Integer = 1000)

Approximately calculate the **expected value of the expected value decision** (`EEV`) of the two-stage `stochasticmodel` to the given `confidence` level, over the scenario distribution induced by `sampler`.

Supply a capable `solver` to solve the intermediate problems. `N` is the number of scenarios to sample in order to determine the EVP decision and `Ñ` is the number of samples in the out-of-sample evaluation of the EVP decision.

See also: [`EVP`](@ref), [`EV`](@ref)
"""
function EEV(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, N::Integer = 100, Ñ::Integer = 1000)
    sp = sample(stochasticmodel, sampler, N)
    x̄ = EVP_decision(sp; solver = internal_solver(solver))
    return evaluate_decision(stochasticmodel, x̄, sampler; solver = internal_solver(solver), confidence = confidence, Ñ = Ñ)
end
"""
    VSS(stochasticprogram::TwoStageStochasticProgram; solver = JuMP.UnsetSolver())

Calculate the **value of the stochastic solution** (`VSS`) of the two-stage `stochasticprogram`.

In other words, calculate the gap between `EEV` and `VRP`. Optionally, supply a capable `solver` to solve the intermediate problems. The default behaviour is to rely on any previously set solver.
"""
function VSS(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver())
    # Solve EVP and determine EEV
    eev = EEV(stochasticprogram; solver = solver)
    # Calculate VRP
    vrp = VRP(stochasticprogram; solver = solver)
    # Return VSS = VRP-EEV
    return abs(vrp-eev)
end
"""
    VSS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver = JuMP.UnsetSolver(), confidence = 0.95, Ñ::Integer = 1000)

Approximately calculate the **value of the stochastic solution** (`VSS`) of the two-stage `stochasticmodel` to the given `confidence` level, over the scenario distribution induced by `sampler`.

In other words, calculate confidence intervals around `EEV` and `VRP`. If they do not overlap, the VSS is statistically significant, and a confidence interval is calculated and returned. Optionally, supply a capable `solver` to solve the intermediate problems. Otherwise, any previously set solver will be used. `Ñ` is the number of samples in the out-of-sample evaluation of EEV.

See also: [`VRP`](@ref), [`EEV`](@ref), [`EVPI`](@ref)
"""
function VSS(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, Ñ::Integer = 1000, tol::AbstractFloat = 1e-1, kwargs...)
    # Condidence level
    α = (1-confidence)/2
    # Calculate confidence interval around VRP
    ss = optimize!(stochasticmodel, sampler; solver = solver, confidence = 1-α, Ñ = Ñ, tol = tol, kwargs...)
    vrp = confidence_interval(ss)
    # Calculate confidence interval around EEV
    eev = EEV(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = ss.N, Ñ = Ñ)
    try
        vss = ConfidenceInterval(lower(vrp) - upper(eev), upper(vrp) - lower(eev), confidence)
        lower(vss) >= -tol || error()
        return vss
    catch
        @warn "VSS is not statistically significant to the chosen confidence level and tolerance"
        return ConfidenceInterval(-Inf, Inf, 1.0)
    end
end
# ========================== #


text:
using Knet: cuda10

function cuda10src(f, j=f, ex="$f(xi,s)"; BLK=256, THR=256)
"""
__global__ void _$(f)_32_10(int n, float *x, float s, float *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    float xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_32_10(int n, float *x, float s, float *y) {
    _$(f)_32_10<<<$BLK,$THR>>>(n,x,s,y);
  }
}
__global__ void _$(f)_64_10(int n, double *x, double s, double *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    double xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_64_10(int n, double *x, double s, double *y) {
    _$(f)_64_10<<<$BLK,$THR>>>(n,x,s,y);
  }
}
"""
end

for a in cuda10
    isa(a,Tuple) || (a=(a,))
    print(cuda10src(a...))
end

ext: jl


text:
using DataDeps

register(DataDep(
  "hayes-roth",
  "http://archive.ics.uci.edu/ml/datasets/Hayes-Roth",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/hayes-roth/hayes-roth.data",
  "7d7ca47e6d9076863fef2810c80a44a28fe7beffbc17473d4a4948100cffb2dc",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=2:5,
        id_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
import Base: show, getindex, length

type OrderBook <: AbstractTimeSeries

    timestamp::Vector{DateTime{ISOCalendar,UTC}}
    values::Matrix{ASCIIString}
    colnames::Vector{ASCIIString}
    ticker::Ticker

    function OrderBook(timestamp::Vector{DateTime{ISOCalendar,UTC}},
                       values::Matrix{ASCIIString},
                       colnames::Vector{ASCIIString},
                       ticker::Ticker)

                       nrow, ncol = size(values, 1), size(values, 2)
                       nrow != size(timestamp, 1) ? error("values must match length of timestamp"):
                       ncol != size(colnames,1) ? error("column names must match width of array"):
                       timestamp != unique(timestamp) ? error("there are duplicate dates"):
                       ~(flipud(timestamp) == sort(timestamp) || timestamp == sort(timestamp)) ? error("dates are mangled"):
                       flipud(timestamp) == sort(timestamp) ?
                       new(flipud(timestamp), flipud(values), colnames):
                       new(timestamp, values, colnames)
    end
end

const orderbookbidvalues   = ["100" "123.12" "bid" "limit" "pending" "" "2.33"]
const orderbookoffervalues = ["100" "123.12" "bid" "limit" "pending" "" "2.33"]
const orderbooksellvalues  = ["100" "123.12" "sell" "market" "pending" "" "2.33"]
const orderbookcovervalues = ["100" "123.12" "sell" "market" "pending" "" "2.33"]
const orderbookcolnames    = ["Qty","Price","Side","Order", "Status", "Status Time", "Fees"]
const orderbookticker      = Ticker("ticker")

OrderBook(t::DateTime{ISOCalendar,UTC}, v::Matrix{ASCIIString}, c::Vector{ASCIIString}, tick::Ticker) = OrderBook([t], v, c, tick)
OrderBook() = OrderBook([date(1980,1,3)], orderbookbidvalues, orderbookcolnames, orderbookticker)

function add!(ob::OrderBook, entry::OrderBook)
    ob.ticker !== entry.ticker ?
    error("ticker values don't match") :
    OrderBook(vcat(ob.timestamp, entry.timestamp), vcat(ob.values, entry.values), orderbookcolnames, ob.ticker)
end

# merge
function merge(ob1::OrderBook, ob2::OrderBook)
    if ob1.ticker !== ob2.ticker
        error("ticker values don't match")
    else
        dt = sort(vcat(ob1.timestamp, ob2.timestamp))
        vals = fill("",length(dt),length(orderbookcolnames))
        for d in 1:length(dt)
            if ob1[dt[d]] != nothing
                vals[d,:] = ob1[dt[d]].values
            elseif ob2[dt[d]] != nothing
                vals[d,:] = ob2[dt[d]].values
            end
        end
    end
    OrderBook(dt, vals, orderbookcolnames)
end

function fillorderbook(s::TimeArray{Bool,1}, timeseries::FinancialTimeSeries{Float64,2})
    if ob1.ticker !== ob2.ticker
        error("ticker values don't match")
    else

        op, hi, lo, cl = timeseries["Open"], timeseries["High"], timeseries["Low"], timeseries["Close"]

        tt = lag(s)              # 495 row of bools, the next day
        t  = discretesignal(tt)  # 78 rows of first true and first false, as floats though

        # entries day after signal, with a bid of signal day's mid price
        entrydates = findwhen(t.==1)
        entries    = OrderBook(entrydates, repmat(orderbookbidvalues, length(entrydates)), orderbookcolnames, timeseries.instrument.ticker)
        bidsignal  = findwhen(discretesignal(s).==1)
        entryprice  = op[entrydates].values .+ .1  # slippage should NOT be here but in the fill algo below
        for i in 1:length(entries)
            entries.values[i,2] = string(round(entryprice[i],2))
        end

        # exits day after signal, with .1 slippage
        exitdates  = findwhen(t.==0)
        exits      = OrderBook(exitdates, repmat(orderbooksellvalues, length(exitdates)), orderbookcolnames, timeseries.instrument..ticker)
        exitprice  = op[exitdates].values .+ .1  # slippage should NOT be here but in the fill algo below
        for i in 1:length(exits)
            exits.values[i,2] = string(round(exitprice[i],2))
        end
    end

    res = merge(entries,exits)
    res.values[1,5] = "open"
    res
end

###### iterator protocol #########

# start(ob::OrderBook)   = 1
# next(ob::OrderBook,i)  = ((ob.timestamp[i],ob.values[i,:]),i+1)
# done(ob::OrderBook,i)  = (i > length(ob))
# isempty(ob::OrderBook) = (length(ob) == 0)

###### show #####################

function show(io::IO, ta::OrderBook)
  # variables
  nrow          = size(ta.values, 1)
  ncol          = size(ta.values, 2)
  spacetime     = strwidth(string(ta.timestamp[1])) + 3
  firstcolwidth = strwidth(ta.colnames[1])
  colwidth      = Int[]
      for m in 1:ncol
          push!(colwidth, max(strwidth(ta.colnames[m]), maximum([strwidth(t) for t in ta.values[:,m]])))
      end

  # summary line
  print(io,@sprintf("%dx%d %s %s to %s", nrow, ncol, typeof(ta), string(ta.timestamp[1]), string(ta.timestamp[nrow])))
  println(io,"")
  println(io,"")

  # row label line
   print_with_color(:magenta, io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

   for p in 2:length(colwidth)
     print_with_color(:magenta, io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
   end
   println(io, "")

  # timestamp and values line
     if nrow > 7
         for i in 1:4
             print(io, ta.timestamp[i], " | ")
         for j in 1:ncol
            ta.values[i,j] == "open" || ta.values[i,j] == "bid" || ta.values[i,j] == "cover" ?
              print_with_color(:green, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "closed" || ta.values[i,j] == "offer" || ta.values[i,j] == "sell" ?
              print_with_color(:red, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "pending" ?
              print_with_color(:yellow, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            print_with_color(:blue, io, rpad(ta.values[i,j], colwidth[j] + 2, " "))
         end
         println(io,"")
         end
         println(io,'\u22EE')
         for i in nrow-3:nrow
             print(io, ta.timestamp[i], " | ")
         for j in 1:ncol
            ta.values[i,j] == "open" || ta.values[i,j] == "bid" || ta.values[i,j] == "cover" ?
              print_with_color(:green, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "closed" || ta.values[i,j] == "offer" || ta.values[i,j] == "sell" ?
              print_with_color(:red, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "pending" ?
              print_with_color(:yellow, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            print_with_color(:blue, io, rpad(ta.values[i,j], colwidth[j] + 2, " "))
         end
         println(io,"")
         end
     else
        for i in 1:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            ta.values[i,j] == "open" || ta.values[i,j] == "bid" || ta.values[i,j] == "cover" ?
              print_with_color(:green, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "closed" || ta.values[i,j] == "offer" || ta.values[i,j] == "sell" ?
              print_with_color(:red, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            ta.values[i,j] == "pending" ?
              print_with_color(:yellow, io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            print_with_color(:blue, io, rpad(ta.values[i,j], colwidth[j] + 2, " "))
        end
        println(io,"")
        end
    end
end

###### getindex #################

# single row
function getindex(b::OrderBook, n::Int)
    OrderBook(b.timestamp[n], b.values[n,:], b.colnames, b.ticker)
end

# range of rows
function getindex(b::OrderBook, r::Range1{Int})
    OrderBook(b.timestamp[r], b.values[r,:], b.colnames, b.ticker)
end

# array of rows
function getindex(b::OrderBook, a::Array{Int})
    OrderBook(b.timestamp[a], b.values[a,:], b.colnames, b.ticker)
end

# single column by name
function getindex(b::OrderBook, s::ASCIIString)
    n = findfirst(b.colnames, s)
    OrderBook(b.timestamp, b.values[:, n], ASCIIString[s], b.ticker)
end

# array of columns by name
function getindex(b::OrderBook, args::ASCIIString...)
    ns = [findfirst(b.colnames, a) for a in args]
    OrderBook(b.timestamp, b.values[:,ns], ASCIIString[a for a in args], b.ticker)
end

# single date
function getindex(b::OrderBook, d::DateTime{ISOCalendar,UTC})
   for i in 1:length(b)
     if [d] == b[i].timestamp
       return b[i]
     else
       nothing
     end
   end
 end

# range of dates
function getindex(b::OrderBook, dates::Array{DateTime{ISOCalendar,UTC}})
  counter = Int[]
#  counter = int(zeros(length(dates)))
  for i in 1:length(dates)
    if findfirst(b.timestamp, dates[i]) != 0
      #counter[i] = findfirst(b.timestamp, dates[i])
      push!(counter, findfirst(b.timestamp, dates[i]))
    end
  end
  b[counter]
end

function getindex(b::OrderBook, r::DateTimeRange{ISOCalendar})
    b[[r]]
end

# day of week
# getindex{T,N}(b::OrderBook{T,N}, d::DAYOFWEEK) = b[dayofweek(b.timestamp) .== d]


text:
"""
    Atomix.@atomic

A superset of `Base.@atomic` supporting atomic operations on array elements.
Atomic operations on fields dispatches to `Base.@atomic`.
"""
:(@atomic)

macro atomic(ex)
    ans = atomic_impl(QuoteNode(:sequentially_consistent), ex)
    ans === nothing || return ans
    esc(:($Base.@atomic($ex)))
end

macro atomic(order, ex)
    ans = atomic_impl(order, ex)
    ans === nothing || return ans
    esc(:($Base.@atomic($order, $ex)))
end

macro atomic(a1, op, a2)
    ans = atomic_impl(QuoteNode(:sequentially_consistent), a1, op, a2)
    ans === nothing || return ans
    esc(:($Base.@atomic($a1, $op, $a2)))
end

macro atomic(order, a1, op, a2)
    ans = atomic_impl(order, a1, op, a2)
    ans === nothing || return ans
    esc(:($Base.@atomic($order, $a1, $op, $a2)))
end

function ref_expr(ex)
    @nospecialize
    isexpr(ex, :ref) || return nothing
    m = esc(ex.args[1])
    indices = map(esc, ex.args[2:end])
    return :(referenceable($m)[$(indices...)])
end

function order_expr(order)
    @nospecialize
    if order isa QuoteNode
        llvm_ordering_from_juila(order.value)
    else
        :(llvm_ordering_from_juila($(esc(order))))
    end
end

function atomic_impl(order, ex)
    @nospecialize
    if ex isa Expr
        if (ref = ref_expr(ex)) !== nothing
            return :(Atomix.get($ref, $(order_expr(order))))
        elseif isexpr(ex, :call, 3)
            return atomic_impl(order, ex.args[2], ex.args[1], ex.args[3])
        elseif ex.head === :(=)
            l, r = ex.args[1], esc(ex.args[2])
            if (ref = ref_expr(l)) !== nothing
                return :(Atomix.set!($ref, $r, $(order_expr(order))))
            end
        elseif length(ex.args) == 2
            shead = string(ex.head)
            if endswith(shead, '=')
                op = Symbol(shead[1:prevind(shead, end)])
                ans = atomic_impl(order, ex.args[1], op, ex.args[2])
                ans === nothing || return :($ans[2])
            end
        end
    end
    return nothing
end

function atomic_impl(order, a1, op, a2)
    @nospecialize
    ref = ref_expr(a1)
    ref === nothing && return nothing
    :(Atomix.modify!($ref, $(esc(op)), $(esc(a2)), $(order_expr(order))))
end

"""
    Atomix.@atomicswap

A superset of `Base.@atomicswap` supporting atomic operations on array elements.
Atomic operations on fields dispatches to `Base.@atomicswap`.
"""
:(@atomicswap)

macro atomicswap(ex)
    ans = atomicswap_impl(QuoteNode(:sequentially_consistent), ex)
    ans === nothing || return ans
    esc(:($Base.@atomicswap($ex)))
end

macro atomicswap(order, ex)
    ans = atomicswap_impl(order, ex)
    ans === nothing || return ans
    esc(:($Base.@atomicswap($order, $ex)))
end

function atomicswap_impl(order, ex)
    @nospecialize
    isexpr(ex, :(=), 2) || return nothing
    ref = ref_expr(ex.args[1])
    ref === nothing && return nothing
    val = esc(ex.args[2])
    return :(first(Atomix.swap!($ref, $val, $(order_expr(order)))))
end

"""
    Atomix.@atomicreplace

A superset of `Base.@atomicreplace` supporting atomic operations on array
elements.  Atomic operations on fields dispatches to `Base.@atomicreplace`.
"""
:(@atomicreplace)

macro atomicreplace(ex, old_new)
    order = QuoteNode(:sequentially_consistent)
    ans = atomicreplace_impl(order, order, ex, old_new)
    ans === nothing || return ans
    esc(:($Base.@atomicreplace($ex, $old_new)))
end

macro atomicreplace(order, ex, old_new)
    ans = atomicreplace_impl(order, order, ex, old_new)
    ans === nothing || return ans
    esc(:($Base.@atomicreplace($order, $ex, $old_new)))
end

macro atomicreplace(success_order, fail_order, ex, old_new)
    ans = atomicreplace_impl(success_order, fail_order, ex, old_new)
    ans === nothing || return ans
    esc(:($Base.@atomicreplace($success_order, $fail_order, $ex, $old_new)))
end

function atomicreplace_impl(success_order, fail_order, ex, old_new)
    @nospecialize
    ref = ref_expr(ex)
    ref === nothing && return nothing
    so = order_expr(success_order)
    fo = order_expr(fail_order)
    if isexpr(old_new, :call, 3) && old_new.args[1] === :(=>)
        exp, rep = esc(old_new.args[2]), esc(old_new.args[3])
        return :(Atomix.replace!($ref, $exp, $rep, $so, $fo))
    else
        old_new = esc(old_new)
        return :(Atomix.replace!($ref, $old_new::Pair..., $so, $fo))
    end
end


text:
function upload_release(user, repo, token, tag, path)
    ghr() do ghr_path
        run(`$ghr_path -replace -u $(user) -r $(repo) -t $(token) $(tag) $(path)`)
    end
end

# Well, to be more precise, last non patch
function last_major_version()
    path = basedir("..", "Project.toml")
    version = VersionNumber(TOML.parse(String(read(path)))["version"])
    return "v" * string(VersionNumber(version.major, version.minor))
end

function upload_reference_images(path=basedir("recorded"), tag=last_major_version(); name="refimages")
    mktempdir() do dir
        tarfile = joinpath(dir, "$(name).tar")
        Tar.create(path, tarfile)
        upload_release("JuliaPlots", "Makie.jl", ENV["GITHUB_TOKEN"], tag, tarfile)
    end
end

function download_refimages(tag=last_major_version(); name="refimages")
    url = "https://github.com/JuliaPlots/Makie.jl/releases/download/$(tag)/$(name).tar"
    images_tar = basedir("$(name).tar")
    images = basedir(name)
    isfile(images_tar) && rm(images_tar)
    isdir(images) && rm(images, recursive=true, force=true)
    Base.download(url, images_tar)
    Tar.extract(images_tar, images)
    return images
end


text:
"""
Read any of the citation network datasets “Cora”, “CiteSeer” and “PubMed”
from the “Revisiting Semi-Supervised Learning with Graph Embeddings” paper.
Nodes represent documents and edges represent citation links.

Data collected from
https://github.com/kimiyoung/planetoid/raw/master/data
"""
function read_planetoid_data(DEPNAME; dir=nothing, reverse_edges=true)
    name = lowercase(DEPNAME)

    x  = read_planetoid_file(DEPNAME, "ind.$(name).x", dir)
    y  = read_planetoid_file(DEPNAME, "ind.$(name).y", dir)
    allx  = read_planetoid_file(DEPNAME, "ind.$(name).allx", dir)
    ally  = read_planetoid_file(DEPNAME, "ind.$(name).ally", dir)
    tx  = read_planetoid_file(DEPNAME, "ind.$(name).tx", dir)
    ty  = read_planetoid_file(DEPNAME, "ind.$(name).ty", dir)
    graph = read_planetoid_file(DEPNAME, "ind.$(name).graph", dir)
    test_indices = read_planetoid_file(DEPNAME, "ind.$(name).test.index", dir)

    ntrain = size(x, 2)
    train_indices = 1:ntrain
    val_indices = ntrain+1:ntrain+500
    sorted_test_index = sort(test_indices)

    if name == "citeseer"
        # There are some isolated nodes in the Citeseer graph, resulting in
        # none consecutive test indices. We need to identify them and add them
        # as zero vectors to `tx` and `ty`.
        len_test_indices = (maximum(test_indices) - minimum(test_indices)) + 1

        tx_ext = zeros(size(tx,1), len_test_indices)
        tx_ext[:, sorted_test_index .- minimum(test_indices) .+ 1] .= tx
        ty_ext = zeros(len_test_indices)
        ty_ext[sorted_test_index .- minimum(test_indices) .+ 1] = ty

        tx, ty = tx_ext, ty_ext
    end
    x = hcat(allx, tx)
    y = vcat(ally, ty)
    x[:, test_indices] = x[:, sorted_test_index]
    y[test_indices] = y[sorted_test_index]
    test_indices = size(allx,2)+1:size(x,2)
    num_nodes = size(x, 2)

    ## Use masks instead of indices?
    # train_mask = falses(num_nodes)
    # train_mask[train_indices] .= 1
    # val_mask = falses(num_nodes)
    # val_mask[val_indices] .= 1
    # test_mask = falses(num_nodes)
    # test_mask[test_indices] .= 1

    adj_list = [Int[] for i=1:num_nodes]
    for (i, neigs) in pairs(graph) # graph is dictionay representing the adjacency list
        neigs = unique(neigs) # remove duplicated edges
        neigs = filter(x -> x!=i, neigs)# remove self-loops
        append!(adj_list[i+1], neigs .+ 1) # convert to 1-indexed
    end
    if reverse_edges
        for (i, neigs) in enumerate(adj_list)
            for j in neigs
                i ∉ adj_list[j] && push!(adj_list[j], i)
            end
        end
    end

    node_data = (features=x, targets=y,
                train_indices, val_indices, test_indices)

    metadata = Dict{String,Any}(
        "num_classes" => length(unique(y)),
        "classes" => sort(unique(y))
    )

    edge_index = adjlist2edgeindex(adj_list)

    g = Graph(; num_nodes,
                num_edges = length(edge_index[1]),
                edge_index,
                node_data)

    return metadata, g
end

function read_pickle_file(filename, name)
    out = Pickle.npyload(filename)
    if name == "graph"
        return out
    end
    if out isa SparseMatrixCSC
        return Matrix(out)
    end
    return out
end

function read_planetoid_file(DEPNAME, name, dir)
    filename = datafile(DEPNAME, name, dir)
    if endswith(name, "test.index")
        out = 1 .+ vec(readdlm(filename, Int))
    else
        out = read_pickle_file(filename, name)
        if out isa Matrix
            out = collect(out')
        end
    end
    if endswith(name, "y")
        out = map(y->y[1], argmax(out, dims=1)) |> vec
    end
    return out
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Julia 1.2.0
#     language: julia
#     name: julia-1.2
# ---

# ## Load data

# +
using Plots, InvertedIndices, StatsPlots, CSV, LinearAlgebra, StatsBase, Distributions, DataFrames
using Mamba:Mamba, quantile, mean, mm

chain1 = CSV.read("output/Matsumoto_MCMCsamples.csv",DataFrame)
nclassmat=CSV.read("../data/NumberOfClasses.csv",DataFrame)[:,2:7]|>Matrix
nstudentsmat=CSV.read("../data/NumberOfClasses.csv",DataFrame)[:,[8:13...]]|>Matrix;

# -

# Q-Q plot
using RCall
RCall.rcall_p(:options, rcalljl_options=Dict(:width => 1000, :height => 1000))
chainmat=Matrix(chain)[:,1:end-1]
@rput chainmat
R"library(royston);par(cex=2);royston.test(chainmat)"
# -

## Class/grade distributuion

gr(fontfamily="Helvetica")
using Random;Random.seed!(2020)
pl=scatter(ceil.(nclassmat).+(rand(29,6).-0.5)./4,nstudentsmat./ceil.(nclassmat), grid=false,
    color=randperm(29),markersize=4,markerstrokewidth=0,legend=false,tick_direction=:out,
    size=(400,300),xguidefontsize=8,yguidefontsize=8)
xlabel!("number of classes per grade")
ylabel!("class size")
plot!(pl,[10.0],linetype=:hline,color=:gray,linestyle=:dash)

## Class/grade against attack ratio
using InvertedIndices
attackrate=[7.9, 11.5, 16.2, 33.5, 50.0, 21.0, 14.8, 23.6, 19.4, 15.7, 10.1, 37.0, 13.5, 16.0, 16.0, 23.1, 18.5, 21.2, 12.5, 40.2, 0.0, 8.2, 29.1, 27.2, 12.6, 30.6, 19.7, 29.9, 21.6]
meanclasssize=[31.833333333333332, 32.19166666666666, 6.166666666666667, 27.833333333333332, 4.666666666666667, 30.0, 28.875, 31.708333333333332, 30.758333333333336, 25.416666666666668, 33.625, 25.666666666666668, 29.0, 29.041666666666668, 28.166666666666668, 33.44444444444445, 31.26388888888889, 21.583333333333332, 30.805555555555554, 21.166666666666668, 5.833333333333333, 29.333333333333332, 32.23055555555556, 25.916666666666668, 30.333333333333332, 29.972222222222225, 30.222222222222218, 30.708333333333332, 36.333333333333336]
meannclasses=[2.0, 4.833333333333333, 1.0, 1.0, 0.5, 2.3333333333333335, 3.8333333333333335, 4.166666666666667, 4.833333333333333, 2.0, 4.5, 1.0, 2.6666666666666665, 3.8333333333333335, 2.0, 3.0, 3.5, 2.0, 2.3333333333333335, 1.0, 0.8333333333333334, 2.0, 4.666666666666667, 2.1666666666666665, 2.0, 2.5, 2.8333333333333335, 3.8333333333333335, 2.0]
attackrate, meanclasssize, meannclasses = getindex.((attackrate, meanclasssize, meannclasses),Ref(Not([3,5,21])))
@rimport stats as st
@show st.cor_test(meanclasssize,attackrate)
@show st.cor_test(meannclasses,attackrate)
plot(scatter.((meanclasssize,meannclasses),Ref(attackrate))...,ylim=(0,50),xlim=[(15,40) (0.5,5)],labels=["r = -0.39 (-0.68,-0.01)" "r = -0.30 (-0.62, 0.10)"],legend=(0.4,0.95),
    xlabel=["mean class size" "mean # classes per grade"], ylabel=["attack ratio %" "attack ratio %"], markercolor=:black,size=(250,500),layout=(2,1))

# Plotting Rs
## Functions
function posterior(chain, parameter::Symbol)
    if parameter==:β0 parameter=:βs;fn=β0
    elseif parameter==:rcom0 parameter=:rcom;fn=rcom0
    elseif parameter in (:suscoef,:infcoef) fn=(slice->slice[parameter].|>exp)
        else fn = (slice->slice[parameter]) end
    v=fn(chain)
    (names=parameter ,samples=v)#[getindex.(v,x) for x in 1:size(first(v),1)])
end

function β0(slice)
    (slice[:logβs]).|>exp end
function rcom0(slice)
    (slice[:logrcom]).|>exp end
bclog(x)=log.(x)
function R_nm(chain,nclasses,sizerange, meansuscovar=0,meaninfcovar=0)
    βs=posterior(chain,:β0).samples
    βsm=βs[:,1:3];βsm[:,3].=vec(mean(βs[:,3:4],dims=2))
    βsm.*=exp.(sum((posterior(chain,:suscoef).samples).*meansuscovar',dims=2)+sum((posterior(chain,:infcoef).samples).*meaninfcovar',dims=2))
    if nclasses==1 βsm[:,2].=0 end
    γs=posterior(chain,:γs).samples
    δs=posterior(chain,:δs).samples
    [quantile.((βsm./(classsize/30).^γs./(nclasses./[3,3,3]').^δs).*[5*nclasses*classsize (nclasses-1)*classsize classsize] |>eachcol,Ref([0.5,0.025,0.975])) for classsize in sizerange]
end
function R_nm(parm::Vector{<:AbstractFloat},nclasses,sizerange,meansuscovar=0,meaninfcovar=0)
    βs=parm[1:4]
    βsm=βs[1:3];βsm[3]=mean(βs[3:4])
    suscoef=parm[12:15];infcoef=parm[16:20]
    if nclasses==1 βsm[2]=0 end
    γs=parm[5:7]
    δs=parm[8:10]
    [(βsm./(classsize/30).^γs./(nclasses./[3,3,3]).^δs).*[5*nclasses*classsize,(nclasses-1)*classsize,classsize] for classsize in sizerange]
end
function R_reduction(chain,nmbefore,nmafter,coefafter=ones(3),coefbefore=ones(3))
    βs=posterior(chain,:β0).samples
    βsm=βs[:,1:3];βsm[:,3].=vec(mean(βs[:,3:4],dims=2))
    if nmbefore[2]==1 coefbefore[2]=0 end
    if nmafter[2]==1 coefafter[2]=0 end
    γs=posterior(chain,:γs).samples
    δs=posterior(chain,:δs).samples
    popbefore=[5*prod(nmbefore) prod(nmbefore)-nmbefore[1] nmbefore[1]]
    popafter=[5*prod(nmafter) prod(nmafter)-nmafter[1] nmafter[1]]
    quantile(sum(popafter.*coefafter'.*βsm./(nmafter[1]/30).^γs./(nmafter[2]./[3,3,3]').^δs,dims=2)./sum((popbefore.*coefbefore'.*βsm./(nmbefore[1]/30).^γs./(nmbefore[2]./[3,3,3]').^δs),dims=2) |>vec, [0.5,0.025,0.975])
end
function Rtotal_nm(chain,nclasses,sizerange, meansuscovar=0,meaninfcovar=0)
    βs=posterior(chain,:β0).samples
    βsm=βs[:,1:3];βsm[:,3].=vec(mean(βs[:,3:4],dims=2))
    βsm.*=exp.(sum((posterior(chain,:suscoef).samples).*meansuscovar',dims=2)+sum((posterior(chain,:infcoef).samples).*meaninfcovar',dims=2))
    if nclasses==1 βsm[:,2].=0 end
    γs=posterior(chain,:γs).samples
    δs=posterior(chain,:δs).samples
    [quantile.(sum((βsm./(classsize/30).^γs./(nclasses./[3,3,3]').^δs).*[5*nclasses*classsize (nclasses-1)*classsize classsize],dims=2) |>eachcol,Ref([0.5,0.025,0.975])) for classsize in sizerange]
end

# ## Posteriors

parnames=["βs[".*string.(1:4).*"]";"rcom";"γs[".*string.(1:3).*"]";"δs[".*string.(1:3).*"]";"suscoef[".*string.(1:4).*"]";"infcoef[".*string.(1:5).*"]";"loglikelihood"]
@time posteriors=posterior.(Ref(chain1),[:β0,:rcom0,:γs,:δs,:suscoef,:infcoef,:loglikelihood])
posteriors=reduce(vcat,posteriors)
hcat(parnames,reduce(vcat,(x->quantile.(x.samples|>eachcol,Ref([0.5,0.025,0.975]))).(posteriors)))

# Plotting school reproduction number

@time Rnms=(R_nm(chain1,nclasses,1:40) for nclasses in 1:6)
medRnms=[reduce(hcat,(x->first.(x)).(Rnm))' for Rnm in Rnms];
uppRnms=[reduce(hcat,(x->last.(x)).(Rnm))' for Rnm in Rnms];
lowRnms=[reduce(hcat,(x->getindex.(x,2)).(Rnm))' for Rnm in Rnms];
function split5(medRnm)
    [fill(medRnm[1],5)./5;(medRnm)[2:3]]
end
function medR(medRnms,n,m)
    (getindex.(medRnms,n,:)[m])|>transpose
end

gr(fontfamily="Helvetica",foreground_color_legend = nothing,background_color_legend = nothing, titlefontsize=9,grid=false, tick_direction=:out,)
# median
barsbyclass=groupedbar.(Ref(1:5),[reduce(vcat,medR.(Ref(medRnms),n,1:5)) for n in [20,30,40]], bar_position = :stack, bar_width=0.9,
    ylim=(0,1),color=reverse([1 :skyblue 9],dims=2),linealpha=0.8)
xlabel!.(barsbyclass,"number of classes per grade")
ylabel!.(barsbyclass,"median (stacked)")
plot(barsbyclass...,layout=(1,3),size=(800,200),xguidefontsize=8,yguidefontsize=8,grid=false,
    label=["classmate" "grademate" "schoolmate" fill("",1,4)], legend=[:none#=(0.75,1.09)=# fill(:none,1,3)],legendfontsize=6,
    top_margin=4mm,bottom_margin=5mm,left_margin=2mm,title="class size: ".* string.([20 30 40]),tick_direction=:out,format=:svg)

# by proximity
medv=[reduce(vcat,reverse.(medR.(Ref(medRnms),n,1:5),dims=2)) for n in [20,30,40]]
uerr=[reduce(vcat,reverse.(medR.(Ref(uppRnms),n,1:5),dims=2)) for n in [20,30,40]].-medv
lerr=medv.-[reduce(vcat,reverse.(medR.(Ref(lowRnms),n,1:5),dims=2)) for n in [20,30,40]]
yerr=((x,y)->Tuple.(vcat.(x,y))).(lerr,uerr)
barsbyclass=[groupedbar(1:5,medv[n], bar_width=0.9,
    ylim=(0,1),color=[1 :skyblue 9],linealpha=0.8,yerrors=yerr[n]) for n in 1:3]
xlabel!.(barsbyclass,"number of classes per grade")
ylabel!.(barsbyclass,"contribution by proximity")
plot(barsbyclass...,layout=(1,3),size=(800,200),xguidefontsize=8,yguidefontsize=8,grid=false,
    label=["classmate" "grademate" "schoolmate" fill("",1,4)], legend=[:none#=(0.75,1.09)=# fill(:none,1,3)],legendfontsize=6,
    top_margin=4mm,bottom_margin=5mm,left_margin=2mm,title="class size: ".* string.([20 30 40]),tick_direction=:out,format=:svg)

# overall Rs
@time Rtot=(Rtotal_nm(chain1,nclasses,1:40) for nclasses in 1:6)
medRtot=[reduce(hcat,(x->first.(x)).(Rnm))' for Rnm in Rtot];
uppRtot=[reduce(hcat,(x->last.(x)).(Rnm))' for Rnm in Rtot];
lowRtot=[reduce(hcat,(x->getindex.(x,2)).(Rnm))' for Rnm in Rtot];
medvtot=[reduce(vcat,reverse.(medR.(Ref(medRtot),n,1:5),dims=2)) for n in [20,30,40]]
uerrtot=[reduce(vcat,reverse.(medR.(Ref(uppRtot),n,1:5),dims=2)) for n in [20,30,40]].-medvtot
lerrtot=medvtot.-[reduce(vcat,reverse.(medR.(Ref(lowRtot),n,1:5),dims=2)) for n in [20,30,40]]
yerrtot=((x,y)->Tuple.(vcat.(x,y))).(lerr,uerrtot)
barsbyclass=[groupedbar(1:5,medvtot[n], bar_width=0.9,
    ylim=(0,1.5),color=[6],linealpha=0.8,yerrors=yerrtot[n]) for n in 1:3]
xlabel!.(barsbyclass,"number of classes per grade")
ylabel!.(barsbyclass,"school reproduction number")
plot(barsbyclass...,layout=(1,3),size=(800,200),xguidefontsize=8,yguidefontsize=8,grid=false,
    label=["classmate" "grademate" "schoolmate" fill("",1,4)], legend=[:none#=(0.75,1.09)=# fill(:none,1,3)],legendfontsize=6,
    top_margin=4mm,bottom_margin=5mm,left_margin=2mm,title="class size: ".* string.([20 30 40]),tick_direction=:out,format=:svg)
Rreductions=reduce(hcat,R_reduction.(Ref(chain1),Ref((40,2)),repeat([(40,2),(20,4),(20,2),(40,1)],3),repeat([ones(3),[0.5,0.5,1],[0.1,0.1,1]],inner=4)))
colors=[:black,:royalblue,3]

plot([Rreductions[1,:] fill(-1.,12,2)],
line=0,marker=([:circle#=,:pent,:hex,:oct=#],6),markerstrokewidth=0,ylim=(0,2),color=[repeat(colors,inner=4) repeat(colors[2:3]|>permutedims,12)], bottom_margin=10mm,
xticks=(1:12,repeat(["no change","split class","staggered attendance\n(within class)", "staggered attendance\n(between class)"],3)),xrotation=60,
label=string.([100 50 10]).*"% outside-class interaction", legend=(0.67,1))
ylabel!("relative change in Rs",yguidefontsize=8)
plot!([1.0],linetype=:hline,color=:gray,linestyle=:dash,label="")
ohlc!(OHLC[Tuple([NaN;Rreductions[[3,2],c];NaN]) for c in 1:size(Rreductions,2)],color=repeat(colors,inner=4),width=2,label="")

# ## Sensitivity analysis

include("SchoolOutbreak.jl")

# +
# Sensitivity analysis
<<<<<<< HEAD:src/outputs.jl
# estim: vector of parameter estimates for each sensitivity analysis
=======
# estim (input): vector of parameter estimates for each sensitivity analysis
>>>>>>> 6620de9142409cbb8d8e2c43ef41ee73e2ceec84:src/outputs_inference.jl
medRnms=(reduce(hcat,R_nm(estim,nclasses,1:40))' for nclasses in 1:6)
barsbyclass=groupedbar.(Ref(1:5),[reduce(vcat,medR.(Ref(medRnms),n,1:5)) for n in [20,30,40]], bar_position = :stack, bar_width=0.9,
    ylim=(0,1),color=reverse([1 :skyblue 9],dims=2))
xlabel!.(barsbyclass,"number of classes per grade")
ylabel!.(barsbyclass,"school reproduction number")
plot(barsbyclass...,layout=(1,3),size=(800,200),xguidefontsize=8,yguidefontsize=8,grid=false,
    label=["classmate" "grademate" "schoolmate" fill("",1,4)], legend=[:none#=(0.75,1.07)=# fill(:none,1,3)],legendfontsize=6,
    top_margin=4mm,bottom_margin=5mm,title="class size: ".* string.([20 30 40]),tick_direction=:out,format=:svg)
# -


text:
function quicksort!(A,i=1,j=length(A))
    if j > i
        pivot = A[rand(i:j)] # random element of A
        left, right = i, j
        while left <= right
            while A[left] < pivot
                left += 1
            end
            while A[right] > pivot
                right -= 1
            end
            if left <= right
                A[left], A[right] = A[right], A[left]
                left += 1
                right -= 1
            end
        end
        quicksort!(A,i,right)
        quicksort!(A,left,j)
    end
    return A
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
	MultipleRepresentation(extractors::Tuple)

	crate a `ProductNode` where each item is the json part
	processed by all extractors in the order

"""
struct MultipleRepresentation{E<:Tuple}
	extractors::E
end

(m::MultipleRepresentation)(x) = TreeNode(map(e -> e(x), m.extractors))

extractsmatrix(s::MultipleRepresentation) = false

function Base.show(io::IO, m::MultipleRepresentation; pad = [], key::String="")
	c = COLORS[(length(pad)%length(COLORS))+1]
	key *=": "
	paddedprint(io,"$(key)MultiRepresentation\n", color = c, pad = pad)
	if m.extractors != nothing
		printtuple(io, m.extractors, c = c, pad = pad)
	end
end

function printtuple(io, d::Tuple; c, pad, last::Bool = true)
	for i in 1:length(d)
  		p = (i == length(d)) ? "  └───" : "  ├───"
  		show(io, d[i], pad = vcat(replacebyspaces(pad), (c,p)) , key = " ")
 	end
end

Base.hash(e::MultipleRepresentation, h::UInt) = hash(e.extractors, h)
Base.:(==)(e1::MultipleRepresentation, e2::MultipleRepresentation) = e1.extractors == e2.extractors


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function image_testing(fname)
    width, height = 4000, 4000
    margin = 500
    Drawing(width, height, fname)
    origin()
    background("grey25")
    setline(5)
    sethue("green")
    image = readpng(dirname(@__FILE__) * "/stackoverflow.png")
    w = image.width
    h = image.height
    pagetiles = Tiler(width, height, 7, 9)
    tw = pagetiles.tilewidth/2
    for (pos, n) in pagetiles
        circle(pos, tw, :stroke)
        circle(pos, tw, :clip)
        gsave()
        translate(pos)
        scale(25, 25)
        rotate(rand(0.0:pi/8:2pi))
        placeimage(image, O, centered=true)
        grestore()
        clipreset()
    end
    @test finish() == true
end

fname = "test-image.png"
image_testing(fname)
println("...finished test: output in $(fname)")


text:
#using Roots

"""
    Datatypes for IC Parameters and Solution
"""

mutable struct SodCRSolution

    x::Array{Float64,1}
    rho::Array{Float64,1}
    rho3::Float64
    rho4::Float64
    P_tot::Array{Float64,1}
    P_th::Array{Float64,1}
    P_cr::Array{Float64,1}
    P_cr_p::Array{Float64,1}
    P_cr_e::Array{Float64,1}
    P34_tot::Float64
    P3_th::Float64
    P3_cr::Float64
    P4_th::Float64
    P4_cr::Float64
    U_tot::Array{Float64,1}
    U_th::Array{Float64,1}
    E_cr::Array{Float64,1}
    E_cr_p::Array{Float64,1}
    E_cr_e::Array{Float64,1}
    v::Array{Float64,1}
    v34::Float64
    vt::Float64
    vs::Float64
    xs::Float64
    xr::Float64
    Mach::Float64

    function SodCRSolution(x::Array{Float64,1})
        N = length(x)
        new(x,
            zeros(N),   # rho
            0.0,        # rho3
            0.0,        # rho4
            zeros(N),   # P_tot
            zeros(N),   # P_th
            zeros(N),   # P_cr
            zeros(N),   # P_cr_p
            zeros(N),   # P_cr_e
            0.0,        # P34_tot
            0.0,        # P3_th
            0.0,        # P3_cr
            0.0,        # P4_th
            0.0,        # P4_cr
            zeros(N),   # U_tot
            zeros(N),   # U_th
            zeros(N),   # E_cr
            zeros(N),   # E_cr_p
            zeros(N),   # E_cr_e
            zeros(N),   # v
            0.0,        # v34
            0.0,        # vt
            0.0,        # vs
            0.0,        # xs
            0.0,        # xr
            0.0)        # Mach
    end
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
function scoping(ast)
    rec = scoping
    @match ast begin
        :([$(frees...)]($(args...), ) -> begin $(stmts...) end) =>
            let stmts = map(rec, stmts),
                arw   = :(($(args...), ) -> begin $(stmts...) end)
                Expr(:scope, (), Tuple(frees), (), arw)
            end
        Expr(hd, tl...) => Expr(hd, map(rec, tl)...)
        a => a
    end
end


text:
__precompile__()
module Reduce
using ForceImport, SyntaxTree, LinearAlgebra
!isdefined(Reduce,:linefilter!) && (linefilter! = SyntaxTree.linefilter)

#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

include(joinpath(@__DIR__,"../deps/svn.jl"))

struct PSL <: Base.AbstractPipe
    input::Pipe
    output::Pipe
    process::Base.Process
    function PSL()
        # Setup pipes and reduce process
        input = Pipe()
        output = Pipe()
        rsl = red()
        dirf = joinpath(@__DIR__,"..","deps")
        if !Sys.iswindows()
            try
                process = _spawn(red(), input, output)
            catch
                process = _spawn(redsys(dirf), input, output)
            end
        else
            process = _spawn(redsys(dirf), input, output)
        end
        # Close the unneeded ends of Pipes
        close(input.out)
        close(output.in)
        return new(input, output, process)
    end
end

Base.kill(rs::PSL) = kill(rs.process)
Base.process_exited(rs::PSL) = process_exited(rs.process)

export error, ReduceError
import Base: error

struct ReduceError <: Exception
    errstr::String
end

Base.showerror(io::IO, err::ReduceError) = print(io,"Reduce: "*chomp(err.errstr))

function ReduceCheck(output) # check for REDUCE errors
    occursin(r"(([*]{5})|([+]{3}) )|( ?  \(Y or N\))",output) && throw(ReduceError(output))
end

function ReduceWarn(output) # check for REDUCE warnings
    if occursin(r"[*]{3}",output)
        @info "REDUCE: $(chomp(chomp(output)))"
        join(split(output,r"[*]{3}.*\n"))
    else
        output
    end
end

function PipeClogged(tf::Bool,c::Int,info::String)
    @warn "Reduce pipe clogged by $info, $(tf ? "success" : "failure") after $c tries"
end

clear(rs::PSL) = (write(rs.input,";\n"); readavailable(rs.output))
clears = (()->(c=true; return (tf=c)->(c≠tf && (c=tf); return c)))()

const EOT = Char(4) # end of transmission character
EOTstr = "symbolic write(int2id $(Int(EOT)))"

function Base.write(rs::PSL, input::String)
    clears() && clear(rs)
    write(rs.input,"$input; $EOTstr;\n")
end

const SOS = "[0-9]+: " # REDUCE terminal prompt
const RES = Regex("\n($EOT\n$SOS)|(\n$SOS\n$EOT)|(\n$SOS$EOT\n)|($EOT\n)")

function Base.read(rs::PSL) # get result and strip prompts/EOT char
    out = String(readuntil(rs.output,EOT))*String(readavailable(rs.output))
    Sys.iswindows() && (out = replace(out,r"\r" => ""))
    out = replace(replace(out,r"\$\n\n" => "\n\n"),RES=>"")
    out = replace(out,Regex(SOS) => "") |> chomp |> chomp |> join
    ReduceCheck(out)
    return ReduceWarn(out)
end

readsp(rs::PSL) = split(read(rs),"\n\n\n")

include("rexpr.jl") # load RExpr features
include("parser.jl") # load parser generator
include("repl.jl") # load repl features
include("switch.jl") # load switch operators

module Algebra
using Reduce, LinearAlgebra
import AbstractTensors
import AbstractTensors: conj, inv, dot, PROD, SUM, -, /, √
import AbstractTensors: sqrt, abs, exp, expm1, log, log1p, sin, cos, sinh, cosh, ^
const *,+ = PROD,SUM

export expm1, log1p, √
for T ∈ (RExpr,Expr,Symbol,String)
    @eval begin
        Base.expm1(z::$T) = exp(z)-1
        Base.log1p(z::$T) = 1+log(z)
        dot(a::$T,b) = conj(a)*b
        dot(a,b::$T) = conj(a)*b
        dot(a::$T,b::$T) = conj(a)*b
    end
end
dot(a::Expr,b::Symbol) = conj(a)*b
dot(a::Symbol,b::Expr) = conj(a)*b

include("unary.jl") # load unary operators
include("args.jl") # load calculus operators

function init_subtype(name)
    Expr(:block,[:(Base.$i(r::$name...)=$i(RExpr.(r)...)|>$name) for i ∈ [alg;iops]]...) |> eval
    Expr(:block,[:($i(r::$name...)=$i(RExpr.(r)...)|>$name) for i ∈ [calculus;cnan;cmat]]...) |> eval
    Expr(:block,[:(Base.$i(r::$name)=$i(RExpr(r))|>$name) for i ∈ [sbas;[:length]]]...) |> eval
    Expr(:block,[:($i(r::$name)=$i(RExpr(r))|>$name) for i ∈ [sfun;snan;snum;scom;sint;sran;smat]]...) |> eval
    :(Base.promote(r::$name,x) = (r,$name(x))) |> eval
    :(Base.promote(x,r::$name) = (r,$name(x))) |> eval
    :(Base.signbit(r::$name) = false) |> eval
end

const variables = [
    :root_multiplicities,
    :requirements,
    :assumptions,
    :low_pow,
    :high_pow,
    :euler_gamma,
    :khinchin,
]

for var ∈ [variables;[:ws]]
    :($var() = rcall(RExpr($(string(var)))) |> parse) |> eval
end

function bypass(op::Symbol,T::ExprSymbol,args)
    1 ∈ args && (@eval $op(x::X) where X<:$T = Base.$op(x))
    if 2 ∈ args
        @eval begin
            $op(a::A,b::B) where {A<:ExprSymbol,B<:$T} = Base.$op(a,b)
            $op(a::A,b::B) where {A<:$T,B<:ExprSymbol} = Base.$op(a,b)
        end
    end
end

import AbstractTensors: TensorAlgebra, value
for op ∈ (:exp,:reverse)
    bypass(op,:TensorAlgebra,(1,))
end
for op ∈ (:+,:*,:/)
    bypass(op,:TensorAlgebra,(2,))
end
@eval bypass(:-,:TensorAlgebra,(1,2))
rank(x::X) where X<:TensorAlgebra = LinearAlgebra.rank(x)

@doc """
    ws()

Several mechanisms are available for saving and retrieving previously evaluated expressions. The simplest of these refers to the last algebraic expression simplified. When an assignment of an algebraic expression is made, or an expression is evaluated at the top level, (i.e., not inside a compound statement or procedure) the results of the evaluation are automatically saved in a variable `ws` that we shall refer to as the workspace. (More precisely, the expression is assigned to the variable `ws` that is then available for further manipulation.)

*Example:* If we evaluate the expression `(x+y)^2` at the top level and next wish to differentiate it with respect to `y`, we can simply say
``
julia> Algebra.df(Algebra.ws(),:y)
```
to get the desired answer.
""" Reduce.Algebra.ws

@doc """
    root_multiplicities()

Solution multiplicities are stored in the global variable `root_multiplicities` rather than the solution list. The value of this variable is a list of the multiplicities of the solutions for the last call of `solve`. For example,
```Julia
julia> Algebra.solve(:(x^2==2x-1),:x); Algebra.root_multiplicities()
```
gives the results
```Julia
(:(x = 1),)

(2,)
```
""" Reduce.Algebra.root_multiplicities

@doc """
    requirements()

The proper design of a variable sequence supplied as a second argument to `solve` is important for the structure of the solution of an equation system. Any unknown in the system not in this list is considered totally free. E.g. the call
```Julia
Algebra.solve((:(x==2z),:(z==2y)),(:z,))
```
produces an empty list as a result because there is no function `z = z(x,y)` which fulfills both equations for arbitrary `x` and `y` values. In such a case the share variable `requirements` displays a set of restrictions for the parameters of the system:

```Julia
julia> Algebra.requirements()
(:(x - 4y),)
```

The non-existence of a formal solution is caused by a contradiction which disappears only if the parameters of the initial system are set such that all members of the requirements list take the value zero. For a linear system the set is complete: a solution of the requirements list makes the initial system solvable. E.g. in the above case a substitution ``x = 4y`` makes the equation set consistent. For a non-linear system only one inconsistency is detected. If such a system has more than one inconsistency, you must reduce them one after the other. 1 The set shows you also the dependency among the parameters: here one of ``x`` and ``y`` is free and a formal solution of the system can be computed by adding it to the variable list of `solve`. The requirement set is not unique – there may be other such sets.
""" Reduce.Algebra.requirements

@doc """
    assumptions()

A system with parameters may have a formal solution, e.g.
```Julia
julia> Algebra.solve((:(x==a*z+1),:(0==b*z-y)),(:z,:x))
(:(z = y // b), :(x = (a * y + b) // b))
```
which is not valid for all possible values of the parameters. The variable `assumptions` contains then a list of restrictions: the solutions are valid only as long as none of these expressions vanishes. Any zero of one of them represents a special case that is not covered by the formal solution. In the above case the value is
```Julia
julia> Algebra.assumptions()
(:b,)
```
which excludes formally the case ``b = 0``; obviously this special parameter value makes the system singular. The set of assumptions is complete for both, linear and non–linear systems.
""" Reduce.Algebra.assumptions

@doc """
    euler_gamma()

Euler's constant, also available as -\$\\psi(1)\$.
""" Reduce.Algebra.euler_gamma

@doc """
    khinchin()

Khinchin's constant, defined as

\$ \\prod_{n=1}^\\infty \\left( 1 + \\frac{1}{n(n+2)} \\right)^{\\log_2 n}. \$
""" Reduce.Algebra.khinchin

end

export Algebra, @force, @vars

macro vars(syms...)
    Expr(:block,[Expr(:(=),esc(sym),QuoteNode(sym)) for sym ∈ syms]...)
end

Base.write(rs::PSL,r::RExpr) = write(rs,convert(String,r))

import Base: zero, one

for T ∈ [:Any,:Expr,:Symbol]
    @eval begin
        zero(::Type{$T}) = 0
        zero(::$T) = 0
        one(::Type{$T}) = 1
        one(::$T) = 1
    end
end

zero(::Type{RExpr}) = R"0"
zero(::RExpr) = R"0"
one(::Type{RExpr}) = R"1"
one(::RExpr) = R"1"

import LinearAlgebra: transpose, adjoint

transpose(r::ExprSymbol) = r
transpose(r::RExpr) = r
adjoint(r::ExprSymbol) = Algebra.conj(r)
adjoint(r::RExpr) = Algebra.conj(r)

## Setup

const offlist = [:nat,:latex,:exp]

export load_package, @load_package

"""
    load_package(::Symbol)

Loads the specified package into REDUCE

## Examples
```julia-repl
julia> load_package(:rlfi)
```
"""
function load_package(pkg::Union{String,Symbol},pkgs...)
    "load_package $pkg" |> rcall
    for extra in pkgs
        load_package(extra)
    end
    return nothing
end
function load_package(pkgs::Union{Array{String,1},Array{Symbol,1}})
    for pkg in pkgs
        load_package(pkg)
    end
end

macro load_package(pkg...)
    load_package(pkg...)
end

"""
    Reduce.Reset()

Kills the REDUCE process and starts a new instance.

## Examples
```julia-repl
julia> Reduce.Reset()
Reduce (Free PSL version, revision 4015),  5-May-2017 ...
```
"""
Reset() = (kill(rs); Load())
__init__() = (Load(); repl_init(); atexit(() -> kill(rs)))

# Server setup

function Load()
    global rs = PSL()

    offs = ""
    for o in offlist
        o != :nat && (offs = offs*"off $o; ")
    end
    write(rs.input,"off nat; $EOTstr;\n")
    banner = readuntil(rs.output,EOT) |> String
    readavailable(rs.output)
    rcsl = occursin(" CSL ",banner)
    if Sys.iswindows()
        banner = replace(banner,r"\r" => "")
    else
        ReduceCheck(banner)
    end
    global banner = split(String(banner),'\n')[rcsl ? 1 : end-3]
    load_package(:rlfi)
    offs |> RExpr |> rcall
    rcall(R"on savestructr")
    return nothing
end

global preload = false
try
    global preload
    (ENV["REDPRE"] ≠ "0") && (preload = true)
catch
end
preload && include("precomp.jl")

end # module


text:
#Stochastic collocation towards a polynomial and the lognormal density
using Polynomials
using Roots
import AQFED.Math: normcdf, normpdf, norminv
import AQFED.Black: blackScholesFormula, blackScholesVega, impliedVolatility
using AQFED.Bachelier
using LeastSquaresOptim
struct LogPolynomial
    p::AbstractPolynomial
    σ::Number
end
struct IsotonicLogCollocation
    p1::AbstractPolynomial
    p2::AbstractPolynomial
    σ::Number
    forward::Number
end

function LogPolynomial(iso::IsotonicLogCollocation; minSlope = 0.0)
    q = degree(iso.p1) + degree(iso.p2)
    p = integrate(iso.p1^2 + Polynomials.Polynomial([0.0, 1.0])*iso.p2^2 + minSlope * iso.forward)
    theoForward = hermiteIntegral(LogPolynomial(p,iso.σ))
    c =  iso.forward / theoForward
    p = p*c
    return LogPolynomial(p,iso.σ)
end

(p::LogPolynomial)(x::S) where {S} = p.p(x)

function solvePositiveStrike(p::LogPolynomial, strike::Number; useNumerical=true)::Number
    if useNumerical

        return find_zero(
            x -> p(x) - strike,
            (1e-7,1e5),
            Bisection()
        )
    else
        r = roots(pk)
        for ri in r
            if (abs(imag(ri)) < eps(strike)) && real(ri) >= 0
                return real(ri)
            end
        end
        throw(DomainError(strike, "no real roots at the given strike"))
    end
end

function priceEuropean(p::LogPolynomial, isCall::Bool, strike::Number, forward::Number, discountDf::Number)::Number
    ck = solvePositiveStrike(p, strike)
    if ck < 0
        throw(DomainError(ck, "expected a positive value corresponding to strike"))
    end
    logck = log(ck)/p.σ
    valuef = hermiteIntegralBounded(p, logck)
    valuek = normcdf(-logck)
    callPrice = valuef - strike * valuek
    putPrice = -(forward - strike) + callPrice
    if isCall
        return callPrice * discountDf
    else
        return putPrice * discountDf
    end
end

function density(p::LogPolynomial, strike::Number)::Number

    ck = solvePositiveStrike(p, strike)
    if ck < 0
        throw(DomainError(ck, "expected a positive value corresponding to strike"))
    end
    logck = log(ck)/p.σ
    #Phi(log(ck)/p.sigma) '
    dp = derivative(p.p)(ck)
    an = normpdf(logck) / ( dp * ck*p.σ )
    #num = (-2*priceEuropean(p, true, strike, 0.0, 1.0) +priceEuropean(p, true, strike+1e-4, 0.0, 1.0) +priceEuropean(p, true, strike-1e-4, 0.0, 1.0) )/(1e-8)
    return an
end

function rawMoment(p::LogPolynomial, moment::Int)::Number
    if moment == 0
        return 1.0
    end
    q = p.p
    for i = 2:moment
        q *= p.p
    end
    return hermiteIntegral(LogPolynomial(q,p.σ))
end

#return mean, standard dev, skew, kurtosis of the collocation
function stats(p::LogPolynomial)
    μ = hermiteIntegral(p)
    μ2 = hermiteIntegral(LogPolynomial((p.p - μ)^2,p.σ))
    μ3 = hermiteIntegral(LogPolynomial((p.p - μ)^3,p.σ))
    μ4 = hermiteIntegral(LogPolynomial((p.p - μ)^4,p.σ))

    skew = μ3 / μ2^1.5
    kurtosis = μ4 / μ2^2
    return μ, sqrt(μ2), skew, kurtosis
end

function hermiteIntegral(p::LogPolynomial)::Number
    c = coeffs(p.p)
    value = sum(c[i]*exp(((i-1)*p.σ)^2 / 2) for i=1:length(c))
    return value
end

function hermiteIntegralBounded(p::LogPolynomial, logck::Number)::Number
    c = coeffs(p.p)
    return sum(c[i]*exp(((i-1)*p.σ)^2 / 2)*normcdf(-logck+(i-1)*p.σ) for i=1:length(c))
end

function makeLogXFromUndiscountedPrices(strikesf::Vector{T}, pricesf::Vector{T}, σ::T) where {T}
    return makeXFromUndiscountedPrices(strikesf, pricesf, s -> exp(σ*norminv(1+s)))
end

function makeIsotonicLogCollocationGuess(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    deg = 3, #1 = Black (trivial, smoother if least squares iterations small). otherwise 3 is good.
) where{T}
    b = fitLogBlack(strikes, callPrices, weights, τ, forward, discountDf)
    if deg >= 3
        strikesf, pricesf, weightsf = filterConvexPrices(strikes, callPrices ./ discountDf, weights, forward)
        strikesf, pif, xif = makeLogXFromUndiscountedPrices(strikesf, pricesf, b.σ)
        # cubic = Polynomials.fit(xif, strikesf, 3; weights = weightsf) #FIXME origin must be = 0
        # if !isCubicMonotone(cubic)
        #     cubic = LeastSquaresCubicMurrayForward(xif, strikesf, weightsf, forward)
        # end
        # if deg == 3
        #     return IsotonicLogCollocation(cubic, b.σ, forward)
        # else
        #     return fitLogMonotonic(xif, strikesf, weightsf, forward, cubic, deg = deg)
        # end
        return fitLogMonotonic(xif, strikesf, weightsf, forward, b, deg = deg)
    else
        return b
    end
end

function makeIsotonicLogCollocation(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    deg = 3, #degree of collocation. 5 is usually best from a stability perspective.
    degGuess = 3, #1 = Bachelier (trivial, smoother if least squares iterations small). otherwise 3 is good.
    minSlope = 1e-4, penalty=0.0
)::Tuple{IsotonicLogCollocation,Number} where {T} #return collocation and error measure
    isoc = makeIsotonicLogCollocationGuess(strikes, callPrices, weights, τ, forward, discountDf, deg = degGuess)
    #optimize towards actual prices
    isoc, m = fit(isoc, strikes, callPrices, weights, forward, discountDf, deg = deg, minSlope=minSlope,penalty=penalty)
    return isoc, m
end

function fitLogBlack(strikes, prices, weights, τ, forward, discountDf)
    m = length(strikes)
    i = findfirst(x -> x > forward, strikes)
    if i == nothing
        i = m
    elseif i == 1
        i = 2
    end
    price = (prices[i] * (forward - strikes[i-1]) + prices[i-1] * (strikes[i] - forward)) / (strikes[i] - strikes[i-1])
    strike = forward
    vol = impliedVolatility(true, price, forward, strike, τ, discountDf)
    σ = vol*sqrt(τ)
    isoc = IsotonicLogCollocation(Polynomials.Polynomial([sqrt(forward*exp(-σ^2 / 2))]), Polynomials.Polynomial([0.0]), σ, forward)
    return isoc
end

Base.length(p::LogPolynomial) = Base.length(p.p)
Base.iterate(p::LogPolynomial) = Base.iterate(p.p)
Base.iterate(p::LogPolynomial, state) = Base.iterate(p.p, state)
Base.broadcastable(p::LogPolynomial) = Ref(p)
Polynomials.degree(p::LogPolynomial) = Polynomials.degree(p.p)

isotonicDegree(deg::Int) = isodd(deg) ? trunc(Int, (deg-1)/2) : trunc(Int, deg / 2)

function fit(isoc::IsotonicLogCollocation, strikes, prices, weights, forward, discountDf; deg = 3, minSlope = 1e-6, penalty=0.0)
    q = isotonicDegree(deg)
    c0length = 2*q
    if isodd(deg)
        q+=1
        c0length = 2*q-1
    end
    iter = 0
    function obj(ct)
        c = ct
        p1 = Polynomials.Polynomial(c[1:q])
        p2 = Polynomials.Polynomial(c[q+1:c0length])
        isoc = IsotonicLogCollocation(p1, p2, isoc.σ, forward)
        p = LogPolynomial(isoc)
        iter += 1
        if penalty > 0
            ip = hermiteIntegral(LogPolynomial(derivative(p.p,2)^2,p.σ)) #does not seem appropriate here
            pvalue = penalty * ip
            return @. weights * sqrt((priceEuropean(p, true, strikes, forward, discountDf) - prices)^2 + pvalue^2)
        else
            return @. weights * (priceEuropean(p, true, strikes, forward, discountDf) - prices)
        end
    end
    c0 = zeros(Float64, c0length)
    c1 = coeffs(isoc.p1)
    for i = 1:min(q, length(c1))
        c0[i] = c1[i]
    end
    c2 = coeffs(isoc.p2)
    for i = 1:min(c0length-q, length(c2))
        c0[q+i] = c2[i]
    end
    for i = 1:length(c0)
        if c0[i] == 0
            c0[i] = eps(forward)  #zero would break automatic differentiation
        end
    end
    fit = optimize(obj, c0, LevenbergMarquardt(); show_trace = false, autodiff = :forward,iterations = deg * 300) #autodiff breaks without Halley. Would need custom

    println(iter, " fit ", fit, obj(fit.minimizer))
    c0 = fit.minimizer
    measure = fit.ssr
    return IsotonicLogCollocation(Polynomials.Polynomial(c0[1:q]), Polynomials.Polynomial(c0[q+1:c0length]), isoc.σ, forward), measure
end

function IsotonicLogCollocation(cguess::LogPolynomial, forward::Number)
    cubic = cguess.p
    if degree(cubic) == 1
        return IsotonicLogCollocation(Polynomials.Polynomial([sqrt(cguess.p[1])]), Polynomials.Polynomial([0.0]), cguess.σ, forward)
    elseif degree(cubic) > 3
        throw(DomainError(degree(cubic), "expected a cubic"))
    end
    e = (cubic[1] - cubic[2]^2 / (3 * cubic[3]))
    e = max(e, forward * 1e-6)
    p2 = Polynomials.Polynomial([sqrt(e)])
    p1 = Polynomials.Polynomial([cubic[2] / sqrt(3 * cubic[3]), sqrt(3 * cubic[3])])
    return IsotonicLogCollocation(p1, p2, cguess.σ, forward)
end

function fitLogMonotonic(xif, strikesf, w1, forward, cubic; deg = 3)
    q = isotonicDegree(deg)
    c0length = 2*q
    if isodd(deg)
        q+=1
        c0length = 2*q-1
    end
    iter = 0
    function obj(c)
        #c = vcat(0.0, ct)
        p1 = Polynomials.Polynomial(c[1:q])
        p2 = Polynomials.Polynomial(c[q+1:c0length])
        isoc = IsotonicLogCollocation(p1, p2, cubic.σ, forward)
        p = LogPolynomial(isoc)
        iter += 1
        v = @. w1 * (p(xif) - strikesf)
        #println(iter, " ", c, " ", p.p, sum(v.^2))
        return v
    end

    isocubic = typeof(cubic) == IsotonicLogCollocation ? cubic : IsotonicLogCollocation(cubic, forward)
    c1 = coeffs(isocubic.p1)
    c2 = coeffs(isocubic.p2)
    c0 = zeros(Float64,c0length)
    c0[1:min(q, length(c1))] = c1
    c0[q+1:min(c0length, q + length(c2))] = c2
    # println("c0 ",c0," isocub ",isocubic)
    xamax = xif[end]
    for i = 1:length(c0)
        if c0[i] == 0
            c0[i] = eps(forward) / (xamax)^(i - 1)  #zero would break automatic differentiation
        end
    end
    fit = optimize(obj, c0, LevenbergMarquardt(); show_trace = false, autodiff = :forward)
    c0 = fit.minimizer
    isoc = IsotonicLogCollocation(Polynomials.Polynomial(c0[1:q]), Polynomials.Polynomial(c0[q+1:c0length]), cubic.σ, forward)
    println(iter, " fitMonotonic ", LogPolynomial(isoc), " ", fit)
    return isoc
end


text:
function Base.write(prob,eqs,sol;dn::String="",fn::String)
    mkpath(dn)
    NPZ.npzwrite(dn*fn*".npz",merge(dumpscalars(prob,sol),
                                    dumpfields(prob,sol),
                                    dumpstats(prob,eqs,sol)))
                                    # dumpcoeffs(prob,eqs,sol)))
end

Base.write(prob,eqs::Vector{AbstractEquations},sols;dn::String,labels::Vector{String}=label(eqs)) = foreach(x->write(prob,x[1],x[2],dn=dn,fn=x[3]),zip(eqs,sols,labels))

tonpz(u) = reshape(cat(u...,dims=length(size(u[1]))),size(u[1])...,length(u))

function dumpcoeffs(prob,eqs::CE2,sol)
    x = zeros(eqs,prob.d)
    x.x[2] .= fcoeffs2(prob,eqs)
    F = forcingspectrum(prob.d,x)
    Dict("t"=>sol.t,"F"=>F)
end

function dumpcoeffs(prob,eqs,sol)
    x = zeros(eqs,prob.d)
    x .= fcoeffs2(prob,eqs)
    F = forcingspectrum(prob.d,x)
    Dict("t"=>sol.t,"F"=>F)
end

function dumpscalars(prob,sol;t0=50.0)
    Et = energy.(Ref(prob.d),sol.u)
    Zt = enstrophy.(Ref(prob.d),sol.u)
    Emt = zonalenergy.(Ref(prob.d),sol.u) |> tonpz
    Em0nt = modalenergy.(Ref(prob.d),sol.u,m=0) |> tonpz
    Em1nt = modalenergy.(Ref(prob.d),sol.u,m=1) |> tonpz
    Em2nt = modalenergy.(Ref(prob.d),sol.u,m=2) |> tonpz
    Em0ntav = modalenergy.(Ref(prob.d),sol.u,m=0) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Em1ntav = modalenergy.(Ref(prob.d),sol.u,m=1) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Em2ntav = modalenergy.(Ref(prob.d),sol.u,m=2) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Emtav = zonalenergy.(Ref(prob.d),sol.u) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Dict("t"=>sol.t,"Et"=>Et,"Emt"=>Emt,"Zt"=>Zt,"Emtav"=>Emtav,"Em0nt"=>Em0nt,"Em1nt"=>Em1nt,"Em2nt"=>Em2nt,
        "Em0ntav"=>Em0ntav,"Em1ntav"=>Em1ntav,"Em2ntav"=>Em2ntav)
end

function dumpfields(prob,sol;t0=50.0)
    Emn = energyspectrum.(Ref(prob.d),sol.u) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Fyym1 = secondcumulant.(Ref(prob.d),sol.u,m=1) |> tonpz
    Fyym2 = secondcumulant.(Ref(prob.d),sol.u,m=2) |> tonpz
    Fyym3 = secondcumulant.(Ref(prob.d),sol.u,m=3) |> tonpz
    Vxyav = vorticity.(Ref(prob.d),sol.u) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Vxy = vorticity.(Ref(prob.d),sol.u) |> tonpz
    Uxy = xvelocity.(Ref(prob.d),sol.u) |> tonpz
    Vyt = zonalvorticity.(Ref(prob.d),sol.u) |> tonpz
    Vytav = zonalvorticity.(Ref(prob.d),sol.u) |> x->timeaverage(sol.t,x,t0=t0) |> tonpz
    Dict("t"=>sol.t,"Emn"=>Emn,"Vxy"=>Vxy,"Vxyav"=>Vxyav,"Uxy"=>Uxy,"Vyt"=>Vyt,"Vytav"=>Vytav,"Fyym1"=>Fyym1,"Fyym2"=>Fyym2,"Fyym3"=>Fyym3)
end

dumpstats(prob,eqs,sol) = Dict("empty"=>0)
dumpstats(prob,eqs::CE2,sol) = Dict("mEVs"=> modaleigvals.(Ref(prob.d),sol.u) |> tonpz)

function dumpadjacency(lx::T,ly::T,nx::Int,ny::Int;fs::String,Λ::Int=nx-1) where {T <: AbstractFloat}
    A,C = adjacency(lx,ly,nx,ny,Λ=Λ)
    d = Dict("A"=>A,"C"=>C)
    NPZ.npzwrite(fs*".npz",d)
end


text:
#
println("Dl) Test of the PhotoDoubleIonization module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-PhotoDoubleIonization.sum")
setDefaults("method: continuum, Galerkin")           ## setDefaults("method: continuum, Galerkin")  "method: continuum, asymptotic Coulomb"
                                                     ## setDefaults("method: normalization, Ong-Russek")
setDefaults("method: normalization, pure sine")      ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
setDefaults("unit: rate", "a.u.")

grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)

if  false
    # Green function for single-photon 2p^2 photoionization of neon
    name             = "Double Auger decay of 1s photoionized neon"
    refConfigs       = [Configuration("1s^2 2s^2 2p^6")]
    levelSymmetries  = [LevelSymmetry(0, Basics.minus), LevelSymmetry(1, Basics.minus)]
    greenSettings    = GreenSettings(3, [0, 1], 0.01, true, LevelSelection())
    greenRep         = Representation(name, Nuclear.Model(10.), Radial.Grid(true), refConfigs,
                                      GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron(), levelSymmetries, 3, greenSettings) )
    greenOut         = generate(greenRep, output=true)
    doubleGreen      = greenOut["Green channels"]

elseif true
    # Single-photon 2p^2 photoionization of neon
    doubleSettings   = PhotoDoubleIonization.Settings(EmMultipole[E1], UseGauge[UseCoulomb, UseBabushkin], [50.0, 70.], doubleGreen, 2,
                                                      false, true, 2, LineSelection(true, indexPairs=[(1,0)]))

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(10.),
                            initialConfigs  =[Configuration("1s^2 2s^2 2p^6")],
                            finalConfigs    =[Configuration("1s^2 2s^2 2p^4")],
                            process = PhotoDouble(),  processSettings = doubleSettings )

    wb = perform(wa)

end
setDefaults("print summary: close", "")


text:
import FastTransforms
import QuadratureRules: chebyshev_points, shift_nodes

@testset "$(rpad("Gauß-Chebyshev",80))" begin

    for s in 2:10
        @test GaussChebyshevQuadrature(s) == ChebyshevQuadrature(s, 1)
        @test GaussChebyshevQuadrature(s) == ChebyshevQuadrature(Float64, s, Val(1))

        @test sum(weights(GaussChebyshevQuadrature(s))) ≈ 1

        c = shift_nodes(reverse(FastTransforms.chebyshevpoints(Float64, s, Val(1))))

        @test chebyshev_points(Float64, s, Val(1)) ≈ c
        @test nodes(GaussChebyshevQuadrature(s))   ≈ c

        # TODO test weights

    end

end


text:
csv_dir = Pkg.dir("D4M") * "/examples/1Intro/2EdgeArt/Edge.csv"
print(csv_dir)

E = ReadCSV(csv_dir);                       # Read CSV file into associative array.

Ev = logical( E[:, StartWith("V,")] );        # Get vertices and convert to numbers.

Av = sqIn(Ev);                                 # Compute vertex adjacency graph.
print(Av);                               # Show in tabular form.

Ae = sqOut(Ev);                                # Compute edge adjacency graph.
print(Ae);                               # Show in tabular form.


text:
# Run all unit tests
using LinearLeastSquares
using Base.Test

include("helpers.jl")
include("add_subtract.jl")
include("multiply_divide.jl")
include("stack.jl")
include("diag.jl")
include("repmat.jl")
include("index.jl")
include("reshape.jl")
include("transpose.jl")
include("sum_squares.jl")
include("compound.jl")
include("constraints.jl")
include("problems.jl")


text:
using .RasterDataSources

using .RasterDataSources: RasterDataSource

const RDS = RasterDataSources

const LayerItr = Union{AbstractArray,Tuple}

"""
    geoarray(T::Type{<:RasterDataSource}, [layer]; kw...) => AbstractArray

Load a `RasterDataSource` as an `AbstractGeoArray`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoArray`.

# Keywords

- `month`: For `Climate` datasets
- `date`: For `Weather` datasets
- `res`: For datasets with multiple resolutions

Normal `geoarray` keywords are passed to the constructor.
"""
geoarray(T::Type{<:RasterDataSource}; kw...) = geoarray(T, first(RDS.layers(T)); kw...)
function geoarray(T::Type{<:RasterDataSource}, layer; kw...)
    rds_kw, gd_kw = _filterkw(kw)
    filename = getraster(T, layer; rds_kw...)
    geoarray(filename; name=_layerkey(T, layer), _sourcekw(T)..., gd_kw...)
end

"""
    stack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractGeoStack

Load a `RasterDataSource` as an `AbstractGeoStack`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoStack`.

# Keywords

- `month`: For `Climate` datasets
- `date`: For `Weather` datasets
- `res`: For datasets with multiple resolutions

Normal `stack` keywords are passed to the constructor.
"""
stack(T::Type{<:RasterDataSource}; kw...) = stack(T, RDS.layers(T); kw...)
stack(T::Type{<:RasterDataSource}, layer::Symbol; kw...) = stack(T, (layer,); kw...)
function stack(T::Type{<:RasterDataSource}, layers::LayerItr; childkwargs=(), kw...)
    rds_kw, gd_kw = _filterkw(kw)
    filenames = map(l -> getraster(T, l; rds_kw...), layers)
    stack(filenames; keys=_layerkey(T, layers), childkwargs=(; _sourcekw(T)..., childkwargs...), gd_kw...)
end

"""
    series(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractGeoSeries

Load a `RasterDataSource` as an `AbstractGeoSeries`. `T`, `args` are
are passed to `getraster`, while `kw` args are for both `getraster` and
`AbstractGeoSeries`.

# Keywords

- `month`: For `Climate` datasets
- `date`: For `Weather` datasets
- `res`: For datasets with multiple resolutions

Normal `stack` keywords are passed to the constructor.
"""
series(T::Type{<:RasterDataSource}; kw...) = series(T, RDS.layers(T); kw...)
series(T::Type{<:RasterDataSource}, layer::Symbol; kw...) = series(T, (layer,); kw...)
# Int month time-series
function series(T::Type{WorldClim{Climate}}, layers::LayerItr;
    res=RDS.defres(T), month=1:12, window=(), kw...
)
    timedim = Ti(month; mode=Sampled(span=Regular(1), sampling=Intervals(Start())))
    stacks = [stack(T, layers; res=res, month=m, window=window) for m in month]
    GeoSeries(stacks, timedim; kw...)
end
# DateTime time-series
function series(T::Type{<:Union{WorldClim{Weather},ALWB,AWAP}}, layers::LayerItr; date, window=(), kw...)
    step = _seriesstep(T)
    dates = RDS._date_sequence(date, step)
    timedim = Ti(dates; mode=Sampled(Ordered(), Regular(step), Intervals(Start())))
    stacks = [stack(T, layers; date=d, window=window) for d in dates]
    GeoSeries(stacks, timedim; kw...)
end

_sourcekw(T) = ()
_sourcekw(T::Type{AWAP}) = (crs=EPSG(4326),)

_layerkey(T::Type{<:RasterDataSource}, keys::LayerItr) = map(k -> _layerkey(T, k), keys)
_layerkey(T::Type{<:Union{CHELSA{BioClim},WorldClim{BioClim}}}, key::Int) = string("BIO", key)
_layerkey(T::Type{<:RasterDataSource}, key) = Symbol(key)

_seriesstep(T::Type{<:ALWB{M,P}}) where {M,P} = P(1)
_seriesstep(T::Type{<:WorldClim{<:Weather}}) = Month(1)
_seriesstep(T::Type{<:AWAP}) = Day(1)

function _filterkw(kw)
    rds = []; gd = []
    for p in kw
        dest = first(p) in (:date, :month, :res) ? rds : gd
        push!(dest, p)
    end
    rds, gd
end


text:
module Filters

export Deflate, Shuffle, Fletcher32, Szip, NBit, ScaleOffset, BloscFilter

import ..HDF5: Properties, h5doc, API

abstract type Filter end

"""
    FILTERS

Maps filter id to filter type.
"""
const FILTERS = Dict{API.H5Z_filter_t, Any}()

"""
    filterid(::F)
    filterid(F)

The internal filter id of a filter of type `F`.
"""
filterid(::F) where {F<:Filter} = filterid(F)

struct UnknownFilter <: Filter
    filter_id::API.H5Z_filter_t
    flags::Cuint
    data::Vector{Cuint}
    name::String
    config::Cuint
end

"""
    FilterPipeline(plist::DatasetCreateProperties)

The filter pipeline associated with `plist`. Acts like a `AbstractVector{Filter}`,
supporting the following operations:

- `length(pipeline)`: the number of filters.
- `pipeline[i]` to return the `i`th filter.
- `pipeline[FilterType]` to return a filter of type `FilterType`
- `push!(pipline, filter)` to add an extra filter to the pipeline.
- `append!(pipeline, filters)` to add multiple filters to the pipeline.
- `delete!(pipeline, FilterType)` to remove a filter of type `FilterType` from the pipeline.
- `empty!(pipeline)` to remove all filters from the pipeline.
"""
struct FilterPipeline{P<:Properties} <: AbstractVector{Filter}
    plist::P
end

function Base.length(f::FilterPipeline)
    API.h5p_get_nfilters(f.plist)
end
Base.size(f::FilterPipeline) = (length(f),)

function Base.getindex(f::FilterPipeline, i::Integer)
    id = API.h5p_get_filter(f.plist, i-1, C_NULL, C_NULL, C_NULL, 0, C_NULL, C_NULL)
    F = get(FILTERS, id, UnknownFilter)
    return getindex(f, F, i)
end

function Base.getindex(f::FilterPipeline, ::Type{UnknownFilter}, i::Integer, cd_values::Vector{Cuint} = Cuint[])
    flags = Ref{Cuint}()
    cd_nelmts = Ref{Csize_t}(length(cd_values))
    namebuf = Array{UInt8}(undef, 256)
    config = Ref{Cuint}()
    id = API.h5p_get_filter(f.plist, i-1, flags, cd_nelmts, cd_values, length(namebuf), namebuf, config)
    if cd_nelmts[] < length(cd_values)
        resize!(cd_values, cd_nelmts[])
        return getindex(f, UnknownFilter, i, cd_values)
    end
    resize!(namebuf, findfirst(isequal(0), namebuf)-1)
    resize!(cd_values, cd_nelmts[])
    return UnknownFilter(id, flags[], cd_values, String(namebuf), config[])
end

function Base.getindex(f::FilterPipeline, ::Type{F}, i::Integer) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        id = API.h5p_get_filter(f.plist, i-1, C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    @assert id == filterid(F)
    return ref[]
end
function Base.getindex(f::FilterPipeline, ::Type{F}) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        API.h5p_get_filter_by_id(f.plist, filterid(F), C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    return ref[]
end

function Base.empty!(filters::FilterPipeline)
    API.h5p_remove_filter(filters.plist, API.H5Z_FILTER_ALL)
    return filters
end
function Base.delete!(filters::FilterPipeline, ::Type{F}) where {F<:Filter}
    API.h5p_remove_filter(filters.plist, filterid(F))
    return filters
end
function Base.append!(filters::FilterPipeline, extra)
    for filter in extra
        push!(filters, filter)
    end
    return filters
end

include("builtin.jl")
include("blosc.jl")

end # module


text:
function base_sdp(solver, seed)

    Random.seed!(seed)
    if Base.libblas_name == "libmkl_rt"
        model = Model()
    else
        model = Model(solver=solver)
    end

    n = 2

    if Base.libblas_name == "libmkl_rt"
        @variable(model, X[1:n, 1:n], PSD)
    else
        @variable(model, X[1:n, 1:n], SDP)
    end
    @objective(model, Min, -3X[1,1]-4X[2,2])
    @constraint(model, 2X[1,1]+1X[2,2] <= 4)
    @constraint(model, 1X[1,1]+2X[2,2] <= 4)

    if Base.libblas_name == "libmkl_rt"
        JuMP.attach(model, solver)
    end
    teste = JuMP.solve(model)

    if Base.libblas_name == "libmkl_rt"
        @show XX = getvalue2.(X)
    else
        @show XX = getvalue.(X)
    end

end
function base_sdp2(solver, seed)

    Random.seed!(seed)
    if Base.libblas_name == "libmkl_rt"
        model = Model()
    else
        model = Model(solver=solver)
    end

    n = 4

    if Base.libblas_name == "libmkl_rt"
        @variable(model, X[1:n, 1:n], PSD)
    else
        @variable(model, X[1:n, 1:n], SDP)
    end
    @objective(model, Min, -3X[1,1]-4X[2,2])
    @constraint(model, 2X[1,1]+1X[2,2]+X[3,3] == 4)
    @constraint(model, 1X[1,1]+2X[2,2]+X[4,4] == 4)

    if Base.libblas_name == "libmkl_rt"
        JuMP.attach(model, solver)
    end
    teste = JuMP.solve(model)

    if Base.libblas_name == "libmkl_rt"
        @show XX = getvalue2.(X)
    else
        @show XX = getvalue.(X)
    end

end

getvalue2(var::JuMP.Variable) = (m=var.m;m.solverinstance.primal[m.solverinstance.varmap[m.variabletosolvervariable[var.instanceindex]]])