text:
module ML

using Random, Statistics
using LinearAlgebra, Flux, Zygote

import Base:
    length, reverse, iterate

export model, train!, validate, preprocess, update_dimension

# NOTE: all data is assumed to be dᵢ x N shaped here!
#       this is to allow for fast matrix multiplication through the network
#       it might be worthwhile to keep both memory layouts stored...

# ------------------------------------------------------------------------
# globals

σ₀(x)  = x + (sqrt(x^2 +1)-1)/2

# ------------------------------------------------------------------------
# types

# Iterator
"""
    struct LayerIterator
        width     :: Array{Int}
        dropout   :: Set{Int}
        normalize :: Set{Int}
        σᵢ        :: Function
        σₒ        :: Function
        σ         :: Function
    end

An iterator used to generate dense latent layers within a neural network.
`width` denotes the widths of each layer; the length of this array immediately determines the depth.
`dropout` denotes the layers, as given by `width` that are followed by a dropout layer.
`normalize` denotes the layers, as given by `width` that are followed by a batch normalization layer.
`σᵢ`, `σₒ`, `σ` is the activation energy on the first, last, and intermediate layers respectively.
"""
struct LayerIterator
    width     :: Array{Int}
    dropout   :: Set{Int}
    normalize :: Set{Int}
    σᵢ        :: Function # activation on input layers
    σₒ        :: Function # activation on output layers
    σ         :: Function # activation on latent layers
end

length(it::LayerIterator)  = length(it.width) + length(it.dropout) + length(it.normalize)
reverse(it::LayerIterator) = LayerIterator(
                                reverse(it.width),
                                Set(length(it.width) - i - 1 for i in it.dropout),
                                Set(length(it.width) - i - 1 for i in it.normalize),
                                it.σᵢ,
                                it.σᵢ, # intentional -> want to make output = input
                                it.σ,
                             )

function iterate(it::LayerIterator)
    w₁ = it.width[1]
    w₂ = it.width[2]
    f  = Dense(w₁, w₂, it.σᵢ)

    return f, (
        index     = 2,
        dropout   = 1 ∈ it.dropout,
        normalize = 1 ∈ it.normalize,
    )
end

function iterate(it::LayerIterator, state)
    return if state.dropout
               Dropout(0.5), (
                   index     = state.index,
                   dropout   = false,
                   normalize = state.normalize,
               )
           elseif state.normalize
               BatchNorm(it.width[state.index]), (
                   index     = state.index,
                   dropout   = false,
                   normalize = false,
               )
           elseif state.index < length(it.width)
                w₁ = it.width[state.index]
                w₂ = it.width[state.index+1]

                i  = state.index+1
                f  = Dense(w₁, w₂, i == length(it.width) ? it.σₒ : it.σ)

                f, (
                     index     = i,
                     dropout   = (i-1) ∈ it.dropout,
                     normalize = (i-1) ∈ it.normalize,
                )
           else
               nothing
           end
end

# ------------------------------------------------------------------------
# functions

"""
    model(dᵢ, dₒ; Ws=Int[], normalizes=Int[], dropouts=Int[], σ=elu)

Initialize an autoencoding neural network with input dimension `dᵢ` and latent layers `dₒ`.
`Ws` specifies both the width and depth of the encoder layer - the width of each layer is given as an entry in the array while the length specifies the depth.
`normalizes` and `dropouts` denote which layers are followed by batch normalization and dropout specifically.
The decoder layer is given the mirror symmetric architecture.
"""
function model(dᵢ, dₒ; Ws=Int[], normalizes=Int[], dropouts=Int[], σ=elu)
    # check for obvious errors here
    length(dropouts)   > 0 && length(Ws) < maximum(dropouts)   ≤ 0 && error("invalid dropout layer position")
    length(normalizes) > 0 && length(Ws) < maximum(normalizes) ≤ 0 && error("invalid normalization layer position")

    layers = LayerIterator(
                [dᵢ; Ws; dₒ],
                Set(dropouts),
                Set(normalizes),
                σ₀, σ₀, σ
             )

    F   = Chain(layers...)
    F¯¹ = Chain(reverse(layers)...)
    𝕀   = Chain(F, F¯¹)

    return (
        pullback=F,
        pushforward=F¯¹,
        identity=𝕀
    )
end

"""
    update_dimension(model, dₒ; ϵ = 1e-6)

Add a colection of new neurons in the encoding layer to encode in the encoding layer to increase dimensions to `dₒ`.
Model weights for the initial dimensions are kept the same.
"""
function update_dimension(model, dₒ; ϵ = 1e-6)
    F, F¯¹, 𝕀 = model

    lᵢ = F[end]
    lₒ = F¯¹[1]

    Wᵢ, bᵢ = params(lᵢ)
    Wₒ, bₒ = params(lₒ)

    size(Wᵢ,1) == dₒ && return nothing
    size(Wᵢ,1) >  dₒ && error("can not reduce dimensionality of model")

    δ  = dₒ - size(Wᵢ, 1)

    W̄ᵢ = vcat(Wᵢ, ϵ*randn(δ, size(Wᵢ,2)))
    b̄ᵢ = vcat(bᵢ, ϵ*randn(δ))
    l̄ᵢ = Dense(W̄ᵢ, b̄ᵢ, lᵢ.σ)

    W̄ₒ = hcat(Wₒ, ϵ*randn(size(Wₒ,1), δ))
    b̄ₒ = bₒ
    l̄ₒ = Dense(W̄ₒ, b̄ₒ, lₒ.σ)

    F   = Chain( (i < length(F) ? f : l̄ᵢ for (i,f) ∈ enumerate(F))...)
    F¯¹ = Chain( (i > 1 ? f : l̄ₒ for (i,f) ∈ enumerate(F¯¹))...)
    𝕀   = Chain(F, F¯¹)

    return (
        pullback=F,
        pushforward=F¯¹,
        identity=𝕀
    )
end

# data batching
"""
    batch(data, n)

Randomly partition `data` into groups of size `n`.
"""
function batch(data, n)
    N = size(data,2)

    lo(i) = (i-1)*n + 1
    hi(i) = min((i)*n, N)

    ι = randperm(N)

    return (data[:,ι[lo(i):hi(i)]] for i in 1:ceil(Int, N/n)),
           (ι[lo(i):hi(i)] for i in 1:ceil(Int, N/n))
end

"""
    validate(data, len)

Reserve `len` samples from `data` during training process to allow for model validation.
"""
function validate(data, len)
    ι = randperm(size(data,2))
    return (
        valid = data[:,ι[1:len]],
        train = data[:,ι[len+1:end]],
    ),(
        valid = ι[1:len],
        train = ι[len+1:end],
    )
end

function noop(epoch) end

# data training
"""
    train!(model, data, index, loss; B=64, η=1e-3, N=100, log=noop)

Trains autoencoder `model` on `data` by minimizing `loss`.
`index` stores the underlying indices of data used for training.
Will mutate the underlying parameters of `model`.
Optional parameters include:
  1. `B` denotes the batch size to be used.
  2. `N` denotes the number of epochs.
  3. `η` denotes the learning rate.
"""
function train!(model, data, index, loss; B=64, η=1e-3, N=100, log=noop)
    Θ   = params(model.identity)
    opt = ADAM(η)

    for n ∈ 1:N
        X, I = batch(data, B)
        for (i,x) ∈ zip(I,X)
            E, backpropagate = pullback(Θ) do
                loss(x, index[i], false)
            end

            isnan(E) && @goto done

            ∇Θ = backpropagate(1f0)
            Flux.Optimise.update!(opt, Θ, ∇Θ)
        end

        log(n)
    end
    @label done
end

# ------------------------------------------------------------------------
# tests

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
"""
    eigenvector_centrality(g)

Compute the eigenvector centrality for the graph `g`.

Eigenvector centrality computes the centrality for a node based on the
centrality of its neighbors. The eigenvector centrality for node `i` is
the \$i^{th}\$ element of \$\\mathbf{x}\$ in the equation
``
    \\mathbf{Ax} = λ \\mathbf{x}
``
where \$\\mathbf{A}\$ is the adjacency matrix of the graph `g`
with eigenvalue λ.

By virtue of the Perron–Frobenius theorem, there is a unique and positive
solution if λ is the largest eigenvalue associated with the
eigenvector of the adjacency matrix \$\\mathbf{A}\$.

### References

- Phillip Bonacich: Power and Centrality: A Family of Measures.
    American Journal of Sociology 92(5):1170–1182, 1986
    http://www.leonidzhukov.net/hse/2014/socialnetworks/papers/Bonacich-Centrality.pdf
- Mark E. J. Newman: Networks: An Introduction.
       Oxford University Press, USA, 2010, pp. 169.
"""
eigenvector_centrality(g::AbstractGraph) = abs.(vec(eigs(adjacency_matrix(g), nev=1)[2]))::Vector{Float64}

ext: jl


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

#-

using MLJ

# ## Getting some data:

using RDatasets
iris = dataset("datasets", "iris"); # a DataFrame
scrambled = shuffle(1:size(iris, 1))
X = iris[scrambled, 1:4];
y = iris[scrambled, 5];

first(X, 4)

#-

y[1:5]

# ## Basic fit and predict:

