text:
struct CutMesh
    mesh::Any
    cellsign::Any
    cutmeshnodeids::Any
    ncells::Any
    numnodes::Any
    function CutMesh(mesh, cellsign::Vector{Int}, cutmeshnodeids::Matrix{Int})
        ncells = number_of_cells(mesh)
        nummeshnodes = number_of_nodes(mesh)
        @assert length(cellsign) == ncells
        @assert size(cutmeshnodeids) == (2, nummeshnodes)

        numnodes = maximum(cutmeshnodeids)
        new(mesh, cellsign, cutmeshnodeids, ncells, numnodes)
    end
end

function Base.show(io::IO, cutmesh::CutMesh)
    dim = dimension(cutmesh)
    x0 = reference_corner(cutmesh)
    meshwidths = mesh_widths(cutmesh)
    numnodes = number_of_nodes(cutmesh)
    ncells = number_of_cells(cutmesh)
    str =
        "CutMesh\n\tDimension : $dim\n\tCorner : $x0\n\tWidth : $meshwidths" *
        "\n\tNum. Nodes : $numnodes\n\tNum. Cells : $ncells"
    print(io, str)
end

function CutMesh(mesh, levelset; tol = 1e-4, perturbation = 1e-2)
    dx = minimum(element_size(mesh))
    cellsign = cell_sign!(levelset, tol, perturbation * dx)

    posactivenodeids = active_node_ids(mesh, +1, cellsign)
    negactivenodeids = active_node_ids(mesh, -1, cellsign)

    totalnumnodes = number_of_nodes(mesh)
    cutmeshnodeids =
        cut_mesh_nodeids(posactivenodeids, negactivenodeids, totalnumnodes)
    return CutMesh(mesh, cellsign, cutmeshnodeids)
end

function cell_sign(cutmesh::CutMesh, cellid)
    return cutmesh.cellsign[cellid]
end

function dimension(cutmesh::CutMesh)
    return dimension(background_mesh(cutmesh))
end

function number_of_nodes(cutmesh::CutMesh)
    return cutmesh.numnodes
end

function reference_corner(cutmesh::CutMesh)
    return reference_corner(background_mesh(cutmesh))
end

function number_of_cells(cutmesh::CutMesh)
    return cutmesh.ncells
end

function background_mesh(cutmesh::CutMesh)
    return cutmesh.mesh
end

function cell_connectivity(cutmesh::CutMesh, faceid, cellid)
    return cell_connectivity(background_mesh(cutmesh), faceid, cellid)
end

function cell_map(cutmesh::CutMesh, cellid)
    return cell_map(background_mesh(cutmesh), cellid)
end

function cell_map(cutmesh::CutMesh, cellsign, cellid)
    return cell_map(cutmesh, cellid)
end

function solution_cell_id(cutmesh::CutMesh, s, cellid)
    return cellid
end

function nodal_connectivity(cutmesh::CutMesh, s, cellid)
    row = cell_sign_to_row(s)
    ncells = cutmesh.ncells
    @assert 1 <= cellid <= ncells
    cs = cell_sign(cutmesh, cellid)
    @assert cs == s ||
            cs == 0 ||
            error(
                "Requested nodeids for sign $s not consistent with cell sign $cs for cellid $cellid",
            )

    ids = nodal_connectivity(background_mesh(cutmesh), cellid)
    nodeids = cutmesh.cutmeshnodeids[row, ids]
    return nodeids
end

function nodal_coordinates(cutmesh::CutMesh)
    return nodal_coordinates(background_mesh(cutmesh))
end

function perturb_lagrange_levelset_coefficients!(levelset, cellid, perturbation)
    newcoeffs = coefficients(levelset, cellid) .+ perturbation
    update_coefficients!(levelset, cellid, newcoeffs)
    load_coefficients!(levelset,cellid)
end

function perturb_hermite_levelset_coefficients!(levelset, cellid, perturbation)
    dim = dimension(levelset)
    interpidx = hermite_interpolation_indices(dim)
    cellcoeffs = coefficients(levelset, cellid)
    perturbvec = zeros(length(cellcoeffs))
    perturbvec[interpidx] .= perturbation
    newcoeffs = cellcoeffs + perturbvec
    update_coefficients!(levelset, cellid, newcoeffs)
    load_coefficients!(levelset,cellid)
end

function cell_sign!(levelset, tol, perturbation)
    ncells = number_of_cells(background_mesh(levelset))
    cellsign = zeros(Int, ncells)
    xL, xR = [-1.0, -1.0], [1.0, 1.0]
    basistype = typeof(basis(interpolater(levelset)))

    for cellid = 1:ncells
        load_coefficients!(levelset, cellid)

        s = sign(interpolater(levelset), xL, xR, tol = tol)
        if (s == +1 || s == 0 || s == -1)
            cellsign[cellid] = s
        else
            # @warn "Perturbing levelset function by perturbation = $perturbation"
            # newcoeffs = coefficients(levelset, cellid) .+ perturbation
            # update_coefficients!(levelset, cellid, newcoeffs)
            # load_coefficients!(levelset, cellid)

            if basistype <: LagrangeTensorProductBasis
                perturb_lagrange_levelset_coefficients!(
                    levelset,
                    cellid,
                    perturbation,
                )
            elseif basistype <: HermiteTensorProductBasis
                perturb_hermite_levelset_coefficients!(
                    levelset,
                    cellid,
                    perturbation,
                )
            else
                error("Unexpected basistype $basistype")
            end

            s = sign(interpolater(levelset), xL, xR, tol = tol)
            if (s == +1 || s == 0 || s == -1)
                cellsign[cellid] = s
            else
                error(
                    "Could not determine cell sign after perturbation = $perturbation",
                )
            end
        end
    end
    return cellsign
end

function active_node_ids(mesh, s, cellsign)
    @assert (s == -1 || s == +1)
    activenodeids = Int[]
    numcells = number_of_cells(mesh)
    @assert length(cellsign) == numcells

    for cellid = 1:numcells
        if cellsign[cellid] == s || cellsign[cellid] == 0
            append!(activenodeids, nodal_connectivity(mesh, cellid))
        end
    end
    return activenodeids
end

function cut_mesh_nodeids(posactivenodeids, negactivenodeids, totalnumnodes)
    cutmeshnodeids = zeros(Int, 2, totalnumnodes)
    counter = 1
    for nodeid in posactivenodeids
        cutmeshnodeids[1, nodeid] = counter
        counter += 1
    end
    for nodeid in negactivenodeids
        cutmeshnodeids[2, nodeid] = counter
        counter += 1
    end
    return cutmeshnodeids
end


text:
using DiffGPFlux;
using Random: randperm;
using Plots;
using Measurements;

f(x) = 2.0*x-3.0
X = reshape(Array(range(0.0, 10.0, length=100)), 1, 100)
y = dropdims(f.(X), dims=1)+0.1*randn(100)
index = randperm(100)
train_index = index[1:20]
test_index = index[21:100]
train_X, train_y = X[:, train_index], y[train_index]
test_X, test_y = X[:, test_index], y[test_index]

θ₀ = seed_duals(rand(2), Float64, 3)
# θ₀ = rand(2)
γ₀ = [θ₀[1], θ₀[2]]
β = 0.1
# c = [θ₀[3]]
# μ = ConstMean(c)
# ∇ₓμ(x) = zero(x)
GP = GaussProcess(γ₀, β, linear, ∇ₓlinear)
θ = [GP.γ]
indexes = [[1, 2]]

opt = ADAM(0.01)

plt = plot(1, marker=2)
@gif for i in 1:350
    nll = negloglik(GP, train_X, train_y)
    push!(plt, nll.value)
    θ̄_array = [g for g in nll.partials]
    θ̄ = grads(θ̄_array, indexes)
    update!(opt, θ, θ̄, 3)
    end every 5

display(GP.γ)
# display(GP.μ.c)

