text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # Cutting stock

# This example solves the cutting stock problem (sometimes also called the
# cutting rod problem) using a column-generation technique. It is based on
# [https://doi.org/10.5281/zenodo.3329388](https://doi.org/10.5281/zenodo.3329388).

# Intuitively, this problem is about cutting large rolls of paper into smaller
# pieces. There is an exact demand of pieces to meet, and all rolls have the
# same size. The goal is to meet the demand while maximising the profits (each
# paper roll has a fixed cost, each sold piece allows earning some money),
# which is roughly equivalent to using the smallest amount of rolls
# to cut (or, equivalently, to minimise the amount of paper waste).

# This function takes five parameters:

#   * `maxwidth`: the maximum width of a roll (or length of a rod)
#   * `widths`: an array of the requested widths
#   * `rollcost`: the cost of a complete roll
#   * `demand`: the demand, in number of pieces, for each width
#   * `prices`: the selling price for each width

# Mathematically, this problem might be formulated with two variables:

#   * `x[i, j] ∈ ℕ`: the number of times the width `i` is cut out of the roll `j`
#   * `y[j] ∈ 𝔹`: whether the roll `j` is used

# Several constraints are needed:

#   * the demand must be satisfied, for each width `i`:
#     ∑j x[i, j] = demand[i]
#   * the roll size cannot be exceed, for each roll `j` that is used:
#     ∑i x[i, j] width[i] ≤ maxwidth y[j]

# If you want to implement this naïve model, you will need an upper bound on the
# number of rolls to use: the simplest one is to consider that each required
# width is cut from its own roll, i.e. `j` varies from 1 to ∑i demand[i].

# This example prefers a more advanced technique to solve this problem:
# column generation.

# It considers a different set of variables: patterns of width to cut a roll.
# The decisions then become the number of times each pattern is used (i.e. the
# number of rolls that are cut following this pattern).

# The intelligence comes from the way these patterns are chosen: not all of them
# are considered, but only the "interesting" ones, within the master problem.

# A "pricing" problem is used to decide whether a new pattern should be
# generated or not (it is implemented in the function `solve_pricing`).
# "Interesting" means, for a pattern, that the optimal solution may use this
# cutting pattern.

# In more detail, the solving process is the following. First, a series of dumb
# patterns are generated (just one width per roll, repeated until the roll is
# completely cut). Then, the master problem is solved with these first patterns
# and its dual solution is passed on to the pricing problem. The latter decides
# if there is a new pattern to include in the formulation or not; if so,
# it returns it to the master problem. The master is solved again, the new dual
# variables are given to the pricing problem, until there is no more pattern to
# generate from the pricing problem: all "interesting" patterns have been
# generated, and the master can take its optimal decision.

# In the implementation, the variables deciding how many times a pattern is
# chosen are called `θ`.

# For more information on column-generation techniques applied on the cutting
# stock problem, you can see:

#   * [Integer programming column generation strategies for the cutting stock
#     problem and its variants](https://tel.archives-ouvertes.fr/tel-00011657/document)
#   * [Tackling the cutting stock problem](https://doi.org/10.5281/zenodo.3329388)

# This example uses the following packages:

using JuMP
import GLPK
import SparseArrays
import Test  #src

# The function `solve_pricing` implements the pricing problem for the function
# `example_cutting_stock`.
#
# It takes, as input, the dual solution from the master problem and the cutting
# stock instance.
#
# It outputs either a new cutting pattern, or `nothing` if no pattern could
# improve the current cost.

function solve_pricing(
    dual_demand_satisfaction, maxwidth, widths, rollcost, demand, prices
)
    reduced_costs = dual_demand_satisfaction + prices
    n = length(reduced_costs)
    ## The actual pricing model.
    submodel = Model(GLPK.Optimizer)
    set_silent(submodel)
    @variable(submodel, xs[1:n] >= 0, Int)
    @constraint(submodel, sum(xs .* widths) <= maxwidth)
    @objective(submodel, Max, sum(xs .* reduced_costs))
    optimize!(submodel)
    new_pattern = round.(Int, value.(xs))
    net_cost = rollcost - sum(new_pattern .* (dual_demand_satisfaction .+ prices))
    ## If the net cost of this new pattern is nonnegative, no more patterns to add.
    return net_cost >= 0 ? nothing : new_pattern
end

function example_cutting_stock(; max_gen_cols::Int = 5_000)
    maxwidth = 100.0
    rollcost = 500.0
    prices = [
        167.0, 197.0, 281.0, 212.0, 225.0, 111.0, 93.0, 129.0, 108.0, 106.0,
        55.0, 85.0, 66.0, 44.0, 47.0, 15.0, 24.0, 13.0, 16.0, 14.0,
    ]
    widths = [
        75.0, 75.0, 75.0, 75.0, 75.0, 53.8, 53.0, 51.0, 50.2, 32.2,
        30.8, 29.8, 20.1, 16.2, 14.5, 11.0, 8.6, 8.2, 6.6, 5.1,
    ]
    demand = [
        38, 44, 30, 41, 36, 33, 36, 41, 35, 37,
        44, 49, 37, 36, 42, 33, 47, 35, 49, 42,
    ]
    nwidths = length(prices)
    n = length(widths)
    ncols = length(widths)
    ## Initial set of patterns (stored in a sparse matrix: a pattern won't
    ## include many different cuts).
    patterns = SparseArrays.spzeros(UInt16, n, ncols)
    for i = 1:n
        patterns[i, i] = min(
            floor(Int, maxwidth / widths[i]), round(Int, demand[i])
        )
    end
    ## Write the master problem with this "reduced" set of patterns.
    ## Not yet integer variables: otherwise, the dual values may make no sense
    ## (actually, GLPK will yell at you if you're trying to get duals for
    ## integer problems).
    m = Model(GLPK.Optimizer)
    set_silent(m)
    @variable(m, θ[1:ncols] >= 0)
    @objective(
        m,
        Min,
        sum(
            θ[p] * (rollcost - sum(patterns[j, p] * prices[j] for j = 1:n))
            for p = 1:ncols
        )
    )
    @constraint(
        m,
        demand_satisfaction[j=1:n],
        sum(patterns[j, p] * θ[p] for p = 1:ncols) >= demand[j]
    )
    ## First solve of the master problem.
    optimize!(m)
    if termination_status(m) != MOI.OPTIMAL
        warn("Master not optimal ($ncols patterns so far)")
    end
    ## Then, generate new patterns, based on the dual information.
    while ncols - n <= max_gen_cols ## Generate at most max_gen_cols columns.
        if ! has_duals(m)
            break
        end
        new_pattern = solve_pricing(
            dual.(demand_satisfaction),
            maxwidth,
            widths,
            rollcost,
            demand,
            prices,
        )
        ## No new pattern to add to the formulation: done!
        if new_pattern === nothing
            break
        end
        ## Otherwise, add the new pattern to the master problem, recompute the
        ## duals, and go on waltzing one more time with the pricing problem.
        ncols += 1
        patterns = hcat(patterns, new_pattern)
        ## One new variable.
        new_var = @variable(m, [ncols], base_name = "θ", lower_bound = 0)
        push!(θ, new_var[ncols])
        ## Update the objective function.
        set_objective_coefficient(
            m,
            θ[ncols],
            rollcost - sum(patterns[j, ncols] * prices[j] for j = 1:n)
        )
        ## Update the constraint number j if the new pattern impacts this production.
        for j = 1:n
            if new_pattern[j] > 0
                set_normalized_coefficient(
                    demand_satisfaction[j], new_var[ncols], new_pattern[j]
                )
            end
        end
        ## Solve the new master problem to update the dual variables.
        optimize!(m)
        if termination_status(m) != MOI.OPTIMAL
            @warn("Master not optimal ($ncols patterns so far)")
        end
    end
    ## Finally, impose the master variables to be integer and resolve.
    ## To be exact, at each node in the branch-and-bound tree, we would need to
    ## restart the column generation process (just in case a new column would be
    ## interesting to add). This way, we only get an upper bound (a feasible
    ## solution).
    set_integer.(θ)
    optimize!(m)
    if termination_status(m) != MOI.OPTIMAL
        @warn("Final master not optimal ($ncols patterns)")
        return
    end
    Test.@test objective_value(m) ≈ 78374.0 atol = 1e-3  #src
    println("Final solution:")
    for i = 1:length(θ)
        if value(θ[i]) > 0.5
            println("$(round(Int, value(θ[i]))) units of pattern $(i)")
        end
    end
    return
end

example_cutting_stock()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

module DataLakeStoreAccountManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_AccountApi.jl")
include("api_FirewallRulesApi.jl")
include("api_LocationsApi.jl")
include("api_OperationsApi.jl")
include("api_TrustedIdProvidersApi.jl")

# export models
export convert, CapabilityInformation
export convert, CheckNameAvailabilityParameters
export convert, DataLakeStoreAccountListResult
export convert, DataLakeStoreAccountPropertiesBasic
export convert, DataLakeStoreAccountUpdateParameters
export convert, DataLakeStoreFirewallRuleListResult
export convert, DataLakeStoreTrustedIdProviderListResult
export convert, EncryptionConfig
export convert, EncryptionIdentity
export convert, ErrorDetails
export convert, FirewallRuleProperties
export convert, KeyVaultMetaInfo
export convert, NameAvailabilityInformation
export convert, Operation
export convert, OperationDisplay
export convert, OperationListResult
export convert, Resource
export convert, SubResource
export convert, TrustedIdProviderProperties
export convert, UpdateDataLakeStoreAccountProperties
export convert, UpdateEncryptionConfig
export convert, UpdateFirewallRuleParameters
export convert, UpdateFirewallRuleProperties
export convert, UpdateKeyVaultMetaInfo
export convert, UpdateTrustedIdProviderParameters
export convert, UpdateTrustedIdProviderProperties
export convert, DataLakeStoreAccount
export convert, DataLakeStoreAccountBasic
export convert, DataLakeStoreAccountProperties
export convert, FirewallRule
export convert, TrustedIdProvider

# export operations
export convert, AccountApi, FirewallRulesApi, LocationsApi, OperationsApi, TrustedIdProvidersApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end

ext: jl


text:
export unbiased_indmax, unbiased_indmin

#=
Agents often select an action that has the best experiential performance.
For example, they pick the action with the minimal expected cost or the
maximum expected reward. And, often, a vector stores the expectations.
This representation makes it tempting to use `Base.indmax`. However,
that function introduces a computational artifact. Given a tie, it
selects the highest index. `unbiased_indmax` and `unbiased_indmin`
operate without this bias. Given a tie, the both collect equal indices.
Then, a random selection over those indices gets returned. For now,
this requires two passes over the collection. Future reimplementations
could do away with such a limitation, but it's good for what I need
right now.
=#