@load SVC()
classifier_ = SVC()
classifier = machine(classifier_, X, y)
fit!(classifier)
ŷ = predict(classifier, X) # or some Xnew

#-

# ## Evaluating the model:
evaluate!(classifier,
          resampling=Holdout(fraction_train=0.8),
          measure=misclassification_rate)

# ## Adding dimension reduction:

@load PCA
dim_reducer_ = PCA()
dim_reducer = machine(dim_reducer_, X)
fit!(dim_reducer)
Xsmall = transform(dim_reducer, X);

first(Xsmall, 3)

#-

classifier = machine(classifier_, Xsmall, y)
fit!(classifier)
ŷ = predict(classifier, Xsmall)

# ## Building a composite model:

# ### Method 1: Compact syntax (but not generalizable):

# (not implemented at time of talk)

## composite_ = @pipeline dim_reducer_ classifier_

## composite = machine(composite_, X, y)
## evaluate!(composite, measure=misclassification_rate)

# ### Method 2: Re-interpret unstreamlined code:

Xraw = X;
yraw = y;

X = source(Xraw)
y = source(yraw)

dim_reducer = machine(dim_reducer_, X)
Xsmall = transform(dim_reducer, X)

classifier = machine(classifier_, Xsmall, y)
ŷ = predict(classifier, Xsmall)

#-

fit!(ŷ)

#-

ŷ(rows=3:4)

#-

dim_reducer_.ncomp = 1  # maximum output dimension
fit!(ŷ)

#-

ŷ(rows=3:4)

#  Changing classifier hyperparameter does not retrigger retraining of
#  upstream dimension reducer:

classifier_.gamma = 0.1
fit!(ŷ)

#-

ŷ(rows=3:4)

# Predicting on new data (`Xraw` in `source(Xraw)` is substituted for `Xnew`):

Xnew = (SepalLength = [4.0, 5.2],
        SepalWidth = [3.2, 3.0],
        PetalLength = [1.2, 1.5],
        PetalWidth = [0.1, 0.4],)
ŷ(Xnew)

# #### Exporting network as stand-alone reusable model:

composite_ = @from_network Composite(pca=dim_reducer_, svc=classifier_) <= (X, y, ŷ)
params(composite_)

#-

composite = machine(composite_, Xraw, yraw)
evaluate!(composite, measure=misclassification_rate)

# ## Evaluating a "self-tuning" random forest (nested resampling):

task = load_boston()
models(task)

#-

# ### Evaluating a single tree:

@load DecisionTreeRegressor # load code

tree_ = DecisionTreeRegressor(n_subfeatures=3)
tree = machine(tree_, task)
evaluate!(tree,
          resampling=Holdout(fraction_train=0.7),
          measure=[rms, mav])

# ### Use ensembling wrapper to create a random forest:

forest_ = EnsembleModel(atom=tree_, n=10)

# ### Wrapping in a tuning strategy creates a "self_tuning" random forest:

r1 = range(forest_, :bagging_fraction, lower=0.4, upper=1.0);
r2 = range(forest_, :(atom.n_subfeatures), lower=1, upper=12)

self_tuning_forest_ = TunedModel(model=forest_,
                          tuning=Grid(),
                          resampling=CV(),
                          ranges=[r1,r2],
                          measure=rms)

# ### Evaluate the self_tuning_forest (nested resampling):

self_tuning_forest = machine(self_tuning_forest_, task)

evaluate!(self_tuning_forest,
          resampling=CV(),
          measure=[rms,rmslp1])


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()

## ============================================================================
@info "load packages"
using CovidSurvey
using Random, Dates, Turing
using Base.Threads
using StatsPlots
using PrettyTables
using DataFrames
using BSON

using Memoization
using ReverseDiff
setadbackend(:reversediff)
Turing.setrdcache(true)

plotlyjs()
@info "number of threads available: $(nthreads())"
## ============================================================================

ps = (
    warmup = 100,
    model = "hospit",
    steps  = 100,
    seed   = nothing,
    observ = "2021-03-01",#"2021-02-06",#"2021-03-25"
    chains = 1,
    preds = "CS,MG,MA", #"CF,CC,CR,CS"
    semipara=true,
    rwstep = 7,
    include= false,
    std = false,
)
## ============================================================================
# load data
@info "load data"

data_params = (
      observationsend  = ps.observ
    , predictors        = ps.preds|> CovidSurvey.parse_predictors
    # , cases_start       = ps.cases
    # , hospitmodel       = Regional.HospitInit1(obs_stop="2020-07-01")
    , casemodel         = Regional.CaseInit(obs_start="2020-06-01")
    , seromodel         = Regional.SeroInit(delay=0, std=1.)
    , rwstep            = ps.rwstep
    , epidemicstart    = 20
    , numimpute        = 6
    , include_early_dynamic = ps.include
    , link              = KLogistic(4.5)
    , invlink           = KLogit(4.5)
    , predictive        = false
    , covariates_kwargs = Dict(
      :semiparametric => ps.semipara,
      :fname          => [projectdir("data/contacts/dk/", "averaged_contact_rates_region=$(r).csv" ) for r in 1:Regional.nregions],#[projectdir("data", "mean_contact_rates_region=$r.csv" ) for r in 1:Regional.nregions],
      :shift          => -1,
      :conditions     => :date => x->x>Date("2020-11-10"), # >= 1000 cases / day
      :datecol        => "dates",
      :aggregation    => ps.rwstep,
      :standardize    => ps.std,
      :normalize      => true,
      :mobility       => [projectdir("data/mobility/mobility_region=$(Regional.regions[i]).csv") for i in 1:Regional.nregions]
      )
    )

data = Regional.load_data(; data_params... )
d = turing_data = data.turing_data;

## ============================================================================
@info "prepare model"

model = Regional.name2model[ps.model]
# model = Regional.model_hospit2
#!isnothing(ps.seed) && Random.seed!(ps.seed);
@info ps
m = model(turing_data; vectorize=true)
m();
## ============================================================================
@info "start sampling"
Turing.emptyrdcache()

@time chain = let #460s; 366s
    thinning = 1
    if ps.chains > 1
        sample(m, NUTS(ps.warmup, 0.9; max_depth=5), MCMCThreads(), ps.steps + ps.warmup, ps.chains; progress=true, thinning)
    else
        sample(m, NUTS(ps.warmup, 0.9; max_depth=5), ps.steps + ps.warmup; progress=true, thinning) #; max_depth=15
    end
end #112.942085

divervences = sum( chain[ps.warmup+1:end, :numerical_error, :] )
@show divervences

## =========================================================================
# MODE ESTIMATION WITH OPTI
using Optim
using Underscores
using Accessors

function initvals(m, chain)
    varinfo = Turing.VarInfo(m)

    param = Dict()
    for key in keys(varinfo.metadata)
        val = @_ group(chain, key) |>
            DataFrame |>
            DataFrames.select(__, Not([:iteration, :chain])) |>
            Array |>
            mean(__, dims=1) |>
            vec
        if key == :latent_Rts_z
            val = reshape(val, (length(val)÷5, 5))
        end

        param[key] = length(val) == 1 ? first(val) : val
    end
    return NamedTuple(param)
end

init_tuple = initvals(
    pp.model(pp.data.turing_data),
    Turing.MCMCChains.get_sections(pp.chain, :parameters)
)

@set θ0.latent_Rts_z = similar(θ0.latent_Rts_z)

function optim_args(model, init_tuple)
    varinfo = Turing.VarInfo(model)
    model(varinfo, Turing.SampleFromPrior(), Turing.PriorContext(θ0))
    init_vec = varinfo[Turing.SampleFromPrior()]

    ctx = Turing.OptimizationContext(DynamicPPL.DefaultContext())
    f = Turing.OptimLogDensity(model, ctx)
    return f, init_vec
end

m = pp.model(pp.data.turing_data)
f, init_vec = optim_args(
      m
    , init_tuple
    )
f(init_vec)
optimize(m, MAP(), init_vec, ConjugateGradient(), Optim.Options(iterations=1000, allow_f_increases=true, show_trace=true))
##
using Accessors, NamedTupleTools
delete(pp.data, :num_total)
d.num_total_days

m = pp.model(pp.data.turing_data)
varinfo = Turing.VarInfo(m)
m(varinfo, Turing.SampleFromPrior(), Turing.PriorContext())
initθ = varinfo[Turing.SampleFromPrior()]

chain = sample(model, NUTS(), 100, init_theta=initθ)

##

x = NamedTuple(name => mode.values.array[i] for (name, i) in pairs(mode.values.dicts[1]))
DynamicPPL.setval_and_resample!(varinfo, x)
model(varinfo)
## ============================================================================
# improve initial parameterschain
function initparams(m)
    param_mod = deepcopy(m)
    varinfo = DynamicPPL.VarInfo(param_mod)
    # varinfo.metadata.R0s.ranges
    f = turing_data.invlink
    param_mod(varinfo, Turing.SampleFromPrior(),
        Turing.PriorContext((
            R0 = 1.,
            #latent_Rts_z = [MvNormal(turing_data.num_rt_steps[m], 0.01) .+ f(1) for m in 1:turing_data.num_regions],
        ))
    )
    varinfo[Turing.SampleFromPrior()]
end