γ_values = [GP.γ[1].value, GP.γ[2].value]
# c_values = [GP.μ.c[1].value]
# non_dual_μ = ConstMean(c_values)
γ_values = [1.0, 1.0]
non_dual_GP = GaussProcess(γ_values, β, linear, ∇ₓlinear)
y_predict, σ_predict = predict(non_dual_GP, train_X, train_y, X)
plot(train_X', train_y, seriestype=:scatter)
plot!(X', y_predict.±σ_predict, seriestype=:scatter)
plot!(X', y, color=:red, lw=3)

ext: jl

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
### Linear Interpolation
struct LinearInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  LinearInterpolation{FT}(u,t) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t)
end

function LinearInterpolation(u,t)
  u, t = munge_data(u, t)
  LinearInterpolation{true}(u,t)
end

### Quadratic Interpolation
struct QuadraticInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  QuadraticInterpolation{FT}(u,t) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t)
end

function QuadraticInterpolation(u,t)
  u, t = munge_data(u, t)
  QuadraticInterpolation{true}(u,t)
end

### Lagrange Interpolation
struct LagrangeInterpolation{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  n::Int
  LagrangeInterpolation{FT}(u,t,n) where FT = new{typeof(u),typeof(t),FT,eltype(u)}(u,t,n)
end

function LagrangeInterpolation(u,t,n=nothing)
  u, t = munge_data(u, t)
  if n == nothing
    n = length(t) - 1 # degree
  end
  if n != length(t) - 1
    error("Currently only n=length(t) - 1 is supported")
  end
  LagrangeInterpolation{true}(u,t,n)
end

### Akima Interpolation
struct AkimaInterpolation{uType,tType,bType,cType,dType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  b::bType
  c::cType
  d::dType
  AkimaInterpolation{FT}(u,t,b,c,d) where FT = new{typeof(u),typeof(t),typeof(b),typeof(c),
                                                     typeof(d),FT,eltype(u)}(u,t,b,c,d)
end

function AkimaInterpolation(u, t)
  u, t = munge_data(u, t)
  n = length(t)
  dt = diff(t)
  m = Array{eltype(u)}(undef, n + 3)
  m[3:end-2] = diff(u) ./ dt
  m[2] = 2m[3] - m[4]
  m[1] = 2m[2] - m[3]
  m[end-1] = 2m[end-2] - m[end-3]
  m[end] = 2m[end-1] - m[end-2]

  b = 0.5 .* (m[4:end] .+ m[1:end-3])
  dm = abs.(diff(m))
  f1 = dm[3:n+2]
  f2 = dm[1:n]
  f12 = f1 + f2
  ind = findall(f12 .> 1e-9 * maximum(f12))
  b[ind] = (f1[ind] .* m[ind .+ 1] .+
            f2[ind] .* m[ind .+ 2]) ./ f12[ind]
  c = (3.0 .* m[3:end-2] .- 2.0 .* b[1:end-1] .- b[2:end]) ./ dt
  d = (b[1:end-1] .+ b[2:end] .- 2.0 .* m[3:end-2]) ./ dt.^2

  AkimaInterpolation{true}(u,t,b,c,d)
end

### ConstantInterpolation Interpolation
struct ConstantInterpolation{uType,tType,dirType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  dir::Symbol # indicates if value to the $dir should be used for the interpolation
  ConstantInterpolation{FT}(u,t,dir) where FT = new{typeof(u),typeof(t),typeof(dir),FT,eltype(u)}(u,t,dir)
end

function ConstantInterpolation(u,t;dir=:left)
  u, t = munge_data(u, t)
  ConstantInterpolation{true}(u,t,dir)
end

Base.@deprecate_binding ZeroSpline ConstantInterpolation

### QuadraticSpline Interpolation
struct QuadraticSpline{uType,tType,tAType,dType,zType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  tA::tAType
  d::dType
  z::zType
  QuadraticSpline{FT}(u,t,tA,d,z) where FT = new{typeof(u),typeof(t),typeof(tA),
                                                  typeof(d),typeof(z),FT,eltype(u)}(u,t,tA,d,z)
end

function QuadraticSpline(u,t)
  u, t = munge_data(u, t)
  s = length(t)
  dl = ones(eltype(t),s-1)
  d_tmp = ones(eltype(t),s)
  du = zeros(eltype(t),s-1)
  tA = Tridiagonal(dl,d_tmp,du)
  d = map(i -> i == 1 ? 0 : 2//1 * (u[i] - u[i-1])/(t[i] - t[i-1]), 1:s)
  z = tA\d
  QuadraticSpline{true}(u,t,tA,d,z)
end

# Cubic Spline Interpolation
struct CubicSpline{uType,tType,hType,zType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  h::hType
  z::zType
  CubicSpline{FT}(u,t,h,z) where FT = new{typeof(u),typeof(t),typeof(h),typeof(z),FT,eltype(u)}(u,t,h,z)
end

function CubicSpline(u,t)
  u, t = munge_data(u, t)
  n = length(t) - 1
  h = vcat(0, map(k -> t[k+1] - t[k], 1:length(t)-1), 0)
  dl = h[2:n+1]
  d_tmp = 2 .* (h[1:n+1] .+ h[2:n+2])
  du = h[2:n+1]
  tA = LinearAlgebra.Tridiagonal(dl,d_tmp,du)
  d = map(i -> i == 1 || i == n + 1 ? 0 : 6(u[i+1] - u[i]) / h[i+1] - 6(u[i] - u[i-1]) / h[i], 1:n+1)
  z = tA\d
  CubicSpline{true}(u,t,h[1:n+1],z)
end

### BSpline Curve Interpolation
struct BSplineInterpolation{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  d::Int    # degree
  p::pType  # params vector
  k::kType  # knot vector
  c::cType  # control points
  pVecType::Symbol
  knotVecType::Symbol
  BSplineInterpolation{FT}(u,t,d,p,k,c,pVecType,knotVecType) where FT =  new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,p,k,c,pVecType,knotVecType)
end

function BSplineInterpolation(u,t,d,pVecType,knotVecType)
  u, t = munge_data(u, t)
  n = length(t)
  s = zero(eltype(u))
  p = zero(t)
  k = zeros(eltype(t),n+d+1)
  l = zeros(eltype(u),n-1)
  p[1] = zero(eltype(t))
  p[end] = one(eltype(t))

  for i = 2:n
    s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2)
    l[i-1] = s
  end
  if pVecType == :Uniform
    for i = 2:(n-1)
      p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1)
    end
  elseif pVecType == :ArcLen
    for i = 2:(n-1)
      p[i] = p[1] + l[i-1]/s * (p[end]-p[1])
    end
  end

  lidx = 1
  ridx = length(k)
  while lidx <= (d+1) && ridx >= (length(k)-d)
    k[lidx] = p[1]
    k[ridx] = p[end]
    lidx += 1
    ridx -= 1
  end

  ps = zeros(eltype(t),n-2)
  s = zero(eltype(t))
  for i = 2:n-1
    s += p[i]
    ps[i-1] = s
  end

  if knotVecType == :Uniform
    # uniformly spaced knot vector
    # this method is not recommended because, if it is used with the chord length method for global interpolation,
    # the system of linear equations would be singular.
    for i = (d+2):n
      k[i] = k[1] + (i-d-1)//(n-d) * (k[end]-k[1])
    end
  elseif knotVecType == :Average
    # average spaced knot vector
    idx = 1
    if d+2 <= n
      k[d+2] = 1//d * ps[d]
    end
    for i = (d+3):n
      k[i] = 1//d * (ps[idx+d] - ps[idx])
      idx += 1
    end
  end
  # control points
  N = spline_coefficients(n,d,k,p)
  c = vec(N\u[:,:])
  BSplineInterpolation{true}(u,t,d,p,k,c,pVecType,knotVecType)
end

### BSpline Curve Approx
struct BSplineApprox{uType,tType,pType,kType,cType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  d::Int    # degree
  h::Int    # number of control points (n => h >= d >= 1)
  p::pType  # params vector
  k::kType  # knot vector
  c::cType  # control points
  pVecType::Symbol
  knotVecType::Symbol
  BSplineApprox{FT}(u,t,d,h,p,k,c,pVecType,knotVecType) where FT =  new{typeof(u),typeof(t),typeof(p),typeof(k),typeof(c),FT,eltype(u)}(u,t,d,h,p,k,c,pVecType,knotVecType)
end

function BSplineApprox(u,t,d,h,pVecType,knotVecType)
  u, t = munge_data(u, t)
  n = length(t)
  s = zero(eltype(u))
  p = zero(t)
  k = zeros(eltype(t),h+d+1)
  l = zeros(eltype(u),n-1)
  p[1] = zero(eltype(t))
  p[end] = one(eltype(t))

  for i = 2:n
    s += √((t[i] - t[i-1])^2 + (u[i] - u[i-1])^2)
    l[i-1] = s
  end
  if pVecType == :Uniform
    for i = 2:(n-1)
      p[i] = p[1] + (i-1)*(p[end]-p[1])/(n-1)
    end
  elseif pVecType == :ArcLen
    for i = 2:(n-1)
      p[i] = p[1] + l[i-1]/s * (p[end]-p[1])
    end
  end

  lidx = 1
  ridx = length(k)
  while lidx <= (d+1) && ridx >= (length(k)-d)
    k[lidx] = p[1]
    k[ridx] = p[end]
    lidx += 1
    ridx -= 1
  end

  ps = zeros(eltype(t),n-2)
  s = zero(eltype(t))
  for i = 2:n-1
    s += p[i]
    ps[i-1] = s
  end

  if knotVecType == :Uniform
    # uniformly spaced knot vector
    # this method is not recommended because, if it is used with the chord length method for global interpolation,
    # the system of linear equations would be singular.
    for i = (d+2):h
      k[i] = k[1] + (i-d-1)//(h-d) * (k[end]-k[1])
    end
  elseif knotVecType == :Average
    # NOTE: verify that average method can be applied when size of k is less than size of p
    # average spaced knot vector
    idx = 1
    if d+2 <= h
      k[d+2] = 1//d * ps[d]
    end
    for i = (d+3):h
      k[i] = 1//d * (ps[idx+d] - ps[idx])
      idx += 1
    end
  end
  # control points
  c = zeros(eltype(u), h)
  c[1] = u[1]
  c[end] = u[end]
  q = zeros(eltype(u), n)
  N = zeros(eltype(t), n, h)
  for i = 1:n
    N[i, :] .= spline_coefficients(h,d,k,p[i])
  end
  for k = 2:n-1
    q[k] = u[k] - N[k, 1]*u[1] - N[k, h]*u[end]
  end
  Q = Matrix{eltype(u)}(undef, h-2, 1)
  for i = 2:h-1
    s = 0.0
    for k = 2:n-1
      s += N[k, i]*q[k]
    end
    Q[i-1] = s
  end
  N = N[2:end-1,2:h-1]
  M = transpose(N) * N
  P = M\Q
  c[2:end-1] .= vec(P)
  BSplineApprox{true}(u,t,d,h,p,k,c,pVecType,knotVecType)
end

### Curvefit
struct CurvefitCache{uType,tType,mType,p0Type,ubType,lbType,algType,pminType,FT,T} <: AbstractInterpolation{FT,T}
  u::uType
  t::tType
  m::mType        # model type
  p0::p0Type      # intial params
  ub::ubType      # upper bound of params
  lb::lbType      # lower bound of params
  alg::algType    # alg to optimize cost function
  pmin::pminType  # optimized params
  CurvefitCache{FT}(u,t,m,p0,ub,lb,alg,pmin) where FT = new{typeof(u),typeof(t),typeof(m),
                                      typeof(p0),typeof(ub),typeof(lb),
                                      typeof(alg),typeof(pmin),FT,eltype(u)}(u,t,m,p0,ub,lb,alg,pmin)
end

function Curvefit(u, t, model, p0, alg, box=false, lb=nothing, ub=nothing)
  u, t = munge_data(u, t)
  errfun(t,u,p) = sum(abs2.(u .- model(t,p)))
  if box == false
      mfit = optimize(p->errfun(t,u,p), p0, alg)
  else
      if lb == nothing || ub == nothing
          error("lower or upper bound should not be nothing")
      end
      od = OnceDifferentiable(p->errfun(t,u,p), p0, autodiff=:finite)
      mfit = optimize(od, lb, ub, p0, Fminbox(alg))
  end
  pmin = Optim.minimizer(mfit)
  CurvefitCache{true}(u,t,model,p0,ub,lb,alg,pmin)
end


text:
module OnlineCCA

using DelimitedFiles:
    writedlm, readdlm
using Statistics:
    mean, var, cor
using LinearAlgebra:
    Diagonal, lu!, qr!, svd, dot, norm, eigvecs, diag
using Random:
    randperm
using ProgressMeter:
  Progress, next!
using ArgParse:
    ArgParseSettings, parse_args, @add_arg_table
using StatsBase:
  percentile
using DataFrames:
  DataFrame
using GLM:
  glm, coef, IdentityLink, @formula
using Distributions:
  Gamma, ccdf, Chisq
using CodecZstd:
  ZstdCompressorStream, ZstdDecompressorStream

export output, common_parse_commandline, csv2bin, sumr, filtering, sgd, gd, rsgd, horst, orthiter, oocmcca, nipals

include("Utils.jl")
include("csv2bin.jl")
include("sumr.jl")
include("sgd.jl")
include("gd.jl")
include("rsgd.jl")
include("orthiter.jl")
include("horst.jl")
include("oocmcca.jl")
include("nipals.jl")

end


text:
"""
    Crash

Collection of crash methods used to generate initial decisions in structured algorithms.

...
# Available crash methods
- [`None`](@ref)
- [`EVP`](@ref)
- [`Scenario`](@ref)
- [`Custom`](@ref)
...

## Examples

The following solves a stochastic program `sp` created in `StochasticPrograms.jl` using an L-shaped algorithm with trust-region and Clp as an `lpsolver` and by generating an initial decision with the `EVP` crash.

```jldoctest
julia> optimize!(sp, solver = LShapedSolver(GLPKSolverLP(), crash=Crash.EVP(), regularize = TrustRegion()))
L-Shaped Gap  Time: 0:00:00 (8 iterations)
  Objective:       -855.8333333333339
  Gap:             0.0
  Number of cuts:  4
  Iterations:      8
:Optimal
```
"""
module Crash

using StochasticPrograms
using StochasticPrograms: AbstractCrash
using MathOptInterface

const MOI = MathOptInterface

"""
    None

Randomize the initial decision (default).

"""
struct None <: AbstractCrash end

function (::None)(stochasticprogram::StochasticProgram)
    return rand(num_decisions(stochasticprogram))
end

function (::None)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 0; optimizer = optimizer)
    return rand(num_decisions(sp))
end

"""
    EVP

Solve the expected value problem corresponding to the stochastic program and use the expected value solution as initial decision.

"""
struct EVP <: AbstractCrash end

function (::EVP)(stochasticprogram::StochasticProgram)
    return expected_value_decision(stochasticprogram)
end

function (crash::EVP)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    FeasiblePoint

Generate a feasible first-stage decision as initial decision.

"""
struct FeasiblePoint <: AbstractCrash end

function (::FeasiblePoint)(stochasticprogram::StochasticProgram)
    # Generate first-stage
    m = stage_one_model(stochasticprogram, optimizer = master_optimizer(stochasticprogram))
    # Solve feasibility problem
    @objective(m, MOI.FEASIBILITY_SENSE, 0)
    optimize!(m)
    # Return feasible point
    return JuMP.value.(all_decision_variables(m, 1))
end

function (crash::FeasiblePoint)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    Scenario

Solve the wait-and-see problem corresponding a supplied scenario and use the optimal solution as initial decision.

"""
struct Scenario{S <: AbstractScenario} <: AbstractCrash
    scenario::S

    function Scenario(scenario::S) where S <: AbstractScenario
        return new{S}(scenario)
    end
end

function (crash::Scenario)(stochasticprogram::StochasticProgram)
    return wait_and_see_decision(stochasticprogram, crash.scenario)
end

function (crash::Scenario)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    Custom(x₀)

Use the user-supplied `x₀` as initial decision.

"""
struct Custom{T <: AbstractFloat} <: AbstractCrash
    x₀::Vector{T}

    function Custom(x₀::Vector{T}) where T <: AbstractFloat
        return new{T}(x₀)
    end
end

function (crash::Custom)(stochasticprogram::StochasticProgram)
    return crash.x₀[1:num_decisions(stochasticprogram)]
end

function (crash::Custom)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 0; optimizer = optimizer)
    return crash.x₀[1:num_decisions(sp)]
end

end


text:
"""
    eval_constraint(S::Array{Int,2}, c::CAcon)
    eval_constraint(S::Array{Int,2}, c::BRcon)
    eval_constraint(S::Array{Int,2}, c::GAcon)
    eval_constraint(S::Array{Int,2}, c::FAcon)
    eval_constraint(S::Array{Int,2}, c::SEcon)

evaluate the penalty for a constraint on a schedule.

# Arguments

- `S::Array{Int,2}`: The schedule to be evaluated on. Assuming teams `i` and `j`
    meet in slot `s` we have `S[i,s] == j`, if team `i` plays at home and
    `S[i,s] == -j` if team `i` plays away.
- `c::T where T <: Constraint`: The constraint to be evaluated.

# Returns

- `Int`: the offset of the constraint on the given schedule.
"""
function eval_constraint end

function eval_constraint(S::Array{Int,2}, c::CAcon)
    count = 0
    for i in c.teams1
        for s in c.slots
            match = S[i, s]
            if c.mode == MODE_HA || sign(match) == c.mode
                if abs(match) in c.teams2
                    count += 1
                end
            end
        end
    end
    return max(count - c.ub, 0)
end

function eval_constraint(S::Array{Int,2}, c::BRcon)
    count = 0
    for i in c.teams
        for s in c.slots
            s == 1 && continue
            if c.mode == MODE_HA || sign(S[i, s]) == c.mode
                if sign(S[i, s]) == sign(S[i, s-1])
                    count += 1
                end
            end
        end
    end
    return max(count - c.ub, 0)
end

function eval_constraint(S::Array{Int,2}, c::GAcon)
    count = 0
    for s in c.slots
        for i = 1:size(S, 1)
            if sign(S[i, s]) == 1 && (i, S[i, s]) in c.meetings
                count += 1
            end
        end
    end
    return max(max(c.lb - count, 0), max(count - c.ub, 0))
end

function eval_constraint(S::Array{Int,2}, c::FAcon)
    hgames = zeros(Int64, size(S, 1))
    maxdiff = zeros(Int64, size(S, 1), size(S, 1))
    lastslot = 0
    for s in c.slots
        for sl = lastslot+1:s
            for i = 1:size(S, 1)
                if sign(S[i, sl]) == 1
                    hgames[i] += 1
                end
            end
        end
        for i in c.teams
            for j in c.teams
                j <= i && continue
                if maxdiff[i, j] < abs(hgames[i] - hgames[j])
                    maxdiff[i, j] = abs(hgames[i] - hgames[j])
                end
            end
        end
        lastslot = s
    end
    count = 0
    for i in c.teams
        for j in c.teams
            j <= i && continue
            count += max(maxdiff[i, j] - c.ub, 0)
        end
    end
    return count
end

function eval_constraint(S::Array{Int,2}, c::SEcon)
    count = 0
    games = Array{Int64}(undef, size(S, 1), size(S, 1))
    for s = 1:size(S, 2)
        for i = 1:size(S, 1)
            if sign(S[i, s]) == 1
                games[i, S[i, s]] = s
            end
        end
    end
    for i in c.teams
        for j in c.teams
            j <= i && continue
            if abs(games[i, j] - games[j, i]) < c.lb
                count += c.lb - abs(games[i, j] - games[j, i]) + 1
            end
        end
    end
    return count
end

"""
    validate_pattern_constraint(rr2::RR2, pat::Int64, con::CAcon)
    validate_pattern_constraint(rr2::RR2, pat::Int64, con::BRcon)

Checks if a pattern constraint `con` is valid for a given pattern `pat`.
"""
function validate_pattern_constraint end

function validate_pattern_constraint(rr2::RR2, pat::Int64, con::CAcon)
    @assert is_pattern_constraint(rr2, con)

    if con.mode == MODE_H
        return count_ones(pat & con.slots.binrepr) <= con.ub
    else # mode_A
        return count_ones(flipbits(pat, 1, rr2.nslots) & con.slots.binrepr) <=
               con.ub
    end
end

function validate_pattern_constraint(rr2::RR2, pat::Int64, con::BRcon)
    @assert is_pattern_constraint(rr2, con)

    breakcount = 0
    if con.mode == MODE_H || con.mode == MODE_HA
        breakcount += count_ones((pat & pat << 1) & con.slots.binrepr)
    end
    if con.mode == MODE_A || con.mode == MODE_HA
        pat = flipbits(pat, 1, rr2.nslots)
        breakcount += count_ones((pat & pat << 1) & con.slots.binrepr)
    end
    return breakcount <= con.ub
end

"""
Groups `cons` into pattern and patternset constraints and adds them to
`patterncons` or `patternSetCons` respectivly.
"""
function group_constraints!(
    patternCons::Vector{Vector{Constraint}},
    patternSetCons::Vector{Constraint},
    cons::Vector{BRcon},
)
    for brc in cons
        if is_pattern_constraint(brc)
            team = trailing_zeros(brc.teams.binrepr) + 1
            push!(patternCons[team], brc)
        else
            push!(patternSetCons, brc)
        end
    end
end

"""
Groups `cons` into pattern, patternset and general constraints and adds them to
`patterncons`, `patternSetCons` or `genCons` respectivly.
"""
function group_constraints!(
    patternCons::Vector{Vector{Constraint}},
    patternSetCons::Vector{Constraint},
    genCons::Vector{Constraint},
    rr2::RR2,
    cons::Vector{CAcon},
)
    for cac in cons
        if is_pattern_constraint(rr2, cac)
            team = trailing_zeros(cac.teams1.binrepr) + 1
            push!(patternCons[team], cac)
        elseif is_pattern_set_constraint(rr2, cac)
            push!(patternSetCons, cac)
        else
            push!(genCons, cac)
        end
    end
end

"""
Adds `cons` to `genCons` and adds a relaxation into two capacity constraints to
`patternCons` or `patternSetCons` respectively.
"""
function group_constraints!(
    patternCons::Vector{Vector{Constraint}},
    patternSetCons::Vector{Constraint},
    genCons::Vector{Constraint},
    rr2::RR2,
    cons::Vector{GAcon},
)
    for gac in cons
        # special case in which GA con can be relaxed to two CA cons
        if gac.lb > 0
            teams1, teams2 = unique.(zip(gac.meetings...))

            # teams1 must play at least lb home games in slots
            cacon1 = CAcon(
                IdxSet(teams1),
                IdxSet(2^(rr2.nteams) - 1),
                gac.slots,
                length(teams1) * length(gac.slots) - gac.lb,
                MODE_A,
                0,
            )
            # teams2 must play at least lb away games in slots
            cacon2 = CAcon(
                IdxSet(teams2),
                IdxSet(2^(rr2.nteams) - 1),
                gac.slots,
                length(teams2) * length(gac.slots) - gac.lb,
                MODE_H,
                0,
            )
            length(teams1) == 1 ? push!(patternCons[only(teams1)], cacon1) :
                push!(patternSetCons, cacon1)
            length(teams2) == 1 ? push!(patternCons[only(teams2)], cacon2) :
                push!(patternSetCons, cacon2)
        end
        push!(genCons, gac)
    end
end

"""
    group_constraints(inst::Instance)

Groups all constraints of `inst` into pattern constraints, pattern-set
constraints and general constraints.
"""
function group_constraints(inst::Instance; f = x -> true)
    nteams = inst.rr2.nteams
    patternCons = [Vector{Constraint}() for i = 1:nteams]
    patternSetCons = Vector{Constraint}()
    genCons = Vector{Constraint}()
    group_constraints!(
        patternCons,
        patternSetCons,
        filter(f, inst.constrBR),
    )
    group_constraints!(
        patternCons,
        patternSetCons,
        genCons,
        inst.rr2,
        filter(f, inst.constrCA),
    )
    group_constraints!(
        patternCons,
        patternSetCons,
        genCons,
        inst.rr2,
        filter(f, inst.constrGA),
    )
    return patternCons, patternSetCons, genCons
end

"""
    conentries(rr2::RR2, cons::Vector{Constraint})

Generates a `BitArray{3}` where the index `(c,i,s)` indicates whether `cons[c]`
constrains games played by team `i` in slot `s`.
"""
function conentries(rr2::RR2, cons::Vector{Constraint})
    A = falses(rr2.nteams, rr2.nslots, length(cons))
    for (i, c) in enumerate(cons)
        constrained_entries!(view(A, :, :, i), rr2, c)
    end
    return A
end

"""
    constrained_entries(M::AbstractArray, rr2::RR2, con::T) where T <:Constraint

Compute all (team, slot) entries that are constrained by `con` and stores them
in `M` where M is either a `BitMatrix` or a `SubArray` of a BitMatrix.
"""
function constrained_entries! end

@inline function constrained_entries!(
    M::AbstractArray,
    rr2::RR2,
    con::Union{CAcon,GAcon},
)
    for i in teams(con)
        for s in slots(con)
            M[i, s] = 1
        end
    end
end

@inline function constrained_entries!(M::AbstractArray, rr2::RR2, con::BRcon)
    for i in teams(con)
        for s in slots(con)
            M[i, s] = 1
            if s > 1
                M[i, s-1] = 1
            end
        end
    end
end

@inline function constrained_entries!(M::AbstractArray, rr2::RR2, con::FAcon)
    maxslot = max(slots(con)...)
    for i in teams(con)
        M[i, 1:maxslot] .= 1
    end
end

@inline function constrained_entries!(M::AbstractArray, rr2::RR2, con::SEcon)
    for i in teams(con)
        M[i, :] .= 1
    end
end


text:
@testset "discreteshape" begin
    img = generate_surface(SynthSphere(), 0.5, [0,0,1], radius = 50, background = false)
    Z, p, q = retrieve_surface(DiscreteShape(), img, smoothness=1000)
    @test Z[38,38] ≈ 0.0003182796154911338
    @test Z[1,1] ≈ 2.0548315482140964e-6
    @test isapprox(Z[76,76], 1.2021784812891086e-19; atol = 2e-15)
    @test Z[45,45] ≈ 0.00040018054317404224
    @test Z[50,18] ≈ 0.00014000747445970892
    @test Z[11, 68] ≈ 0.00018127096811786046
    @test Z[42,42] ≈ 0.0003883024349324123
end


text:
# --------------------------------------------------------- #
# Stuff to construct basis matrices of complete polynomials #
# --------------------------------------------------------- #

using Base.Cartesian: @nloops, @nexprs

"""
Construct basis matrix for complete polynomial of degree `d`, given
input data `z`. `z` is assumed to be the degree 1 realization of each
variable. For example, if variables are `q`, `r`, and `s`, then `z`
should be `z = [q r s]`

Output is a basis matrix. In our example, with `d` set to 2 we would have

TODO: update docstring to properly give order of terms

```julia
out = [ones(size(z,1)) q r s q.*r q.*s r.*s q.^2 q.*r q.*s r.^2 r.*s s.^2]
```
"""
:complete_polynomial

immutable Degree{N} end

function n_complete(n::Int, D::Int)
    out = 1
    for d=1:D
        tmp = 1
        for j=0:d-1
            tmp *= (n+j)
        end
        out += div(tmp, factorial(d))
    end
    out
end

@generated function complete_polynomial!{N}(z::Matrix, d::Degree{N},
                                            out::Matrix)
    complete_polynomial_impl!(z, d, out)
end

function complete_polynomial_impl!{T,N}(z::Type{Matrix{T}}, ::Type{Degree{N}},
                                        out::Type{Matrix{T}})
    quote
        nobs, nvar = size(z)
        if size(out) != (nobs, n_complete(nvar, $N))
            error("z, out not compatible")
        end

        # reset first column to ones
        @inbounds for i=1:nobs
            out[i, 1] = one($T)
        end

        # set next nvar columns to input matrix
        @inbounds for n=2:nvar+1, i=1:nobs
            out[i, n] = z[i, n-1]
        end

        ix = nvar+1
        @nloops($N, # number of loops
                i,  # counter
                d->((d == $N ? 1 : i_{d+1}) : nvar),  # ranges
                d->(d == $N ? nothing :
                    (begin
                        ix += 1
                        @inbounds @simd for r=1:nobs
                            tmp = one($T)
                            @nexprs $N-d+1 j->(tmp *= z[r, i_{$N-j+1}])
                            out[r, ix]=tmp
                        end
                    end)),  # preexpr
                Expr(:block, :nothing)  # bodyexpr
                )
        out
    end
end

function complete_polynomial{T}(z::Matrix{T}, d::Int)
    nobs, nvar = size(z)
    out = Array(T, nobs, n_complete(nvar, d))
    complete_polynomial!(z, Degree{d}(), out)::Matrix{T}
end

function complete_polynomial!{T}(z::Matrix{T}, d::Int, out::Matrix{T})
    complete_polynomial!(z, Degree{d}(), out)::Matrix{T}
end


text:
struct GetIndex{T,N} <: ArrayOp{T,N}
    input::ArrayOp
    idx::Tuple
end

GetIndex(input::ArrayOp, idx::Tuple) =
    GetIndex{eltype(input), ndims(input)}(input, idx)

function stage(ctx, gidx::GetIndex)
    inp = cached_stage(ctx, gidx.input)

    dmn = domain(inp)
    idxs = [if isa(gidx.idx[i], Colon)
        indexes(dmn)[i]
    else
        gidx.idx[i]
    end for i in 1:length(gidx.idx)]

    # Figure out output dimension
    view(inp, ArrayDomain(idxs))
end

function size(x::GetIndex)
    map(a -> a[2] isa Colon ?
        size(x.input, a[1]) : length(a[2]),
        enumerate(x.idx)) |> Tuple
end

struct GetIndexScalar <: Computation
    input::ArrayOp
    idx::Tuple
end

function stage(ctx, gidx::GetIndexScalar)
    inp = cached_stage(ctx, gidx.input)
    s = view(inp, ArrayDomain(gidx.idx))
    delayed(identity)(collect(s)[1])
end

Base.getindex(c::ArrayOp, idx::ArrayDomain) = GetIndex(c, indexes(idx))
Base.getindex(c::ArrayOp, idx...)           = GetIndex(c, idx)
Base.getindex(c::ArrayOp, idx::Integer...)  = compute(GetIndexScalar(c, idx))


text:
using FileIO
include("julia2ims.jl")

field = ["s3-2" "s3-5" "s3-2" "s3-5" "s3-6"	"s4-9"	"s4-6"	"s4-4"	"s5-4"	"s5-2"	"s5-6" "s7-4"	"s8-9"	"s8-8"	"s11-3"	"s11-7"	"s13-3"	"s14-2"	"s16-5"	"s16-6" "s17-2"	"s20-1"	"s21-2"	"s23-4"	"s25-1"	"s30-4"	"s31-6"	"s32-2"	"s34-4" "s36-1"	"s36-2"	"s36-4" ];

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200314";
output_dir = "$(data_dir)_visualization"

try
    mkdir("$output_dir")
catch
end

function max_projection(cell)
    println("Processing $cell")
    #raw_img = load(File(format"TIFF", "$data_dir/HE7-11-1-80uw-const_3_$pos.ome.tiff"));
    raw_img = loadims2("$(data_dir)_result/output/d14$cell.ims");
    img_size = size(raw_img);
    img_max = zeros(N0f16, img_size[1], img_size[2], Int(img_size[3]/20));
    for i in 1:Int(img_size[3]/20)
        img_max[:, :, i] = maximum(raw_img[:, :, (i-1)*20+1:i*20], dims=3);
    end
    save("$output_dir/$cell.tiff", img_max);
end

for cell in field
    max_projection(cell)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    deformspace(domain, localaniso, metric; kwargs)
    deformspace(domain, localaniso, metric, refvariogram; kwargs)
	deformspace(samples, domain, localaniso, metric; kwargs)
    deformspace(samples, domain, localaniso, metric, refvariogram; kwargs)
    deformspace(graphobject, metric=GraphDistance; kwargs)
	kwargs = (anchors=1500, maxoutdim=10, weights=nothing)

Spatial deformation is a method that transforms coordinates to an isotropic
space in high dimension where the local anisotropies information are embedded
into the new distances. Traditional estimations and simulations can be performed
in this new data configuration. It's necessary to inform the `domain` object
(and the `samples` if applied). The local anisotropies are passed via `localaniso`
and `metric` define how the local anisotropies are used to calculate the new
data distances. Available metrics to calculate distance between two points are:

* `AnisoDistance`   - averaged anisotropic distance
* `KernelVariogram` - non-stationary variogram kernel estimator
* `GraphDistance`   - geodesic distances of a informed graph. More details
  of the graph construction in [`graph`](@ref) docstring.

A reference variogram `refvariogram` is necessary if metric is `KernelVariogram`.
Additional keyword arguments are the number of `anchors` to perform landmark MDS
if the number of data is too big for a traditional MDS. `maxoutdim` is the
maximum number of dimensions to retain after MDS. `weights` are the declustering
weights of the domain data (+ samples) in order to help selecting declustered
anchor points for the transformations.
"""
function deformspace(obj::SpatialData, lpar::LocalAnisotropy, metric::Type{<:LocalMetric};
	anchors=1500, maxoutdim=10, weights=nothing)
	D = LocalGeoData(obj,lpar)
	deformspace(D, metric, anchors=anchors, maxoutdim=maxoutdim, weights=weights)
end

function deformspace(obj::SpatialData, lpar::LocalAnisotropy, metric::Type{<:LocalMetric},
	refvario::Variogram; anchors=1500, maxoutdim=10, weights=nothing)
	D = LocalGeoData(obj,lpar,refvario)
	deformspace(D, metric, anchors=anchors, maxoutdim=maxoutdim, weights=weights)
end

function deformspace(hd::SpatialData, obj::SpatialData, lpar::LocalAnisotropy,
	metric::Type{<:LocalMetric}; anchors=1500, maxoutdim=10, weights=nothing)
	D = LocalGeoData(hd,obj,lpar)
	deformspace(D, metric, anchors=anchors, maxoutdim=maxoutdim, weights=weights)
end

function deformspace(hd::SpatialData, obj::SpatialData, lpar::LocalAnisotropy,
	metric::Type{<:LocalMetric}, refvario::Variogram; anchors=1500, maxoutdim=10, weights=nothing)
	D = LocalGeoData(hd,obj,lpar,refvario)
	deformspace(D, metric, anchors=anchors, maxoutdim=maxoutdim, weights=weights)
end

# metric = KernelVariogram, AnisoDistance, GraphDistance
function deformspace(D::LocalGeoData, metric::Type{<:LocalMetric}=GraphDistance;
	anchors=1500, maxoutdim=10, weights=nothing)

	#@assert graph exists if GraphDistance
	dim, n = ndims(D), nall(D)
	n < anchors && (anchors = n)

	ianchors = setanchors(n, anchors, weights)
	ADM = Array{Float64}(undef, (anchors, anchors))
	dissmatrix!(ADM, D, metric, ianchors)

	if n > anchors
		iothers = setdiff(1:n, ianchors)
		atcoords, M1, M3 = anchors_mds(ADM, maxoutdim)
		dim = size(atcoords, 1)
		otcoords = Array{Float64}(undef, (dim, length(iothers)))
		Threads.@threads for i in 1:length(iothers)
			otcoords[:,i] .= triangulation(D,metric,iothers[i],ianchors,M1,M3)
		end
		tcoords = Array{Float64}(undef, (dim,n))
		tcoords[:,ianchors] .= atcoords
		tcoords[:,iothers] .= otcoords
	else
		M = fit(MDS, ADM, maxoutdim=maxoutdim, distances=true)
    	tcoords = transform(M)
		#println("Explained variance: $(sum(sortλ[1:maxdim])/sum(sortλ[λ .> 0]))")
	end

	outobj(D, tcoords)
end

function setanchors(n,anchors,weights)
	n == anchors && (return collect(1:n))
	# maybe need to convert GeoStatsBase weights to StatsBase format
	args = weights == nothing ? (1:n, anchors) : (1:n, weights, anchors)
	ianchors = sample(args..., replace=false)
	sort!(ianchors)
end

function dissmatrix!(ADM, D::LocalGeoData, metric::Type{<:LocalMetric}, iax::Vector{Int})
	n = size(ADM,1)
	for (i,ia) in enumerate(iax)
		dcols = colwise(D, metric, ia, iax[i:n])
		for (d,j) in zip(dcols,i:n)
			ADM[i,j] = ADM[j,i] = d
		end
	end
end

function anchors_mds(ADM, maxoutdim)
	nx = size(ADM,1)
	G = dmat2gram(ADM)
    F = eigen(Symmetric(G))
	λ = F.values
	maxdim = minimum([maxoutdim, nx-1, sum(λ .> 0)])
	sorti = sortperm(F.values,rev=true)
	sortλ = λ[sorti]
    EM = (F.vectors[:,sorti])[:,1:maxdim]
	#println("Explained variance: $(sum(sortλ[1:maxdim])/sum(sortλ[sortλ .> 0]))")
    sq_eigenvals = sortλ[1:maxdim].^0.5
    AM = Diagonal(sq_eigenvals)
    atcoords = permutedims(EM*AM)

	# to use later for another points allocations
	M1 = permutedims(EM)
	M1 ./= reshape(sq_eigenvals,maxdim,1)
	M3 = Array{Float64}(undef,nx) #Array{Float64}(undef,(nx,1))
	mean!(M3,ADM.^2)

	atcoords, M1, M3
end

function triangulation(D::LocalGeoData, metric::Type{<:LocalMetric},i,j,M1,M3)
	M2 = colwise(D, metric, i, j) .^ 2
	-0.5*M1*(M2-M3)
end

function outobj(D, coord)
	dom = values(obj(D))
	out = if sdata(D)
		n, hn = nvals(D), nall(D)
		h1 = n+1
		dc   = view(coord,:,1:n)
		dobj = dom isa Domain ? PointSet(dc) : georef(dom, dc)
		georef(values(sobj(D)), view(coord,:,h1:hn)), dobj
	else
		dobj = dom isa Domain ? PointSet(coord) : georef(dom, coord)
		dobj
	end
	out
end

"""
    to_3d(geodata)

The spatial deformation method often return an output in dimensions higher
than 3-D. This function create a new spatial object retaining only the
first three dimensions in order to help plotting of the data.
"""
function to_3d(s)
	dom = domain(s)
	c = reduce(hcat,[coordinates(centroid(dom,x))[1:3] for x in 1:nelements(dom)])
	georef(values(s),PointSet(c))
end


text:
# Core traits and accessor functions

"""
`alpha(p)` extracts the alpha component of a color. For a color
without an alpha channel, it will always return 1.
"""
alpha(c::TransparentColor) = c.alpha
alpha(c::Color)   = oneunit(eltype(c))
alpha(c::RGB24)   = N0f8(1)
alpha(c::ARGB32)  = N0f8((c.color & 0xff000000)>>24, 0)
alpha(c::AGray32) = N0f8((c.color & 0xff000000)>>24, 0)

"`red(c)` returns the red component of an `AbstractRGB` opaque or transparent color."
red(c::AbstractRGB   ) = c.r
red(c::TransparentRGB) = c.r
red(c::RGB24)  = N0f8((c.color & 0x00ff0000)>>16, 0)
red(c::ARGB32) = N0f8((c.color & 0x00ff0000)>>16, 0)

"`green(c)` returns the green component of an `AbstractRGB` opaque or transparent color."
green(c::AbstractRGB   ) = c.g
green(c::TransparentRGB) = c.g
green(c::RGB24)  = N0f8((c.color & 0x0000ff00)>>8, 0)
green(c::ARGB32) = N0f8((c.color & 0x0000ff00)>>8, 0)

"`blue(c)` returns the blue component of an `AbstractRGB` opaque or transparent color."
blue(c::AbstractRGB   ) = c.b
blue(c::TransparentRGB) = c.b
blue(c::RGB24)  = N0f8(c.color & 0x000000ff, 0)
blue(c::ARGB32) = N0f8(c.color & 0x000000ff, 0)

"`gray(c)` returns the gray component of a grayscale opaque or transparent color."
gray(c::Gray)    = c.val
gray(c::TransparentGray) = c.val
gray(c::Gray24)  = N0f8(c.color & 0x000000ff, 0)
gray(c::AGray32) = N0f8(c.color & 0x000000ff, 0)
gray(x::Number)  = x

Base.real(g::Gray) = gray(g)

"""
`chroma(c)` returns the chroma of a `Lab`, `Luv` or their variants color.
!!! note
    The other color types (e.g. `RGB`, `HSV` or `YCbCr`) are not supported
    because their definitions of *chroma* are not clear. Colorfulness, chroma
    and saturation are defined as distinct aspects by the CIE.
"""
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{Lab,DIN99,DIN99o,DIN99d}} = sqrt(c.a^2 + c.b^2)
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Luv} = sqrt(c.u^2 + c.v^2)
chroma(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{LCHab,LCHuv}} = c.c

"""
`hue(c)` returns the hue in degrees. This function does not guarantee that the
return value is in [0, 360].
"""
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{HSV,HSL,HSI,LCHab,LCHuv}} = c.h
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Union{Lab,DIN99,DIN99o,DIN99d}} =
    (h = atand(c.b, c.a); h < 0 ? h + 360 : h)
hue(c::Union{C,AlphaColor{C},ColorAlpha{C}}) where {C<:Luv} =
    (h = atand(c.v, c.u); h < 0 ? h + 360 : h)

# Extract the first, second, and third arguments as you'd
# pass them to the constructor
"""
`comp1(c)` extracts the first component you'd pass to the constructor
of the corresponding object.  For most color types without an alpha
channel, this is just the first field, but for types like `BGR` that
reverse the internal storage order this provides the value that you'd
use to reconstruct the color.

Specifically, for any `Color{T,3}`,

    c == typeof(c)(comp1(c), comp2(c), comp3(c))

returns true.
"""
comp1(c::AbstractRGB) = red(c)
comp1(c::Union{AlphaColor{C},ColorAlpha{C}}) where {C<:AbstractRGB} = red(c)
comp1(c::Union{Color,ColorAlpha}) = getfield(c, 1)
comp1(c::AlphaColor) = getfield(c, 2)
comp1(c::AGray32) = gray(c)

"`comp2(c)` extracts the second constructor argument (see `comp1`)."
comp2(c::AbstractRGB) = green(c)
comp2(c::Union{AlphaColor{C},ColorAlpha{C}}) where {C<:AbstractRGB} = green(c)
comp2(c::Union{Color,ColorAlpha}) = getfield(c, 2)
comp2(c::AlphaColor) = getfield(c, 3)

"`comp3(c)` extracts the third constructor argument (see `comp1`)."
comp3(c::AbstractRGB) = blue(c)
comp3(c::Union{AlphaColor{C},ColorAlpha{C}}) where {C<:AbstractRGB} = blue(c)
comp3(c::Union{Color,ColorAlpha}) = getfield(c, 3)
comp3(c::AlphaColor) = getfield(c, 4)

"`color(c)` extracts the opaque color component from a Colorant (e.g., omits the alpha channel, if present)."
color(c::Color) = c
color(c::TransparentColor{C,T,4}) where {C,T} = C(comp1(c), comp2(c), comp3(c))
color(c::TransparentColor{C,T,2}) where {C,T} = C(comp1(c))

# Some of these traits exploit a nice trick: for subtypes, walk up the
# type hierarchy until we get to a stage where we can define the
# function in general

# recurse up the type hierarchy until you get to Colorant{T,N} for
# specific T,N.
to_top(::Type{Colorant{T,N}}) where {T,N} = Colorant{T,N}
@pure to_top(::Type{C}) where {C<:Colorant} = to_top(supertype(C))

to_top(c::Colorant) = to_top(typeof(c))

# eltype(RGB{Float32}) -> Float32
eltype(::Type{Colorant{T}}) where {T          }   = T
eltype(::Type{Colorant{T,N}}) where {T,N        } = T
@pure eltype(::Type{C}) where {C<:Colorant} = eltype(supertype(C))

eltype(c::Colorant) = eltype(typeof(c))

# eltypes_supported(Colorant{T<:X}) -> T<:X (pre 0.6) or X (post 0.6)
@pure eltypes_supported(::Type{C}) where {C<:Colorant} =
    Base.parameter_upper_bound(base_colorant_type(C), 1)

eltypes_supported(c::Colorant) = eltypes_supported(typeof(c))

@pure issupported(::Type{C}, ::Type{T}) where {C<:Colorant,T} = T <: eltypes_supported(C)

"""
`color_type(c)` or `color_type(C)` (`c` being a color instance and `C`
being the type) returns the type of the Color object (without
alpha channel).  This, and related functions like `base_color_type`,
`base_colorant_type`, and `ccolor` are useful for manipulating types for
writing generic code.

For example,

    color_type(RGB)          == RGB
    color_type(RGB{Float32}) == RGB{Float32}
    color_type(ARGB{N0f8})     == RGB{N0f8}
"""
color_type(::Type{TransparentColor})        = Color
color_type(::Type{C}) where {C<:Color} = C
color_type(c::Colorant) = color_type(typeof(c))

# Return the number of components in the color
# Note this is different from div(sizeof(c), sizeof(eltype(c))) (e.g., XRGB)
length(c::Colorant) = length(typeof(c))

length(::Type{C}) where C<:(Colorant{T,N} where T) where N = N
# This definition should be unnecessary, but julia currently incorrectly
# dispatches to the first definition below, even when the second is
# applicable.
_color_type(::Type{TC}) where TC<:(TransparentColor{C, T, N} where T where N) where C = C
color_type(::Type{TC}) where TC<:TransparentColor =
    isa(TC, UnionAll) ? Base.parameter_upper_bound(TC, 1) : _color_type(TC)
color_type(::Type{TC}) where TC<:(TransparentColor{C, T, N} where T where N) where C = C

@pure color_type(::Type{C}) where {C<:AlphaColor} = color_type(supertype(C))
@pure color_type(::Type{C}) where {C<:ColorAlpha} = color_type(supertype(C))

"""
`base_color_type` is similar to `color_type`, except it "strips off" the
element type.  For example,

    color_type(RGB{N0f8})     == RGB{N0f8}
    base_color_type(RGB{N0f8}) == RGB

This can be very handy if you want to switch element types. For example:

    c64 = base_color_type(c){Float64}(color(c))

converts `c` into a `Float64` representation (potentially discarding
any alpha-channel information).
"""
base_color_type(::Type{C}) where {C<:Colorant} = base_colorant_type(color_type(C))

base_color_type(c::Colorant) = base_color_type(typeof(c))
base_color_type(x::Number)   = Gray

@pure basetype(T) = Base.typename(T).wrapper
base_colorant_type(::Type{C}) where {C<:Colorant} = basetype(C)

"""
`base_colorant_type` is similar to `base_color_type`, but it preserves the
"alpha" portion of the type.

For example,

    base_color_type(ARGB{N0f8})  == RGB
    base_colorant_type(ARGB{N0f8})  == ARGB

If you just want to switch element types, this is the safest default
and the easiest to use:

    c64 = base_colorant_type(c){Float64}(c)
"""
base_colorant_type(c::Colorant) = base_colorant_type(typeof(c))

"""
    floattype(::Type{T}) where T<:Colorant

Promote storage data type of colorant `T` to `AbstractFloat` while keep the
`base_colorant_type` the same.

!!! info

    Non-parametric colorants will be promote to corresponding parametric
    colorants. For example, `floattype(RGB24) == RGB{Float32}`
"""
floattype(::Type{T}) where T <: Colorant =
    base_colorant_type(T){floattype(eltype(T))} # 1 parameter
# 0 parameter
floattype(::Type{RGB24}) = RGB{Float32}
floattype(::Type{Gray24}) = Gray{Float32}
floattype(::Type{ARGB32}) = ARGB{Float32}
floattype(::Type{AGray32}) = AGray{Float32}

colorant_string(::Type{Union{}}) = "Union{}"
colorant_string(::Type{C}) where {C<:Colorant} = string(nameof(C))
function colorant_string_with_eltype(::Type{C}) where {C<:Colorant}
    io = IOBuffer()
    colorant_string_with_eltype(io, C)
    String(take!(io))
end
colorant_string_with_eltype(io::IO, ::Type{Union{}}) = show(io, Union{})
function colorant_string_with_eltype(io::IO, ::Type{C}) where {C<:Colorant}
    print(io, colorant_string(C), '{')
    showcoloranttype(io, eltype(C))
    print(io, '}')
end
# Nonparametric types
colorant_string_with_eltype(io::IO, ::Type{Gray24})  = print(io, "Gray24")
colorant_string_with_eltype(io::IO, ::Type{AGray32}) = print(io, "AGray32")
colorant_string_with_eltype(io::IO, ::Type{RGB24})   = print(io, "RGB24")
colorant_string_with_eltype(io::IO, ::Type{ARGB32})  = print(io, "ARGB32")

showcoloranttype(io, ::Type{Union{}}) = show(io, Union{})
showcoloranttype(io, ::Type{T}) where {T<:FixedPoint} = FixedPointNumbers.showtype(io, T)
showcoloranttype(io, ::Type{T}) where {T} = show(io, T)

"""
 `ccolor` ("concrete color") helps write flexible methods. The idea is
that users may write `convert(HSV, c)` or even `convert(Array{HSV},
A)` without specifying the element type explicitly (e.g.,
`convert(Array{HSV{Float32}}, A)`). `ccolor` implements the logic "choose the
user's eltype if specified, otherwise retain the eltype of the source
object." However, when the source object has FixedPoint element type,
and the destination only supports AbstractFloat, we choose Float32.

Usage:

    ccolor(desttype, srctype) -> concrete desttype

Example:

    convert{C<:Colorant}(::Type{C}, p::Colorant) = cnvt(ccolor(C,typeof(p)), p)

where `cnvt` is the function that performs explicit conversion.
"""
ccolor(::Type{Colorant   }, ::Type{Csrc}) where {   Csrc<:Colorant} = Csrc
ccolor(::Type{Colorant{T}}, ::Type{Csrc}) where {T, Csrc<:Colorant} = base_colorant_type(Csrc){T}
ccolor(::Type{Colorant{T,3}}, ::Type{Csrc}) where {T, Csrc<:Color3  } = Csrc
ccolor(::Type{Colorant{T,3}}, ::Type{Csrc}) where {T, Csrc<:Transparent3} = base_color_type(Csrc)
ccolor(::Type{Color   }, ::Type{Csrc}) where {   Csrc<:Colorant} = color_type(Csrc)
ccolor(::Type{Color{T}}, ::Type{Csrc}) where {T, Csrc<:Colorant} = base_color_type(Csrc){T}

ccolor(::Type{TransparentColor}, ::Type{Csrc}) where {Csrc<:Color} =
          error("Ambiguous storage order, choose AlphaColor or ColorAlpha")
ccolor(
::Type{TransparentColor{C    }}, ::Type{Csrc}) where {C<:Color,    Csrc<:Color} =
           error("Ambiguous storage order, choose AlphaColor or ColorAlpha")
ccolor(
::Type{TransparentColor{C,T  }}, ::Type{Csrc}) where {C<:Color,T,  Csrc<:Color} =
           error("Ambiguous storage order, choose AlphaColor or ColorAlpha")
ccolor(
::Type{TransparentColor{C,T,N}}, ::Type{Csrc}) where {C<:Color,T,N,Csrc<:Color} =
           error("Ambiguous storage order, choose AlphaColor or ColorAlpha")

ccolor(::Type{TransparentColor}, ::Type{Csrc}) where {Csrc<:TransparentColor} = Csrc

ccolor(::Type{AlphaColor}, ::Type{Csrc}) where {Csrc<:Colorant} = alphacolor(Csrc)
ccolor(
::Type{AlphaColor{C    }}, ::Type{Csrc}) where {C<:Color,    Csrc<:Colorant} = ccolor(alphacolor(C), Csrc)
ccolor(
::Type{AlphaColor{C,T  }}, ::Type{Csrc}) where {C<:Color,T,  Csrc<:Colorant} = ccolor(alphacolor(C){T}, Csrc)
ccolor(
::Type{AlphaColor{C,T,N}}, ::Type{Csrc}) where {C<:Color,T,N,Csrc<:Colorant} = ccolor(alphacolor(C){T}, Csrc)

ccolor(::Type{ColorAlpha}, ::Type{Csrc}) where {Csrc<:Colorant} = coloralpha(Csrc)
ccolor(
::Type{ColorAlpha{C    }}, ::Type{Csrc}) where {C<:Color,    Csrc<:Colorant} = ccolor(coloralpha(C), Csrc)
ccolor(
::Type{ColorAlpha{C,T  }}, ::Type{Csrc}) where {C<:Color,T,  Csrc<:Colorant} = ccolor(coloralpha(C){T}, Csrc)
ccolor(
::Type{ColorAlpha{C,T,N}}, ::Type{Csrc}) where {C<:Color,T,N,Csrc<:Colorant} = ccolor(coloralpha(C){T}, Csrc)

ccolor(::Type{AbstractRGB},    ::Type{Csrc}) where {  Csrc<:AbstractRGB} = Csrc
ccolor(::Type{AbstractRGB{T}}, ::Type{Csrc}) where {T,Csrc<:AbstractRGB} = base_colorant_type(Csrc){T}

# Generic concrete types
ccolor(::Type{Cdest}, ::Type{Csrc}) where {Cdest<:Colorant,Csrc<:Colorant} = _ccolor(Cdest, Csrc, pick_eltype(Cdest, eltype(Cdest), eltype(Csrc)))
ccolor(::Type{Cdest}, ::Type{T}) where {Cdest<:AbstractGray,T<:Number} = _ccolor(Cdest, Gray, pick_eltype(Cdest, eltype(Cdest), T))

_ccolor(::Type{Cdest}, ::Type{Csrc}, ::Type{T}) where {Cdest,Csrc,T<:Number} =
    isconcretetype(T) ? base_colorant_type(Cdest){T} :
                    base_colorant_type(Cdest){S} where S<:T

_ccolor(          ::Type{Cdest}, ::Type{Csrc}, ::Any) where {Cdest,Csrc}     = Cdest

# Specific concrete types
ccolor(::Type{RGB24},   ::Type{Csrc}) where {Csrc<:Colorant} = RGB24
ccolor(::Type{ARGB32},  ::Type{Csrc}) where {Csrc<:Colorant} = ARGB32
ccolor(::Type{Gray24},  ::Type{Csrc}) where {Csrc<:Colorant} = Gray24
ccolor(::Type{AGray32}, ::Type{Csrc}) where {Csrc<:Colorant} = AGray32

pick_eltype(::Type{C}, ::Type{T1}, ::Type{T2}) where {C,T1<:Number,T2<:Number} = T1
pick_eltype(::Type{C}, ::Any, ::Any) where {C} = eltypes_supported(C)
pick_eltype(::Type{C}, ::Any, ::Type{T2}) where {C,T2<:Number} = issupported(C, T2) ? T2 : eltype_default(C)

### Equality
function ==(c1::AbstractRGB, c2::AbstractRGB)
    red(c1) == red(c2) && green(c1) == green(c2) && blue(c1) == blue(c2)
end
==(c1::HSV, c2::HSV) = c1.h == c2.h && c1.s == c2.s && c1.v == c2.v
==(c1::HSI, c2::HSI) = c1.h == c2.h && c1.s == c2.s && c1.i == c2.i
==(c1::HSL, c2::HSL) = c1.h == c2.h && c1.s == c2.s && c1.l == c2.l
==(c1::XYZ, c2::XYZ) = c1.x == c2.x && c1.y == c2.y && c1.z == c2.z
==(c1::xyY, c2::xyY) = c1.x == c2.x && c1.y == c2.y && c1.Y == c2.Y
==(c1::Lab, c2::Lab) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::Luv, c2::Luv) = c1.l == c2.l && c1.u == c2.u && c1.v == c2.v
==(c1::LCHab, c2::LCHab) = c1.l == c2.l && c1.c == c2.c && c1.h == c2.h
==(c1::LCHuv, c2::LCHuv) = c1.l == c2.l && c1.c == c2.c && c1.h == c2.h
==(c1::DIN99, c2::DIN99) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::DIN99d, c2::DIN99d) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::DIN99o, c2::DIN99o) = c1.l == c2.l && c1.a == c2.a && c1.b == c2.b
==(c1::LMS, c2::LMS) = c1.l == c2.l && c1.m == c2.m && c1.s == c2.s
==(c1::YIQ, c2::YIQ) = c1.y == c2.y && c1.i == c2.i && c1.q == c2.q
==(c1::YCbCr, c2::YCbCr) = c1.y == c2.y && c1.cb == c2.cb && c1.cr == c2.cr