function unbiased_indmax(xs)
    isempty(xs) && throw(ArgumentError("collection must be non-empty"))

    i = start(xs)
    max_i = i
    max_x, i = next(xs, i)
    tie_flag = false
    while !done(xs, i)
        last_i = i
        x, i = next(xs, i)

        if x > max_x
            max_x = x
            max_i = last_i
            tie_flag = false
        elseif x == max_x
            tie_flag = true
        end
    end

    if tie_flag
        return rand(collect(filter(x -> x[2] == max_x, enumerate(xs))))[1]
    else
        max_i
    end
end

function unbiased_indmin(xs)
    isempty(xs) && throw(ArgumentError("collection must be non-empty"))

    i = start(xs)
    min_i = i
    min_x, i = next(xs, i)
    tie_flag = false
    while !done(xs, i)
        last_i = i
        x, i = next(xs, i)

        if x < min_x
            min_x = x
            min_i = last_i
            tie_flag = false
        elseif x == min_x
            tie_flag = true
        end
    end

    if tie_flag
        rand(collect(filter(x -> x[2] == min_x, enumerate(xs))))[1]
    else
        min_i
    end
end


text:
"""
    RandomSequences(n::Int, sequence_length)

Indicates that resampling should be performed on discrete,
continuous subsets of a dataset. The lengths of each of the
`n` subsets is indicated by `sequence_length`, which should
be an integer (fixed sequence length) or an iterable of
integers (sequences may have different lengths).
"""
struct RandomSequences <: AbstractUncertainDataResampling
    n::Int
    sequence_length
end

export RandomSequences


text:
using VehicleCalculate
using Test

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

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


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

mutable struct VirtualMachineCaptureResultProperties <: SwaggerModel
    output::Any # spec type: Union{ Nothing, Any } # spec name: output

    function VirtualMachineCaptureResultProperties(;output=nothing)
        o = new()
        validate_property(VirtualMachineCaptureResultProperties, Symbol("output"), output)
        setfield!(o, Symbol("output"), output)
        o
    end
end # type VirtualMachineCaptureResultProperties

const _property_map_VirtualMachineCaptureResultProperties = Dict{Symbol,Symbol}(Symbol("output")=>Symbol("output"))
const _property_types_VirtualMachineCaptureResultProperties = Dict{Symbol,String}(Symbol("output")=>"Any")
Base.propertynames(::Type{ VirtualMachineCaptureResultProperties }) = collect(keys(_property_map_VirtualMachineCaptureResultProperties))
Swagger.property_type(::Type{ VirtualMachineCaptureResultProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineCaptureResultProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineCaptureResultProperties }, property_name::Symbol) =  _property_map_VirtualMachineCaptureResultProperties[property_name]

function check_required(o::VirtualMachineCaptureResultProperties)
    true
end

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


text:
export alloc_opt!, barrier_noop!, gc_invariant_verifier!, lower_exc_handlers!,
       combine_mul_add!, multi_versioning!, propagate_julia_addrsp!, lower_ptls!,
       lower_simdloop!, late_lower_gc_frame!, final_lower_gc!, remove_julia_addrspaces!

alloc_opt!(pm::PassManager) =
    API.LLVMAddAllocOptPass(ref(pm))

barrier_noop!(pm::PassManager) =
    API.LLVMAddBarrierNoopPass(ref(pm))

gc_invariant_verifier!(pm::PassManager, strong::Bool=false) =
    API.LLVMAddGCInvariantVerifierPass(ref(pm), convert(Bool, strong))

lower_exc_handlers!(pm::PassManager) =
    API.LLVMAddLowerExcHandlersPass(ref(pm))

combine_mul_add!(pm::PassManager) =
    API.LLVMAddCombineMulAddPass(ref(pm))

multi_versioning!(pm::PassManager) =
    API.LLVMAddMultiVersioningPass(ref(pm))

propagate_julia_addrsp!(pm::PassManager) =
    API.LLVMAddPropagateJuliaAddrspaces(ref(pm))

lower_ptls!(pm::PassManager, imaging_mode::Bool=false) =
    API.LLVMAddLowerPTLSPass(ref(pm), convert(Bool, imaging_mode))

lower_simdloop!(pm::PassManager) =
    API.LLVMAddLowerSimdLoopPass(ref(pm))

late_lower_gc_frame!(pm::PassManager) =
    API.LLVMAddLateLowerGCFramePass(ref(pm))

if VERSION >= v"1.3.0-DEV.95"
    final_lower_gc!(pm::PassManager) =
        API.LLVMAddFinalLowerGCPass(ref(pm))
else
    final_lower_gc!(pm::PassManager) = nothing
end

if VERSION >= v"1.5.0-DEV.802"
    remove_julia_addrspaces!(pm::PassManager) = API.LLVMAddRemoveJuliaAddrspacesPass(LLVM.ref(pm))
else
    remove_julia_addrspaces!(pm::PassManager) = nothing
end


text:
export
ruleSPNonlinearEOutNG,
ruleSPNonlinearEOutNGX,
ruleSPNonlinearEIn1GG,
ruleSPNonlinearEInGX,
ruleMNonlinearEInGX

"""
Concatenate a vector of vectors and return with original dimensions (for splitting)
"""
function concatenate(xs::Vector{Vector{Float64}})
    ds = [length(x_k) for x_k in xs] # Extract dimensions
    x = vcat(xs...)

    return (x, ds)
end

"""
Return local linearization of g around expansion point x_hat
"""
function localLinearization(V::Type{Univariate}, g::Function, x_hat::Float64)
    a = ForwardDiff.derivative(g, x_hat)
    b = g(x_hat) - a*x_hat

    return (a, b)
end

function localLinearization(V::Type{Multivariate}, g::Function, x_hat::Vector{Float64})
    A = ForwardDiff.jacobian(g, x_hat)
    b = g(x_hat) - A*x_hat

    return (A, b)
end

function localLinearization(V::Type{Univariate}, g::Function, x_hat::Vector{Float64})
    g_unpacked(x::Vector) = g(x...)
    A = ForwardDiff.gradient(g_unpacked, x_hat)'
    b = g(x_hat...) - A*x_hat

    return (A, b)
end

function localLinearization(V::Type{Multivariate}, g::Function, x_hat::Vector{Vector{Float64}})
    (x_cat, ds) = concatenate(x_hat)
    g_unpacked(x::Vector) = g(split(x, ds)...)
    A = ForwardDiff.jacobian(g_unpacked, x_cat)
    b = g(x_hat...) - A*x_cat

    return (A, b)
end

#-----------------------
# Extended Update Rules
#-----------------------

# Forward rule
function ruleSPNonlinearEOutNG(g::Function,
                               msg_out::Nothing,
                               msg_in1::Message{<:Gaussian, V}) where V<:VariateType

    (m_in1, V_in1) = unsafeMeanCov(msg_in1.dist)
    (A, b) = localLinearization(V, g, m_in1)

    return Message(GaussianMeanVariance, A*m_in1 + b, A*V_in1*A') # Automatically determine VariateType
end

# Multi-argument forward rule
function ruleSPNonlinearEOutNGX(g::Function, # Needs to be in front of Vararg
                                msg_out::Nothing,
                                msgs_in::Vararg{Message{<:Gaussian, V}}) where V<:VariateType

    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearization(V, g, ms_fw_in)
    (m_fw_in, V_fw_in, _) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)

    return Message(GaussianMeanVariance, A*m_fw_in + b, A*V_fw_in*A') # Automatically determine VariateType
end

# Backward rule with given inverse
function ruleSPNonlinearEIn1GG(g::Function,
                               g_inv::Function,
                               msg_out::Message{<:Gaussian, V},
                               msg_in1::Nothing) where V<:VariateType

    (m_out, V_out) = unsafeMeanCov(msg_out.dist)
    (A, b) = localLinearization(V, g_inv, m_out)

    return Message(GaussianMeanVariance, A*m_out + b, A*V_out*A') # Automatically determine VariateType
end

# Multi-argument backward rule with given inverse
function ruleSPNonlinearEInGX(g::Function, # Needs to be in front of Vararg
                              g_inv::Function,
                              msg_out::Message{<:Gaussian},
                              msgs_in::Vararg{Union{Message{<:Gaussian, V}, Nothing}}) where V<:VariateType

    (ms, Vs) = collectStatistics(msg_out, msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearization(V, g_inv, ms)
    (mc, Vc) = concatenateGaussianMV(ms, Vs)

    return Message(V, GaussianMeanVariance, m=A*mc, v=A*Vc*A')
end

# Backward rule with unknown inverse
function ruleSPNonlinearEIn1GG(g::Function,
                               msg_out::Message{<:Gaussian},
                               msg_in1::Message{<:Gaussian, V}) where V<:VariateType

    m_in1 = unsafeMean(msg_in1.dist)
    d_out = convert(ProbabilityDistribution{V, GaussianMeanPrecision}, msg_out.dist)
    m_out = d_out.params[:m]
    W_out = d_out.params[:w]
    (A, b) = localLinearization(V, g, m_in1)

    return Message(V, GaussianWeightedMeanPrecision, xi=A'*W_out*(m_out - b), w=A'*W_out*A)
end

# Multi-argument backward rule with unknown inverse
function ruleSPNonlinearEInGX(g::Function,
                              inx::Int64, # Index of inbound interface inx
                              msg_out::Message{<:Gaussian},
                              msgs_in::Vararg{Message{<:Gaussian, V}}) where V<:VariateType

    # Approximate joint inbounds
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearization(V, g, ms_fw_in)

    (m_fw_in, V_fw_in, ds) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)
    m_fw_out = A*m_fw_in + b
    V_fw_out = A*V_fw_in*A'
    C_fw = V_fw_in*A'

    # RTS Smoother
    (m_bw_out, V_bw_out) = unsafeMeanCov(msg_out.dist)
    (m_in, V_in) = smoothRTS(m_fw_out, V_fw_out, C_fw, m_fw_in, V_fw_in, m_bw_out, V_bw_out)

    # Marginalize joint belief on in's
    (m_inx, V_inx) = marginalizeGaussianMV(V, m_in, V_in, ds, inx) # Marginalization is overloaded on VariateType V
    W_inx = cholinv(V_inx) # Convert to canonical statistics
    xi_inx = W_inx*m_inx

    # Divide marginal on inx by forward message
    (xi_fw_inx, W_fw_inx) = unsafeWeightedMeanPrecision(msgs_in[inx].dist)
    xi_bw_inx = xi_inx - xi_fw_inx
    W_bw_inx = W_inx - W_fw_inx # Note: subtraction might lead to posdef inconsistencies

    return Message(V, GaussianWeightedMeanPrecision, xi=xi_bw_inx, w=W_bw_inx)
end

function ruleMNonlinearEInGX(g::Function,
                             msg_out::Message{<:Gaussian},
                             msgs_in::Vararg{Message{<:Gaussian, V}}) where V<:VariateType

    # Approximate joint inbounds
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearization(V, g, ms_fw_in)

    (m_fw_in, V_fw_in, ds) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)
    m_fw_out = A*m_fw_in + b
    V_fw_out = A*V_fw_in*A'
    C_fw = V_fw_in*A'

    # RTS Smoother
    (m_bw_out, V_bw_out) = unsafeMeanCov(msg_out.dist)
    (m_in, V_in) = smoothRTS(m_fw_out, V_fw_out, C_fw, m_fw_in, V_fw_in, m_bw_out, V_bw_out)

    return ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=m_in, v=V_in)
end


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking
using CSV, DataFrames
using StanSample, MonteCarloMeasurements

ProjDir = @__DIR__

df = DataFrame!(CSV.File(joinpath(ProjDir, "sports.csv"), delim=','));
scale!(df, [:C, :FL, :TM, :WUE, :PGP])
println()

m5 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] WUE; // Outcome
 vector[N] TM; // Predictor
 vector[N] PGP; // Predictor
 vector[N] FL; // Predictor
 vector[N] C; // Predictor
}