chain = sample(m, NUTS(ps.warmup, 0.95; max_depth=5), ps.steps + ps.warmup; progress=true, discard_adapt = false)
initp = initparams(m)
chain = sample(m, NUTS(ps.warmup, 0.95; max_depth=5), ps.steps + ps.warmup; progress=true, init_params = initp)
plot(chain)
chain = sample(m, Turing.NUTS(1000, 0.8), 1000; init_params = init_theta, discard_adapt = false)

chains_params = Turing.MCMCChains.get_sections(chain, :parameters)
m_pred = model(turing_data, predict = true, context=:prediction)
gq = Turing.generated_quantities(m_pred, chains_params)
generated_posterior = Regional.posterior(model, turing_data, chains_params)

## ==========================================================================
@info "save"
fdir = projectdir("reports/tmp", "$(today())")
fname = Regional.save_results(fdir, ps, data, chain)
## ==========================================================================
@info "post processing"
using InlineStrings
using WeakRefStrings
Regional.postprocessing(fname; plot_results=false)
##
# fname = "/home/and/dev/CovidSurvey/reports/regional_effect_estimation/mobility/apple_totalcontacts/2021-09-05/CHAIN_chains=10_model=hospit_observ=2021-02-01_preds=CT,MA_regional=true_rwstep=7_semipar=true_std=false_steps=500_warmup=2000.jls"
fname = "/home/and/dev/CovidSurvey/reports/regional_prediction_end_of_peak/2021-09-07/CHAIN_chains=10_model=hospit_observ=2021-03-01_regional=true_rwstep=1_semipar=true_std=false_steps=500_warmup=2000.jls"
# fname = "/home/and/dev/CovidSurvey/reports/tmp/2021-08-31/CHAIN_cases=2020-06-01_chains=1_include=true_model=hospit_observ=2021-02-01_preds=CT,MNRES_rwstep=7_semipara=false_std=true_steps=10_warmup=10.jls"
fdir, ps, ignores = Regional.parse_fname(fname)
pp = Regional.PostProcessing2(fdir, ps, ignores, [], fname)
gp = Regional.generate_posterior(pp)
# gp = BSON.load(fname) |> NamedTuple;

## ===========================================================================
#  ARVIZ
# chain = deepcopy(chain)
turing_data_predictive = Regional.predictive(d)

predictive_model = model(turing_data_predictive, false, false)
posterior_predictive = predict(predictive_model, chain)

prior = sample(predictive_model, Prior(), 1000; progress=false)
prior_predictive = predict(predictive_model, prior)

loglikelihoods = Turing.pointwise_loglikelihoods(
    param_mod, MCMCChains.get_sections(chain, :parameters)
)

data_params_dct = convert(Dict, data_params)
data_params_dct[:predictive] = true
data_params = namedtuple(data_params_dct)
data = Regional.load_data(; data_params... )
turing_data = data.turing_data;

predictive_model = model1(turing_data, false)
posterior_predictive = predict(predictive_model, chain);

prior = sample(m, Prior(), 1000; progress=false)
prior_predictive = predict(param_mod_predict, prior):

loglikelihoods = Turing.pointwise_loglikelihoods(
    m, MCMCChains.get_sections(chain, :parameters)
)

using LinearAlgebra
# Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
ynames = string.(keys(posterior_predictive))
loglikelihoods_vals = getindex.(Ref(loglikelihoods), ynames)
# Reshape into `(nchains, nsamples, size(y)...)`
loglikelihoods_arr = permutedims(cat(loglikelihoods_vals...; dims=3), (2, 1, 3))

using LinearAlgebra
# Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
ynames = string.(keys(posterior_predictive))
loglikelihoods_vals = getindex.(Ref(loglikelihoods), ynames)
# Reshape into `(nchains, nsamples, size(y)...)`
loglikelihoods_arr = permutedims(cat(loglikelihoods_vals...; dims=3), (2, 1, 3))

## ==========================================================================
@info "post processing"
fdir, ps, ignores = Regional.parse_fname(fname)
p = Regional.PostProcessing2(fdir, ps, ignores, [], fname)
Regional.savechain(p)
## ==========================================================================
@info "plot chain"
plot_results = true
Regional.plot_chains(p; plot_results)

## ==========================================================================
@info "meanplot"
Regional.plot_means(p; plot_results)
## ==========================================================================
# perform diagnostics
p = Regional.skip_warmup(p)
Regional.diagnostics(p)
## ==========================================================================
@info "make predictions"
gp = Regional.generate_posterior(p);
## ==========================================================================
@info "plot regions"
Regional.plot_regions(p, gp; plot_results=true)

recipe = Regional.RegionPlottingRecipe(p.data, gp, "capital")
plot(recipe)
## ==========================================================================
@info "plot predictors"
# pgfplotsx()
# default(titlefontsize = 20, legendfontsize = 18, labelfontsize = 18, guidefontsize = 18, tickfontsize = 12, framestyle = :zerolines, yminorgrid = true)
Regional.plot_effects(p, gp; plot_results , grouped = false, effect_on_Rt = 0.)
Regional.plot_effects(p, gp; plot_results, grouped = false, effect_on_Rt = 2.)
Regional.plot_effects(p, gp; plot_results, grouped = false, effect_on_Rt = -0.5)

Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = 0.)
Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = 2.)
Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = -0.5)

## -----------------------------------------------------------------------------
@info "store reproduction number"
Regional.save_rt(p, gp)
## ===========================================================================
# ArviZ

using ArviZ
import PyPlot as plt

c = MCMCChains.get_sections(chain, :parameters)
n = filter(
    x-> !occursin("latent_Rts", x) &&
        !occursin("effect", x)     &&
        !occursin("ys", x)     &&
        !occursin("R0", x), String.(names(c)))