for T in (RGB24, ARGB32, Gray24, AGray32)
    @eval begin
        reinterpret(::Type{UInt32}, x::$T) = x.color
        reinterpret(::Type{$T}, x::UInt32) = $T(x, Val{true})
    end
end
==(x::AbstractGray, y::AbstractGray) = gray(x) == gray(y)
==(x::Number, y::AbstractGray) = x == gray(y)
==(x::AbstractGray, y::Number) = ==(y, x)

function ==(x::TransparentColor, y::TransparentColor)
    color(x) == color(y) && alpha(x) == alpha(y)
end

struct BoolTuple end
@inline BoolTuple(args::Bool...) = (args...,)

function _isapprox(a::Colorant, b::Colorant; kwargs...)
    componentapprox(x, y) = isapprox(x, y; kwargs...)
    all(ColorTypes._mapc(BoolTuple, componentapprox, a, b))
end
isapprox(a::C, b::C; kwargs...) where {C<:Colorant} =
    _isapprox(a, b; kwargs...)
isapprox(a::Colorant, b::Colorant; kwargs...) =
    _isapprox(base_colorant_type(a), base_colorant_type(b), a, b; kwargs...)
_isapprox(::Type{C}, ::Type{C}, a, b; kwargs...) where {C<:Colorant} =
    _isapprox(a, b; kwargs...)