parameters {
 real a;     // Intercept
 real bPGP;  // Slope (regression coefficients)
 real bTM;   // Slope (regression coefficients)
 real bFL;   // Slope (regression coefficients)
 real bC;    // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;                 // mu is a vector
  a ~ normal(0, 0.2);           //Priors
  bPGP ~ normal(0, 0.5);
  bTM ~ normal(0, 0.5);
  bFL ~ normal(0, 0.5);
  bC ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bC * C + bPGP * PGP + bTM * TM + bFL * FL;
  WUE ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5s = SampleModel("m5", m5);

# Input data for cmdstan

m5_data = Dict("N" => size(df, 1), "PGP" => df[:, :PGP_s],
    "TM" => df[:, :TM_s], "WUE" => df[:, :WUE_s],
    "FL" => df[:, :FL_s], "C" => df[:, :C]);

# Sample using StanSample

rc = stan_sample(m5s, data=m5_data);

if success(rc)

  # Describe the draws

  dfa5 = read_samples(m5s; output_format=:dataframe)

  p5 = Particles(dfa5)

end


text:
using Test, bad

# Simons minimax design
# (Optimal Two-Stage Designs for Phase II Clinical Trials, Richard Simon)
# r1/n1 = 4/18; r/n = 10/33
function simon(r1::Int, n1::Int, r::Int, n::Int)
    nn2 = zeros(Int, n1 + 1)
    nn2[(r1 + 2):end] .= n - n1
    cc = repeat([Inf], n1 + 1)
    cc[(r1 + 2):end] .= r
    cc2 = cc .- collect(0:n1)
    Design(nn2, cc2)
end

pnull, palt  = .2, .4
α, β         = .05, .2
p            = Beta(1, 1)

design_simon = simon(4, 18, 10, 33)
essnull      = SampleSize(p | pnull)
@test essnull(design_simon) ≈ 22.3 atol = .1 # precision given in Simons table

problem = Problem(
        MiniMaxSampleSize(.2, p | pnull),
        Power(p | pnull) <= α,
        Power(p | palt)  >= 1 - β;
        # again, just removing as many as possible heuristics from the marginal fesible space
        n1values = collect(10:35),
        nmax = 35,
        n2mincontinueabs = 5,
        n2ton1fctrs      = (.9, 5.0),
        curtail_stage_one_buffer = 5
    )
design = optimise(problem; verbosity = 0)

# show that its actually better than the group sequential one ;)
@test maximum(n(design)) < maximum(n(design_simon))


text:
## deferred loading
##-----------------

require("Gadfly")
## require("Winston")

#####################
## Gadfly plotting ##
#####################

function gdfPlot(tn::AbstractTimenum, args...)
    nams = names(tn)
    tn = narowrm(tn)
    vals = asArr(tn, Float64, NaN) # NA to NaN
    dats = dat2num(tn)
    df = composeDataFrame([dats vals], [:Idx; nams])
    stackedData = stack(df, nams)
    Gadfly.plot(stackedData, x="Idx", y="value",
         color="variable", Gadfly.Geom.line, args...)
end

######################
## Winston plotting ##
######################

## ## no legend, but fast
## function wstPlot(tn::AbstractTimenum; title::String = "",
##                  xlabel::String = "time",
##                  ylabel::String = "value")
##     vals = asArr(tn, Float64, NaN) # transform to floating array
##     dats = dat2num(idx(tn))
##     Winston.plot(dats, vals)
##     Winston.title(title)
##     Winston.xlabel(xlabel)
##     Winston.ylabel(ylabel)
## end

###############
## histogram ##
###############

## function wstHist(x::Vector, nBins::Int = 20; title = "", xlabel = "", ylabel = "")
##     p = Winston.FramedPlot(title = title, xlabel = xlabel, ylabel = ylabel)
##     Winston.add(p, Winston.Histogram(hist(x[!isnan(x)], nBins)...))
## end

## function wstHist(da::DataArray, nBins::Int = 20; title = "", xlabel = "", ylabel = "")
##     x = asArr(da, Float64, NaN)
##     wstHist(x, nBins, title = title, xlabel = xlabel, ylabel = ylabel)
## end

## import Gadfly.plot
## function plot(tm::Timematr, settings...)
##     ## plot Timematr object
##     ## - multi-line plot
##     ## - index must be dates convertible to strings
##     ## - plot only some maximum number of randomly picked paths
##     ## - plot some statistics

##     dats = datsAsStrings(tm)
##     (nObs, nAss) = size(tm)

##     ## maximum number of paths plotted
##     maxPaths = 10

##     ## get values for paths
##     if nAss > maxPaths
##         ## get randomly drawn paths

##         chosenAssets = randperm(nAss)[1:maxPaths]
##         randPaths = core(tm)[:, chosenAssets]
##         nams = names(tm)[chosenAssets]

##         ## get statistics
##         allVals = core(tm)
##         averagePath = mean(allVals, 2)
##         lowestEndingPath = indmin(allVals[end, :])
##         lowPath = allVals[:, lowestEndingPath]
##         highestEndingPath = indmax(allVals[end, :])
##         highPath = allVals[:, highestEndingPath]

##         vals = [averagePath lowPath highPath randPaths]
##         valsDf = composeDataFrame(vals,
##                                   [:mean, :low, :high,
##                                    nams])
##     else
##         valsDf = copy(tm.vals)
##     end

##     nams = names(valsDf)
##     valsDf[:Idx] = dats

##     stackedData = stack(valsDf, nams)

##     plot(stackedData, x="Idx", y="value",
##          color="variable", Geom.line, settings...)
## end


text:
# include("sei3r.jl")
# include("evolution.jl")

using XLSX
using Plots
using DiseaseOutbreak
# TODO: make numerical tests to replace plot tests
using Test

# default parameters for RCAN model
d0 = getParams(5.0,6.0,0.5,0.1,0.1,0.15,0.06,6.0,8.0,0.4,SEI3R())
# social distancing reducing transmission probs by 1/2
d1 = getParams(5.0,6.0,0.5/2,0.1/2,0.1/2,0.15,0.06,6.0,8.0,0.4,SEI3R())

N = 1000.0
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,300)
plotEvolution(N,states0,d0)
# savefig("dist300.pdf")
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,60)
(states1,deltaStates1) = evolve(N,states0[:,60],d1,241)
plotEvolution(N,[states0 states1[:,2:end]],d0)
# savefig("dist60.pdf")
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,30)
(states1,deltaStates1) = evolve(N,states0[:,30],d1,271)
plotEvolution(N,[states0 states1[:,2:end]],d0)
# savefig("dist30.pdf")
state1 = initialize(1.0/N,d1)
(states1,deltaStates1) = evolve(N,state1,d1,300)
plotEvolution(N,states1,d1)
# savefig("dist0.pdf")

sheet = XLSX.readdata(joinpath("../memphis-covid-data",
             "covid_status_updates/",
             "COVID Status April 23.xlsx"),
             "Status","A1:C49")
using CSV
enc = CSV.read(joinpath("../memphis-covid-data",
               "mobility_data","unacast-shelby.csv"))
memCases = convert(Vector{Int64},sheet[2:49,2])
memPop = 1400000

fit = fitCaseModel(memCases[13:end],memPop,0.2,0.5,10.0,d0)
plot(memCases[13:end],legend=:bottomright,label="cumulative cases")
prd = predictCases(fit,60)
plot!(prd,label="fit")

length(memCases[13:end])
encMetric = reverse(enc.encountersMetric[1:36])
x = [ones(36) encMetric]
fit1 = fitCaseModel((memCases[13:end]),x[1:35,:],memPop,0.3,10.0,d0)
plot(memCases[13:end],legend=:bottomright,label="cumulative cases")
prd = predictCases(fit,60)
plot!(prd,label="fit0")

prd1 = predictCases(fit1,x[1:35,:])
plot!(prd1,label="fit1")

x2 = ones(35,2);
x2[1:13,2] = zeros(13)
fit2 = fitCaseModel(memCases[13:end],x2,memPop,0.01,10.0,d0)
prd2 = predictCases(fit2,x2)
plot!(prd2,label="fit2")

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

tn = CSV.read("/Users/sen/covid-memphis/tn.csv")

tnCases = reverse(tn.ConfirmedCases[1:42])
tnPop = 6800000
fittn = fitCaseModel(tnCases,tnPop,0.01,0.5,100.0,d0)
plot(tnCases,legend=:bottomright,yaxis=:log,label="cumulative cases",color=:black)
prdtn = predictCases(fittn,30)
plot!(prdtn,label="fit0",yaxis=:log)
xtn = ones(length(tnCases)-1,2)
xtn[1:10,2] = zeros(10)

fittn1 = fitCaseModel(tnCases,xtn,tnPop,0.01,500.0,d0)
xtn1 = ones(60,2)
xtn1[1:10,2] = zeros(10)
prdtn1 = predictCases(fittn1,xtn1)
plot!(prdtn1,label="fit1",yaxis=:log)

plot(diff(tnCases),legend=:bottomright)
plot!(diff(prdtn1))

fittn1.fit.minimizer

7*invlogit(0.895-2.50)
7*invlogit(0.895)

mob = hcat(ones(41),reverse(enc.encountersMetric[4:44]) )

fittn2 = fitCaseModel(tnCases,mob,tnPop,0.001,1000.0,d0)
prdtn2 = predictCases(fittn2,mob)
plot!(prdtn2,label="mob",yaxis=:log)


text:
function MOI.get(::Optimizer, ::MOI.ObjectiveFunctionType)
    return MOI.ScalarAffineFunction{Float64}
end
function MOI.get(m::Optimizer, ::MOI.ObjectiveFunction{F}) where F
    obj = MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    return convert(F, obj)