plot_autocorr(c; var_names=n);
plt.gcf()
##
idata = from_mcmcchains(
    p.chain;
    #coords=Dict("school" => schools),
    #dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)
plot_pair(
    idata;
    #coords=Dict("school" => ["Choate", "Deerfield", "Phillips Andover"]),
    divergences=true,
);

fn = fdir*"tmp.png"
plt.gcf().savefig(fn)
run(`firefox $(fn)`, wait=false)
##

# function posterior(model, turing_data, chain)
m_pred = model(turing_data, predict = true, context=:prediction)
gq = Turing.generated_quantities(m_pred, chain)
cat(gq...; dims=3)

gp = vectup2tupvec(reshape(gq, length(gq)))
groupregions(gp)
# end

turing_data_predictive = Regional.predictive(d)
predictive_model = model(turing_data_predictive,
    turing_data_predictive.hospits,
    vectorize = true,
    predict = false)
posterior_predictive = Turing.predict(predictive_model, chain)

predictive_model()

m = 1
s = d.hospitmodel.starts[m]
e = d.hospitmodel.stops[m]
ts= s:e
length(ts)

prior = sample(predictive_model, Prior(), 100; progress=false)
prior_predictive = predict(predictive_model, prior)

param_mod = model(turing_data, false, false)
loglikelihoods = Turing.pointwise_loglikelihoods(
    param_mod, MCMCChains.get_sections(chain, :parameters)
)

using LinearAlgebra
# Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
ynames = string.(keys(posterior_predictive))
loglikelihoods_vals = getindex.(Ref(loglikelihoods), ynames)
# Reshape into `(nchains, nsamples, size(y)...)`
loglikelihoods_arr = permutedims(cat(loglikelihoods_vals...; dims=3), (2, 1, 3))

idata = from_mcmcchains(
    turing_chns;
    posterior_predictive=posterior_predictive,
    log_likelihood=Dict("y" => loglikelihoods_arr),
    prior=prior,
    prior_predictive=prior_predictive,
    observed_data=Dict("y" => y),
    coords=Dict("school" => schools),
    dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)
ArviZ.loo(idata)

using MCMCChains
MCMCChains.header(chain, section = :parameters)
MCMCChains.sections(chain)

idata = from_mcmcchains(
    turing_chns,
    coords = Dict("school" => schools),
    dims = Dict(
        "y" => ["school"],
        "σ" => ["school"],
        "θ" => ["school"],
    ),
    library = "Turing",
)

plot_pair(
    idata;
    coords=Dict("school" => ["Choate", "Deerfield", "Phillips Andover"]),
    divergences=true,
);
plt.gcf()

plot_rank(post; var_names=["μ", "τ"]);
gcf()

plot_density(
    [idata.posterior_predictive, idata.prior_predictive];
    data_labels=["Post-pred", "Prior-pred"],
    var_names=["y"],
)
gcf()

##
TV = Vector{Float64}
V = Float64
vectorize = true
θ = turing_data
@unpack num_observations, num_total_days, num_regions, num_rt_steps, hospitmodel, seromodel, num_covariates, rwscale, include_early_dynamic = θ
# If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
num_time_steps = num_observations

############# 2.) time varying reproduction number

R0_std     = truncated(Normal(0.5, 0.2), 0.01, Inf) |> rand
R0         = truncated(Normal(3.5, 1.), 1., 4.5) |> rand
R0s        = filldist(truncated(Normal(R0, R0_std), 1., 5.), num_regions) |> rand
# R0         = truncated(Normal(3.0, 0.5), 0., 6.)
# R0s        = filldist(truncated(Normal(R0, 0.2), 0., 6.), num_regions) #CovidSurvey.CustomMvNormal1(R0, 0.2, 1., 4.5, num_regions) #filldist(truncated(Normal(R0, 0.2), 1., 4.5), num_regions)
σ_rt       = truncated(Normal(0.3*rwscale, .02*rwscale), 0, .5*rwscale) |> rand #N(0.32,0.5)

init_rw = if include_early_dynamic
        R0s
    else
        # R1  = truncated(Normal(0.8, 0.1), 0., 2.)
        # R1s = filldist(truncated(Normal(R1, 0.1), 0., 2.), num_regions)

        R1_std    = truncated(Normal(0.1, .05), 0.05, Inf) |> rand
        R1s       = filldist(truncated(Normal(0.8, R1_std), 0., 1.5), num_regions) |> rand
        R1s
end

latent_Rts = if include_early_dynamic
    latent_Rts_z = [TV(undef, num_rt_steps[m]) for m in 1:num_regions]
    latent_Rts = [TV(undef, num_rt_steps[m]) for m in 1:num_regions]
    for m in 1:num_regions
        latent_Rts_z[m] = CovidSurvey.RandomWalk(num_rt_steps[m]) |> rand
        @. latent_Rts[m] = latent_Rts_z[m] * σ_rt + invlink(init_rw[m])
    end
    latent_Rts
else
    n = first( num_rt_steps )
    latent_Rts_z = filldist( CovidSurvey.RandomWalk( n ), num_regions) |> rand
    CovidSurvey.rescale( latent_Rts_z, σ_rt, θ.invlink.(init_rw) )
end
# latent_Rts_z = [TV(undef, num_rt_steps[m]) for m in 1:num_regions]
# latent_Rts = [TV(undef, num_rt_steps[m]) for m in 1:num_regions] |> rand
# for m in 1:num_regions
# 	latent_Rts_z[m] = RandomWalk(num_rt_steps[m])
# 	latent_Rts[m] = latent_Rts_z[m] * σ_rt .+ invlink(init_rw[m])
# end

Rts = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
if num_covariates > 0
    grouped_effects = filldist( Laplace(0,0.2), num_covariates) |> rand #Exponential(0.2)
    effect_std = filldist( GammaMeanCv(0.1, 0.5), num_covariates) |> rand
    effects_z = filldist( MvNormal( num_covariates, 1.), num_regions) |> rand
    effects = [ effects_z[:,m] .* effect_std .+ grouped_effects for m in 1:num_regions]

    # grouped_effect = filldist( Exponential(0.2), num_covariates)
    # effect_std = filldist( GammaMeanCv(0.1, 0.5), num_covariates)
    # pooled_effects_distr_raw = Normal.( grouped_effect, effect_std )
    # pooled_effects_distr = arraydist( truncated.( pooled_effects_distr_raw, zero(V), V(Inf) ) )
    # effects = filldist( pooled_effects_distr, num_regions )

    if θ.semiparametric
        CovidSurvey.semiparametric!(Rts, θ, false, latent_Rts, R0s, σ_rt, effects)
    else
        CovidSurvey.mixed!(Rts, θ, false, latent_Rts, R0s, σ_rt, effects)
    end
else
    CovidSurvey.random_walks!(Rts, θ, false, latent_Rts, R0s, σ_rt)
end
m = 1
plot(latent_Rts[:,m])
plot(d.dates[m], Rts[m])
############ 3.) infection dynamics
τ  = Exponential(20) |> rand# `Exponential` has inverse parameterization of the one in Stan
T  = typeof(τ)
# ys = filldist(truncated(Normal(τ, 30),T(1),T(500)), num_regions)
ys = filldist(Exponential(τ), num_regions) |> rand

newly_infecteds      = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
cumulative_infecteds = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
# effective_Rts        = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]

CovidSurvey.infections!(newly_infecteds, cumulative_infecteds, θ, ys, Rts)
plot(d.dates[m], newly_infecteds[m], yticks=:native)
# infections!(newly_infecteds, cumulative_infecteds, effective_Rts, θ, τ, ys, Rts)
########### 4.) derive observables
μ_i2h = truncated(Normal(11.7, .5), 8, 14) |> rand
σ_i2h = truncated(Normal(6.5,0.5), 1, 14) |> rand
ihr   = truncated(Normal(1.8/100,0.5/100), 0.5/100, 5/100) |> rand #GammaMeanCv(0.025, 0.1)
ϕ_h   = truncated(Normal(50, 10), 20, Inf) |> rand

expected_daily_hospit    = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
hospit_observation_model = CovidSurvey.SimpleObsModel2(hospitmodel, μ_i2h, σ_i2h, ihr, ϕ_h, expected_daily_hospit, vectorize)
# hospit_observation_model = CovidSurvey.SimpleObsModel(hospitmodel, μ_i2h, ihr, ϕ_h, expected_daily_hospit, vectorize)
sero_observation_model   = CovidSurvey.SimpleSeroObsModel(seromodel, cumulative_infecteds)
CovidSurvey.expected!(hospit_observation_model, newly_infecteds)
plot(d.dates[m], expected_daily_hospit[m], yticks=:native)

########### 4.) compare model to observations
## 4.1) observe hospitalizations

ℓ  = zero(V)
ℓ += CovidSurvey.logpdf(hospit_observation_model, θ.hospits)
ℓ += CovidSurvey.logpdf(sero_observation_model, θ.sero)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
using FinEtools
using FinEtools.TetRemeshingModule
using FinEtools
using FinEtools.VoxelBoxModule
using FinEtools.VoxelTetMeshingModule
using FinEtools.MeshExportModule
using NIfTI

File = "Head_256x256x126-256x256x252"
file = niread(File * ".nii")
bs = voxel_size(file.header) .* [size(file,1)-100, size(file,2)-100, size(file,3)-100]
vb = VoxelBoxVolume(file.raw[101:end, 101:end, 101:end, 1], bs)
# vtkexport(File * ".vtk", vb)

im = ImageMesher(vb, zero(eltype(vb.data)), eltype(vb.data)[255])
mesh!(im)
println("Mesh size: initial = $(size(im.t,1))")
im.elementsizeweightfunctions = [ElementSizeWeightFunction(20.0, vec([112.0, 192.5, 26.5]), 45.0), ElementSizeWeightFunction(20.0, vec([157.0, 160.5, 0.5]), 45.0), ElementSizeWeightFunction(20.0, vec([56.0, 160.5, 0.5]), 45.0)]
for i = 1:5
    mesh!(im, 1.2)
    println("Mesh size: intermediate = $(size(im.t,1))")
end

fens = FENodeSet(im.v)
fes = FESetT4(im.t)
setlabel!(fes, im.tmid)

File = "Head_mesh_3.vtk"
vtkexportmesh(File, fens, fes)
@async run(`"paraview.exe" $File`)

ne = NASTRANExporter("Head_mesh_3.nas")
BEGIN_BULK(ne)
for i = 1:count(fens)
    GRID(ne, i, fens.xyz[i, :])
end
for i = 1:count(fes)
    CTETRA(ne, i, 1, fes.conn[i, :])
end
PSOLID(ne, 1, 1)
MAT1(ne, 1, 20.0e9, 0.3, 1800.0)
ENDDATA(ne)
close(ne)


text:
function LayerTrace(gml,target,tileReg;omitThemes::Array=[],colour::String="#aaaaaa",CropOutOfBounds::Bool=false)
    iKeep = findall([length(findall(descr.==target))>0 for descr in gml.summary.descr])
    if !isempty(omitThemes)
        iKeep2 = vcat([findall([length(findall(theme.==omitTheme))==0 for theme in gml.summary.theme]) for omitTheme in omitThemes]...)
        iKeep = intersect(iKeep,iKeep2)
    end

    xs = vcat(gml.summary.eastings[iKeep]...)
    ys = vcat(gml.summary.northings[iKeep]...)

    if CropOutOfBounds == true
        ob_up = findall([sum(y.>tileReg.northings_max)==length(y) for y in ys])
        ob_down = findall([sum(y.<tileReg.northings_min)==length(y) for y in ys])
        ob_left = findall([sum(x.<tileReg.eastings_min)==length(x) for x in xs])
        ob_right = findall([sum(x.>tileReg.eastings_max)==length(x) for x in xs])

        iOB = unique(vcat(ob_up,ob_down,ob_left,ob_right))
        xs = xs[Not(iOB)]
        ys = ys[Not(iOB)]

        ob_up = [findall(y.>tileReg.northings_max) for y in ys]
        ob_down = [findall(y.<tileReg.northings_min) for y in ys]
        ob_left = [findall(x.<tileReg.eastings_min) for x in xs]
        ob_right = [findall(x.>tileReg.eastings_max) for x in xs]

        iRm = [unique(vcat(ob_up[i],ob_down[i],ob_left[i],ob_right[i])) for i in 1:length(xs)]
        [xs[i] = xs[i][Not(iRm[i])] for i in 1:length(xs)]
        [ys[i] = ys[i][Not(iRm[i])] for i in 1:length(xs)]

    end

    xs = vcat([vcat(pts,NaN) for pts in xs]...)
    ys = vcat([vcat(pts,NaN) for pts in ys]...)

    trace = scatter(x=xs,y=ys,name=gml.tile,mode="lines",line=attr(width=0.2,color=colour),hoverinfo="skip",hovertemplate=nothing)
    return trace
end

function SelectTile(env,project,ui,tileRef)

    # Index of current tile in register
    iTile = findfirst(project.tileRegister.tileRefs.==tileRef)

    # "loading" pop-up
    w_popup = Window(Dict(:title=>"Loading...",:width=>500,:height=>150))
    body!(w_popup,"Loading $(tileRef) ...")
    sleep(0.1)

    # Clear blank initialising trace
    iRm = findall([tr[:name] for tr in ui.p0.plot.data].=="blank")
    [deletetraces!(ui.p0,i) for i in iRm]

    # Check for pre-existing parsed GML
    GmlParsedDir = readdir(joinpath(env.paths[:projects],".GmlParsed"))
    if sum(occursin.(tileRef,GmlParsedDir)).==1
        # Load data for new tile
        data_dict = readJSON(".GmlParsed",tileRef)
        data = DataFrame(data_dict)

        gml = ParaDwell.GML(tileRef,"",[""],[[""]],data)
    else
        # Load data for new tile
        gml = ParaDwell.ProcessGMLs(env.paths[:OS],tileRef)
        json_string = JSON.json(gml.summary)
        writeJSON(json_string,".GmlParsed",tileRef)
    end

    if !haskey(project.dat,"gml")
        push!(project.dat,"gml"=>[gml])
    else
        push!(project.dat["gml"],gml)
    end

    traces = []

    # Plot tile boundary
    x_boundary_OSEN = [project.tileRegister.eastings_min[iTile];project.tileRegister.eastings_max[iTile];project.tileRegister.eastings_max[iTile];project.tileRegister.eastings_min[iTile];project.tileRegister.eastings_min[iTile]]
    y_boundary_OSEN = [project.tileRegister.northings_min[iTile];project.tileRegister.northings_min[iTile];project.tileRegister.northings_max[iTile];project.tileRegister.northings_max[iTile];project.tileRegister.northings_min[iTile]]
    tmp = [ParaDwell.OSENtoLLA(x_boundary_OSEN[i],y_boundary_OSEN[i]) for i in 1:length(y_boundary_OSEN)]
    x_boundary_LLA = [i[2] for i in tmp]
    y_boundary_LLA = [i[1] for i in tmp]
    push!(traces,scatter(
        x=x_boundary_OSEN,
        y=y_boundary_OSEN,
        mode="lines",
        line=attr(color="#aaaaaaaa",width=1),
        name=tileRef,
        hoverinfo="skip",hovertemplate=nothing
    ))
    traces_boundary_macroLLA = scatter(
        x=x_boundary_LLA,
        y=y_boundary_LLA,
        mode="lines",
        line=attr(color="#aaaaaaaa",width=1),
        name=tileRef,
        hoverinfo="skip",hovertemplate=nothing
    )

    # Plot map layers
    push!(traces,LayerTrace(gml,"Natural Environment",project.tileRegister[iTile,:];colour="#aaaaaa",CropOutOfBounds=true))
    push!(traces,LayerTrace(gml,"Road Or Track",project.tileRegister[iTile,:];colour="#ff9f43",CropOutOfBounds=true))
    push!(traces,LayerTrace(gml,"Building",project.tileRegister[iTile,:];omitThemes=["Land"],colour="#1e3799"))

    # Add all new traces
    [addtraces!(ui.p0,trace) for trace in traces] # Main UI
    addtraces!(ui.plt_macro,traces_boundary_macroLLA) # Macro-map UI
    sleep(0.1)
    close(w_popup)

    return project
end

function ClickTile(env,project,ui)

    on(ui.plt_macro["click"]) do dat

        # Determine selected point
        X = dat["points"][1]["x"]
        Y = dat["points"][1]["y"]

        # Find selected tile
        iTile  = intersect(findall(project.tileRegister.Lon_mid.==X),findall(project.tileRegister.Lat_mid.==Y))[1]
        tileRef = project.tileRegister.tileRefs[iTile]

        # Print to console
        println("TileRef selected: $tileRef")

        # Gather and plot data for selected tile
        project = SelectTile(env,project,ui,tileRef)
    end
end


text:
# A Julia implementation of QuickCheck, a specification-based tester
#
# QuickCheck was originally written for Haskell by Koen Claessen and John Hughes
# http://www.cse.chalmers.se/~rjmh/QuickCheck/

module QuickCheck

export property
export condproperty
export quantproperty

function lambda_arg_types(f::Function)
    if !isa(f.code, LambdaStaticData)
        error("You must supply either an anonymous function with typed arguments or an array of argument types.")
    end
    [eval(var.args[2]) for var in Base.uncompressed_ast(f.code).args[1]]
end

# Simple properties
function property(prop::Function, typs::Vector, ntests)
    arggens = [size -> generator(typ, size) for typ in typs]
    quantproperty(prop, typs, ntests, arggens...)
end
property(prop::Function, typs::Vector) = property(prop, typs, 100)
property(prop::Function, ntests) = property(prop, lambda_arg_types(prop), ntests)
property(prop::Function) = property(prop, 100)

# Conditional properties
function condproperty(prop::Function, typs::Vector, ntests, maxtests, argconds...)
    arggens = [size -> generator(typ, size) for typ in typs]
    check_property(prop, arggens, argconds, ntests, maxtests)
end
condproperty(prop::Function, args...) = condproperty(prop, lambda_arg_types(prop), args...)

# Quantified properties (custom generators)
function quantproperty(prop::Function, typs::Vector, ntests, arggens...)
    argconds = [(_...)->true for t in typs]
    check_property(prop, arggens, argconds, ntests, ntests)
end
quantproperty(prop::Function, args...) = quantproperty(prop, lambda_arg_types(prop), args...)

function check_property(prop::Function, arggens, argconds, ntests, maxtests)
    totalTests = 0
    for i in 1:ntests
        goodargs = false
        args = {}
        while !goodargs
            totalTests += 1
            if totalTests > maxtests
                println("Arguments exhausted after $i tests.")
                return
            end
            args = [arggen(div(i,2)+3) for arggen in arggens]
            goodargs = all([cond(args...) for cond in argconds])
        end
        if !prop(args...)
            error("Falsifiable, after $i tests:\n$args")
        end
    end
    println("OK, passed $ntests tests.")
end

# Default generators for primitive types
generator{T<:Unsigned}(::Type{T}, size) = convert(T, rand(1:size))
generator{T<:Signed}(::Type{T}, size) = convert(T, rand(-size:size))
generator{T<:FloatingPoint}(::Type{T}, size) = convert(T, (rand()-0.5).*size)
# This won't generate interesting UTF-8, but doing that is a Hard Problem
generator{T<:String}(::Type{T}, size) = convert(T, randstring(size))

generator(::Type{Any}, size) = error("Property variables cannot by typed Any.")

# Generator for array types
function generator{T,n}(::Type{Array{T,n}}, size)
    dims = [rand(1:size) for i in 1:n]
    reshape([generator(T, size) for x in 1:prod(dims)], dims...)
end

# Generator for composite types
function generator{C}(::Type{C}, size)
    if C.types == ()
        error("No generator defined for type $C.")
    end
    C([generator(T, size) for T in C.types]...)
end

end


text:
# Testing continuous univariate distributions

using Distributions
using Base.Test

### load reference data
#
#   Note
#   -------
#   To generate the reference data:
#   (1) make sure that python, numpy, and scipy are installed in your system
#   (2) enter the sub-directory test
#   (3) run: python discrete_ref.py > discrete_ref.csv
#
#   For most cases, you don't have. You only need to run this when you
#   implement a new distribution and want to add new test cases, then
#   you should add the new test cases to discrete_ref.py and run this
#   procedure to update the reference data.
#
n_tsamples = 100

# additional distributions that have no direct counterparts in scipy
println("    -----")

for distr in [
    Biweight(),
    Biweight(1,3),
    Epanechnikov(),
    Epanechnikov(1,3),
    Frechet(0.5, 1.0),
    Frechet(3.0, 1.0),
    Frechet(20.0, 1.0),
    Frechet(120.0, 1.0),
    Frechet(0.5, 2.0),
    Frechet(3.0, 2.0),
    GeneralizedPareto(),
    GeneralizedPareto(1.0, 1.0),
    GeneralizedPareto(1.0, 1.0, 1.0),
    GeneralizedPareto(0.1, 2.0, 0.0),
    GeneralizedPareto(0.0, 0.5, 0.0),
    GeneralizedPareto(-1.5, 0.5, 2.0),
    InverseGaussian(1.0, 1.0),
    InverseGaussian(2.0, 7.0),
    Levy(0.0, 1.0),
    Levy(2.0, 8.0),
    Levy(3.0, 3.0),
    LogNormal(0.0, 1.0),
    LogNormal(0.0, 2.0),
    LogNormal(3.0, 0.5),
    LogNormal(3.0, 1.0),
    LogNormal(3.0, 2.0),
    NoncentralBeta(2,2,0),
    NoncentralBeta(2,6,5),
    NoncentralChisq(2,2),
    NoncentralChisq(2,5),
    NoncentralF(2,2,2),
    NoncentralF(8,10,5),
    NoncentralT(2,2),
    NoncentralT(10,2),
    Triweight(),
    Triweight(1,3),
]
    println("    testing $(distr)")
    test_distr(distr, n_tsamples)
end

# for distr in [
#     VonMises(0.0, 1.0),
#     VonMises(0.5, 1.0),
#     VonMises(0.5, 2.0) ]

#     println("    testing $(distr)")
#     test_samples(distr, n_tsamples)
# end


text:
include("splinesParallel.jl")

knots = [0,0,0,0,1,2,3,4,4,4,4]
ncontrols = 8
degree = 2
basis = pbsplineBasis(degree, knots, ncontrols)
point = [0.7,0.2,3.0]
basis(point)


text:
using MATLAB

type BaseLurieAdCertificate
  lf::LurieFormAd

  B::AbstractMatrix
  C::AbstractMatrix
  R::AbstractMatrix

  D::AbstractMatrix
  E::AbstractMatrix

  function BaseLurieAdCertificate(lf::LurieFormAd, R)
    A = create_A(lf)
    B, C = create_BC(lf, A, R)
    D, E = create_DE(lf, lf.SQ, lf.SQJ, R)
    new(lf, B, C, R, D, E)
  end
end

function create_A(lf::LurieFormAd)
  @import_fields lf.sys
  E = sparse([1:np;1:np], [from; to], [ones(np); -ones(np)], np, nb)

  A = [ spzeros(npq, npv+ npq)    speye(npq) #voltage pertubation limits
        spzeros(np, npq)   0*E[:,pv]  E[:,pq] #voltage discrepancy limits
        E[:,pv]  E[:,pq]   spzeros(np, npq) ] #angular discrepancy limits
  return A
end

# assuming δs = R*u
function create_BC(lf::LurieFormAd, A::AbstractMatrix, R::AbstractMatrix)
  B =  A*(lf.J⁻¹)*R
  C = -A*(lf.J⁻¹)*lf.S
  return (B, C)
end

function create_DE(lf::LurieFormAd, SQ::AbstractMatrix, SQJ::AbstractMatrix, R::AbstractMatrix)
  D = SQJ*(lf.J⁻¹)*R
  E = SQ - SQJ*(lf.J⁻¹)*lf.S
  return (D, E)
end

function create_box_R(lf::LurieFormAd, r::Real)
  @import_fields lf.sys

  # Find loads with nonzero value
  nzpq = find(b->(p⁰[pq[b]]!=0)||(q⁰[pq[b]]!=0), 1:npq)
  rnzpq = r*ones(length(nzpq))

  return create_loadability_R(lf, nzpq, rnzpq)
end

create_loadability_R(lf::LurieFormAd) = create_loadability_R(lf, lf.s⁰)
create_loadability_R(lf::LurieFormAd, ds::Vector) = reshape(ds, length(ds), 1)

function create_loadability_R(lf::LurieFormAd, nzpq::Vector, r::Vector)
  @import_fields lf.sys
  v⁰ = lf.sys.v⁰

  R = zeros(npv + 2*npq, length(nzpq))
  for b in 1:length(nzpq)
    R[npv + nzpq[b], b] = r[b]*p⁰[pq[nzpq[b]]]/(v⁰[pq[nzpq[b]]]^2)
    R[npv + npq + nzpq[b], b] = r[b]*q⁰[pq[nzpq[b]]]/(v⁰[pq[nzpq[b]]]^2)
  end

  return R
end

function create_RX_branch(cert::BaseLurieAdCertificate)
    sys = cert.lf.sys
    np = cert.lf.sys.np
    RX_branch = zeros(np,1)*im
    for p in 1:np
      RX_branch[p] = -1/sys.Y[sys.from[p],sys.to[p]]
    end
    return real(RX_branch), imag(RX_branch), RX_branch

end

function create_Imax(lf::LurieFormAd, rate::Real)
  YF = lf.sys.Yᶠ' - lf.sys.Yᵗ'
  V⁰ = lf.sys.v⁰ .* exp(lf.sys.ϑ⁰*im)
  Iₑ = YF*V⁰

  return abs(Iₑ)*rate
end

function get_Qlims(cert::BaseLurieAdCertificate, pm::GenericPowerModel)
  q⁰ = cert.lf.sys.q⁰
  net = cert.lf.sys.buses
  pv = cert.lf.sys.pv
  sl = cert.lf.sys.sl
  npv = length(pv)
  Qmin = []
  Qmax = []
  genbus = []
  for (gen, gen_data) in pm.data["gen"]
    Qmin = [Qmin; gen_data["qmin"]]
    Qmax = [Qmax; gen_data["qmax"]]
    genbus = [genbus; gen_data["gen_bus"]]
  end
   ngen = length(genbus)
   pos = []
  for k in 1:ngen
      pos = [pos; find(net -> net == genbus[k], net)]
  end
  pos[find(pos -> pos == sl[1], pos)] = 1e+7

  #reorder
  pos_in_pm = sortperm(pos)
  Qmin = Qmin[pos_in_pm]
  Qmax = Qmax[pos_in_pm]

  #allowed deviations
  dQmin = 1e+7*ones(npv)
  dQmax = 1e+7*ones(npv)

  for k in 1:npv
    if Qmin[k] <= q⁰[pv[k]]
       dQmin[k] = q⁰[pv[k]] - Qmin[k]
    end
    if q⁰[pv[k]]<= Qmax[k]
       dQmax[k] = -q⁰[pv[k]] + Qmax[k]
     end
  end
  return [dQmax dQmin]
end

function cal_box(R, Δ₁u, Δ₁v, npv, buspJ, basep, v⁰)
    resB = zeros(2, 2)
    res = R*Δ₁u[:,1]
    for k in 1:length(buspJ)
        resB[:, k] = -res[[npv+buspJ[1],npv+buspJ[2]]]
    end
    for k in 1:length(buspJ)
      resB[k, 1] = (v⁰[buspJ[k]] - Δ₁v[buspJ[k],2])^2 * (resB[k, 1] + basep[k]/(v⁰[buspJ[k]])^2) #max
      resB[k, 2] = (v⁰[buspJ[k]] - Δ₁v[buspJ[k],2])^2 * (-resB[k, 2] + basep[k]/(v⁰[buspJ[k]])^2) #min
    end

    return resB
end

function polytope_plot(casename,busp,nump,resB,basep,nb,ratio,thermalrate,dV,stepsize,adaptive,limitQQ,limitIQ,Δ₁u)
    if Δ₁u[1,1] > 0
        @mput casename busp nump resB basep nb ratio thermalrate dV stepsize adaptive limitQQ limitIQ
        eval_string("
                    close all
                    addpath(genpath('/Users/hunghtd/Dropbox (Personal)/Journal/Feasibility region/matpower6.0b2'));

                    resR = constrainedCPF(casename, busp, nump, thermalrate, dV, stepsize, adaptive, limitQQ, limitIQ, 2*pi);

                    polytope_plot(resB, resR, basep, busp)

                    ")
    end

end

function add_Δ₁y_constraints_LP(cert::BaseLurieAdCertificate, m::Model, thermalrate::Real, T, K, dV::Real, limitIQ)
  @import_fields cert.lf.sys

  Δ₁y = getvariable(m, :Δ₁y)

  #current limits
  Imax = create_Imax(cert.lf, 2)
  R, X, Z = create_RX_branch(cert)
  Z = abs(Z)
  v⁰ = cert.lf.sys.v⁰
  ϑ⁰ = cert.lf.sys.ϑ⁰
  vt⁰ = v⁰[to]
  ρₑ⁰ = log(v⁰[from]./v⁰[to])
  Θₑ⁰ = ϑ⁰[from] - ϑ⁰[to]
  npvpq = npv+npq

  @variables m begin
    Δ₁v[1:npq,±] >= 0(start = 0.0)
    Δ₁ρ[1:np,±] >= 0(start = 0.0)
    Δ₁Θ[1:np,±] >= 0(start = 0.0)

    Δ₁ρᵐ[1:np] >= 0, (start = 0.0) #max of   Δ₁ρ\^±
    Δ₁Θᵐ[1:np] >= 0, (start = 0.0)
  end

  @constraints m begin
    [b in 1:npq,   s in ±], Δ₁y[b,s] == Δ₁v[b,s]
    [b in 1:npq,   s in ±], Δ₁y[b,s] <= dV
    # [p in 1:np,    s in ±], Δ₁y[npq + p,s] == Δ₁ρ[p,s]
    [p in 1:np,    s in ±], Δ₁y[npq + p,s] == Δ₁ρᵐ[p]
    # [p in 1:np,    s in ±], Δ₁y[npq + p,s] <= dV
    # [p in 1:np,    s in ±], Δ₁y[npq + np + p,s] == Δ₁Θ[p,s]
    [p in 1:np,    s in ±], Δ₁y[npq + np + p,s] == Δ₁Θᵐ[p]
    # [p in 1:np,    s in ±], Δ₁y[npq + np + p,s] <= T

    [b in 1:npq,   s in ±], Δ₁Θ[b,s] <= Δ₁Θᵐ[b]
    [p in 1:np           ], Δ₁Θᵐ[p] <= T

    [p in 1:np,   s in ±], Δ₁ρ[p,s] <= Δ₁ρᵐ[p]
    [p in 1:np          ], Δ₁ρᵐ[p] <= K
  end

  if limitIQ > 0 #if current limits are considered
      for p in 1:np, s = ±
          t = to[p]
          s̄ = negate(s)
          istobuspq = length(find(pq -> pq == t,pq)) > 0
          dVmaxload = Z[p]*Imax[p]/(vt⁰[p]+dV) #max voltage drop %
          dVmaxgen = Z[p]*Imax[p]/vt⁰[p]

          if istobuspq #to bus is a pq bus
              @constraints m begin
                (1+2*Δ₁ρᵐ[p])*exp(2*ρₑ⁰[p]) - 2*exp(ρₑ⁰[p])*cos(T+Θₑ⁰[p])*(T*(1+Δ₁ρᵐ[p]+ρₑ⁰[p]) + Δ₁Θᵐ[p]*(1+K+ρₑ⁰[p]) - T*(1+K+ρₑ⁰[p])) - 1.2 <= dVmaxload^2
              end
          else
              @constraints m begin
                (1+2*Δ₁ρᵐ[p])*exp(2*ρₑ⁰[p]) - 2*exp(ρₑ⁰[p])*cos(T+Θₑ⁰[p])*(T*(1+Δ₁ρᵐ[p]+ρₑ⁰[p]) + Δ₁Θᵐ[p]*(1+K+ρₑ⁰[p]) - T*(1+K+ρₑ⁰[p])) - 1.2 <= dVmaxgen^2
              end
          end
        end
    end

end

function add_Δ₂f_constraints_LP(cert::BaseLurieAdCertificate, m::Model, Qlim::AbstractMatrix, T, K, limitQQ)
  @import_fields cert.lf.sys
  coshcos = range_coshcos(cert.lf)
  sinhcos = range_sinhcos(cert.lf)
  coshsin = range_coshsin(cert.lf)
  sinhsin = range_sinhsin(cert.lf)

  Δ₂f = getvariable(m, :Δ₂f);
  Δ₁ρ = getvariable(m, :Δ₁ρ);
  Δ₁Θ = getvariable(m, :Δ₁Θ);
  Δ₁ρᵐ = getvariable(m, :Δ₁ρᵐ);
  Δ₁Θᵐ = getvariable(m, :Δ₁Θᵐ);
  Δ₁u = getvariable(m, :Δ₁u);

  #constraints on reactive powers
  pv = cert.lf.sys.pv
  v⁰ = cert.lf.sys.v⁰
  npv = length(pv)

  SQ = cert.lf.SQ #last rows of the full matrix S
  SQ⁺, SQ⁻ = split_matrix(SQ)

  B = cert.B
  B⁺,B⁻ = split_matrix(B)

  D = cert.D
  D⁺, D⁻ = split_matrix(D)

  E = cert.E
  E⁺, E⁻ = split_matrix(E)

  nu = size(B,2)
  nf = size(cert.C,2)

  @variables m begin
   Δ₁cos[1:np,⁻] >= 0.0, (start = 0.0)
   Δ₂cos[1:np,⁻] >= 0.0, (start = 0.0)
   Δ₁sin[1:np,±] >= 0.0, (start = 0.0)
   Δ₂sin[1:np,±] >= 0.0, (start = 0.0)
   Δ₁cosh[1:np,⁺] >= 0.0, (start = 0.0)
   Δ₂cosh[1:np,⁺] >= 0.0, (start = 0.0)
   Δ₁sinh[1:np,±] >= 0.0, (start = 0.0)
   Δ₂sinh[1:np,±] >= 0.0, (start = 0.0)
  end

  for p in 1:np
    f = from[p]
    t = to[p]
    @constraints m begin
          Δ₁cos[p,⁻] == (1-cos(T))*Δ₁Θᵐ[p]/T
          Δ₂cos[p,⁻] == (1-cos(T))*Δ₁Θᵐ[p]/T

          Δ₁cosh[p,⁺] == (cosh(K)-1)*Δ₁ρᵐ[p]/K
          Δ₂cosh[p,⁺] == (cosh(K)-1)*Δ₁ρᵐ[p]/K

          Δ₂f[coshcos[p],⁺] == Δ₂cosh[p,⁺]

          Δ₂f[coshcos[p],⁻] >= (cosh(K)-1)*(1-cos(T))*Δ₁Θᵐ[p]/T + Δ₂cos[p,⁻]
          Δ₂f[coshcos[p],⁻] >= (cosh(K)-1)*(1-cos(T))*Δ₁ρᵐ[p]/K + Δ₂cos[p,⁻]
    end

    for s in ±
      s̄ = negate(s)
      @constraints m begin
            Δ₁sin[p,s] == Δ₁Θ[p,s]
            Δ₂sin[p,s] == (T-sin(T))*Δ₁Θ[p,s̄]/T
            # Δ₁sin[p,s] == Δ₁Θᵐ[p]
            # Δ₂sin[p,s] == (T-sin(T))*Δ₁Θᵐ[p]/T

            Δ₁sinh[p,s] == sinh(K)*Δ₁ρ[p,s]/K
            Δ₂sinh[p,s] == (sinh(K)-K)*Δ₁ρ[p,s]/K
            # Δ₁sinh[p,s] == sinh(K)*Δ₁ρᵐ[p]/K
            # Δ₂sinh[p,s] == (sinh(K)-K)*Δ₁ρᵐ[p]/K

            Δ₂f[coshsin[p],s] >= (cosh(K)-1)*Δ₁Θ[p,s] + Δ₂sin[p,⁺]
            # Δ₂f[coshsin[p],s] >= (cosh(K)-1)*Δ₁Θᵐ[p] + Δ₂sin[p,⁺]
            Δ₂f[coshsin[p],s] >= (cosh(K)-1)*Δ₁ρᵐ[p]*T/K + Δ₂sin[p,⁺]

            Δ₂f[sinhcos[p],s] >= sinh(K)*(1-cos(T))*Δ₁Θᵐ[p]/T + Δ₂sinh[p,s]
            Δ₂f[sinhcos[p],s] >= sinh(K)*(1-cos(T))*Δ₁ρ[p,s̄]/K + Δ₂sinh[p,s]
            # Δ₂f[sinhcos[p],s] >= sinh(K)*(1-cos(T))*Δ₁ρᵐ[p]/K + Δ₂sinh[p,s]

            Δ₂f[sinhsin[p],s] >= sinh(K)*Δ₁Θᵐ[p]
            Δ₂f[sinhsin[p],s] >= sinh(K)*Δ₁ρᵐ[p]*T/K
      end
    end

    if limitQQ > 0
        for k in 1:npv, s in ±
          s̄ = negate(s)
          @constraints m begin
            Qlim[k,s]/(v⁰[pv[k]]^2) >=
                sum(D⁺[k,j]*Δ₁u[j,s] + D⁻[k,j]*Δ₁u[j,s̄] for j=1:nu) +
                sum(E⁺[k,j]*Δ₂f[j,s] + E⁻[k,j]*Δ₂f[j,s̄] for j=1:nf)
          end
        end
    end

  end
end

export BaseLurieAdCertificate, create_A, create_BC
export create_box_R, create_loadability_R
export create_RX_branch, create_Imax, get_Qlims, cal_box, polytope_plot
export add_Δ₁y_constraints_LP, add_Δ₂f_constraints_LP

max_stars_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_issues_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_forks_repo_name: hunghtd/PowerFlowAnalysis-tutorial


text:
# Base 26 number where alphabets are numbers.

"""
```
    LabelNumeral{T<:Integer}
```
Wrapper around an `Integer` type that provides the following caabilities:

1. Prefix - like A-1, A-2 etc...
2. Lower case or upper case conversions
3. show and print options.
4. Mathematical operators like `+, -, <=, ==, >, isless, max and min`

The wrapped `struct` should implement the following methods:
```
1. T(::String)
2. T(::Int)
3. Base.hash(::T)
4. Base.convert{S <: Integer}(::Type{S}, num::T)  <-- to covert to standard numeral types
```
"""
mutable struct LabelNumeral{T<:Integer} <: Integer
    val::T
    prefix::String
    caselower::Bool
end

"""
```
    LabelNumeral{T <: Integer}(::T; prefix="", caselower=false)
    LabelNumeral{T <: Integer}(::Type{T}, i::Integer; prefix="", caselower=false)
    LabelNumeral{T <: Integer}(::Type{T}, s::String; prefix="", caselower=false)
```
Example:
```
julia> using RomanNumerals

julia> a = LabelNumeral(rn"XXIV"; prefix="A-", caselower=true)
A-xxiv

julia> a = LabelNumeral(rn"XXIV"; prefix="A-")
A-XXIV
```
Constructors for LabelNumeral
"""
LabelNumeral(t::T; prefix="", caselower=false) where T <: Integer =
    LabelNumeral(t, prefix, caselower)
LabelNumeral(t::Type{T}, i::Integer; prefix="", caselower=false) where T <: Integer =
    LabelNumeral(t(i), prefix, caselower)
LabelNumeral(t::Type{T}, str::String; prefix="", caselower=false) where T <: Integer =
    LabelNumeral(parse(T, str), prefix, caselower)
LabelNumeral{T}(i::S) where {T <: Integer, S <: Integer} = LabelNumeral(T, Int(i))

"""
```
    LabelNumeral(s::String; prefix="", caselower=false) --> LabelNumeral{Int}
```
Specialized constructor to return a LabelNumeral{Int} type.

Example:
```
julia> a = LabelNumeral("23"; prefix="A-", caselower=true)
A-23
```
The `prefix` does not get affected by the `caselower` parameter.
"""
LabelNumeral(str::String; prefix="", caselower=false) =
    LabelNumeral(parse(Int, str), prefix, caselower)

# Standard functions
# Conversion + promotion
Base.convert(::Type{Bool}, num::LabelNumeral) = true
Base.convert(::Type{BigInt}, num::LabelNumeral) = BigInt(Int(num))
Base.convert(::Type{T}, num::LabelNumeral) where T <: Integer = convert(T, num.val)
Base.convert(::Type{LabelNumeral{T}}, num::Int) where T <: Integer = LabelNumeral(T, num)
Base.promote_rule(::Type{LabelNumeral{S}}, ::Type{T}) where {T <: Integer, S <: Integer} = T
Base.Int(num::LabelNumeral) = convert(Int, num.val)

Base.string(num::LabelNumeral) = begin
    sval = num.val.str
    if num.caselower
        sval = lowercase(sval)
    end
    num.prefix*sval
end

# IO
Base.print(io::IO, num::LabelNumeral) = print(io, Base.string(num))
Base.show(io::IO, num::LabelNumeral) = write(io, Base.string(num))

Base.hash(num::LabelNumeral) = xor(hash(num.prefix), hash(num.val))

import Base: ==, isless, <=, <, >,
    +, -, max, min

# Equality operators
==(n1::LabelNumeral, n2::LabelNumeral) = Int(n1) == Int(n2)

# Comparisons
isless(n1::LabelNumeral, n2::LabelNumeral) = Int(n1) < Int(n2)
<(n1::LabelNumeral, n2::LabelNumeral) = Int(n1) < Int(n2)
>(n1::LabelNumeral, n2::LabelNumeral) = Int(n1) > Int(n2)
<=(n1::LabelNumeral, n2::LabelNumeral) = Int(n1) <= Int(n2)

## Arithmetic
# Multiple argument operators
for op in [:+, :-,:max, :min]
    @eval ($op)(n1::LabelNumeral{T}, n2::LabelNumeral{T},
        ns::LabelNumeral{T}...) where {T <: Integer} =
        $(op)(Int(n1), Int(n2), map(Int, ns)...) |> LabelNumeral{T}
end


text:
"""
    Colorbar(range, tick, scale, margin, colors)

Return a `Colorbar` object containing the data that defines the colorbar
associated to a plot.

The fields contained in a `Legend` object are:

* **`range`**: a Tuple{Float64, Float64} with the range of the color scale
    represented in the bar.
* **`tick`**: a `Float64` that gives the distance between tick marks drawn as
    guide next to the bar.
* **`scale`**: an `Int` code used by `GR.setscale` to define the
    scale of the bar (basically, if it is presented as linear or log scale).
* **`margin`**: a `Float64` with the size of the extra margin between the main
    plot frame and the bar.
* **`colors`**: an `Int` indicating the number of different grades in the color scale.
"""
struct Colorbar
    range::Tuple{Float64, Float64}
    tick::Float64
    scale::Int
    margin::Float64
    colors::Int
end

const EMPTYCOLORBAR = Colorbar(NULLPAIR, 0.0, 0, 0.0, 0)
Colorbar() = EMPTYCOLORBAR

"""
    Colorbar(axes::Axes [, colors=256])

Return a `Colorbar` defined by an `Axes` object that is used to calculate
its different properties, depending on the kind of axis and the range of the
`c` axis. If the `c` axis is not defined, this will return an empty `Colorbar`.
"""
function Colorbar(axes, colors=256)
    range = axes.ranges[:c]
    if !all(isfinite.(range))
        return EMPTYCOLORBAR
    end
    # Set the colorbar scale taking into account the scale of the main axes
    axscale = get(axes.options, :scale, 0)
    tick = get(axes.options, :clog, false) ? 2.0 : 0.5 * GR.tick(range...)
    if get(axes.options, :yflip, false)
        scale = axscale & ~GR.OPTION_FLIP_Y & ~GR.OPTION_FLIP_X
    else
        scale = axscale & ~GR.OPTION_FLIP_X
    end
    # Increase margin in 3-D axes (except if seen as a 2-D plane) and in polar axes
    if all(axes.perspective .≠ 0) || axes.kind == :polar
        margin = 0.05
    else
        margin = 0.0
    end
    Colorbar(range, tick, scale, margin, colors)
end

####################
## `draw` methods ##
####################

function draw(cb::Colorbar, crange=cb.range)
    cb == EMPTYCOLORBAR && return nothing
    zmin, zmax = crange
    # Redefine ticks if the range does not coincide with the default
    tick = (crange == cb.range) ? cb.tick : 0.5 * GR.tick(zmin, zmax)
    mainvp = GR.inqviewport()
    _, charheight = _tickcharheight(mainvp)
    GR.savestate()
    GR.setviewport(mainvp[2] + 0.02 + cb.margin, mainvp[2] + 0.05 + cb.margin, mainvp[3], mainvp[4])
    GR.setscale(cb.scale)
    GR.setwindow(0, 1, zmin, zmax)
    # Draw grade of colors
    l = round.(Int32, range(1000, stop=1255, length=cb.colors))
    h = 0.5 * (zmax - zmin) / (cb.colors - 1)
    GR.cellarray(0, 1, zmax, zmin, 1, cb.colors, l)
    GR.cellarray(0, 1, zmax + h, zmin - h, 1, cb.colors, l)
    # Draw ruled box
    GR.setlinecolorind(1)
    GR.setcharheight(charheight)
    GR.axes(0, tick, 1, zmin, 0, 1, 0.005)
    GR.restorestate()
    return nothing
end


text:
@marginalrule BIFM(:in_zprev_znext) (m_out::MultivariateNormalDistributionsFamily, m_in::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, m_znext::MultivariateNormalDistributionsFamily, meta::BIFMMeta{Float64}) = begin
    # Note: this rules is equal to the rule for the :in_zprev marginal. However, it yield identical performance with respect to the RTS smoother (aside from the additional prior on the last z)

    # fetch statistics from meta data
    A = getA(meta)
    B = getB(meta)
    C = getC(meta)
    ξ_ztilde = getξztilde(meta)
    Λ_ztilde = getΛztilde(meta)

    # # extract parameters from messages
    ξ_zprev_marginal, Λ_zprev_marginal = weightedmean_precision(m_zprev)
    ξ_in, Λ_in       = weightedmean_precision(m_in)
    ξ_out, Λ_out     = weightedmean_precision(m_out)

    # calculate outgoing message of zprev
    ξ_zprev_message = A' * ξ_ztilde
    Λ_zprev_message = A' * Λ_ztilde * A

    # calculate message towards znext from y
    dist1 = MvNormalWeightedMeanPrecision(C'*ξ_out, C'*Λ_out*C)

    # calculate message from z towards the addition node
    dist2 = prod(ProdAnalytical(), dist1, m_znext)
    ξ2, Λ2 = weightedmean_precision(dist2)

    # calculate joint message from the addition node
    BA = hcat(B, A)'
    ξ3 = BA*ξ2
    Λ3 = BA*Λ2*BA'

    # create a joint message of the input messages
    ξ4 = vcat(ξ_in, (ξ_zprev_marginal - ξ_zprev_message))
    Λ4 = cat(Λ_in, (Λ_zprev_marginal - Λ_zprev_message); dims=(1,2))

    # Actual return type depends on meta object as well, so we explicitly cast the result here
    # Should be noop if type matches
    T = promote_type(eltype(m_out), eltype(m_in), eltype(m_zprev), eltype(m_znext))

    # return joint marginal
    left  = convert(MvNormalWeightedMeanPrecision{T}, MvNormalWeightedMeanPrecision(ξ3, Λ3))
    right = convert(MvNormalWeightedMeanPrecision{T}, MvNormalWeightedMeanPrecision(ξ4, Λ4))

    return prod(ProdAnalytical(), left, right)

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GLPKBuilder"
version = v"4.64"

# Collection of sources required to build GLPKBuilder
sources = [
    "http://ftpmirror.gnu.org/gnu/glpk/glpk-4.64.tar.gz" =>
    "539267f40ea3e09c3b76a31c8747f559e8a097ec0cda8f1a3778eec3e4c3cc24",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd glpk-4.64/
cd $WORKSPACE/srcdir
cd glpk-4.64/
export LDFLAGS="-L${prefix}/lib"
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then export CPPFLAGS="-I${prefix}/include -D__WOE__=1";   else export CPPFLAGS="-I${prefix}/include"; fi
./configure --prefix=$prefix  --host=${target} --with-gmp
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl, :eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libglpk", :libglpk)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: juan-pablo-vielma/GLPKBuilder

max_issues_repo_name: juan-pablo-vielma/GLPKBuilder

max_forks_repo_name: juan-pablo-vielma/GLPKBuilder