_isapprox(::Type{<:AbstractRGB}, ::Type{<:AbstractRGB}, a, b; kwargs...) =
    _isapprox(RGB(a), RGB(b); kwargs...)
_isapprox(::Type{<:AbstractGray}, ::Type{<:AbstractGray}, a, b; kwargs...) =
    isapprox(gray(a), gray(b); kwargs...)
_isapprox(TA::Type, TB::Type, a, b; kwargs...) = false
isapprox(x::Number, y::AbstractGray; kwargs...) =
    isapprox(x, gray(y); kwargs...)
isapprox(x::AbstractGray, y::Number; kwargs...) =
    isapprox(y, x; kwargs...)

zero(::Type{C}) where {C<:Gray} = C(0)
oneunit(::Type{C}) where {C<:Gray} = C(1)

function Base.one(::Type{C}) where {C<:Gray}
    Base.depwarn("one($C) will soon switch to returning 1; you might need to switch to `oneunit`", :one)
    C(1)
end

Base.broadcastable(x::Colorant) = Ref(x)


text:
include("../src/PushRelabelMaxFlow.jl")

function clique_expansion_penalties(k)
    """
    Symmetric clique splitting penalty for hyperedge of size k.
    Only gives first r = floor(k/2) penalties
    """
    @assert(k > 1)
    r = floor(Int64,(k)/2)
    I = collect(0:k)
    w = I .* (k .- I)
    w = w[1:r+1]
    return w