end
function MOI.get(m::Optimizer,
                 ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    cis = MOI.get(m, MOI.ListOfVariableIndices())
    cols = columns(m, cis).values
    coeffs = getclist(m.task, cols)
    constant = getcfix(m.task)
    @assert length(coeffs) == length(cis)
    terms = MOI.ScalarAffineTerm{Float64}[
        MOI.ScalarAffineTerm(coeffs[i], cis[i]) for i in 1:length(cis)]
    # TODO add matrix terms
    return MOI.ScalarAffineFunction(terms, constant)
end

const ObjF = Union{MOI.VariableIndex, MOI.ScalarAffineFunction{Float64}}
MOI.supports(::Optimizer,::MOI.ObjectiveFunction{<:ObjF})  = true
MOI.supports(::Optimizer,::MOI.ObjectiveSense) = true

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction,
                 func::MOI.VariableIndex)
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            convert(MOI.ScalarAffineFunction{Float64}, func))
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction,
                 func::MOI.ScalarAffineFunction{Float64})
    cols, values = split_scalar_matrix(m, MOIU.canonical(func).terms,
                                       (j, ids, coefs) -> putbarcj(m.task, j, ids, coefs))
    c = zeros(Float64, getnumvar(m.task))
    for (col, val) in zip(cols, values)
        c[col] += val
    end
    putclist(m.task, convert(Vector{Int32}, 1:length(c)), c)
    putcfix(m.task,func.constant)
    m.has_objective = true
    return
end

function MOI.modify(m::Optimizer,
                    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
                    change :: MOI.ScalarConstantChange)
    putcfix(m.task,change.new_constant)
    m.has_objective = true
    return
end

function MOI.modify(m::Optimizer,
                    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
                    change :: MOI.ScalarCoefficientChange)
    putcj(m.task, column(m, change.variable).value, change.new_coefficient)
    m.has_objective = true
    return
end


text:
# Copyright 2017 Jose Ortiz-Bejar
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module Nets
export Net, enet, kmnet, dnet, gen_features, KernelClassifier
import KernelMethods.Kernels: sigmoid, gaussian, linear, cauchy
import KernelMethods.Scores: accuracy, recall
import KernelMethods.Supervised: NearNeighborClassifier, NaiveBayesClassifier, optimize!, predict_one, predict_one_proba, LabelEncoder, transform, inverse_transform
using SimilaritySearch: KnnResult, L2Distance, L2SquaredDistance, CosineDistance, DenseCosine, JaccardDistance
using TextModel
using PyCall

@pyimport sklearn.naive_bayes as nb
@pyimport sklearn.model_selection as ms

#using JSON
#using DataStructures

mutable struct Net{ItemType,LabelType}
    data::Vector{ItemType}
    labels::Vector{Int}
    le::LabelEncoder{LabelType}
    references::Vector{Int32}
    partitions::Vector{Int32}
    centers::Vector{ItemType}
    centroids::Vector{ItemType}
    dists::Vector{Float64}
    csigmas::Vector{Float64}
    sigmas::Dict{Int,Float64}
    stats::Dict{LabelType,Float64}
    reftype::Symbol
    distance
    kernel
end

function Net(data::Vector{ItemType},labels::Vector{LabelType}) where {ItemType, LabelType}
    le = LabelEncoder(labels)
    y = transform.(le, labels)

    Net(data,y,le,
        Int32[],Int32[],ItemType[],
        ItemType[],Float64[],Float64[],
        Dict{Int,Float64}(), Dict{LabelType,Float64}(),:centroids,L2SquaredDistance(),gaussian)
end

function cosine(x1,x2)::Float64
    xc1=DenseCosine(x1)
    xc2=DenseCosine(x2)
    d=CosineDistance()
    return d(xc1,xc2)
end

function maxmin(data,centers,ind,index::KnnResult,distance,partitions)::Tuple{Int64,Float64}
    c=last(centers)
    if length(index)==0
        for i in ind
            if i!=c
                push!(index,i,Inf)
            end
        end
    end
    nindex=KnnResult(length(index))
    for fn in index
        dist=distance(data[fn.objID],data[c])
        #push!(lK[fn.objID],dist)
        dist = if (dist<fn.dist) dist else fn.dist end
        partitions[fn.objID] = if (dist<fn.dist) c else partitions[fn.objID] end
        if fn.objID!=c
            push!(nindex,fn.objID,dist)
        end
    end
    index.k=nindex.k
    index.pool=nindex.pool
    fn=pop!(index)
    return fn.objID,fn.dist
end

function get_centroids(data::Vector{T}, partitions::Vector{Int})::Vector{T} where T
    centers=[j for j in Set(partitions)]
    sort!(centers)
    centroids=Vector{T}(length(centers))
    for (ic,c) in enumerate(centers)
        ind=[i for (i,v) in enumerate(partitions) if v==c]
        centroids[ic]=mean(data[ind])
    end
    return centroids
end

# Epsilon Network using farthest first traversal Algorithm

function enet(N::Net,num_of_centers::Int; distance=L2SquaredDistance(),
              per_class=false,reftype=:centroids, kernel=linear)
    N.distance=distance
    N.kernel=kernel
    n=length(N.data)
    partitions=[0 for i in 1:n]
    gcenters,dists,sigmas=Vector{Int}(0),Vector{Float64}(num_of_centers-1),Dict{Int,Float64}()
    indices=[[i for i in  1:n]]
    for ind in indices
        centers=Vector{Int}(0)
        s=rand(1:length(ind))
        push!(centers,ind[s])
        #ll=N.labels[ind[s]]
        index=KnnResult(length(ind))
        partitions[ind[s]]=ind[s]
        k=1
        while  k<=num_of_centers-1 && k<=length(ind)
            fnid,d=maxmin(N.data,centers,ind,index,distance,partitions)
            push!(centers,fnid)
            dists[k]=d
            partitions[fnid]=fnid
            k+=1
        end
        sigmas[0]=minimum(dists)
        gcenters=vcat(gcenters,centers)
    end
    N.references,N.partitions,N.dists,N.sigmas=gcenters,partitions,dists,sigmas
    N.centers,N.centroids=N.data[gcenters],get_centroids(N.data,partitions)
    N.csigmas,N.stats=get_csigmas(N.data,N.centroids,N.partitions,distance=N.distance)
    N.reftype=reftype
end

# KMeans ++ seeding Algorithm

function kmpp(N::Net,num_of_centers::Int)::Vector{Int}
    n=length(N.data)
    s=rand(1:n)
    centers, d = Vector{Int}(num_of_centers), L2SquaredDistance()
    centers[1]=s
    D=[d(N.data[j],N.data[s]) for j in 1:n]
    for i in 1:num_of_centers-1
        cp=cumsum(D/sum(D))
        r=rand()
        sl=[j for j in 1:length(cp) if cp[j]>=r]
        s=sl[1]
        centers[i+1]=s
        for j in 1:n
            dist=d(N.data[j],N.data[s])
            if dist<D[j]
                D[j]=dist
            end
        end
    end
    centers
end

#Assign Elementes to thier nearest centroid

function assign(data,centroids,partitions;distance=L2SquaredDistance())
    d=distance
    for i in 1:length(data)
        partitions[i]=sortperm([d(data[i],c) for c in centroids])[1]
    end
end

#Distances for each element to its nearest cluster centroid

function get_distances(data,centroids,partitions;distance=L2SquaredDistance())::Vector{Float64}
    dists=Vector{Float64}(length(centroids))
    for i in 1:length(centroids)
        ind=[j for (j,l) in enumerate(partitions) if l==i]
        if length(ind)>0
            X=data[ind]
            dd=[distance(centroids[i],x) for x in X]
            dists[i]=maximum(dd)
        end
    end
    sort!(dists)
    return dists
end

#Calculated the sigma for each ball

function get_csigmas(data,centroids,partitions;distance=L2SquaredDistance())::Tuple{Vector{Float64},Dict{String,Float64}}
    stats=Dict("SSE"=>0.0,"BSS"=>0.0)
    refs=[j for j in Set(partitions)]
    sort!(refs)
    csigmas=Vector{Float64}(length(refs))
    df=distance
    m=mean(data)
    for (ii,i) in enumerate(refs)
        ind=[j for (j,l) in enumerate(partitions) if l==i]
        #if length(ind)>0
        X=data[ind]
        dd=[df(data[i],x) for x in X]
        csigmas[ii]=max(0,maximum(dd))
        stats["SSE"]+=sum(dd)
        stats["BSS"]+=length(X)*(sum(mean(X)-m))^2
        #end
    end
    return csigmas,stats
end

#Feature generator using kmeans centroids

function kmnet(N::Net,num_of_centers::Int; max_iter=1000,kernel=linear,distance=L2SquaredDistance(),reftype=:centroids)
    n=length(N.data)
    #lK,partitions=[[] for i in 1:n],[0 for i in 1:n],[0 for i in 1:n]
    partitions=[0 for i in 1:n]
    dists=Vector{Float64}
    init=kmpp(N,num_of_centers)
    centroids=N.data[init]
    i,aux=1,Vector{Float64}(length(centroids))
    while centroids != aux && i < max_iter
        i=i+1
        aux = centroids
        assign(N.data,centroids,partitions)
        centroids=get_centroids(N.data,partitions)
    end
    N.distance=distance
    dists=get_distances(N.data,centroids,partitions,distance=N.distance)
    N.partitions,N.dists=partitions,dists
    N.centroids,N.sigmas[0]=centroids,maximum(N.dists)
    N.csigmas,N.stats=get_csigmas(N.data,N.centroids,N.partitions,distance=N.distance)
    N.sigmas[0]=maximum(N.csigmas)
    N.reftype=:centroids
    N.kernel=kernel
end

#Feature generator using naive algorithm for density net

function dnet(N::Net,num_of_elements::Int64; distance=L2SquaredDistance(),kernel=linear,reftype=:centroids)
    n,d,k=length(N.data),distance,num_of_elements
    partitions,references=[0 for i in 1:n],Vector{Int}(0)
    pk=1
    dists,sigmas=Vector{Float64},Dict{Int,Float64}()
    while 0 in partitions
        pending=[j for (j,v) in enumerate(partitions) if partitions[j]==0]
        s=rand(pending)
        partitions[s]=s
        pending=[j for (j,v) in enumerate(partitions) if partitions[j]==0]
        push!(references,s)
        pc=sortperm([d(N.data[j],N.data[s]) for j in pending])
        if length(pc)>=k
            partitions[pending[pc[1:k]]]=s
        else
            partitions[pending[pc]]=s
        end
    end
    N.references,N.partitions=references,partitions
    N.centers,N.centroids=N.data[references],get_centroids(N.data,partitions)
    N.csigmas,N.stats=get_csigmas(N.data,N.centroids,N.partitions,distance=N.distance)
    N.sigmas[0]=maximum(N.csigmas)
    N.reftype=:centroids
    N.distance=distance
    N.kernel=kernel
end

#Generates feature espace using cluster centroids or centers

function gen_features(Xo::Vector{T},N::Net)::Vector{Vector{Float64}} where T
    n=length(Xo)
    sigmas,Xr=N.csigmas,Vector{T}(n)
    Xm = N.reftype==:centroids || length(N.centers)==0 ? N.centroids : N.centers
    nf=length(Xm[1])
    kernel=N.kernel
    for i in 1:n
        xd=Vector{Float64}(nf)
        for j in 1:nf
            xd[j]=kernel(Xo[i],Xm[j],sigma=sigmas[j],distance=N.distance)
        end
        Xr[i]=xd
    end
    Xr