end

function clique_expansion_EdgesW(order::Vector{Int64})
    EdgesW = Vector{Vector{Float64}}()
    for k in order
        weight = 1/(k-1)
        push!(EdgesW,weight*clique_expansion_penalties(k))
    end
    return EdgesW
end

function sqrt_EdgesW(order::Vector{Int64})
    EdgesW = Vector{Vector{Float64}}()
    for k in order
        r = floor(Int64,(k)/2)
        w = sqrt.(collect(0:1:r))
        # w = min.(w,delta)
        push!(EdgesW,w)
    end
    return EdgesW
end

function generic_EdgesW(order::Vector{Int64},fun)
    # fun is a generic function
    EdgesW = Vector{Vector{Float64}}()
    for k in order
        r = floor(Int64,(k)/2)
        w = fun(collect(0:1:r))
        push!(EdgesW,w)
    end
    return EdgesW
end

function clique_expansion_EdgesW(Edges::Vector{Vector{Int64}},weighted=true)
    EdgesW = Vector{Vector{Float64}}()
    for e in Edges
        k = length(e)
        if weighted
            weight = 1/k
        else
            weight = 1
        end
        push!(EdgesW,weight*clique_expansion_penalties(k))
    end
    return EdgesW
end

function delta_linear_EdgesW(Edges::Vector{Vector{Int64}},delta::Float64)
    EdgesW = Vector{Vector{Float64}}()
    for e in Edges
        k = length(e)
        r = floor(Int64,(k)/2)
        w = collect(0:1:r)
        w = min.(w,delta)
        push!(EdgesW,w)
    end
    return EdgesW