end

function traintest(N; op_function=recall, runs=3, folds=0, trainratio=0.7, testratio=0.3)
    clf, avg = nb.GaussianNB(), Vector{Float64}(runs)
    skf = ms.ShuffleSplit(n_splits=runs, train_size=trainratio, test_size=testratio)
    X=gen_features(N.data,N)
    y=N.labels
    skf[:get_n_splits](X,y)
    for (ei,vi) in skf[:split](X,y)
        ei,vi=ei+1,vi+1
        xt,xv=X[ei],X[vi]
        yt,yv=y[ei],y[vi]
        clf[:fit](xt,yt)
        y_pred=clf[:predict](xv)
        push!(avg,op_function(yv,y_pred))
    end
    #println("========== ",length(N.centroids) ,"  ",avg/folds)
    #@show typeof(clf)
    clf[:fit](X,y)
    return clf,mean(avg)
end

#function transductive(){
#    continue
#}

function predict_test(xt,y,xv,desc,cl)::Vector{Int64}
    y_pred=[]
    if contains(desc,"KNN")
        #@show typeof(xt), typeof(y), typeof(cl.X.dist), Symbol(cl.X.dist)
        cln=NearNeighborClassifier(xt,y,cl.X.dist,cl.k,cl.weight)
        y_pred=[predict_one(cln,x)[1] for x in xv]
    else
        cln=nb.GaussianNB()
        cln[:fit](xt,y)
        y_pred=cln[:predict](xv)
    end
    y_pred
end

HammingDistance(x1,x2)::Float64 = length(x1)-sum(x1.==x2)

L2Squared = L2SquaredDistance()

function KlusterClassifier(Xe, Ye; op_function=recall,
                        K=[4, 8, 16, 32],
                        kernels=[:gaussian, :sigmoid, :linear, :cauchy],
                        runs=3,
                        trainratio=0.6,
                        testratio=0.4,
                        folds=0,
                        top_k=32,
                        threshold=0.03,
                        distances=[:cosine, :L2Squared],
                        nets=[:enet, :kmnet, :dnet], nsplits=3)::Vector{Tuple{Tuple{Any,Net},Float64,String}}
    top=Vector{Tuple{Float64,String}}(0)
    DNNC=Dict()

    for (k, nettype, reftype, kernel, distancek) in zip(K, nets, [:centers, :centroids], kernels, distances)
        if (distancek==:L2Squared || reftype==:centers) && nettype=="kmeans"
            continue
        else
            N=Net(Xe,Ye)
            eval(nettype)(N, k, kernel=eval(kernel), distance = eval(distancek), reftype=reftype)
            X=gen_features(N.data, N)
        end
        for distance in distances
            nnc = NearNeighborClassifier(X,Ye, eval(distance))
            opval,_tmp=optimize!(nnc, op_function,runs=runs, trainratio=trainratio,
                                    testratio=testratio,folds=folds)[1]
            kknn,w = _tmp
            key="$nettype/$kernel/$k/KNN$kknn/$reftype/$distance/$w"
            push!(top,(opval,key))
            DNNC[key]=(nnc,N)
        end
        key="$nettype/$kernel/$k/NaiveBayes/$reftype/NA"
        nbc,opval=traintest(N,op_function=op_function,folds=folds,
                            trainratio=trainratio,testratio=testratio)
        push!(top,(opval,key))
        DNNC[key]=(nbc,N)
    end

    sort!(top, rev=true)
    top=top[1:min(12, length(top))]
    # if top_k>0
    #     top=ctop[1:k]
    # else
    #     top=[t for t in top if (ctop[1][1]-t[1])<threshold]
    # end
    # @show length(top)
    LN=[(DNNC[t[2]],t[1],t[2]) for t in top ]
    #@show typeof(LN)
end

function ensemble_cfft(knc,k::Int64=7;testratio=0.4,distance=HammingDistance)::Vector{Tuple{Tuple{Any,Net},Float64,String}}
    (cl,n),opv,desc=knc[1]
    data=Vector{Vector{String}}(length(knc))
    for i in 1:length(knc)
        kc,opv,desc=knc[i]
        v=split(desc,"/")[1:6]
        data[i]=v
    end
    ind=[i for (i,x) in enumerate(data)]
    partitions,centers,index=[0 for x in ind],[1],KnnResult(k)
    while length(centers)<k && length(centers)<=length(ind)
        oid,dist=maxmin(data,centers,ind,index,distance,partitions)
        push!(centers,oid)
    end
    knc[centers]
end

function ensemble_pfft(knc,k::Int64=7;trainratio=0.6,distance=HammingDistance)::Vector{Tuple{Tuple{Any,Net},Float64,String}}
    (cl,N),opv,desc=knc[1]
    n=length(N.data)
    tn=Int(trunc(n*trainratio));
    data=Vector{Vector{Int}}(length(knc))
    perm=randperm(n)
    ti,vi=perm[1:tn],perm[tn+1:n]
    for i in 1:length(knc)
        (cl,N),opv,desc=knc[i]
        xv=gen_features(N.data[vi],N)
        xt=gen_features(N.data[ti],N)
        y=N.labels[ti]
        data[i]=predict_test(xt,y,xv,desc,cl)
    end
    ind=[i for (i,x) in enumerate(data)]
    partitions,centers,index=[0 for x in ind],[1],KnnResult(k)
    while length(centers)<k && length(centers)<=length(ind)
        oid,dist=maxmin(data,centers,ind,index,distance,partitions)
        push!(centers,oid)
    end
    knc[centers]
end

function predict(knc,X;ensemble_k=1)::Vector{Int64}
    y_t=Vector{Int}(0)
    for i in 1:ensemble_k
        kc,opv,desc=knc[i]
        cl,N=kc
        xv=gen_features(X,N)
        if contains(desc,"KNN")
            y_i=[predict_one(cl,x)[1] for x in xv]
        else
            y_i=cl[:predict](xv)
        end
        y_t = length(y_t)>0 ? hcat(y_t,y_i) : hcat(y_i)
    end
    y_pred=Vector{Int}(length(X))
    for i in 1:length(y_t[:,1])
        y_r=y_t[i,:]
        y_pred[i]=last(sort([(count(x->x==k,y_r),k) for k in unique(y_r)]))[2]
    end
    y_pred
end

function predict_proba(knc,X;ensemble_k=1):Vector{Vector{Float64}}
    kc,opv,desc=knc[1]
    cl,N=kc
    xv=gen_features(X,N)
    if contains(desc,"KNN")
        y_pred=[predict_one_proba(cl,x) for x in xv]
    else
        y_pred=cl[:predict_proba](xv)
    end
    y_pred
end

end


text:
macro define_trait(args)
    def = if @capture(args, trait_Symbol)
        trait == :AbstractTrait && error("Can not redefine the AbstractTrait")
        trait == :NullTrait && error("Can not redefine the NullTrait")
        @q begin
            abstract type $(esc(trait)){$(esc(:T))} <: AbstractTrait{$(esc(:T))} end
        end
    elseif @capture(args, trait_Symbol = block_)
        trait == :AbstractTrait && error("Can not redefine the AbstractTrait")
        trait == :NullTrait && error("Can not redefine the NullTrait")
        code = @q begin
            abstract type $(esc(trait)){$(esc(:T))} <: AbstractTrait{$(esc(:T))} end
        end
        for ex in rmlines(block).args
            if @capture(ex,subtrait_Symbol)
                subtrait == :AbstractTrait && error("Can not redefine the AbstractTrait")
                subtrait == :NullTrait && error("Can not redefine the NullTrait")
                subcode = @q begin
                    abstract type $(esc(subtrait)){$(esc(:T))} <: $(esc(trait)){$(esc(:T))} end
                end
                push!(code.args,unblock(subcode))
            elseif @capture(ex,f_(x__)) || @capture(ex,f_(x__) = body_)
                push!(code.args,esc(:(@define_traitfn($trait,$ex))))
            else
                error("Expected subtrait or traitfn definition, got ", ex)
            end
        end
        code
    elseif @capture(args, subtrait_Symbol <: trait_Symbol)
        subtrait == :AbstractTrait && error("Can not redefine the AbstractTrait")
        subtrait == :NullTrait && error("Can not redefine the NullTrait")
        @q begin
            $(esc(trait)) == NullTrait && error("The NullTrait can not have subtraits")
            !($(esc(trait)) <: AbstractTrait) && error($(esc(trait))," is not a trait")
            abstract type $(esc(subtrait)){$(esc(:T))} <: $(esc(trait)){$(esc(:T))} end
        end
    end
    prettify(def)
end

macro implement_trait(x,trait)
    code = @q begin
        # Sanity checks
        !(isa($(esc(x)),DataType) || isa($(esc(x)),UnionAll)) && error("Can only give traits to types")
        ($(esc(x)) <: AbstractTrait) && error($(esc(x))," is a trait, and can not implement traits")
        $(esc(trait)) == AbstractTrait && error("Cannot give the AbstractTrait")
        $(esc(trait)) == NullTrait && error("Cannot give the NullTrait")
        !($(esc(trait)) <: AbstractTrait) && error($(esc(trait))," is not a trait")
        !isleaftrait($(esc(trait))) && error("Trait must be a leaf trait, ",$(esc(trait))," has subtraits: ",[@sprintf("%s ",t) for t in subtraits($(esc(trait)))]...)
        if implementstrait($(esc(x)),$(esc(trait)))
            @warn @sprintf("%s already has trait %s, ignoring",$(esc(x)),$(esc(trait)))
        elseif implementstrait($(esc(x)),$(esc(supertype))($(esc(trait))))
            @warn @sprintf("%s already has parent trait %s, ignoring",$(esc(x)),$(esc(supertype))($(esc(trait))).body.name)
        else
            $(esc(:TraitDispatch)).traitdispatch(::Type{$(esc(trait)){$(esc(:T))}}) where $(esc(:T)) <: $(esc(x)) = $(esc(trait)){$(esc(:T))}
            super = $(esc(supertype))($(esc(trait)))
            if super != AbstractTrait
                $(esc(:TraitDispatch)).traitdispatch(::Type{super{$(esc(:T))}}) where $(esc(:T)) <: $(esc(x)) = $(esc(trait)){$(esc(:T))}
            end
        end
        nothing
    end
    prettify(code)
end


text:
#using SurfaceGeometry
using Optim

using Parameters

@with_kw struct Erdmanis2016
    C::AbstractFloat = 0.4
    ftol::AbstractFloat = 1e-6
end