end

function gen_hypergraph_cut(Edges,EdgesW,S,n)
    """
    For a set S, evaluate generalized hypergraph cut using splitting functions
    from EdgesW (symmetric and cardinality-based).
    """
    eS = zeros(n)
    eS[S] = ones(length(S))

    # higher-order cut penalties
    cutval=0
    for i = 1:length(Edges)
        e = Edges[i]
        k = length(e)
        w = EdgesW[i]
        t = round(Int64,sum(eS[e]))

        # do this because it's a symmetric splitting function,
        # so w is only half the size
        t = min(t,k-t)
        cutval += w[t+1]
    end

    return cutval

end

# For a set S in a graph with adjacency matrix A, return some information about
# S including its conductance, number of interior edges, volume, and cut.
function set_stats(A::SparseMatrixCSC{Float64,Int64},
    S::Vector{Int64},volA::Float64)

    if volA == 0.0
        volA = sum(A.nzval)
    end

    if length(S) == size(A,1)
        # then we have an indicator vector
        S = findall(x->x!=0,eS)
        AS = A[:,S]
    else
        # then we have a subset
        @assert(minimum(S) >= 1)
        @assert(maximum(S) <= size(A,1))
        AS = A[:,S]
    end

    vol = sum(AS.nzval);
    SAS = AS[S,:]
    edges = sum(SAS.nzval);
    cut = vol-edges

    cond = cut/minimum([vol,volA-vol]);

    return cut, vol, edges, cond

end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
# This file contains the definition for the nodes struct.

mutable struct Nodes
	id
	systemload
	load
	exchange
	bmva
	slack

	function Nodes(nodes_df::DataFrame,
				   load_df::DataFrame,
				   exchange_df::DataFrame,
				   bmva::Number,
				   slack::Number
				   )

		# Initialise dataframe, including all nodes with zeros
		load_full = DataFrame(zeros(length(load_df[!, 1]),
										   length(nodes_df[!, 1])
									)
							 , Symbol.(nodes_df[!, 1])
							 )

		systemload = load_df[!, 1]

		# Calculate the absolute load at a node
		for col in 2:ncol(load_df)
			load_df[!, col] = load_df[!, col] .* load_df[!, 1]
		end

		load_df = load_df[!, 2:end]

		for col in names(load_df)
			load_full[!, col] = load_df[!, col]
		end

		N = parse.(Int,String.(names(load_full)))
		load_dict = Dict(N[col] => load_full[!, col] for col in 1:length(N))

		# Exchange
		exchange_full = DataFrame(zeros(length(load_df[!, 1]),
										   length(nodes_df[!, 1])
									)
							 , Symbol.(nodes_df[!, 1])
							 )

		for col in names(exchange_df)
			exchange_full[!, col] = exchange_df[!, col]
		end

		N_exchange = parse.(Int,String.(names(exchange_full)))
		exchange_dict = Dict(N_exchange[col] => exchange_full[!, col] for col in 1:length(N_exchange))

		return new(N,
				   systemload,
				   load_dict,
				   exchange_dict,
				   bmva,
				   slack
				   )
	end
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
"""
Get a list of all AWS service API definition files from the `awsk-sdk-js` GitHub repository.
"""
function _get_aws_sdk_js_files(repo_name::String, auth::GitHub.OAuth2)
    master_tree = @mock tree(repo_name, "master"; auth=auth)
    apis_sha = [t for t in master_tree.tree if t["path"]=="apis"][1]["sha"]
    files = @mock tree(repo_name, apis_sha)
    files = files.tree

    filter!(f -> endswith(f["path"], ".normal.json"), files)

    return _filter_latest_service_version(files)
end

"""
Return a list of all AWS Services and their latest version.
"""
function _filter_latest_service_version(services::AbstractArray)
    seen_services = Set{String}()
    latest_versions = OrderedDict[]

    for service in reverse(services)
        service_name, _ = _get_service_and_version(service["path"])

        if !(service_name in seen_services)
            push!(seen_services, service_name)
            push!(latest_versions, service)
        end
    end

    return latest_versions
end

"""
Return a string with line breaks added such that lines are wrapped at or before the limit.
"""
function _wraplines(str, limit=92; delim="\n")
    lines = String[]

    while !isempty(str)
        line, str = _splitline(str, limit)
        push!(lines, rstrip(line))  # strip trailing whitespace
    end

    return join(lines, delim)
end

"""
Split the string `str` at or before `limit`.
Prefers splitting the string on whitespace rather than mid-word, when possible.
`limit` is measured in codeunits, which is an upper-bound on the number of characters.
"""
function _splitline(str, limit)
    limit >= 1 || throw(DomainError(limit, "Lines cannot be split before the first char."))
    ncodeunits(str) <= limit && return (str, "")
    limit = _validindex(str, limit)
    first_line = str[1:limit]
    # split on whitespace if possible, else just split when we hit the limit.
    split_point = something(findlast(==(' '), first_line), limit)
    stop = _validindex(first_line, split_point)

    while ispunct(first_line[stop])  # avoid splitting escaped characters.
        stop = prevind(first_line, stop)
    end

    restart = nextind(first_line, stop)

    return (str[1:stop], str[restart:end])
end

"""
Return a valid index into the string `str`, rounding towards the first index of `str` if `i` is not itself a valid index into `str`.
`i` must be within the bounds of `string`. `_validindex(str, i)` only protects against a `StringIndexError`, not a `BoundsError`.
"""
function _validindex(str, limit)
    prev = max(firstindex(str), prevind(str, limit))
    next = nextind(str, prev)

    return next == limit ? limit : prev
end

"""
Convert a function name from CamelCase to snake_case
"""
function _format_name(function_name::String)
    # Replace a string of uppercase characters with themselves prefaced by an underscore
    # [A-Z](?![A-Z]) => Match a single uppercase character that is not followed by another uppercase character
    # |(A-Z]{1,})    => Match 1-Infinite amounts of uppercase characters
    function_name = replace(function_name, r"[A-Z](?![A-Z])|([A-Z]{1,})" => s"_\g<0>")

    # Lowercase the entire string
    function_name = lowercase(function_name)

    # Chop off the leading underscore
    return if startswith(function_name, "_")
        chop(function_name, head=1, tail=0)
    else
        function_name
    end
end

"""
Replace URI parameters with the appropriate syntax for Julia interpolation.

Find all URI parameters, and apply the following replacements:
    * { => \$(
    * } => )
    * - => _
    * + => empty string

Example: "/v1/configurations/{configuration-id}" => "/v1/configurations/\$(configuration_id)"
"""
function _clean_uri(uri::String)
    uri_parameters = eachmatch(r"{.*?}", uri) # Match anything surrounded in "{ }"

    for param in uri_parameters
        match = param.match
        original_match = match

        match = replace(match, '{' => "\$(")  # Replace { with $(
        match = replace(match, '}' => ')')  # Replace } with )
        match = replace(match, '-' => '_')  # Replace hyphens with underscores
        match = replace(match, '+' => "")  # Remove +

        uri = replace(uri, original_match => match)
    end

    return uri
end

"""
Clean up the documentation to make it Julia compiler and human-readable.
    * Remove any HTML tags
    * Remove any dollar signs
    * Remove any backslashes
    * Escape any double-quotes
"""
function _clean_documentation(documentation::String)
    documentation = replace(documentation, r"\<.*?\>" => "")
    documentation = replace(documentation, '$' => "")
    documentation = replace(documentation, '\\' => "")
    documentation = replace(documentation, '"' => "\\\"")

    return documentation
end

"""
Get the `service` and `version` from a filename.
Example filename: `{Service}-{Version}.normal.json`
"""
function _get_service_and_version(filename::String)
    try
        # Remove ".normal.json" suffix
        service_and_version = join(split(filename, '.')[1:end-2],'.')

        service_and_version = split(service_and_version, '-')
        service = join(service_and_version[1:end-3], '-')
        version = join(service_and_version[end-2:end], '-')

        return (service, version)
    catch e
        if e isa BoundsError
            throw(InvalidFileName("$filename is an invalid AWS JSON filename."))
        else
            rethrow()
        end
    end
end

"""
Get the required and optional parameters for a given operation.
"""
function _get_function_parameters(input::String, shapes::AbstractDict{String})
    """
        _get_parameter_name(parameter::String, input_shape::AbstractDict{String, <:Any})

    Find the correct parameter name for making requests. Certain ones have a specific locationName, for Batch requests this locationName is nested one shape deeper.

    # Arguments
    - `parameter::String`: Name of the original parameter
    - `input_shape::AbstractDict{String, <:Any}`: The parameter shape

    # Returns
    - `String`: Either the original parameter name, the locationName for the parameter, or the locationName nested one shape deeper
    """
    function _get_parameter_name(parameter::String, input_shape::AbstractDict{String, <:Any})
        # If the parameter has a locationName, return it
        if haskey(input_shape["members"][parameter], "locationName")
            return input_shape["members"][parameter]["locationName"]
        end

        # Check to see if the nested shape has a locationName
        nested_shape = input_shape["members"][parameter]["shape"]
        nested_member = get(shapes[nested_shape], "member", Dict())

        # If nested_shape[member] exists return locationName (if exists), otherwise return the original parameter name
        return get(nested_member, "locationName", parameter)
    end

    required_parameters = LittleDict{String, Any}()
    optional_parameters = LittleDict{String, Any}()
    input_shape = shapes[input]

    if haskey(input_shape, "required")
        for parameter in input_shape["required"]
            parameter_name = _get_parameter_name(parameter, input_shape)

            # Check if the parameter needs to be in a certain place
            parameter_location = get(input_shape["members"][parameter], "location", "")

            documentation = _clean_documentation(get(input_shape["members"][parameter], "documentation", ""))

            required_parameters[parameter_name] = LittleDict{String, String}(
                "location" => parameter_location,
                "documentation" => documentation,
            )
        end
    end

    if haskey(input_shape, "members")
        for (member_key, member_value) in input_shape["members"]
            parameter_name = get(input_shape["members"][member_key], "locationName", member_key)

            if !haskey(required_parameters, parameter_name)
                documentation = _clean_documentation(get(member_value, "documentation", ""))
                idempotent = get(member_value, "idempotencyToken", false)

                optional_parameters[parameter_name] = LittleDict{String, Union{String, Bool}}(
                    "documentation" => documentation,
                    "idempotent" => idempotent,
                )
            end
        end
    end

    return (sort(required_parameters), sort(optional_parameters))
end


text:
### subroutines to read data

"""
    read_plink_data(xpath) -> x::BEDFile, y::SharedVector

Read the trio of binary PLINK files (BED, BIM, FAM) from `xpath`. The response `y` is loaded from the rightmost column of the FAM file.
Unlike the `BEDFile` constructors, the variable `xpath` follows the PLINK convention in that it should point to the trio of files.
For example,

    read_plink_data("mydata")

would read the trio `mydata.bed`, `mydata.bim`, and `mydata.fam` in the current working directory.
"""
function read_plink_data(
    T      :: Type,
    xpath  :: String;
    pids   :: DenseVector{Int} = procs(),
    header :: Bool = false,
    delim  :: Char = ' '
)

    # load genotype data
    x = BEDFile(T, xpath * ".bed", pids=pids)

    # read the FAM file
    famfile = xpath * ".fam"
    Y = readdlm(famfile, delim, header=header)

    # check that the FAM file has six columns
    p = size(Y,2)
    @assert p == 6 "FAM file does not have six columns, is it formatted correctly?"

    # we cannot know for certain that Y is loaded as floating point
    # to be safe, explicitly convert the phenotype column to the correct type
    ytemp = convert(Vector{T}, Y[:,end])

    # in FAM file, the phenotype is the rightmost column
    # initialize a SharedVector and fill it with the phenotype
    y = SharedArray{T}((x.geno.n,), pids=pids)
    copy!(y, ytemp)

    return x, y
end

# default for previous function is Float64
read_plink_data(xpath::String; pids::DenseVector{Int} = procs(), header::Bool = false, delim::Char = ' ') = read_plink_data(Float64, xpath, pids=pids, header=header, delim=delim)

"""
    read_plink_data(xpath, ypath [,isbin=false, header=false, delim=' ']) -> x::BEDFile, y::SharedVector

Read the trio of binary PLINK files (BED, BIM, FAM) based on the file path `xpath`.
By default, the response `y` is loaded from a delimited text file located at `ypath`.
Set `isbin=true` to load from a _binary_ file.
"""
function read_plink_data(
    T      :: Type,
    xpath  :: String,
    ypath  :: String;
    pids   :: DenseVector{Int} = procs(),
    isbin  :: Bool = false,
    header :: Bool = false,
    delim  :: Char = ' '
)

    # load genotype data
    x = BEDFile(T, xpath * ".bed", pids=pids)

    # load phenotype from file
    # binary files -> make SharedVector directly
    # delimited text -> must wrangle into SharedVector
    if isbin
        y = SharedArray{T}(abspath(ypath), (x.geno.n,), pids=pids)
    else
        ytemp = readdlm(ypath, delim, header=header)
        @assert size(ytemp,2) == 1 "File at ypath must have exactly one column"
        y = SharedArray{T}((x.geno.n,), pids=pids)
        copy!(y, ytemp)
    end

    return x, y
end

# default for previous function is Float64
read_plink_data(xpath::String, ypath::String, pids::DenseVector{Int} = procs(), isbin::Bool = false, header::Bool = false, delim::Char = ' ') = read_plink_data(Float64, xpath, ypath, pids=pids, isbin=isbin, header=header, delim=delim)

"""
    read_plink_data(xpath, covpath, ypath [,isbin=false, header=false, delim=' ']) -> x::BEDFile, y::SharedVector

Read the trio of binary PLINK files (BED, BIM, FAM) based on the file path `xpath`, plus covariates from the file path `covpath`.
By default, the response `y` is loaded from a delimited text file located at `ypath`.
Set `isbin=true` to load from a _binary_ file.
"""
function read_plink_data(
    T         :: Type,
    xpath     :: String,
    covpath   :: String,
    ypath     :: String;
    pids      :: DenseVector{Int} = procs(),
    covheader :: Bool = false,
    covdelim  :: Char = ' ',
    isbin     :: Bool = false,
    yheader   :: Bool = false,
    ydelim    :: Char = ' ',
)

    # load genotype data
    x = BEDFile(T, xpath * ".bed", covpath, pids=pids, header=covheader)

    # load phenotype from file
    # binary files -> make SharedVector directly
    # delimited text -> must wrangle into SharedVector
    if isbin
        y = SharedArray{T}(abspath(ypath), (x.geno.n,), pids=pids)
    else
        ytemp = readdlm(ypath, ydelim, header=yheader)
        @assert size(ytemp,2) == 1 "File at ypath must have exactly one column"
        y = SharedArray{T}((x.geno.n,), pids=pids)
        copy!(y, ytemp)
    end

    return x, y
end

# default for previous function is Float64
read_plink_data(xpath::String, covpath::String, ypath::String, pids::DenseVector{Int} = procs(), covheader::Bool = false, covdelim::Char = ' ', isbin::Bool = false, yheader::Bool = false, ydelim::Char = ' ' ) = read_plink_data(Float64, xpath, covpath, ypath, pids=pids, covheader=covheader, covdelim=covdevlim, isbin=isbin, yheader=yheader, ydelim=ydelim)

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
struct HeteroscedasticLikelihood{T<:Real} <: RegressionLikelihood{T}

end


text:
using LinearAlgebra, SparseArrays, Test

using MultivariatePolynomials
using SumOfSquares

# Taken from JuMP/test/solvers.jl
function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

if try_import(:DynamicPolynomials)
    import DynamicPolynomials.@polyvar
    import DynamicPolynomials.@ncpolyvar
else
    if try_import(:TypedPolynomials)
        import TypedPolynomials.@polyvar
    else
        error("No polynomial implementation installed : Please install TypedPolynomials or DynamicPolynomials")
    end
end

include("certificate.jl")

include("gram_matrix.jl")

include("variable.jl")
include("constraint.jl")

include("Mock/mock_tests.jl")

# Tests needing a solver
# FIXME these tests should be converted to Literate and moved to `examples` or
# converted to be used with `MockOptimizer` and moved to `test/Tests`
include("solvers.jl")
include("sospoly.jl")
include("sosquartic.jl")
include("equalitypolyconstr.jl")


text:
using Gridap.Algebra: residual
using Gridap.Algebra: jacobian
import Gridap.Algebra: NonlinearSolver
import Gridap.Algebra: NonlinearOperator
import Gridap.Algebra: solve!
import Gridap.ODEs.ODETools: solve_step!
import Gridap.ODEs.ODETools: ODESolver
import Gridap.ODEs.ODETools: zero_initial_guess
import Gridap.ODEs.ODETools: residual!
import Gridap.ODEs.ODETools: jacobian!
import Gridap.ODEs.ODETools: solve!
import Gridap.ODEs.ODETools: allocate_residual
import Gridap.ODEs.ODETools: allocate_jacobian

struct OperatorMock <: NonlinearOperator
  odeop
  tf::Float64
  dt::Float64
  u0::AbstractVector
  cache
end

function OperatorMock(odeop::ODEOperator,tf::Real,dt::Real,u0::AbstractVector)
  cache = nothing
  OperatorMock(odeop,tf,dt,u0,cache)
end

function residual!(b::AbstractVector,op::OperatorMock,x::AbstractVector)
  uf = x
  uf_t = (x-op.u0)/op.dt
  residual!(b,op.odeop,op.tf,(uf,uf_t),op.cache)
end

function jacobian!(A::AbstractMatrix,op::OperatorMock,x::AbstractVector)
  uf = x
  uf_t = (x-op.u0)/op.dt
  fill!(A,0.0)
  # jacobian!(A,op.odeop,op.tf,(uf,uf_t),1,1.0,op.cache)
  # jacobians!(A,op.odeop,op.tf,(uf,uf_t),(1/op.dt),op.cache)
  jacobians!(A,op.odeop,op.tf,(uf,uf_t),(1.0,1.0/op.dt),op.cache)
end

function allocate_residual(op::OperatorMock,x::AbstractVector)
  allocate_residual(op.odeop,x,op.cache)
end

function allocate_jacobian(op::OperatorMock,x::AbstractVector)
  allocate_jacobian(op.odeop,x,op.cache)
end

function zero_initial_guess(op::OperatorMock)
  x0 = similar(op.u0)
  fill!(x0,zero(eltype(x0)))
  x0
end

struct NLSolverMock <: NonlinearSolver
end

function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cache::Nothing)
  r = residual(nlop,x)
  J = jacobian(nlop,x)
  dx = inv(Matrix(J))*(-r)
  x.= x.+dx
  cache = (r,J,dx)
end

function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cache)
  r, J, dx = cache
  residual!(r, nlop, x)
  jacobian!(J, nlop, x)
  dx = inv(Matrix(J))*(-r)
  x.= x.+dx
end

struct ODESolverMock <: ODESolver
  nls::NLSolverMock
  dt::Float64
end

function solve_step!(
  uf::AbstractVector,solver::ODESolverMock,op::ODEOperator,u0::AbstractVector,t0::Real, cache) # -> (uF,tF)

  dt = solver.dt
  tf = t0+dt
  if (cache == nothing)
    ode_cache = allocate_cache(op)
  else
    ode_cache, nl_cache = cache
    update_cache!(ode_cache,op,tf)
  end

  nlop = OperatorMock(op,tf,dt,u0,ode_cache)

  if (cache==nothing)
    nl_cache = solve!(uf,solver.nls,nlop)
  else
    nl_cache = solve!(uf,solver.nls,nlop,nl_cache)
  end

  cache = ode_cache, nl_cache

  return (uf, tf, cache)