function F(points,faces,v,zc::Erdmanis2016)

    Cp = zc.C

    s = 0
    for ti in 1:size(faces,2)
        v1,v2,v3 = faces[:,ti]

        a = norm(points[:,v2] - points[:,v1])
        b = norm(points[:,v3] - points[:,v2])
        c = norm(points[:,v1] - points[:,v3])

        xava = dot(points[:,v2]-points[:,v1],v[:,v2]-v[:,v1])
        xbvb = dot(points[:,v3]-points[:,v2],v[:,v3]-v[:,v2])
        xcvc = dot(points[:,v1]-points[:,v3],v[:,v1]-v[:,v3])

        ### This part is responsible for first part in the sum
        s += v2>v1 ? 4 * xava^2 : 0
        s += v3>v2 ? 4 * xbvb^2 : 0
        s += v1>v3 ? 4 * xcvc^2 : 0

        ### This part is needed for ordinary things
        Cdelta = 1/4*sqrt(1 - 2* (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        DCdelta = -A*xava - B*xbvb - C*xcvc

        s += Cp*DCdelta^2/Cdelta^2
    end

    return s
end

function gradF!(points,faces,v,storage,zc::Erdmanis2016)

    Cp = zc.C

    storage[:,:] .= 0

    for i in 1:size(faces,2)
        v1,v2,v3 = faces[:,i]
        a = norm(points[:,v2] - points[:,v1])
        b = norm(points[:,v3] - points[:,v2])
        c = norm(points[:,v1] - points[:,v3])

        Cdelta = 1/4*sqrt(1 - 2 * (a^4 + b^4 + c^4)/(a^2 + b^2 + c^2)^2)
        A = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( a^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        B = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( b^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)
        C = 1/4/Cdelta/(a^2 + b^2 + c^2)^3 * ( c^2*(a^2 + b^2 + c^2) - a^4 - b^4 - c^4)

        xa = points[:,v2] - points[:,v1]
        xb = points[:,v3] - points[:,v2]
        xc = points[:,v1] - points[:,v3]

        xava = dot(points[:,v2]-points[:,v1],v[:,v2]-v[:,v1])
        xbvb = dot(points[:,v3]-points[:,v2],v[:,v3]-v[:,v2])
        xcvc = dot(points[:,v1]-points[:,v3],v[:,v1]-v[:,v3])

        DCdelta = -A*xava - B*xbvb - C*xcvc

        A_ = 2*Cp*DCdelta*A/Cdelta^2
        B_ = 2*Cp*DCdelta*B/Cdelta^2
        C_ = 2*Cp*DCdelta*C/Cdelta^2

        storage[:,v1] += (A_ - 8*xava)*xa - C_*xc
        storage[:,v2] += (B_ - 8*xbvb)*xb - A_*xa
        storage[:,v3] += (C_ - 8*xcvc)*xc - B_*xb
    end
end

function stabilise!(points,faces,n,v,zc::Erdmanis2016; op=:optim)
    stabiliseV3Optim!(points,faces,n,v,zc)
end

function stabiliseV3Optim!(points,faces,n,v,zc::Erdmanis2016)
    #h2 = hvec(points,faces,n)

    function f(x::Vector)
        vv = reshape(x,3,div(length(x),3))
        return F(points,faces,vv,zc)
    end

    function g!(storage::Vector,x::Vector)

        vv = reshape(x,size(points)...)
        gradf = reshape(storage,size(points)...)

        gradF!(points,faces,vv,gradf,zc)

        for i in 1:size(v,2)
            P = eye(3)-n[:,i]*n[:,i]'
            gradf[:,i] = P*gradf[:,i]
        end
    end

    res = optimize(f,g!,v[:],method=Optim.ConjugateGradient(),ftol=zc.ftol)
    v[:,:] = reshape(res.minimum, size(v)...)[:,:]
end

function stabiliseV3!(points,faces,n,v,zc::Erdmanis2016)

    gradF_v = Array{Float64}(undef,size(points)...)
    gradF!(points,faces,v,gradF_v,zc)

    f = Array{Float64}(undef,size(points)...)
    for i in 1:size(points,2)
        f[:,i] = - (eye(3) - n[:,i]*n[:,i]') * gradF_v[:,i]
    end

    Fs = Inf
    eps = 1e-4

    gradF_f = Array{Float64}(undef,size(points)...)
    gradF!(points,faces,f,gradF_f,zc)

    q = 2*dot(gradF_f[:],f[:])
    s = dot(gradF_v[:],f[:]) + dot(gradF_f[:],v[:])
    ksi = - s/q

    gradF_dv = Array{Float64}(undef,size(points)...)
    vp = copy(v)
    v[:,:] = v[:,:] + ksi*f[:,:]

    # step = 1
    for step in 1:150

        gradF!(points,faces,v,gradF_v,zc)
        for i in 1:size(points,2)
            f[:,i] = -(eye(3) - n[:,i]*n[:,i]') * gradF_v[:,i]
        end

        gradF!(points,faces,f,gradF_f,zc)
        gradF!(points,faces,v-vp,gradF_dv,zc)

        S = [dot(gradF_v[:],f[:]) + dot(gradF_f[:],v[:]), dot(gradF_dv[:],v[:]) + dot(gradF_v[:],v[:] - vp[:])]
        Q = [2*dot(gradF_f[:],f[:])   dot(gradF_f[:],v[:]-vp[:])+dot(gradF_dv[:],f[:]); dot(gradF_dv[:],f[:])+dot(gradF_v[:],v[:]-vp[:])   2*dot(gradF_dv[:],v[:]-vp[:])]
        ksi, eta = -Q\S

        v[:,:],vp[:,:] = v[:,:] + ksi*f[:,:]  + eta*(v[:,:] - vp[:,:]), v[:,:]

        Fsnew = F(points,faces,v,zc)
        if (abs(Fs - Fsnew) < eps*Fsnew )
            break
        else
            Fs = Fsnew
        end
    end

    # if step>100
    #     warn("stabilistation used $step steps")
    # end
end


text:
using Test, StaticArrays, Colors

# push!(LOAD_PATH, "./src", "src")
using Solids
import Solids: _FIXED, Vec, Path, Point
import Solids: from_clipper, to_clipper

function test_from_to(T, x)
	return @test from_clipper(T, to_clipper(T, x)) == x
end

@testset "Types" begin #<<<1
@testset "Basic types" begin #<<<2
V = Point(1,2)
@test V.coords === SA[1,2]
end
@testset "Conversion to/from Clipper.jl" begin #<<<2
V = Point(1,2)
test_from_to(Int, 3)
test_from_to(Int, V)
test_from_to(Int, [V, V])
test_from_to(_FIXED, 3)
test_from_to(_FIXED, Point{2,_FIXED}(1,2))
test_from_to(Float64, 3)
end
end

@testset "Handling of objects" begin #<<<1
using Solids: Square, Circle
using Solids: children
using Solids: mult_matrix, translate, scale
using Solids: color
s = Square(1)
@testset "Primitives" begin #<<<2
@test s == Square([1,1])
end
@testset "Operations" begin #<<<2
@test union(s, union()) === s
@test length(children(union(s, union(s)))) == 2
end
@testset "Transforms" begin #<<<2
@test 2s == scale(2, s)
@test scale(2)*s == scale(2, s)
@test scale(2)*[s] == scale(2, s)
@test color("red", s) == color(parse(Colorant, "red"), s)
end
end
@testset "Clipper" begin #<<<1
s = Square(1)
# FIXME
end
@testset "Extrusion" begin #<<<1
using Solids: path_extrude, points
C = points(Circle(3.),(precision=.01,accuracy=1))
c = [Point(20*cos(i),20*sin(i)) for i in 0:.1:π]; c=[c;[Point(0.,-1.)]]
@test (path_extrude(c, C)) != 0
end
@testset "Convex hull" begin #<<<1
using Solids: convex_hull, convex_hull_list
P(x...) = Point(Float64.(x)...)
CH = convex_hull([P(0,0,0),P(0,0,10),P(10,0,0),P(0,10,0),P(1,1,1),P(1,0,0),])
@test Set(CH[1]) == Set([P(0,0,0),P(0,0,10),P(10,0,0),P(0,10,0),])
@test length(CH[2]) ==4
@test convex_hull_list([ P(0., 0), P(1., 1), P(-1., 1), P(.2, 2.), P(0, .8), P(-.2, .8), ]) == [1,2,4,3]
# @test convex_hull_list(Solids.rows(SA[0. 0;1. 1;-1. 1;.2 2.;0 .8;-.2 .8])) == [1,2,4,3]
@test convex_hull_list([
	P(-2.627798062316817, 1.075268817204301),
	P(-0.5030257403564974, -1.720430107526882),
	P(0.7927283156659947, 2.7956989247311825),
	P(0.0, 2.396978520135108),
	P(0.0, 0.03278003249806397),
	]) == [2,5,3,1]
@test convex_hull_list([
	P(-2.150537634408602, 1.3494700327417308),
	P(-0.4301075268817205, -2.097580910437773),
	P(2.3655913978494625, 0.04739817471346019),
	P(0.0, 0.0),
	P(2.3038140933536018, 0.0),
	P(0.0, 0.7294358146330306),
	])== [2,5,3,6,1]
end
@testset "Surfaces" begin #<<<1
using Solids: connected_components, Surface, merge, select_faces
using Solids: nvertices, nfaces
v=[[0,-1],[1,0],[0,1],[-1,0]]
for j in eachindex(v), i in 1:j-1
	@test Solids.circular_lt(v[i], v[j])
end
@test connected_components([:a, :b, :c, :d, :e], [[1,2],[1,3],[4,5]]) ==
	[([:a, :b, :c], [[1,2],[1,3]]),
	 ([:d, :e], [[1,2]]) ]
function pyramid(t=[0,0,0], n=0)
	points = Solids.Point{3}.([ t, t+[2.,0,0], t+[2,2,0], t+[0,2,0], t+[1,1,1]])
	faces = [[4,3,2],[2,3,1],[1,2,5],[2,3,5],[3,4,5],[4,1,5]]
	faces1 = [ f .+ n for f in faces ]
	return Surface(points, faces1)
end
p1 = pyramid()
p2 = pyramid([1,0,0])
p3 = merge(p1, p2)
p1bis = select_faces(1:6, p3)
@test nfaces(p1bis) == nfaces(p1)
@test nvertices(p1bis) == nvertices(p1)
u12 = union(p1, p2)
i12 = intersect(p1, p2)
@test length(u12.points) == 14 && length(u12.faces) == 24
@test length(i12.points) == 8  && length(u12.faces) == 12
end
#>>>1

# vim: noet ts=2 fmr=<<<,>>>


text:
# s-pair, a pair of polynomials
struct SPair
    # first generator as index from the basis array
    poly1::Int
    # second generator -//-
    poly2::Int
    # position of lcm(poly1, poly2) in hashtable
    lcm::Int
    # total degree of lcm
    deg::UInt
end

mutable struct Pairset
    pairs::Vector{SPair}
    # number of filled pairs,
    # Initially zero
    load::Int
end

function initialize_pairset(; initial_size=2^6) # TODO: why 64?
    pairs = Vector{SPair}(undef, initial_size)
    return Pairset(pairs, 0)
end

function Base.isempty(ps::Pairset)
    return ps.load == 0
end

function check_enlarge_pairset!(ps::Pairset, added::Int)
    sz = length(ps.pairs)
    if ps.load + added >= sz
        newsz = max(2 * sz, ps.load + added)
        resize!(ps.pairs, newsz)
    end
end

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

# TODO: get rid of `Tv` dependency
mutable struct Basis
    # vector of polynomials, each polynomial is a vector of monomials,
    # each monomial is represented with it's position in hashtable
    gens::Vector{Vector{Int}}
    # polynomial coefficients
    coeffs::Vector{Vector{UInt64}}

    #= Keeping track of sizes   =#
    #=  ndone <= ntotal <= size =#
    # total allocated size,
    # size == length(gens) is always true
    size::Int
    # number of processed polynomials in `gens`
    # Iitially a zero
    ndone::Int
    # total number of polys filled in `gens`
    # (these will be handled during next update! call)
    # Iitially a zero
    ntotal::Int

    #= Keeping track of redundancy =#
    #= invariant =#
    #= length(lead) == length(nonred) == count(isred) == nlead =#
    # if element of the basis
    # is redundant
    isred::Vector{Int8}
    # positions of non-redundant elements in the basis
    nonred::Vector{Int}
    # division masks of leading monomials of
    # non redundant basis elements
    lead::Vector{UInt32}
    # number of filled elements in lead
    nlead::Int

    # characteristic of ground field
    ch::UInt64
end

function initialize_basis(ring::PolyRing, ngens::Int)
    #=
        always true
        length(gens) == length(coeffs) == length(isred) == size
    =#

    sz     = ngens * 2
    ndone  = 0
    ntotal = 0
    nlead  = 0

    gens   = Vector{Vector{Int}}(undef, sz)
    coeffs = Vector{Vector{UInt64}}(undef, sz)
    isred  = zeros(Int8, sz)
    nonred = Vector{Int}(undef, sz)
    lead   = Vector{UInt32}(undef, sz)

    ch = ring.ch

    Basis(gens, coeffs, sz, ndone, ntotal, isred, nonred, lead, nlead, ch)
end

function check_enlarge_basis!(basis::Basis, added::Int)
    if basis.ndone + added >= basis.size
        basis.size = max(basis.size * 2, basis.ndone + added)
        resize!(basis.gens, basis.size)
        resize!(basis.coeffs, basis.size)
        resize!(basis.isred, basis.size)
        basis.isred[basis.ndone+1:end] .= 0
        resize!(basis.nonred, basis.size)
        resize!(basis.lead, basis.size)
    end
end

# Normalize each element of the input basis
# by dividing it by leading coefficient
function normalize_basis!(basis::Basis)
    cfs = basis.coeffs
    for i in 1:basis.ntotal
        # mul = inv(cfs[i][1])
        mul = uinvmod(cfs[i][1], basis.ch)
        for j in 2:length(cfs[i])
            # cfs[i][j] *= mul
            cfs[i][j] = umultmod(cfs[i][j], mul, basis.ch)
        end
        cfs[i][1] = one(cfs[i][1])
    end
    basis
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

@testset "demand file i/o" begin
	demand = readDemandFile("data/demand/demand_1.csv")

	# check that demand can be written to file and read again
	tempDemandFilename = "temp/demand_1.csv"
	writeDemandFile(tempDemandFilename, demand)
	tempDemand = readDemandFile(tempDemandFilename)
	rm(tempDemandFilename)

	# check that input and output files have same raster filenames
	@test abspath.(demand.rasterFilenames) == abspath.(tempDemand.rasterFilenames)

	# check that demand files were read correctly
	for d in [demand, tempDemand]
		# check demand
		@test d.numRasters == 4
		@test d.numModes == 6
		@test d.numSets == 2

		# check demand modes
		@test [mode.index for mode in d.modes] == [1,2,3,4,5,6]
		@test [mode.rasterIndex for mode in d.modes] == [1,3,3,2,4,4]
		@test [Int(mode.priority) for mode in d.modes] == [1,2,3,1,2,3]
		@test [mode.arrivalRate for mode in d.modes] == [25,15,5,30,20,10]

		# check demand sets
		@test d.modeLookup == [1 2 3; 4 5 6]
		@test demand.setsStartTimes == sort(vcat(0,[0.7:1:10;],[0.3:1:10;]))
		@test demand.setsTimeOrder == [2-isodd(i) for i = 1:length(demand.setsTimeOrder)]
	end
end

@testset "demand mode lookup" begin
	demand = readDemandFile("data/demand/demand_1.csv")

	for t = 0:0.1:1, priority in JEMSS.priorities
		i = JEMSS.getDemandMode!(demand, priority, t).index
		if 0.3 <= mod(t,1) < 0.7
			@test (priority == highPriority && i == 4) ||
				(priority == medPriority && i == 5) ||
				(priority == lowPriority && i == 6)
		else
			@test (priority == highPriority && i == 1) ||
				(priority == medPriority && i == 2) ||
				(priority == lowPriority && i == 3)
		end
	end

	# change in time beyond the last set start time should not cause demand set to change
	for t = demand.setsStartTimes[end] .+ (0:0.1:1)
		JEMSS.updateDemandToTime!(demand, t)
		@test demand.recentSetsStartTimesIndex == length(demand.setsStartTimes)
	end
end


text:
export heaviside
function heaviside(x,c = 0.15,delta = 0.1,epsilon = 0.1*delta)
a = c - delta;
b = c + delta;

y = zeros(eltype(x),size(x));
dy = zeros(eltype(x),size(x));
I = (x .>= (a-epsilon)) .& (x .<= (a+epsilon));
xI = x[I];
height = 1.0/(b-a);
t = (xI .- a .+ epsilon).*((0.25*height)/epsilon)

dy[I] = 2.0*t;
y[I] = (xI .- a .+ epsilon).*t;

I = (x .> a + epsilon) .& (x .<= b - epsilon);
y[I] .= height*(x[I].-a);
dy[I] .= height;
I = (x .> b+epsilon);
y[I] .= 1.0;
dy[I] .= 0.0;
I = (x .> b-epsilon) .& (x .<= b+epsilon);
xI = x[I];
t = (xI.-(b+epsilon)).*(-(0.25*height)/epsilon);
y[I] = (xI.-(b+epsilon)).*t .+ 1.0;
dy[I] = 2.0.*t;
return y,Diagonal(dy[:]);
end

export heaviside!
function heaviside!(x::Union{Array{Float64},Array{Float32}},y::Union{Array{Float64},Array{Float32}},c::Float64 = 0.15,delta::Float64 = 0.1,epsilon::Float64 = 0.1*delta)
a = c - delta;
b = c + delta;
height = 1.0/(b-a);
temp = ((0.25*height)/epsilon);
tmp2 = 0.0;
# y = zeros(size(x));
y[:] .= 0.0;
t1 = a-epsilon;
t2 = a+epsilon;
t3 = b-epsilon;
t4 = b+epsilon;
for k=1:length(x)
	xk = x[k];
	if xk <= t1
		x[k] = 0.0;
	elseif xk >= t4
		x[k] = 1.0;
	elseif xk >= t2 && xk <=t3
		x[k] = height*(xk-a);
		y[k] = height;
	elseif xk >= t1 && xk <=t2
		tmp2 = (xk - t1)*temp;
		x[k] = (xk - t1);
		x[k] *= tmp2;
		y[k] = 2.0;
		y[k] *= tmp2;
	elseif xk >= t3 && xk <=t4
		tmp2 = (t4 - xk)*temp;
		x[k] = (xk-t4).*tmp2 + 1.0;
		y[k] = 2.0*tmp2;
	end
end
return y;
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMIExport

using FMICore: FMU2, FMU2Component, fmi2ModelDescription, fmi2ValueReference, fmi2Component, fmi2ComponentEnvironment, fmi2Status, fmi2EventInfo
using FMICore: fmi2Causality, fmi2CausalityOutput, fmi2CausalityInput
using FMICore: fmi2ScalarVariable, fmi2Variability, fmi2Initial
using FMICore: fmi2ModelDescriptionReal, fmi2ModelDescriptionInteger, fmi2ModelDescriptionBoolean, fmi2ModelDescriptionString, fmi2ModelDescriptionEnumeration
using FMICore: fmi2VariableDependency, fmi2Unknown, fmi2DependencyKind
using FMICore: fmi2VariableNamingConventionStructured
using FMICore: fmi2CausalityToString, fmi2VariabilityToString, fmi2InitialToString
using FMICore: fmi2String, fmi2Boolean, fmi2Integer, fmi2Real, fmi2Enum
using FMICore: fmi2Type, fmi2TypeModelExchange, fmi2TypeCoSimulation
using FMICore: fmi2CallbackFunctions
using FMICore: fmi2StatusOK, fmi2StatusWarning, fmi2StatusError
using FMICore: fmi2True, fmi2False
using FMICore: fmi2ModelDescriptionModelExchange, fmi2ModelDescriptionCoSimulation

include("FMI2_md.jl")
export fmi2CreateModelDescription
export fmi2ModelDescriptionAddEvent
export fmi2GetIndexOfScalarVariable
export fmi2ModelDescriptionAddRealStateAndDerivative
export fmi2ModelDescriptionAddRealState, fmi2ModelDescriptionAddRealDerivative, fmi2ModelDescriptionAddRealInput, fmi2ModelDescriptionAddRealOutput, fmi2ModelDescriptionAddRealParameter, fmi2ModelDescriptionAddEventIndicator
export fmi2ModelDescriptionAddModelVariable, fmi2ModelDescriptionAddModelStructureOutputs, fmi2ModelDescriptionAddModelStructureDerivatives, fmi2ModelDescriptionAddModelStructureInitialUnknowns

export fmi2Create
export fmi2AddRealStateAndDerivative, fmi2AddStateAndDerivative
export fmi2AddRealOutput, fmi2AddOutput
export fmi2AddRealInput, fmi2AddInput
export fmi2AddRealParameter, fmi2AddParameter, fmi2AddEventIndicator

export fmi2SetFctGetTypesPlatform, fmi2SetFctGetVersion
export fmi2SetFctInstantiate, fmi2SetFctFreeInstance, fmi2SetFctSetDebugLogging, fmi2SetFctSetupExperiment, fmi2SetEnterInitializationMode, fmi2SetFctExitInitializationMode
export fmi2SetFctTerminate, fmi2SetFctReset
export fmi2SetFctGetReal, fmi2SetFctGetInteger, fmi2SetFctGetBoolean, fmi2SetFctGetString, fmi2SetFctSetReal, fmi2SetFctSetInteger, fmi2SetFctSetBoolean, fmi2SetFctSetString
export fmi2SetFctSetTime, fmi2SetFctSetContinuousStates, fmi2SetFctEnterEventMode, fmi2SetFctNewDiscreteStates, fmi2SetFctEnterContinuousTimeMode, fmi2SetFctCompletedIntegratorStep
export fmi2SetFctGetDerivatives, fmi2SetFctGetEventIndicators, fmi2SetFctGetContinuousStates, fmi2SetFctGetNominalsOfContinuousStates

include("FMI2_simple.jl")
export fmi2CreateSimple

# export fmi2ComponentStruct

function fmi2SetFctGetTypesPlatform(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2String, ())
    fmu.cGetTypesPlatform = c_fun.ptr
end

function fmi2SetFctGetVersion(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2String, ())
    fmu.cGetVersion = c_fun.ptr
end

function fmi2SetFctInstantiate(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Component, (fmi2String, fmi2Type, fmi2String, fmi2String, Ptr{fmi2CallbackFunctions}, fmi2Boolean, fmi2Boolean))
    fmu.cInstantiate = c_fun.ptr
end

function fmi2SetFctFreeInstance(fmu::FMU2, fun)
    c_fun = @cfunction($fun, Cvoid, (fmi2Component,))
    fmu.cFreeInstance = c_fun.ptr
end

function fmi2SetFctSetDebugLogging(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, fmi2Boolean, Csize_t, Ptr{fmi2String}))
    fmu.cSetDebugLogging = c_fun.ptr
end

function fmi2SetFctSetupExperiment(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, fmi2Boolean, fmi2Real, fmi2Real, fmi2Boolean, fmi2Real))
    fmu.cSetupExperiment = c_fun.ptr
end

function fmi2SetEnterInitializationMode(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cEnterInitializationMode = c_fun.ptr
end

function fmi2SetFctExitInitializationMode(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cExitInitializationMode= c_fun.ptr
end

function fmi2SetFctTerminate(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cTerminate = c_fun.ptr
end

function fmi2SetFctReset(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cReset = c_fun.ptr
end

function fmi2SetFctGetReal(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Real}))
    fmu.cGetReal = c_fun.ptr
end

function fmi2SetFctGetInteger(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Integer}))
    fmu.cGetInteger = c_fun.ptr
end

function fmi2SetFctGetBoolean(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Boolean}))
    fmu.cGetBoolean = c_fun.ptr
end

function fmi2SetFctGetString(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2String}))
    fmu.cGetString = c_fun.ptr
end

function fmi2SetFctSetReal(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Real}))
    fmu.cSetReal = c_fun.ptr
end

function fmi2SetFctSetInteger(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Integer}))
    fmu.cSetInteger = c_fun.ptr
end

function fmi2SetFctSetBoolean(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Boolean}))
    fmu.cSetBoolean = c_fun.ptr
end

function fmi2SetFctSetString(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2String}))
    fmu.cSetString = c_fun.ptr
end

function fmi2SetFctSetTime(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, fmi2Real))
    fmu.cSetTime = c_fun.ptr
end

function fmi2SetFctSetContinuousStates(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2Real}, Csize_t))
    fmu.cSetContinuousStates = c_fun.ptr
end

function fmi2SetFctEnterEventMode(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cEnterEventMode = c_fun.ptr
end

function fmi2SetFctNewDiscreteStates(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2EventInfo}))
    fmu.cNewDiscreteStates = c_fun.ptr
end

function fmi2SetFctEnterContinuousTimeMode(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component,))
    fmu.cEnterContinuousTimeMode = c_fun.ptr
end

function fmi2SetFctCompletedIntegratorStep(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, fmi2Boolean, Ptr{fmi2Boolean}, Ptr{fmi2Boolean}))
    fmu.cCompletedIntegratorStep = c_fun.ptr
end

function fmi2SetFctGetDerivatives(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2Real}, Csize_t))
    fmu.cGetDerivatives = c_fun.ptr
end

function fmi2SetFctGetEventIndicators(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2Real}, Csize_t))
    fmu.cGetEventIndicators = c_fun.ptr
end

function fmi2SetFctGetContinuousStates(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2Real}, Csize_t))
    fmu.cGetContinuousStates = c_fun.ptr
end

function fmi2SetFctGetNominalsOfContinuousStates(fmu::FMU2, fun)
    c_fun = @cfunction($fun, fmi2Status, (fmi2Component, Ptr{fmi2Real}, Csize_t))
    fmu.cGetNominalsOfContinuousStates = c_fun.ptr
end

"""
ToDo
"""
function fmi2Create(instanceName::String="", type=fmi2TypeModelExchange)
    fmu = FMU2()

    fmu.instanceName = instanceName
    fmu.type = type
    fmu.modelDescription = fmi2CreateModelDescription()
    fmu.fmuResourceLocation = pwd()
    #fmu.callbackFunctions
    #fmu.visible = false
    #fmu.loggingOn = false

    return fmu
end

function fmi2AddRealStateAndDerivative(fmu, stateName; stateDescr, derivativeDescr)
    fmi2ModelDescriptionAddRealStateAndDerivative(fmu.modelDescription, stateName; stateDescr=stateDescr, derivativeDescr=derivativeDescr)
end
fmi2AddStateAndDerivative = fmi2AddRealStateAndDerivative

function fmi2AddRealOutput(fmu, name; description)
    fmi2ModelDescriptionAddRealOutput(fmu.modelDescription, name; description=description)
end
fmi2AddOutput = fmi2AddRealOutput

function fmi2AddRealInput(fmu, name; description)
    fmi2ModelDescriptionAddRealInput(fmu.modelDescription, name; description=description)
end
fmi2AddInput = fmi2AddRealInput

function fmi2AddRealParameter(fmu, name; description)
    fmi2ModelDescriptionAddRealParameter(fmu.modelDescription, name; description=description)
end
fmi2AddParameter = fmi2AddRealParameter

function fmi2AddEventIndicator(fmu)
    fmi2ModelDescriptionAddEventIndicator(fmu.modelDescription)
end

"""
ToDo
"""
function fmi2Check(fmu::FMU2)
    return true
end

end # module


text:
import JSON2

function validate_serialization(sys::System; time_series_read_only = false)
    test_dir = mktempdir()
    orig_dir = pwd()
    cd(test_dir)

    try
        path = joinpath(test_dir, "test_system_serialization.json")
        io = open(path, "w")
        @info "Serializing to $path"
        sys_ext = get_ext(sys)
        sys_ext["data"] = 5
        ext_test_bus_name = ""
        try
            IS.prepare_for_serialization!(sys.data, path; force = true)
            bus = collect(get_components(PSY.Bus, sys))[1]
            ext_test_bus_name = PSY.get_name(bus)
            ext = PSY.get_ext(bus)
            ext["test_field"] = 1
            to_json(io, sys)
        finally
            close(io)
        end

        sys2 = System(path; time_series_read_only = time_series_read_only)
        sys_ext2 = get_ext(sys2)
        sys_ext2["data"] != 5 && return false
        bus = PSY.get_component(PSY.Bus, sys2, ext_test_bus_name)
        ext = PSY.get_ext(bus)
        ext["test_field"] != 1 && return false
        return sys2, IS.compare_values(sys, sys2)
    finally
        cd(orig_dir)
    end
end

@testset "Test JSON serialization of RTS data with mutable time series" begin
    sys = create_rts_system()

    # Add an AGC service to cover its special serialization.
    control_area = get_component(Area, sys, "1")
    AGC_service = PSY.AGC(
        name = "AGC_Area1",
        available = true,
        bias = 739.0,
        K_p = 2.5,
        K_i = 0.1,
        K_d = 0.0,
        delta_t = 4,
        area = control_area,
    )
    initial_time = Dates.DateTime("2020-01-01T00:00:00")
    end_time = Dates.DateTime("2020-01-01T23:00:00")
    dates = collect(initial_time:Dates.Hour(1):end_time)
    data = collect(1:24)
    label = "get_activepower"
    contributing_devices = Vector{Device}()
    for g in get_components(
        ThermalStandard,
        sys,
        x -> (x.primemover ∈ [PrimeMovers.ST, PrimeMovers.CC, PrimeMovers.CT]),
    )
        if get_area(get_bus(g)) != control_area
            continue
        end
        ta = TimeSeries.TimeArray(dates, data, [Symbol(get_name(g))])
        forecast = IS.Deterministic(label, ta)
        add_forecast!(sys, g, forecast)

        t = RegulationDevice(g, participation_factor = (up = 1.0, dn = 1.0), droop = 0.04)
        add_component!(sys, t)
        @test isnothing(get_component(ThermalStandard, sys, get_name(g)))
        push!(contributing_devices, t)
    end
    add_service!(sys, AGC_service, contributing_devices)

    sys2, result = validate_serialization(sys; time_series_read_only = false)
    @test result

    # Ensure the forecasts attached to the ThermalStandard got deserialized.
    for rd in get_components(RegulationDevice, sys2)
        @test get_forecast(Deterministic, rd, initial_time, label) isa Deterministic
    end

    clear_forecasts!(sys2)
end

@testset "Test JSON serialization of RTS data with immutable time series" begin
    sys = create_rts_system()
    sys2, result = validate_serialization(sys; time_series_read_only = true)
    @test result
    @test_throws ErrorException clear_forecasts!(sys2)
    # Full error checking is done in IS.
end

@testset "Test JSON serialization of matpower data" begin
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))

    # Add a Probabilistic forecast to get coverage serializing it.
    bus = Bus(nothing)
    bus.name = "Bus1234"
    add_component!(sys, bus)
    tg = RenewableFix(nothing)
    tg.bus = bus
    add_component!(sys, tg)
    tProbabilisticForecast =
        PSY.Probabilistic("scalingfactor", Hour(1), DateTime("01-01-01"), [0.5, 0.5], 24)
    add_forecast!(sys, tg, tProbabilisticForecast)

    _, result = validate_serialization(sys)
    @test result
end

@testset "Test JSON serialization of ACTIVSg2000 data" begin
    path = joinpath(DATA_DIR, "ACTIVSg2000", "ACTIVSg2000.m")
    sys = System(PowerSystems.PowerModelsData(path))
    _, result = validate_serialization(sys)
    @test result
end

@testset "Test JSON serialization of dynamic inverter" begin
    sys = create_system_with_dynamic_inverter()
    _, result = validate_serialization(sys)
    @test result
end

@testset "Test deepcopy of a system" begin
    sys = create_rts_system()
    sys2 = deepcopy(sys)
    clear_forecasts!(sys2)
    @test !isempty(collect(PSY.iterate_forecasts(sys)))
end


text:
push!(LOAD_PATH, "@stdlib")  # NOTE: this is needed to enable this test to run from the Pkg manager
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Test
using Pkg
import ParallelStencil: SUPPORTED_PACKAGES, PKG_CUDA
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "incremental compilation" begin
            Pkg.activate(joinpath(@__DIR__, "test_projects", "Diffusion3D_$(nameof($package))"))
            using $(Symbol("Diffusion3D_$package"))
            @test $(Symbol("Diffusion3D_$package")).diffusion3D()
        end;
    end;
)) end == nothing || true;