end


text:
# Autogenerated wrapper script for Lz4_jll for x86_64-w64-mingw32
export liblz4, lz4, lz4c, lz4cat, unlz4

JLLWrappers.@generate_wrapper_header("Lz4")
JLLWrappers.@declare_library_product(liblz4, "msys-lz4-1.dll")
JLLWrappers.@declare_executable_product(lz4)
JLLWrappers.@declare_executable_product(lz4c)
JLLWrappers.@declare_executable_product(lz4cat)
JLLWrappers.@declare_executable_product(unlz4)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liblz4,
        "bin\\msys-lz4-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lz4,
        "bin\\lz4.exe",
    )

    JLLWrappers.@init_executable_product(
        lz4c,
        "bin\\lz4c.exe",
    )

    JLLWrappers.@init_executable_product(
        lz4cat,
        "bin\\lz4cat.exe",
    )

    JLLWrappers.@init_executable_product(
        unlz4,
        "bin\\unlz4.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Random, IntervalLinearAlgebra, Plots

sizes = Iterators.flatten((2:9, 10:10:90, 100:100:1000)) |> collect

times = zeros(length(sizes), 4)

seed = 42
Random.seed!(seed)

random_interval_matrix(N) = Interval.(randn(N, N) .± abs.(randn(N, N)))

for (j, mult_mode) in enumerate((:fast, :rank1, :slow))

    set_multiplication_mode(mult_mode)
    Random.seed!(seed)
    for (i, N) in enumerate(sizes)
        A = random_interval_matrix(N)
        B = random_interval_matrix(N)
        t = @benchmark $A * $B
        times[i, j] = minimum(t.times)
        @show N
    end

end

# normal multiplication with accurate rounding mode
set_multiplication_mode(:slow)
setrounding(Interval, :accurate)
Random.seed!(seed)
for (i, N) in enumerate(sizes)
    A = random_interval_matrix(N)
    B = random_interval_matrix(N)
    t = @benchmark $A * $B
    times[i, 4] = minimum(t.times)
    @show N
end

# back to default settings
setrounding(Interval, :tight)
set_multiplication_mode(:fast)

# floating point multiplication
times_float = zeros(size(sizes))
for (i, N) in enumerate(sizes)
    A = randn(N, N)
    B = randn(N, N)
    t = @benchmark $A * $B
    times_float[i] = minimum(t.times)
    @show N
end

# plotting
labels = ["fast" "rank1" "slow-tight" "slow-accurate"]
plot(sizes, times/1e6; label=labels, axis=:log, m=:auto, legend=:topleft)
plot!(sizes, times_float/1e6, label="Float64", m=:auto)
xlabel!("size")
ylabel!("time [ms]")
yticks!(10.0 .^(-3:5))

ratios = times ./ times_float
labels = ["fast" "rank1" "slow-tight" "slow-accurate"]
plot(sizes, ratios; label=labels, axis=:log, m=:auto, legend=:topleft)
xlabel!("size")
ylabel!("ratio")
yticks!(10.0 .^(-3:5))

## overestimate benchmarks

radii = 10.0 .^ (-15:15)
sizes = [10, 50, 100, 200]
Nsamples = 10

min_overestimate = zeros(length(radii), length(sizes))
max_overestimate = zeros(length(radii), length(sizes))
mean_overestimate = zeros(length(radii), length(sizes))

# function random_interval_matrix(N, r)
#     Ac = randn(N, N)
#     return Ac .± (r * abs.(Ac))
# end

random_interval_matrix(N, r) = randn(N, N) .± r*abs.(randn(N, N))
Random.seed!(seed)
for (i, N) in enumerate(sizes)
    for (j, r) in enumerate(radii)

        tmp_statistics = zeros(3)
        for _ in 1:Nsamples
            A = random_interval_matrix(N, r)
            B = random_interval_matrix(N, r)

            set_multiplication_mode(:fast)
            Cfast = A * B
            set_multiplication_mode(:slow)
            Cslow = A * B

            ratios = (diam.(Cfast)./diam.(Cslow) .- 1) * 100
            tmp_statistics .+= (minimum(ratios), mean(ratios), maximum(ratios))
        end

        min_overestimate[j, i] = tmp_statistics[1]/Nsamples
        mean_overestimate[j, i] = tmp_statistics[2]/Nsamples
        max_overestimate[j, i] = tmp_statistics[3]/Nsamples
        println("N=$(i)/$(length(sizes)) r=$(j)/$(length(radii))")
    end
end

for (i, N) in enumerate(sizes)
    plot(radii, max_overestimate[:, i], xaxis=:log, m=:o, label="max")
    plot!(radii, min_overestimate[:, i], m=:o, label="min")
    plot!(radii, mean_overestimate[:, i], m=:o, label="mean")
    xlabel!("radius")
    ylabel!("%-overestimate")
    title!("matrix size N=$N")
    savefig("accurate_overestimate_N$(N)_random_radius.pdf")
end


text:
export  dimension

"""
    dimension( vertices :: Vector{Vector{Int64}} ) :: Int64

For the provided vertices (points), finds the dimension of the affine space spanned
by the vertices. This method computes the dimension of a polytope, facet, or
collection of points.

This also accepts matrices and `DeterministicStrategy` types as arguments:

    dimension( vertices :: Vector{Matrix{Int64}} ) :: Int64

    dimension( vertices :: Vector{DeterministicStrategy} ) :: Int64
"""
function dimension(vertices :: Vector{Vector{Int64}}) :: Int64
    vectors = map( v -> v - vertices[1], vertices[2:end])
    rank(hcat(vectors...))
end

function dimension(vertices :: Vector{Matrix{Int64}}) :: Int64
    dimension(map( v -> v[:], vertices))
end

function dimension(vertices :: Vector{DeterministicStrategy}) :: Int64
    dimension(map(v -> v[:], vertices))
end


text:
#!/usr/bin/env julia
####
#### Janzen-Connell Model
####
#### A simple forest ecosystem model to investigate the effect of pathogen-
#### induced density-dependent mortality on tree diversity.
####
#### (c) 2020 Daniel Vedder <daniel.vedder@stud-mail.uni-wuerzburg.de>
####     Licensed under the terms of the MIT license.
####

module jcm

const jcm_version = v"1.0"

using Dates,
    ArgParse,
    Logging,
    Random

const settings = Dict("species" => 16,              # The number of species that will be created
                      "worldsize" => 500,           # The extent from the center of the square world arena in meters
                      "runtime" => 500,             # The number of updates the simulation will run
                      "datafile" => "jcm_data.csv", # The name of the recorded data file
                      "datafreq" => 10,             # How long between data recordings?
                      "pathogens" => false,         # Include pathogens in the simulation?
                      "transmission" => 50,         # Pathogen infection radius
                      "neutral" => false,           # All species have identical trait values?
                      "verbosity" => "Info",        # The log level (Debug, Info, Warn, Error)
                      "seed" => 0)                  # The seed for the RNG (0 -> random)

include("trees.jl")
include("forest.jl")

"""
Set parameters from the commandline.
"""
function parsecommandline()
    s = ArgParseSettings(prog="jcm.jl", version="$jcm_version", add_help=true, add_version=true,
                         description="Investigate the Janzen-Connell effect in a forest model.",
                         epilog="JCM $jcm_version, (c) 2020 Daniel Vedder \n\nEcosystem Modelling Group, University of Würzburg")
    @add_arg_table! s begin
        "--species", "-s"
        help = "the number of tree species to simulate"
        arg_type = Int
        default = settings["species"]
        "--worldsize", "-w"
        help = "the extent from the center of the square world arena in meters"
        arg_type = Int
        default = settings["worldsize"]
        "--runtime", "-t"
        help = "the number of updates the simulation will run"
        arg_type = Int
        default = settings["runtime"]
        "--datafile", "-f"
        help = "the file to which simulation data is written"
        arg_type = String
        default = settings["datafile"]
        "--datafreq", "-d"
        help = "the frequency in updates with which data is recorded"
        arg_type = Int
        default = settings["datafreq"]
        "--pathogens", "-p"
        help = "run a simulation with pathogens"
        action = :store_true
        "--transmission", "-i"
        help = "maximum pathogen infection radius"
        arg_type = Int
        default = settings["transmission"]
        "--neutral", "-n"
        help = "all species have identical trait values"
        action = :store_true
        "--verbosity", "-v"
        help = "set the log level (Debug, Info, Warn, Error)"
        arg_type = String
        default = settings["verbosity"]
        "--seed", "-r"
        help = "set the seed for the RNG (0 -> random seed)"
        arg_type = Int
        default = settings["seed"]
    end
    return parse_args(s)
end

"""
Initialise the random number generator.
"""
function initRNG()
    if settings["seed"] == 0
        settings["seed"] = abs(rand(Random.RandomDevice(), Int32))
    end
    Random.seed!(settings["seed"])
end

"""
Initialise the world with one mature tree from each species at a random location.
"""
function initworld()
    createspecies()
    for n in 1:settings["species"]
        ws = settings["worldsize"]
        xpos = rand(-ws:ws)
        ypos = rand(-ws:ws)
        sp = getspecies(n)
        settings["pathogens"] ? infection = Pathogen(n) : infection = nothing
        tree = Tree(nextid(), sp, convert(Int16, round(sp.max_age/2)),
                    sp.max_size, true, infection, (x=xpos, y=ypos))
        recordindividual(tree)
        planttree!(tree)
    end
end

let updatelog::String = "", update::Int=0
    """
    Set the internal update counter to the current update number.
    (Needed by run())
    """
    global function setupdate(t)
        update = t
    end

    """
    Save the data records of one individual in this turn.
    Called during `grow!()`
    """
    global function recordindividual(tree::Tree)
        update < 0 && return # block unless we've reached a recording point
        datastring = "$update,$(tree.species.id),$(tree.uid),$(tree.age),$(tree.size)"
        datastring *= ",$(tree.mature ? "TRUE" : "FALSE"),$(tree.infection == nothing ? "FALSE" : "TRUE")"
        datastring *= ",$(tree.position.x),$(tree.position.y)"
        updatelog *= datastring * "\n"
    end

    """
    Write out a csv file with all necessary analysis data.
    """
    global function recorddata()
        open(settings["datafile"], "a") do df
            print(df,updatelog)
        end
        updatelog = ""
        update = -1 # block until the next recorded update
    end

    """
    Initialise the data file and the logger.
    """
    global function initrecording()
        global_logger(ConsoleLogger(stdout, eval(
            Meta.parse("Logging.$(settings["verbosity"])"))))
        open(settings["datafile"], "w") do df
            time = Dates.format(Dates.now(), "dd/mm/yyy HH:MM")
            println(df, "# Janzen-Connell Model data file, created $time")
            println(df, "# settings = $settings")
            println(df, "update,species,uid,age,size,mature,infected,x,y")
        end
    end
end

"""
The main simulation function.
"""
function run(noinit::Bool=false)
    merge!(settings, parsecommandline())
    initRNG()
    initrecording()
    if !noinit #intended for use by test functions
        @info "Initialising the world"
        initworld()
    end
    recorddata()
    for u in 1:settings["runtime"]
        @info "UPDATE $u"
        record = u%settings["datafreq"] == 0
        record && setupdate(u)
        @info "Dispersal"
        disperse!()
        @info "Competition"
        compete!()
        if settings["pathogens"]
            @info "Infection"
            infect!()
        end
        @info "Growth"
        grow!() # includes data recording
        record && recorddata()
    end
end

if !isinteractive()
    @time run()
end

end

max_stars_repo_name: veddox/jcm

max_issues_repo_name: veddox/jcm

max_forks_repo_name: veddox/jcm


text:
const TEST_FLOAT32 = false
const TEST_DOUBLE64 = false
const TEST_STATIC_SIZED = false

using Manifolds
using ManifoldsBase
using ManifoldsBase: number_of_coordinates

using LinearAlgebra
using Distributions
using DoubleFloats
using ForwardDiff
using Quaternions
using Random
using ReverseDiff
using StaticArrays
using Statistics
using StatsBase
using Test
using LightGraphs
using SimpleWeightedGraphs

function include_test(path)
    @info "Testing $path"
    @time include(path)  # show basic timing, (this will print a newline at end)
end

function our_base_ambiguities()
    ambigs = Test.detect_ambiguities(Base)
    modules_we_care_about =
        [Base, LinearAlgebra, Manifolds, ManifoldsBase, StaticArrays, Statistics, StatsBase]
    our_ambigs = filter(ambigs) do (m1, m2)
        we_care = m1.module in modules_we_care_about && m2.module in modules_we_care_about
        return we_care && (m1.module === Manifolds || m2.module === Manifolds)
    end
    return our_ambigs
end


text:
let
  state = 0
  global counter
  counter() = state += 1
end;

println(counter())
println(counter())
println(counter())

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox