text:
using SparseArrays
using SuiteSparse.CHOLMOD

"""
    Cholmod <: LSBackend

Use CHOLMOD backend.

Options available:
* `Float64` only
* Augmented system with LDLᵀ factorization
* Normal equations with Cholesky factorization
"""
struct Cholmod <: LSBackend end

# ==============================================================================
#   SparseIndefLinearSolver
# ==============================================================================

"""
    SparseIndefLinearSolver

Linear solver for the 2x2 augmented system with ``A`` sparse.

Uses an LDLᵀ factorization of the quasi-definite augmented system.
"""
mutable struct SparseIndefLinearSolver <: AbstractLinearSolver{Float64}
    m::Int  # Number of rows
    n::Int  # Number of columns

    # TODO: allow for user-provided ordering,
    #   and add flag (default to false) to know whether user ordering should be used

    # Problem data
    A::SparseMatrixCSC{Float64}
    θ::Vector{Float64}
    regP::Vector{Float64}  # primal-dual regularization
    regD::Vector{Float64}  # dual regularization

    # Left-hand side matrix
    S::SparseMatrixCSC{Float64, Int}  # TODO: use `CHOLMOD.Sparse` instead

    # Factorization
    F::CHOLMOD.Factor{Float64}

    # TODO: constructor with initial memory allocation
    function SparseIndefLinearSolver(A::SparseMatrixCSC{Float64, Int})
        m, n = size(A)
        θ = ones(Float64, n)

        S = [
            spdiagm(0 => -θ)  A';
            spzeros(Float64, m, n) spdiagm(0 => ones(m))
        ]

        # Symbolic factorization
        # See code in src/cholmod.jl of SuiteSparse.jl for `ldlt` function
        F = ldlt(Symmetric(S))

        return new(m, n, A, θ, ones(Float64, n), ones(Float64, m), S, F)
    end

end

AbstractLinearSolver(
    ::Cholmod,
    ::AugmentedSystem,
    A::AbstractMatrix{Float64}
) = SparseIndefLinearSolver(sparse(A))

AbstractLinearSolver(
    ::Cholmod,
    ::DefaultSystem,
    A::AbstractMatrix{Float64}
) = SparseIndefLinearSolver(sparse(A))

backend(::SparseIndefLinearSolver) = "CHOLMOD"
linear_system(::SparseIndefLinearSolver) = "Augmented system"

"""
    update_linear_solver!(ls, θ, regP, regD)

Update LDLᵀ factorization of the augmented system.

Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute
    the factorization.
Throws a `PosDefException` if matrix is not quasi-definite.
"""
function update_linear_solver!(
    ls::SparseIndefLinearSolver,
    θ::AbstractVector{Float64},
    regP::AbstractVector{Float64},
    regD::AbstractVector{Float64}
)
    # Sanity checks
    length(θ)  == ls.n || throw(DimensionMismatch(
        "θ has length $(length(θ)) but linear solver is for n=$(ls.n)."
    ))
    length(regP) == ls.n || throw(DimensionMismatch(
        "regP has length $(length(regP)) but linear solver has n=$(ls.n)"
    ))
    length(regD) == ls.m || throw(DimensionMismatch(
        "regD has length $(length(regD)) but linear solver has m=$(ls.m)"
    ))

    # Update diagonal scaling
    ls.θ .= θ
    # Update regularizers
    ls.regP .= regP
    ls.regD .= regD

    # Update S. S is stored as upper-triangular and only its diagonal changes.
    @inbounds for j in 1:ls.n
        k = ls.S.colptr[1+j] - 1
        ls.S.nzval[k] = -ls.θ[j] - regP[j]
    end
    @inbounds for i in 1:ls.m
        k = ls.S.colptr[1+ls.n+i] - 1
        ls.S.nzval[k] = regD[i]
    end

    # `S` is first converted into CHOLMOD's internal data structure,
    # so this line allocates.
    ldlt!(ls.F, Symmetric(ls.S))

    return nothing
end

"""
    solve_augmented_system!(dx, dy, ls, ξp, ξd)

Solve the augmented system, overwriting `dx, dy` with the result.
"""
function solve_augmented_system!(
    dx::Vector{Float64}, dy::Vector{Float64},
    ls::SparseIndefLinearSolver,
    ξp::Vector{Float64}, ξd::Vector{Float64}
)
    m, n = ls.m, ls.n

    # Set-up right-hand side
    ξ = [ξd; ξp]

    # Solve augmented system
    d = ls.F \ ξ

    # Recover dx, dy
    @views dx .= d[1:n]
    @views dy .= d[(n+1):(m+n)]

    # Iterative refinement
    # TODO:
    # * Max number of refine steps
    # * Check for residuals before refining
    # * Check whether residuals did improve
    # resP = ls.A*dx + ls.regD .* dy - ξp
    # resD = - dx .* (ls.θ + ls.regP) + ls.A' * dy - ξd
    # println("\n|resP| = $(norm(resP, Inf))\n|resD| = $(norm(resD, Inf))")

    # ξ1 = [resD; resP]
    # d1 = ls.F \ ξ1

    # # Update search direction
    # @views dx .-= d1[1:n]
    # @views dy .-= d1[(n+1):(m+n)]

    return nothing
end

# ==============================================================================
#   SparsePosDefLinearSolver
# ==============================================================================

"""
    SparsePosDefLinearSolver

Linear solver for the 2x2 augmented system
```
    [-(Θ⁻¹ + Rp)   Aᵀ] [dx] = [xi_d]
    [   A          Rd] [dy]   [xi_p]
```
with ``A`` sparse.

Uses a Cholesky factorization of the positive definite normal equations system
```
(A * ((Θ⁻¹ + Rp)⁻¹ * Aᵀ + Rd) dy = xi_p + A * (θ⁻¹ + Rp)⁻¹ * xi_d
                              dx = (Θ⁻¹ + Rp)⁻¹ * (Aᵀ * dy - xi_d)
```
"""
mutable struct SparsePosDefLinearSolver <: AbstractLinearSolver{Float64}
    m::Int  # Number of rows
    n::Int  # Number of columns

    # TODO: allow for user-provided ordering,
    #   and add flag (default to false) to know whether user ordering should be used

    # Problem data
    A::SparseMatrixCSC{Float64, Int}
    θ::Vector{Float64}   # Diagonal scaling
    # Regularization
    regP::Vector{Float64}  # primal
    regD::Vector{Float64}  # dual

    # Factorization
    F::CHOLMOD.Factor{Float64}

    # Constructor and initial memory allocation
    # TODO: symbolic only + allocation
    function SparsePosDefLinearSolver(A::SparseMatrixCSC{Float64, Int})
        m, n = size(A)
        θ = ones(Float64, n)

        S = A * A' + spdiagm(0 => ones(m))

        # TODO: PSD-ness checks
        F = cholesky(Symmetric(S))

        return new(m, n, A, θ, zeros(Float64, n), ones(Float64, m), F)
    end

end

AbstractLinearSolver(
    ::Cholmod,
    ::NormalEquations,
    A::AbstractMatrix{Float64}
) = SparsePosDefLinearSolver(sparse(A))

backend(::SparsePosDefLinearSolver) = "CHOLMOD - Cholesky"
linear_system(::SparsePosDefLinearSolver) = "Normal equations"

"""
    update_linear_solver!(ls, θ, regP, regD)

Compute normal equation system matrix, and update the factorization.
"""
function update_linear_solver!(
    ls::SparsePosDefLinearSolver,
    θ::AbstractVector{Float64},
    regP::AbstractVector{Float64},
    regD::AbstractVector{Float64}
)

    # Sanity checks
    length(θ) == ls.n || throw(DimensionMismatch(
        "θ has length $(length(θ)) but linear solver is for n=$(ls.n)."
    ))
    length(regP) == ls.n || throw(DimensionMismatch(
        "regP has length $(length(regP)) but linear solver has n=$(ls.n)"
    ))
    length(regD) == ls.m || throw(DimensionMismatch(
        "regD has length $(length(regD)) but linear solver has m=$(ls.m)"
    ))

    ls.θ .= θ
    ls.regP .= regP  # Primal regularization is disabled for normal equations
    ls.regD .= regD

    # Re-compute factorization
    # D = (Θ^{-1} + Rp)^{-1}
    D = Diagonal(one(Float64) ./ (ls.θ .+ ls.regP))
    Rd = spdiagm(0 => ls.regD)
    S = ls.A * D * ls.A' + Rd

    # Update factorization
    cholesky!(ls.F, Symmetric(S), check=false)
    issuccess(ls.F) || throw(PosDefException(0))

    return nothing
end

"""
    solve_augmented_system!(dx, dy, ls, ξp, ξd)

Solve the augmented system, overwriting `dx, dy` with the result.
"""
function solve_augmented_system!(
    dx::Vector{Float64}, dy::Vector{Float64},
    ls::SparsePosDefLinearSolver,
    ξp::Vector{Float64}, ξd::Vector{Float64}
)
    m, n = ls.m, ls.n

    d = one(Float64) ./ (ls.θ .+ ls.regP)
    D = Diagonal(d)

    # Set-up right-hand side
    ξ_ = ξp .+ ls.A * (D * ξd)

    # Solve augmented system
    dy .= (ls.F \ ξ_)

    # Recover dx
    dx .= D * (ls.A' * dy - ξd)

    # TODO: Iterative refinement
    # * Max number of refine steps
    # * Check for residuals before refining
    # * Check whether residuals did improve
    # resP = ls.A * dx + ls.regD .* dy - ξp
    # resD = - D \ dx + ls.A' * dy - ξd
    # println("\n|resP| = $(norm(resP, Inf))\n|resD| = $(norm(resD, Inf))")

    return nothing
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Threads

include("threadingconstructs.jl")
include("atomics.jl")
include("locks.jl")

end

ext: jl


text:
abstract type EBayesShrinkageLocation end

#-----------------------------------------------
# FixedLocation
#-----------------------------------------------

struct FixedLocation{TS} <: EBayesShrinkageLocation
    offset::TS
end
FixedLocation(; offset=0.0) = FixedLocation(offset)

fit(fixed_loc::FixedLocation, ss) = nothing

function predict(fixed_loc::FixedLocation, ss::NormalSamples)
    zeros(ss) .+ fixed_loc.offset
end

# fitted_loc will be nothing
function predict(fixed_loc::FixedLocation, fitted_loc, ss::NormalSamples)
    predict(fixed_loc, ss)
end

function leverage(fixed_loc::FixedLocation, fitted_loc, ss)
    zero(eltype(response(ss)))
end

#-----------------------------------------------
# GrandMeanLocation
#-----------------------------------------------

struct GrandMeanLocation{TS} <: EBayesShrinkageLocation
    offset::TS
end

GrandMeanLocation(; offset=0.0) = GrandMeanLocation(offset)

fit(loc::GrandMeanLocation, ss) = mean( response(ss) .- loc.offset)

function predict(loc::GrandMeanLocation, fitted_loc, ss::NormalSamples)
    zeros(ss) .+ loc.offset .+ fitted_loc
end

function leverage(loc::GrandMeanLocation, fitted_loc, ss)
    1/length(ss)
end

#----------------------------------------------------------
#  RegressionLocation
#----------------------------------------------------------

#struct RegressionLocation{TS, MT <: Supervised} <: EBayesShrinkageLocation
#    offset::TS
#    model::MT
#end

#RegressionLocation(model; offset=0.0) = FixedLocation(offset, model)

#fit(fixed_loc::FixedLocation, ss) = nothing

#function predict(fixed_loc::FixedLocation, ss::NormalSamples)
#    zeros(ss) .+ fixed_loc.offset
#end

# fitted_loc will be nothing
#function predict(fixed_loc::FixedLocation, fitted_loc, ss::NormalSamples)
#    predict(fixed_loc, ss)
#end

# GrandMeanLocation(μs) = GrandMeanLocation(μs, nothing)
# GrandMeanLocation() = GrandMeanLocation(0.0)
# GrandMeanLocation(::Type{AnalyticWeights}) = GrandMeanLocation(0.0, AnalyticWeights)
#
# weights(gm::GrandMeanLocation) = gm.ws
# dof(gm::GrandMeanLocation) = 1
#
# # TODO: Allow varying μ here as well
# function AnalyticWeights(s::NormalSamples)
#     AnalyticWeights(1 ./ s.σ.^2)
# end
#
# function fit(loc::GrandMeanLocation, s::NormalSamples, ::Type{AnalyticWeights})
#     mean(s.Z, AnalyticWeights(s))
# end
#
# function fit(loc::GrandMeanLocation,
#                   s::Union{HomoskedasticNormalSamples, StandardNormalSamples}, ws)
#     mean(s.Zs)
# end

# function dof_residual(ebloc::EBayesShrinkageLocation, fitted_loc, ss::EBayesSamples)
#     length(ss) - dof(ebloc, fitted_loc)
# end
# function dof_residual(ebloc::EBayesShrinkageLocation, ss::EBayesSamples)
#     length(ss) - dof(ebloc)
# end
#
# dof(ebloc::EBayesShrinkageLocation, fitted_loc) = dof(ebloc)
#
# function fit(loc::EBayesShrinkageLocation, ss::EBayesSamples)
#     fit(loc, ss, weights(loc))
# end
# function fit(loc::EBayesShrinkageLocation, ss::EBayesSamples, ws::Nothing)
#     fit(loc)
# end
#
# fit(loc::EBayesShrinkageLocation, X, ss::EBayesSamples) = fit(loc, ss)
#
# fit(loc::EBayesShrinkageLocation) = loc #by default return itself
# function leverage(loc::EBayesShrinkageLocation, ss::EBayesSamples)
#     leverage(loc)
# end
#
#
# weights(loc::EBayesShrinkageLocation) = nothing
#
#
#
# struct FixedLocation{T} <: EBayesShrinkageLocation
#     μs::T
# end
# FixedLocation() = FixedLocation(0.0)
#
# dof(fixed_loc::FixedLocation) = 0
# leverage(loc::FixedLocation) = 0 # don't need this to be a vector.
#
# function predict(loc::FixedLocation, ss::NormalSamples)
#     μs = loc.μs #fitted objects return the fit
#     preds = copy(ss.Z)
#     preds[:] .= μs
#     preds
# end #let's keep this
#
#
# struct SURELocation{T} <: EBayesShrinkageLocation
#     μs::T
# end
#
# SURELocation() = SURELocation(0.0)
#
# struct GrandMeanLocation{T, W} <: EBayesShrinkageLocation
#     μs::T
#     ws::W
# end
#
# GrandMeanLocation(μs) = GrandMeanLocation(μs, nothing)
# GrandMeanLocation() = GrandMeanLocation(0.0)
# GrandMeanLocation(::Type{AnalyticWeights}) = GrandMeanLocation(0.0, AnalyticWeights)
#
# weights(gm::GrandMeanLocation) = gm.ws
# dof(gm::GrandMeanLocation) = 1
#
# # TODO: Allow varying μ here as well
# function AnalyticWeights(s::NormalSamples)
#     AnalyticWeights(1 ./ s.σ.^2)
# end
#
# function fit(loc::GrandMeanLocation, s::NormalSamples, ::Type{AnalyticWeights})
#     mean(s.Z, AnalyticWeights(s))
# end
#
# function fit(loc::GrandMeanLocation,
#                   s::Union{HomoskedasticNormalSamples, StandardNormalSamples}, ws)
#     mean(s.Zs)
# end
#
# struct FittedGrandMeanLocation{T, GM<:GrandMeanLocation} <: EBayesShrinkageLocation
#     zbar::T
#     gm::GM
# end
#
# function fit(loc::GrandMeanLocation, ss::NormalSamples, ::Nothing)
#     # want to zero center (ss.Z - offset)
#     zbar = mean(ss.Z .- loc.μs)
#     FittedGrandMeanLocation(zbar, loc)
# end
#
# dof(fixed_loc::FittedGrandMeanLocation) = 1
# function leverage(loc::FittedGrandMeanLocation, ss::NormalSamples)
#     1/length(ss)
# end
#
# function predict(loc::FittedGrandMeanLocation, ss::NormalSamples)
#     μs = loc.zbar .+ loc.gm.μs #fitted objects return the fit
#     preds = copy(ss.Z)
#     preds[:] .= μs
#     preds
# end
#
#
#
#
# ########### Code for ...
#
# # could have the μs field as well, but don't need to...
# struct RegressionLocation{M , W, F<:StatsModels.FormulaTerm} <: EBayesShrinkageLocation
#     model::M
#     f::F
#     ws::W # M is the model
#     isfitted::Bool
# end
#
# RegressionLocation(model, f; isfitted=false) = RegressionLocation(model, f, nothing, isfitted)
# # maybe have default formula include everything + 1?
# RegressionLocation(model) = RegressionLocation(model, @formula(1~1))
#
# weights(reg::RegressionLocation) = reg.ws
# dof(reg::RegressionLocation, fitted_loc) = dof(fitted_loc)
# leverage(reg::RegressionLocation) = leverage(reg.model)
#
# #fit(loc, X, ss)
#
# function fit(loc::RegressionLocation, X, ss::NormalSamples)
#     if !loc.isfitted
#         fitted_model = fit(loc.model, X, ss, loc.f, loc.ws)
#         loc = RegressionLocation(fitted_model, loc.f, loc.ws, true)
#     end
#     loc
# end
#
# function fit(model, X_df::AbstractDataFrame, ss::NormalSamples, form, ::Nothing)
#     X = modelmatrix(form, X_df)
#     fitted_model = fit(model, X, ss.Z) #default...
#     fitted_model
# end
#
# # TODO: get weights working again...
# #function fit(loc::RegressionLocation, ssc::CovNormalSamples, ::Type{AnalyticWeights})
# #    X = modelmatrix(loc.f,
# #                    ssc.X_df)
#     # TODO: there will be some code duplication here, fix later
# #    fitted_model = fit(loc.model, X, ssc.S.Zs; wts=AnalyticWeights(ssc).values)
# #end
#
# function predict(loc::RegressionLocation, X, ss::NormalSamples) #with NormalSamples?
#     # TODO: go back to same fun for arbitrary Location
#     if !loc.isfitted
#         loc = fit(loc, X, ss)
#     end
#     predict(loc, X)
# end
#
# # for a fitted model
# function predict(loc::RegressionLocation, X) #with NormalSamples?
#     predict(loc.model, X)
# end
#
# function predict(loc::RegressionLocation, ss::NormalSamples) #with NormalSamples?
#     predict(loc.model)
# end
#
# predict(model, X) = predict(model)
#


text:
mutable struct DbId
    iid::String
end

function DbId()
    DbId("")
end

struct FieldName_Type
    fieldname::Symbol
    type::Type
end

const Build_Types = Union{<:AbstractAttribute, <:AbstractRelation}
const Direct_Storables = Union{<:AbstractString, <:Number, Bool, Date, DateTime}

max_stars_repo_name: Humans-of-Julia/TypeDBORM

max_issues_repo_name: Humans-of-Julia/TypeDBORM

max_forks_repo_name: Humans-of-Julia/TypeDBORM


text:
"""
    aperture(xrange, yrange, d, o=0)
    aperture(dom::CartesianDomain2D, d, o)

Create circular aperture in array `xrange × yrange` or given by `dom`ain with diameter `d` and center at `o`.

# Example
```jdoctests
julia> dom1= PhaseBases.CartesianDomain2D(-1:.2:1, -.9:.2:.9)
SampledDomains.CartesianDomain2D(-1.0:0.2:1.0, -0.9:0.2:0.9)

julia> ap, mask = PhaseBases.aperture(dom1, 1);

julia> ap
10×11 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  1  0  0  0  0  0
 0  0  0  1  1  1  1  1  0  0  0
 0  0  0  1  1  1  1  1  0  0  0
 0  0  0  1  1  1  1  1  0  0  0
 0  0  0  1  1  1  1  1  0  0  0
 0  0  0  0  0  1  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0
```
"""
function aperture(xrange::AbstractRange, yrange::AbstractRange, d, o=(0,0))
ap = [ ((xc-o[1])^2 + (yc-o[2])^2) <= d^2 /4 ? 1. : 0. for yc ∈ yrange,  xc ∈ xrange]
# area = +(ap[:]...)
phmask = [ ((xc-o[1])^2 + (yc-o[2])^2) <= d^2 /4 ? 1. : NaN for yc ∈ yrange,  xc ∈ xrange]
return(ap, phmask)
end

aperture(dom::CartesianDomain2D, d, o=(0,0)) = aperture(dom.xrange, dom.yrange, d,o)


text:
function main(P=nothing; p=1.2, nstep=100)

    if P == nothing
        P = Params(p=p, nstep=nstep)
    end
    @printf("BWH - vegetation patterns\n")
    @printf("-------------------------\n\n")
    @show P

   if P.freadinit
        @printf("Reading initial condition from file %s\n", P.initfile)
        a = readdlm(P.initfile)
        b = reshape(a[:,1], P.nx, P.ny)
        w = reshape(a[:,2], P.nx, P.ny)
    else
        @printf("Assigning random initial conditions\n")
        b = rand(P.nx, P.ny)*0.5;
        w = rand(P.nx, P.ny).*0.1.+0.5
    end

    @printf("Integrating:\n")
    b, w, h = integrate( b, w, P, P.nstep, fplot=P.fplot, fsave=P.fsave, nsave=P.nsave)

    u = reshape(cat(b, w, h, dims=3),P.nx*P.ny,3)
    @printf("Writing final file %s\n", P.finalfile)
    writedlm(P.finalfile, u, " ")
    return b, w, h, P
end

max_stars_repo_name: jhardenberg/bwh

max_issues_repo_name: jhardenberg/bwh

max_forks_repo_name: jhardenberg/bwh


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

# Generates a Linear Programming model for computing PageRank of a random graph.

# Sample usage:
# $ julia --project=. generate_pagerank_lp.jl --num_nodes 10000 \
#    --approx_num_edges 30000 --output_filename /tmp/pagerank.10k.mps.gz

# The LP formulation for PageRank is taken from:
# Nesterov, Y., "Subgradient methods for huge-scale optimization problems",
# Mathematical Programming, 2014. https://doi.org/10.1007/s10107-013-0686-4
# http://www.optimization-online.org/DB_FILE/2012/02/3339.pdf (preprint)

# Barabasi Albert preferential attachment graphs are used as a model of a
# random "web" graph.

import ArgParse
import JuMP
import LightGraphs

const SimpleGraph = LightGraphs.SimpleGraph

"""
Populates 'model' with the linear program model for computing the pagerank
vector for 'graph', with the given damping factor. Given an adjacency matrix A,
define S (the stochastic transition matrix) by dividing each column by it's
L_1 norm. The linear programming model is then:
  find x
  s.t. damping_factor * (S*x)[i] + (1 - damping_factor)/num_nodes <= x[i]
       sum_i x[i] = 1
       x >= 0
Note: Pagerank is usually computed for directed graphs. However, this example is
using the barabasi_albert generator, which generates an undirected graph
(SimpleGraph) instead of a directed graph (SimpleDiGraph).
"""
function populate_pagerank_model(
  model::JuMP.Model,
  graph::SimpleGraph,
  damping_factor::Float64 = 0.99,
)
  num_nodes = LightGraphs.nv(graph)
  degrees = [length(LightGraphs.neighbors(graph, i)) for i in 1:num_nodes]
  JuMP.@variable(model, x[i = 1:num_nodes], lower_bound = 0.0)
  for i in 1:num_nodes
    transition_sum = @expression(
      sum(x[j] / degrees[j] for j in LightGraphs.neighbors(graph, i))
    )
    JuMP.@constraint(
      model,
      damping_factor * transition_sum + (1 - damping_factor) / num_nodes <=
      x[i]
    )
  end
  JuMP.@constraint(
    model,
    L1_norm,
    sqrt(num_nodes) * sum(x[i] for i in 1:num_nodes) == sqrt(num_nodes)
  )
  return model
end

"""
Defines parses and args.

# Returns
A dictionary with the values of the command-line arguments.
"""
function parse_command_line()
  arg_parse = ArgParse.ArgParseSettings()

  ArgParse.@add_arg_table! arg_parse begin
    "--num_nodes"
    help = "The number of nodes in the generated graph."
    arg_type = Int
    required = true

    "--approx_num_edges"
    help = "The approximate number of edges in the generated graph."
    arg_type = Int
    required = true

    "--output_filename"
    help = "Filename for the output .mps (or .mps.gz) model."
    arg_type = String
    required = true

    "--damping_factor"
    help = "Damping factor for the pagerank model."
    arg_type = Float64
    default = 0.99

    "--random_seed"
    help = "Seed for the random generator."
    arg_type = Int
    default = 0
  end

  return ArgParse.parse_args(arg_parse)
end

function main()
  parsed_args = parse_command_line()

  filename = parsed_args["output_filename"]
  model = JuMP.Model()
  num_nodes = parsed_args["num_nodes"]
  approx_num_edges = parsed_args["approx_num_edges"]
  degree = round(Int, approx_num_edges / num_nodes)
  graph = LightGraphs.SimpleGraphs.barabasi_albert(
    num_nodes,
    degree,
    seed = parsed_args["random_seed"],
  )
  populate_pagerank_model(model, graph, parsed_args["damping_factor"])
  JuMP.write_to_file(model, parsed_args["output_filename"])
end

main()


text:
import GRAYTR: lerp, quadratic

@testset "Utility functions" begin

@testset "Linear extrapolation/interpolation" begin
    @test lerp(0, 0, 1) == 0
    @test lerp(1, 0, 1) == 1
    @test lerp(0.0, 0, 1) ≈ 0.0
    @test lerp(1.0, 0, 1) ≈ 1.0

    for t = [-0.2, 0.15, 0.5, 1.25]
        @test lerp(t, 0.0, 1.0) ≈ t
        @test lerp(t, 0.0, 2.0) ≈ 2.0 * t
        @test lerp(t, -1.0, 1.0) ≈ 2.0 * t - 1.0
    end
end

@testset "Quadratic solver" begin
    @testset "No solution" begin
        solvable, t1, t2 = quadratic(2.0, 1.0, 1.0)
        @test !solvable
        @test t1 == t2 == Inf
    end

    @testset "Golden ratio" begin
        solvable, t1, t2 = quadratic(1.0, -1.0, -1.0)
        @test solvable
        @test t1 ≈ (1 - sqrt(5)) / 2
        @test t2 ≈ (1 + sqrt(5)) / 2
    end
end

end # outer testset

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
transD_GP.SurfaceRegression.plot_surface_posterior(opt, optlog10λ, cmapmean="RdBu", property="moho depth",
            property_units="km", cmappdf="viridis")
## plot shapefile
ax = gcf().axes
using Shapefile, ArchGDAL, GeoInterface, GeoFormatTypes
t = Shapefile.Table("/home/547/ar0754/Downloads/STE_2016_AUST.shp")
for iaxis in 1:2
    for row in t
        if (row.STE_NAME16 == "Queensland") || (row.STE_NAME16 == "Northern Territory")
            xy = GeoInterface.coordinates(Shapefile.shape(row))
            for (i, points) in enumerate(xy)
                XY = permutedims(reverse(reduce(hcat, points[1]), dims=1))
                reprojXY = Vector{Array{Float64,1}}(undef, size(XY,1))
                for i in 1:size(XY, 1)
                        reprojXY[i] = XY[i,:]
                end
                XY = permutedims(reduce(hcat, ArchGDAL.reproject(reprojXY, EPSG(4283), EPSG(4462))))
                ax[iaxis].plot(XY[:,1], XY[:,2], "-k", linewidth=0.5)
            end
        end
    end
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
#= # for local testing, need this:
using Test
using PhyloNetworks
using PhyloPlots
using CSV
=#

@testset "auxiliary" begin

@testset "setlengths and setgammas" begin
originalstdout = stdout
redirect_stdout(open("/dev/null", "w")) # not portable to Windows
@test_nowarn PhyloNetworks.citation()
redirect_stdout(originalstdout)

str_level1_s = "(((S8,S9),((((S1,S4),(S5)#H1),(#H1,(S6,S7))))#H2),(#H2,S10));" # indviduals S1A S1B S1C go on leaf 1
net = readTopology(str_level1_s)

PhyloNetworks.setlengths!([net.edge[1]], [1.1])
@test net.edge[1].length == 1.1
PhyloNetworks.setlengths!([net.edge[3], net.edge[4]], [3.3, 4.4])
@test net.edge[3].length == 3.3
@test net.edge[4].length == 4.4

PhyloNetworks.setmultiplegammas!([net.edge[18]], [0.25])
@test net.edge[18].gamma == 0.25
@test net.edge[16].gamma == 0.75

@test PhyloNetworks.getlengths([net.edge[1]]) == [net.edge[1].length]
@test PhyloNetworks.getlengths([net.edge[1], net.edge[5]]) == [net.edge[1].length, net.edge[5].length]
end

@testset "hashybridladder" begin
tree = readTopology("(A:3.0,(B:2.0,(C:1.0,D:1.0):1.0):1.0);");
@test !PhyloNetworks.hashybridladder(tree)
PhyloNetworks.addhybridedge!(tree, tree.edge[5], tree.edge[1], true)
PhyloNetworks.addhybridedge!(tree, tree.edge[2], tree.edge[1], true)
@test PhyloNetworks.hashybridladder(tree)
end # of testing hashybridladder

@testset "shrink 2/3 cycles" begin
net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);")
@test !PhyloNetworks.shrink2cycles!(net)
@test !PhyloNetworks.shrink3cycles!(net)
PhyloNetworks.addhybridedge!(net, net.edge[7], net.edge[4], true, 0.1, 0.2)
@test !PhyloNetworks.shrink2cycles!(net)
@test PhyloNetworks.shrink3cycles!(net) # tree case
@test writeTopology(net) == "(((A:2.0,(B:1.0)#H1:0.1::0.9):1.37,(C:0.6,#H1:1.0::0.1):0.9):0.6,D:2.0);"
PhyloNetworks.addhybridedge!(net, net.edge[3], net.edge[6], true, 0.3, 0.4)
@test PhyloNetworks.shrink3cycles!(net) # hybrid case
@test writeTopology(net, round=true, digits=5) == "(((A:2.0,(B:1.0)#H1:0.13191::0.94):1.37,(C:0.6,#H1:1.0::0.06):0.9):0.6,D:2.0);"
net0 = readTopology("((((((a:1)#H1:1::.9)#H2:1::.8)#H3:1::.7,#H3:0.5):1,#H2:1):1,(#H1:1,b:1):1,c:1);")
net = deepcopy(net0) # new 2/3 cycles appear when some are shrunk
@test PhyloNetworks.shrink2cycles!(net)
@test writeTopology(net, round=true) == "((#H1:1.0::0.1,b:1.0):1.0,c:1.0,(a:1.0)#H1:4.48::0.9);"
@test PhyloNetworks.shrink3cycles!(net0)
writeTopology(net0, round=true) == "(c:1.1,a:5.132,b:1.9);"
# non-tree-child network: w shape
net0 = readTopology("((a:1,#H1:.1::.1):1,(((b:.5)#H3:1)#H1:1,(#H3:0.8::.4)#H2:1):1,(#H2:.2::.0,c:1):1);")
PhyloNetworks.deletehybridedge!(net0, net0.edge[2], false,true,false,false)
net = deepcopy(net0) # delete edge 10 (with γ=0) then shrink 2-cycle
PhyloNetworks.deletehybridedge!(net, net.edge[7], false,true,false,false)
@test PhyloNetworks.shrink3cycles!(net)
@test writeTopology(net) == "(a:2.0,c:2.0,b:3.42);"
# shrink 3-cycle, which deletes edge 10 (with γ=0) : same result in the end
@test PhyloNetworks.shrink3cycles!(net0)
@test writeTopology(net) == "(a:2.0,c:2.0,b:3.42);"
end

end # of set of auxiliary test sets


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using JLD2,MAT
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_1.jld2") results_1
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3.jld2") results_3

incidence_epsilon_0_no_control = results_1[1][1][1:20,:,1]
incidence_epsilon_0_25_no_control = results_3[1][1][1:20,:,1]
incidence_epsilon_0_25_CI = results_3[6][1][1:20,:,1]

matwrite("data/selected_incidence.mat",Dict("inc_0_nc"=>incidence_epsilon_0_no_control,
                                            "inc_025_nc"=>incidence_epsilon_0_25_no_control,
                                            "inc_025_CI"=>incidence_epsilon_0_25_CI))


text:
symbolmatch(s::Symbol, st::AS) where {AS <: AbstractString} = ismatch(Regex("$(st)"), string(s))
symbolmatch(s::Symbol, st::Char) = symbolmatch(s, string(st))

landmarksymbol(i::Int) = Symbol("l$i")
landmarknumber(sy::Symbol) = parse(Int, string(sy)[2:end])

function addoccurance!(lmo::Dict{Symbol, Int}, lm::Symbol)
  lmo[lm] = haskey(lmo, lm) ? lmo[lm]+1 : 1
  nothing
end

"""
    $(SIGNATURES)

For multiple measurements::Dict{Int,T} with unique Int signature, add "LANDMARK" `;softtype=Point2` to the fg<:AbstractDFG.  Type T not used in this function.
"""
function findAddNewLandmarks!(fgl::G, lmoccur::Dict{Symbol, Int}, measurements::Dict{Int, T}; min_occurance::Int=2, softtype=Point2) where {G <: AbstractDFG, T}
  # Count occurances of each landmark
  lmsyms = landmarksymbol.(collect(keys(measurements)))
  addoccurance!.(lmoccur, lmsyms)

  # add new landmarks not already in factor graph
  lmsinfgl = ls(fgl)[2]
  for lm in lmsyms
    if lmoccur[lm] >= min_occurance
      if !(lm in lmsinfgl)
        # add a new landmark
        println("adding lm=$(lm)")
        addVariable!(fgl, lm, softtype, tags=[:LANDMARK;])
      end
      # @show measurements[landmarknumber(lm)]
    end
  end
  nothing
end

function findAddLandmFactorsByPose!(fgl::G, psym::Symbol, meas::Dict{Int, T};
                    range_sigma::Float64=1.0, bearing_sigma::Float64=0.05, N::Int=100) where {G <: AbstractDFG, T}
  lmsyms = landmarksymbol.(collect(keys(meas)))
  posels = ls(fgl, psym)
  lmsinfg = ls(fgl)[2]
  lms = intersect(lmsyms, lmsinfg)
  for lm in lms
    btwfcts = intersect(ls(fgl, lm), posels)
    if length(btwfcts) == 0
      z_rb = meas[landmarknumber(lm)]
      println("addFactor! $(psym) -- $(lm), z_rb=$(z_rb)")
      pp = Pose2Point2BearingRange( Normal(z_rb[2],bearing_sigma),Normal(z_rb[1],range_sigma) )
      addFactor!(fgl, [psym; lm], pp)
    end
  end
  nothing
end

function drive!(fgl::G, idx::Int, odos::Dict, meas::Dict, lmoccurl::Dict, Podo;
                N::Int=100, min_occurance::Int=2 ) where G <: AbstractDFG
  #
  prev, X, nextn = getLastPose2D(fgl)
  vp, fp = addOdoFG!(fgl, nextn, odos[idx][1:3], Podo, N=N)

  findAddNewLandmarks!(fgl, lmoccurl, meas[idx], min_occurance=min_occurance)

  # walk backwards
  donels = Symbol[]
  for i in 0:(min_occurance-1)
    @show psym = Symbol("x$(idx-i-1)") # off by one
    if idx-i > 0
      currmeas = meas[idx-i]
      findAddLandmFactorsByPose!(fgl, psym, currmeas, N=N)
    end
  end
  nothing
end

function donextframe!(fgl::G, idx::Int, odos::Dict, meas::Dict, lmoccurl::Dict, Podo; N::Int=100) where G <: AbstractDFG
  drive!(fgl, idx, odos, meas, lmoccurl, Podo, N=N)

  tree = wipeBuildNewTree!(fgl, drawpdf=true)
  inferOverTree!(fgl, tree,N=N)

  drawPosesLandms(fgl)
end


text:
import MySQL, JSON
import Tables: rowtable, columntable
import DelimitedFiles: writedlm
import Random: randstring, seed!
seed!(0)

host = "127.0.0.1"
user = "root"
password = "test"
conn = MySQL.connect(host, user, password, port=3306)

MySQL.execute!(conn, "truncate table testdb.testtb")
MySQL.execute!(conn, "create database if not exists testdb")
create_table = """
create table if not exists testdb.testtb (
    id int,
    name varchar(255),
    rand decimal(10,8)
)
"""
MySQL.execute!(conn, create_table)

for id in 1:1000
    sql = """insert into testdb.testtb values
            ($(id), \"$(randstring())\", $(rand()));"""
    MySQL.execute!(conn, sql)
end

data = MySQL.Query(conn, "select * from testdb.testtb;")
json_data = JSON.json(rowtable(data))
println("Writing a json file")
write("test.json", json_data)

data = MySQL.Query(conn, "select * from testdb.testtb;")
columntable(data)

println("Writing a CSV file")
data = MySQL.Query(conn, "select * from testdb.testtb;")
writedlm("test.csv", data, ',')
MySQL.disconnect(conn)


text:
# These 3 methods are defined on concrete formulations of the branches to avoid ambiguity
construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranch},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranch},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchBounds},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchBounds},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::Union{NetworkModel{CopperPlatePowerModel}, NetworkModel{AreaBalancePowerModel}},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.DCBranch, HVDCDispatch},
    ::NetworkModel{CopperPlatePowerModel},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.DCBranch, HVDCDispatch},
    ::NetworkModel{CopperPlatePowerModel},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.DCBranch, HVDCLossless},
    ::NetworkModel{CopperPlatePowerModel},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.DCBranch, HVDCLossless},
    ::NetworkModel{CopperPlatePowerModel},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::NetworkModel{<:PM.AbstractPowerModel},
) = nothing

construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{<:PSY.ACBranch, StaticBranchUnbounded},
    ::NetworkModel{<:PM.AbstractPowerModel},
) = nothing

# For DC Power only. Implements constraints
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranch},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractActivePowerModel} end

# For DC Power only. Implements constraints
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranch},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractActivePowerModel}
    @debug "construct_device" _group = LOG_GROUP_BRANCH_CONSTRUCTIONS

    devices = get_available_components(B, sys)
    add_constraints!(container, RateLimitConstraint, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

# For DC Power only
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranch},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)

    add_variables!(container, S, devices, StaticBranch())
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranch},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)

    add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)

    add_constraints!(container, RateLimitConstraint, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranchBounds},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)
    add_variables!(container, S, devices, StaticBranchBounds())
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranchBounds},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)

    add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)

    branch_rate_bounds!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranchUnbounded},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)
    add_variables!(container, S, devices, StaticBranchUnbounded())
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranchUnbounded},
    network_model::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: StandardPTDFModel}
    devices = get_available_components(B, sys)

    add_constraints!(container, NetworkFlowConstraint, devices, model, network_model)
    add_constraint_dual!(container, sys, model)
    return
end

# For AC Power only. Implements Bounds on the active power and rating constraints on the aparent power
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranch},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel} end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranch},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    branch_rate_bounds!(container, devices, model, S)

    add_constraints!(container, RateLimitConstraintFromTo, devices, model, S)
    add_constraints!(container, RateLimitConstraintToFrom, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, StaticBranchBounds},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel} end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, StaticBranchBounds},
    ::NetworkModel{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    branch_rate_bounds!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, F},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, F <: HVDCLossless, S <: PM.AbstractPowerModel} end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, F},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, F <: HVDCUnbounded, S <: PM.AbstractPowerModel} end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, F},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, F <: HVDCDispatch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    add_variables!(container, HVDCTotalPowerDeliveredVariable, devices, F())
    add_to_expression!(
        container,
        ActivePowerBalance,
        HVDCTotalPowerDeliveredVariable,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, F},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, F <: HVDCDispatch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    branch_rate_bounds!(container, devices, model, S)
    add_constraints!(container, FlowRateConstraintFromTo, devices, model, S)
    add_constraints!(container, FlowRateConstraintToFrom, devices, model, S)
    add_constraints!(container, HVDCPowerBalance, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

# Repeated method to avoid ambiguity between HVDCUnbounded and HVDCLossless
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCLossless,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)

    add_variables!(container, FlowActivePowerVariable, devices, U())
    add_to_expression!(
        container,
        ActivePowerBalance,
        FlowActivePowerVariable,
        devices,
        model,
        S,
    )
    return
end

# Repeated method to avoid ambiguity between HVDCUnbounded and HVDCLossless
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCUnbounded,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)

    add_variables!(container, FlowActivePowerVariable, devices, U())
    add_to_expression!(
        container,
        ActivePowerBalance,
        FlowActivePowerVariable,
        devices,
        model,
        S,
    )
    return
end

# Repeated method to avoid ambiguity between HVDCUnbounded and HVDCLossless
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCUnbounded,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)

    add_constraints!(container, FlowRateConstraint, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

# Repeated method to avoid ambiguity between HVDCUnbounded and HVDCLossless
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCLossless,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)
    branch_rate_bounds!(container, devices, model, S)
    add_constraints!(container, FlowRateConstraint, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, U <: HVDCLossless, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    branch_rate_bounds!(container, devices, model, S)
    add_constraints!(container, FlowRateConstraintFromTo, devices, model, S)
    add_constraints!(container, FlowRateConstraintToFrom, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {B <: PSY.DCBranch, U <: HVDCUnbounded, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    add_constraints!(container, FlowRateConstraintFromTo, devices, model, S)
    add_constraints!(container, FlowRateConstraintToFrom, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCDispatch,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)

    add_variables!(container, FlowActivePowerVariable, devices, U())
    # Note: Don't add the flow variables to the expressions for PTDF Models in this constructor
    # otherwise the flow calculation will include the flow variables as part of the nodal balances
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{B, U},
    ::NetworkModel{S},
) where {
    B <: PSY.DCBranch,
    U <: HVDCDispatch,
    S <: Union{StandardPTDFModel, PTDFPowerModel},
}
    devices = get_available_components(B, sys)

    add_constraints!(container, FlowRateConstraint, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end


text:
module WoodburyMatrices

using LinearAlgebra
import LinearAlgebra: det, ldiv!, mul!, adjoint
import Base: *, \, convert, copy, show, similar, size

export Woodbury, SymWoodbury, liftFactor

#### Woodbury matrices ####
"""
`W = Woodbury(A, U, C, V)` creates a matrix `W` identical to `A + U*C*V` whose inverse will be calculated using
the Woodbury matrix identity.
"""
struct Woodbury{T,AType,UType,VType,CType,CpType} <: Factorization{T}
    A::AType
    U::UType
    C::CType
    Cp::CpType
    V::VType
    tmpN1::Vector{T}
    tmpN2::Vector{T}
    tmpk1::Vector{T}
    tmpk2::Vector{T}
end

function Woodbury(A, U::AbstractMatrix{T}, C, V::AbstractMatrix{T}) where {T}
    N = size(A, 1)
    k = size(U, 2)
    if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
        throw(DimensionMismatch("Sizes of U ($(size(U))) and/or V ($(size(V))) are inconsistent with A ($(size(A)))"))
    end
    if k > 1
        if size(C, 1) != k || size(C, 2) != k
            throw(DimensionMismatch("C should be $(k)x$(k)"))
        end
    end
    Cp = inv(convert(Matrix, inv(C) .+ V*(A\U)))
    # temporary space for allocation-free solver
    tmpN1 = Array{T,1}(undef, N)
    tmpN2 = Array{T,1}(undef, N)
    tmpk1 = Array{T,1}(undef, k)
    tmpk2 = Array{T,1}(undef, k)

    # Construct the struct based on the types of the copies,
    # not the originals. See: https://github.com/JuliaLang/julia/issues/26294
    # don't copy A, it could be huge.
    Woodbury(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
end

Woodbury(A, U::Vector{T}, C, V::Matrix{T}) where {T} = Woodbury(A, reshape(U, length(U), 1), C, V)

Woodbury(A, U::AbstractVector, C, V::Adjoint) = Woodbury(A, U, C, Matrix(V))

size(W::Woodbury) = size(W.A)
size(W::Woodbury, d) = size(W.A, d)

function show(io::IO, W::Woodbury)
    println(io, "Woodbury factorization:\nA:")
    show(io, MIME("text/plain"), W.A)
    print(io, "\nU:\n")
    Base.print_matrix(IOContext(io, :compact=>true), W.U)
    if isa(W.C, Matrix)
        print(io, "\nC:\n")
        Base.print_matrix(IOContext(io, :compact=>true), W.C)
    else
        print(io, "\nC: ", W.C)
    end
    print(io, "\nV:\n")
    Base.print_matrix(IOContext(io, :compact=>true), W.V)
end

Base.Matrix(W::Woodbury{T}) where {T} = convert(Matrix{T}, W)
convert(::Type{Matrix{T}}, W::Woodbury{T}) where {T} = Matrix(W.A) + W.U*W.C*W.V

copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)

## Woodbury matrix routines ##

*(W::Woodbury, B::StridedVecOrMat)=W.A*B + W.U*(W.C*(W.V*B))

function \(W::Woodbury, R::StridedVecOrMat)
    AinvR = W.A\R
    AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
end

det(W::Woodbury)=det(W.A)*det(W.C)/det(W.Cp)

function ldiv!(W::Woodbury, B::AbstractVector)
    length(B) == size(W, 1) || throw(DimensionMismatch("Vector length $(length(B)) must match matrix size $(size(W,1))"))
    copyto!(W.tmpN1, B)
    Alu = lu(W.A) # Note. This makes an allocation (unless A::LU). Alternative is to destroy W.A.
    ldiv!(Alu, W.tmpN1)
    mul!(W.tmpk1, W.V, W.tmpN1)
    mul!(W.tmpk2, W.Cp, W.tmpk1)
    mul!(W.tmpN2, W.U, W.tmpk2)
    ldiv!(Alu, W.tmpN2)
    for i = 1:length(W.tmpN2)
        @inbounds B[i] = W.tmpN1[i] - W.tmpN2[i]
    end
    B
end

"""
Given a type `T`, an integer `n` and `m` tuples of the form (i, j, v), build
sparse matrices `rows`, `vals`, `cols` such that the product
`out = rows * vals * cols` is equivalent to:

```julia
out = zeros(T, n, n)

for (i, j, v) in args
    out[i, j] = v
end
```

The first two components (`i` and `j`) of each tuple should be integers
whereas the third component should be of type `T`

Example:

```
julia> r, v, c = WoodburyMatrices.sparse_factors(Float64, 3,
                                                 (1, 1, 2.0),
                                                 (2, 2, 3.0),
                                                 (3, 3, 4.0));

julia> r*c*v - Diagonal([2.0, 3.0, 4.0])
3x3 sparse matrix with 0 Float64 entries:
```

"""
function sparse_factors(::Type{T}, n::Int, args::Tuple{Int, Int, Any}...) where {T}
    m = length(args)
    rows = spzeros(T, n, m)
    cols = spzeros(T, m, n)
    vals = zeros(T, m, m)

    ix = 1
    for (i, (row, col, val)) in enumerate(args)
        rows[row, ix] = 1
        cols[ix, col] = 1
        vals[ix, ix] = val
        ix += 1
    end

    rows, vals, cols
end

include("SymWoodburyMatrices.jl")

end


text:
import Base.in
export in, Perm, fenchel

# the set of permutations of 1,...,n
mutable struct Perm <: ContiSet
  dim::Int
end

# check if w is in the Permutation set
function in(w::AbstractArray, p::Perm)
  n = p.dim
  @assert length(w) == n
  return Set(w) == Set(collect(1:n))
end

# compute min w^x st x in p
function fenchel(p::Perm, w::AbstractArray)
  n = p.dim
  @assert length(w) == n
  i = sortperm(w, rev = true) # biggest index of c is first element of i, etc.
  h = zeros(n)
  for ii=1:n
    h[i[ii]] = ii
  end
  return h
end


text:
using DictFiles
try rm("/tmp/test") catch end
for i = 1:1000
    p = spawn(`julia stresstesthelper.jl`)
    sleep(10)
    println("killing ...")
    kill(p)
    debug = readall(`h5debug /tmp/test`)
    println("\n i == $i\n $debug")
    if ismatch(r"Error", debug)
        break
    end
end


text:
include("../src/PollyBenchmarks.jl")
PollyBenchmarks.run()


text:
tagrng(ω, rng::AbstractRNG) = tag(ω, (rng = rng,))
rng(ω::Ω) = Random.GLOBAL_RNG
rng(ω::TaggedΩ) = haskey(ω.tags, :rng) ? ω.tags.rng : Random.GLOBAL_RNG


text:
# Julia native arrays are ones based
x = [1,2,3];
println( "x[1] = ", x[1] );

# Can use 2n and skip "*" in some cases
a=7;
println( "2a = ", 2a );

# Can use unicode, including for some operators
b=4
println( "√b = ",√b);

# Not every unicode as operator works
# in second println cᵈ is seen as a variable
c=2
d=2
println( "cᵈ = ",c^d);
#E println( "cᵈ = ",cᵈ);

# Dot used for element wise operations
# Otherwise interpreted as whole object function
println( "x.^3 = ",x.^3);
#E println( "x^3 = ",x^3);

#### ERROR: LoadError: MethodError:  .....
#    is a very common message from Julia runtime error

# Assignment copies values for scalars, pointers for arrays
a1=3
a2=a1
a1=4
println("a1, a2 = ",a1,", ", a2);
a1=[3]
a2=a1
a1[1]=4
println("a1, a2 = ",a1,", ", a2);

# Dot is used to "broadcast" arrays in assingment (which copies values)
a1=[3]
a2.=a1
a1[1]=4
println("a1, a2 = ",a1,", ", a2);


text:
abstract type Params end

"""
```
keypoint = Keypoint(y, x)
keypoint = Keypoint(feature)
```

A `Keypoint` may be created by passing the coordinates of the point or from a feature.
"""
const Keypoint = CartesianIndex{2}

"""
```
keypoints = Keypoints(boolean_img)
keypoints = Keypoints(features)
```

Creates a `Vector{Keypoint}` of the `true` values in a boolean image or from a list of features.
"""
const Keypoints = Vector{CartesianIndex{2}}

"""
```
feature = Feature(keypoint, orientation = 0.0, scale = 0.0)
```

The `Feature` type has the keypoint, its orientation and its scale.
"""
struct Feature
    keypoint::Keypoint
    orientation::Float64
    scale::Float64
end

"""
```
features = Features(boolean_img)
features = Features(keypoints)
```

Returns a `Vector{Feature}` of features generated from the `true` values in a boolean image or from a
list of keypoints.
"""
const Features = Vector{Feature}

Feature(k::Keypoint) = Feature(k, 0.0, 0.0)

Feature(k::Keypoint, ori::Number) = Feature(k, ori, 0.0)

Features(keypoints::Keypoints) = map(k -> Feature(k), keypoints)

Features(img::AbstractArray) = Features(Keypoints(img))

Keypoint(feature::Feature) = feature.keypoint

function Keypoints(img::AbstractArray)
    I = findall(!iszero, img)
    r = getindex.(I, 1)
    c = getindex.(I, 2)

    map((ri, ci) -> Keypoint(ri, ci), r, c)
end

Keypoints(features::Features) = map(f -> f.keypoint, features)

const OrientationPair = Tuple{Int16, Int16}
const OrientationWeights = Tuple{Float16, Float16}
const SamplePair = Tuple{Float16, Float16}

"""
```
distance = hamming_distance(desc_1, desc_2)
```

Calculates the hamming distance between two descriptors.
"""
hamming_distance(desc_1, desc_2) = mean(xor.(desc_1, desc_2))

"""
```
matches = match_keypoints(keypoints_1, keypoints_2, desc_1, desc_2, threshold = 0.1)
```

Finds matched keypoints using the [`hamming_distance`](@ref) function having distance value less than `threshold`.
"""
function match_keypoints(keypoints_1::Keypoints, keypoints_2::Keypoints, desc_1, desc_2, threshold::Float64 = 0.1)
    smaller = desc_1
    larger = desc_2
    s_key = keypoints_1
    l_key = keypoints_2
    order = false
    if length(desc_1) > length(desc_2)
        smaller = desc_2
        larger = desc_1
        s_key = keypoints_2
        l_key = keypoints_1
        order = true
    end
    hamming_distances = [hamming_distance(s, l) for s in smaller, l in larger]
    matches = Keypoints[]
    for i in 1:length(smaller)
        if any(hamming_distances[i, :] .< threshold)
            id_min = argmin(hamming_distances[i, :])
            push!(matches, order ? [l_key[id_min], s_key[i]] : [s_key[i], l_key[id_min]])
            hamming_distances[:, id_min] .= 1.0
        end
    end
    matches
end

"""
```
grade = grade_matches(keypoints_1, keypoints_2, limit, difference_method)
```
Returns the fraction of keypoint pairs which have
`difference_method(keypoint_1,keypoint_2)` less than `limit`.
"""
function grade_matches(keypoints_1::Keypoints, keypoints_2::Keypoints, limit::Real, diff::Function = (i,j) -> (sqrt( (i[1]-j[1])^2 + (i[2]-j[2])^2 )))
    @assert length(keypoints_1)==length(keypoints_2) "Keypoint lists are of different lengths."
    @assert length(keypoints_1)!=0 "Keypoint list is of size zero."
    mean(map((keypoint_1,keypoint_2)->((diff(keypoint_1,keypoint_2) < limit) ? 1.0 : 0.0), keypoints_1, keypoints_2))
end


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

"""
    IQ(var₁=>param₁, var₂=>param₂, ...)

Image quilting simulation solver as described in Hoffimann et al. 2017.

## Parameters

### Required

* `trainimg` - Training image from which to extract tiles
* `tilesize` - Tuple with tile size for each dimension

### Optional

* `overlap`  - Overlap size (default to (1/6, 1/6, ..., 1/6))
* `path`     - Simulation path (:raster (default), :dilation, or :random)
* `mapping`  - Data mapping method (default to `NearestMapping()`)
* `inactive` - Vector of inactive voxels (i.e. `CartesianIndex`) in the grid
* `soft`     - A vector of `(data,dataTI)` pairs
* `tol`      - Initial relaxation tolerance in (0,1] (default to 0.1)

## Global parameters

### Optional

* `threads`      - Number of threads in FFT (default to number of physical CPU cores)
* `gpu`          - Whether to use the GPU or the CPU (default to false)
* `showprogress` - Whether to show or not the estimated time duration (default to false)
* `rng`          - Random number generator (default to `Random.GLOBAL_RNG`)

## References

* Hoffimann et al 2017. *Stochastic simulation by image quilting of process-based geological models.*
* Hoffimann et al 2015. *Geostatistical modeling of evolving landscapes by means of image quilting.*
"""
@simsolver IQ begin
  @param trainimg
  @param tilesize
  @param overlap       = nothing
  @param path          = :raster
  @param mapping       = NearestMapping()
  @param inactive      = nothing
  @param soft          = []
  @param tol           = .1
  @global threads      = cpucores()
  @global gpu          = false
  @global showprogress = false
  @global rng          = Random.GLOBAL_RNG
end

function preprocess(problem::SimulationProblem, solver::IQ)
  # retrieve problem info
  pdata   = data(problem)
  pdomain = domain(problem)
  simsize = size(pdomain)
  dims    = embeddim(pdomain)

  # result of preprocessing
  preproc = Dict{Symbol,Tuple}()

  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # training image as simple array
      TI = varparams.trainimg
      trainimg = asarray(TI, var)

      # default overlap
      overlap = varparams.overlap ≠ nothing ? varparams.overlap :
                                              ntuple(i->1/6, dims)

      # determine data mappings
      vmapping = if hasdata(problem)
        map(pdata, pdomain, (var,), varparams.mapping)[var]
      else
        Dict()
      end

      # create hard data object
      hdata = Dict{CartesianIndex{dims},Real}()
      for (loc, datloc) in vmapping
        push!(hdata, lin2cart(simsize, loc) => pdata[var][datloc])
      end

      # disable inactive voxels
      shape = Dict{CartesianIndex{dims},Real}()
      if varparams.inactive ≠ nothing
        for icoords in varparams.inactive
          push!(shape, icoords => NaN)
        end
      end

      hard = merge(hdata, shape)

      preproc[var] = (varparams, trainimg, simsize, overlap, hard)
    end
  end

  preproc
end

function solvesingle(::SimulationProblem, covars::NamedTuple, solver::IQ, preproc)
  # random number generator
  rng = solver.rng

  varreal = map(covars.names) do var
    # unpack preprocessed parameters
    par, trainimg, simsize, overlap, hard = preproc[var]

    # run image quilting core function
    reals = iqsim(trainimg, par.tilesize, simsize;
                  overlap=overlap, path=par.path,
                  soft=par.soft, hard=hard, tol=par.tol,
                  threads=solver.threads, gpu=solver.gpu,
                  showprogress=solver.showprogress, rng=rng)

    # flatten result
    var => vec(reals[1])
  end

  Dict(varreal)
end


text:
using TranscodingStreams: splitkwargs

const BUF_SIZE = 16*1024
const LZ4_FOOTER_SIZE = 4

mutable struct LZ4Compressor <: TranscodingStreams.Codec
    ctx::Ref{Ptr{LZ4F_cctx}}
    prefs::Ref{LZ4F_preferences_t}
    header::Memory
    write_header::Bool
end

"""
    LZ4Compressor(; kwargs...)

Creates an LZ4 compression codec.

# Keywords
- `blocksizeid::BlockSizeID=default_size`: `max64KB`, `max256KB`, `max1MB`, or `max4MB` or `default_size`
- `blockmode::BlockMode=block_linked`:  `block_linked` or `block_independent`
- `contentchecksum::Bool=false`: if `true`, frame is terminated with a
    32-bits checksum of decompressed data
- `frametype::FrameType=normal_frame)`:  `normal_frame` or `skippable_frame`
- `contentsize::Integer=0`: Size of uncompressed content (0 for unknown)
- `blockchecksum::Bool=false`: if `true`, each block is followed by a
    checksum of block's compressed data
- `compressionlevel::Integer=0`: compression level (-1..12)
- `autoflush::Bool=false`: always flush if `true`
"""
function LZ4Compressor(; kwargs...)
    x, y = splitkwargs(kwargs, (:compressionlevel, :autoflush))
    ctx = Ref{Ptr{LZ4F_cctx}}(C_NULL)
    frame = LZ4F_frameInfo_t(; y...)
    prefs = Ref(LZ4F_preferences_t(frame; x...))
    return LZ4Compressor(ctx, prefs, Memory(Vector{UInt8}(undef, LZ4F_HEADER_SIZE_MAX)), false)
end

const LZ4CompressorStream{S} = TranscodingStream{LZ4Compressor,S} where S<:IO

"""
    LZ4CompressorStream(stream::IO; kwargs...)

Creates an LZ4 compression stream. See `LZ4Compressor()` and `TranscodingStream()` for arguments.
"""
function LZ4CompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:blocksizeid, :blockmode, :contentchecksum, :blockchecksum, :frametype, :contentsize, :compressionlevel, :autoflush))
    return TranscodingStream(LZ4Compressor(; x...), stream; y...)
end

"""
Returns the expected size of the transcoded data.
"""
function TranscodingStreams.expectedsize(codec::LZ4Compressor, input::Memory)::Int
    LZ4F_compressBound(input.size, codec.prefs) + LZ4F_HEADER_SIZE_MAX + LZ4_FOOTER_SIZE
end

"""
Returns the minimum output size of `process`.
"""
function TranscodingStreams.minoutsize(codec::LZ4Compressor, input::Memory)::Int
    LZ4F_HEADER_SIZE_MAX
end

"""
Initializes the LZ4F Compression Codec.
"""
function TranscodingStreams.initialize(codec::LZ4Compressor)::Nothing
    LZ4F_createCompressionContext(codec.ctx, LZ4F_getVersion())
    nothing
end

"""
Finalizes the LZ4F Compression Codec.
"""
function TranscodingStreams.finalize(codec::LZ4Compressor)::Nothing
    LZ4F_freeCompressionContext(codec.ctx[])
    nothing
end

"""
Starts processing with the codec
Creates the LZ4F header to be written to the output.
"""
function TranscodingStreams.startproc(codec::LZ4Compressor, mode::Symbol, error::Error)::Symbol
    try
        header = Vector{UInt8}(undef, LZ4F_HEADER_SIZE_MAX)
        headerSize = LZ4F_compressBegin(codec.ctx[], header, convert(Csize_t, LZ4F_HEADER_SIZE_MAX), codec.prefs)
        codec.header = Memory(resize!(header, headerSize))
        codec.write_header = true
        :ok
    catch err
        error[] = err
        :error
    end
end

"""
Compresses the data from `input` and writes to `output`.
The LZ4 compression algorithm may simply buffer the input data a full frame can be produced, so `data_written` may be 0.
`flush()` may be used to force `output` to be written.
"""
function TranscodingStreams.process(codec::LZ4Compressor, input::Memory, output::Memory, error::Error)::Tuple{Int,Int,Symbol}
    data_read = 0
    data_written = 0
    if codec.write_header
        if output.size < codec.header.size
            error[] = ErrorException("Output buffer too small for header.")
            return (data_read, data_written, :error)
        end
        unsafe_copyto!(output.ptr, codec.header.ptr, codec.header.size)
        data_written = codec.header.size
        codec.write_header = false
    end

    try
        if input.size == 0
            data_written += LZ4F_compressEnd(codec.ctx[], output.ptr + data_written, output.size - data_written, C_NULL)
            (data_read, data_written, :end)
        else
            data_written += LZ4F_compressUpdate(codec.ctx[], output.ptr + data_written, output.size - data_written, input.ptr, input.size, C_NULL)
            (input.size, data_written, :ok)
        end

    catch err
        error[] = err
        (data_read, data_written, :error)
    end

end

struct LZ4Decompressor <: TranscodingStreams.Codec
    dctx::Ref{Ptr{LZ4F_dctx}}
end

"""
    LZ4Compressor()

Creates an LZ4 decompression codec.
"""
function LZ4Decompressor()
    dctx = Ref{Ptr{LZ4F_dctx}}(C_NULL)
    return LZ4Decompressor(dctx)
end

const LZ4DecompressorStream{S} = TranscodingStream{LZ4Decompressor,S} where S<:IO

"""
    LZ4CompressorStream(stream::IO; kwargs...)

Creates an LZ4 decompression stream. See `TranscodingStream()` for arguments.
"""
function LZ4DecompressorStream(stream::IO; kwargs...)
    return TranscodingStream(LZ4Decompressor(), stream; kwargs...)
end

"""
Initializes the LZ4F Decompression Codec.
"""
function TranscodingStreams.initialize(codec::LZ4Decompressor)::Nothing
    LZ4F_createDecompressionContext(codec.dctx, LZ4F_getVersion())
    nothing
end

"""
Finalizes the LZ4F Decompression Codec.
"""
function TranscodingStreams.finalize(codec::LZ4Decompressor)::Nothing
    LZ4F_freeDecompressionContext(codec.dctx[])
    nothing
end

"""
Deompresses the data from `input` and writes to `output`.
If the input data is not properly formatted this function will throw an error.
"""
function TranscodingStreams.process(codec::LZ4Decompressor, input::Memory, output::Memory, error::Error)::Tuple{Int,Int,Symbol}
    data_read = 0
    data_written = 0

    try
        if input.size == 0
            (data_read, data_written, :end)
        else
            src_size = Ref{Csize_t}(input.size)
            dst_size = Ref{Csize_t}(output.size)
            LZ4F_decompress(codec.dctx[], output.ptr, dst_size, input.ptr, src_size, C_NULL)
            (src_size[], dst_size[], :ok)
        end

    catch err
        if isa(err, LZ4Exception) && err.msg == "ERROR_frameType_unknown"
            codec.dctx[] = C_NULL
        end
        error[] = err
        (data_read, data_written, :error)
    end

end


text:
function dist(coord_01, coord_02)
    num_dims = length(coord_01)
    dist = 0
    for dim = 1 : num_dims
        dist += (coord_01[dim] - coord_02[dim])^2
    end
    sqrt(dist)
end
function str2float(str::AbstractString)
    str = split(str)
    [parse(Float64, str[n]) for n = 1:length(str)]
end

function diag(vec::AbstractArray)
    len = length(vec)
    result = zeros(len, len)
    for i = 1 : len
        result[i, i] = vec[i]
    end
    result
end

function copy_array(goal)
    res = zeros(typeof(goal).parameters[1], size(goal))
    res[:] = goal[:]
    return res
end

function leastsq(x, y; dims=2)
    mat_left = Array(hcat([
        [sum(x.^i) for i = j:j+dims-1] for j = 0:dims-1
    ]...)')

    mat_right = [
        sum(y.* x.^i) for i = 0:dims-1
    ]
    return inv(mat_left) * mat_right
end


text:
function terminator_loop(output::IO, input::IO, ppid::Int)
    mypid = getpid()
    timers = Dict{UInt,Timer}()
    for ln in eachline(input)
        @debug "[$mypid] GOT: $ln"
        sop, sid, stime, label = split(ln; limit = 4, keepempty = true)
        id = parse(UInt, sid)
        time = parse(Float64, stime)
        print(output, "ack $sid\n")
        flush(output)
        if sop == "start"
            timer = Timer(time) do _
                msglabel = label == "" ? "" : " for $label;"
                @error "[$mypid] Timeout ($time)$msglabel terminating process $ppid"
                isactive() = isdir("/proc/$ppid")
                function signalling(sig, n = typemax(Int))
                    @info "[$mypid] Trying to terminate process $ppid with $sig"
                    for _ in 1:n
                        run(`kill -s $sig $ppid`)
                        sleep(0.1)
                        isactive() || return true
                    end
                    return false
                end
                signalling("SIGINT", 10) && return
                signalling("SIGTERM", 10) && return
                signalling("SIGHUP", 100) && return
                signalling("SIGKILL")
            end
            timers[id] = timer
        elseif sop == "stop"
            close(pop!(timers, id))
        else
            error("unknown op: $sop")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    terminator_loop(stdout, stdin, parse(Int, ARGS[1]))
end


text:
module Meta
export quot, is_expr, is_fdef, split_fdef, subs_ex, @expect, @get!

macro expect(pred)
    quote
        $(esc(pred)) ? nothing : error(
          $(string("expected: ", sprint(Base.show_unquoted, pred)", == true")))
    end
end

# ---- Metaprogramming --------------------------------------------------------

quot(ex) = Expr(:quote, ex)

is_expr(ex::Expr, head)          = ex.head === head
is_expr(ex::Expr, heads::Set)    = ex.head in heads
is_expr(ex::Expr, heads::Vector) = ex.head in heads
is_expr(ex,       head)          = false
is_expr(ex,       head, n::Int)  = is_expr(ex, head) && length(ex.args) == n

function subs_ex(subs::Function, ex)
    s = subs(ex)
    s === nothing ? subsubs_ex(subs, ex) : s
end

subsubs_ex(subs::Function, ex) = ex
function subsubs_ex(subs::Function, ex::Expr)
    ex.head === :quote ? ex : Expr(ex.head, Any[subs_ex(subs,a) for a in ex.args]...)
end

function is_fdef(ex::Expr)
    is_expr(ex,:function,2) || (is_expr(ex,:(=),2)&&is_expr(ex.args[1],:call))
end
is_fdef(ex) = false

# Return the signature and body from a named method definition, either syntax.
# E.g. split_fdef( :( f(x) = x^2) ) == (:(f(x)), :(x^2))
function split_fdef(fdef::Expr)
    @expect (fdef.head == :function) || (fdef.head == :(=))
    @expect length(fdef.args) == 2
    signature, body = fdef.args
    @expect is_expr(signature, :call)
    @expect length(signature.args) >= 1
    (signature, body)
end
split_fdef(f::Any) = error("split_fdef: expected function definition, got\n$f")

# ---- Other stuff (todo: move?) ---------------------------------------------

macro get!(d, k, default)
    quote
        d, k = $(esc(d)), $(esc(k))
        haskey(d, k) ? d[k] : (d[k] = $(esc(default)))
    end
end

end # module


text:
# NOTE: there's a `kmeans!` function that updates centers, maybe a candidate
# for the `update` machinery. Same for `kmedoids!`
# NOTE: if the prediction is done on the original array, just the assignment
# should be returned, unclear what's the best way of doing this.

module MLJClusteringInterface

# ===================================================================
# IMPORTS
import Clustering
import MLJModelInterface
import MLJModelInterface: Continuous, Count, Finite, Multiclass, Table, OrderedFactor,
    @mlj_model, metadata_model, metadata_pkg

using Distances

# ===================================================================
## EXPORTS
export KMeans, KMedoids

# ===================================================================
## CONSTANTS
# Define constants for easy referencing of packages
const MMI = MLJModelInterface
const Cl = Clustering

# Definitions of model descriptions for use in model doc-strings.
const KMeansDescription ="""
K-Means algorithm: find K centroids corresponding to K clusters in the data.
"""

const KMedoidsDescription ="""
K-Medoids algorithm: find K centroids corresponding to K clusters in the data.
Unlike K-Means, the centroids are found among data points themselves.
"""

const KMFields ="""
    ## Keywords

    * `k=3`     : number of centroids
    * `metric`  : distance metric to use
"""

const PKG = "MLJClusteringInterface"

####
#### KMeans
####
"""
KMeans(; kwargs...)

$KMeansDescription

$KMFields

See also the
[package documentation](http://juliastats.github.io/Clustering.jl/latest/kmeans.html).
"""

@mlj_model mutable struct KMeans <: MMI.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
end

####
#### KMeans
####

function MMI.fit(model::KMeans, verbosity::Int, X)
    # NOTE: using transpose here to get a LinearAlgebra.Transpose object
    # which Kmeans can handle.
    Xarray = transpose(MMI.matrix(X))
    result = Cl.kmeans(Xarray, model.k; distance=model.metric)
    cluster_labels = MMI.categorical(1:model.k)
    fitresult = (result.centers, cluster_labels) # centers (p x k)
    cache = nothing
    report = (
        assignments=result.assignments, # size n
        cluster_labels=cluster_labels
    )
    return fitresult, cache, report
end

MMI.fitted_params(::KMeans, fitresult) = (centers=fitresult[1],)

function MMI.transform(model::KMeans, fitresult, X)
    # pairwise distance from samples to centers
    X̃ = pairwise(
        model.metric,
        transpose(MMI.matrix(X)),
        fitresult[1],
        dims=2
    )
    return MMI.table(X̃, prototype=X)
end

"""
KMedoids(; kwargs...)

$KMedoidsDescription

$KMFields

See also the
[package documentation](http://juliastats.github.io/Clustering.jl/latest/kmedoids.html).
"""
@mlj_model mutable struct KMedoids <: MMI.Unsupervised
    k::Int = 3::(_ ≥ 2)
    metric::SemiMetric = SqEuclidean()
end

function MMI.fit(model::KMedoids, verbosity::Int, X)
    # NOTE: using transpose=true will materialize the transpose (~ permutedims), KMedoids
    # does not yet accept LinearAlgebra.Transpose
    Xarray = MMI.matrix(X, transpose=true)
    # cost matrix: all the pairwise distances
    cost_array = pairwise(model.metric, Xarray, dims=2) # n x n
    result = Cl.kmedoids(cost_array, model.k)
    cluster_labels = MMI.categorical(1:model.k)
    fitresult = (view(Xarray, :, result.medoids), cluster_labels) # medoids
    cache = nothing
    report = (
        assignments=result.assignments, # size n
        cluster_labels=cluster_labels
    )
    return fitresult, cache, report
end

MMI.fitted_params(::KMedoids, fitresult) = (medoids=fitresult[1],)

function MMI.transform(model::KMedoids, fitresult, X)
    # pairwise distance from samples to medoids
    X̃ = pairwise(
        model.metric,
        MMI.matrix(X, transpose=true),
        fitresult[1], dims=2
    )
    return MMI.table(X̃, prototype=X)
end

####
#### Predict methods
####

function MMI.predict(model::Union{KMeans,KMedoids}, fitresult, Xnew)
    locations, cluster_labels = fitresult
    Xarray = MMI.matrix(Xnew)
    (n, p), k = size(Xarray), model.k
    pred = zeros(Int, n)

    @inbounds for i in 1:n
        minv = Inf
        @inbounds @simd for j in 1:k
            curv = evaluate(
                model.metric, view(Xarray, i, :), view(locations, :, j)
            )
            P = curv < minv
            pred[i] = j * P + pred[i] * !P # if P is true --> j
            minv = curv * P + minv * !P # if P is true --> curvalue
        end
    end
    return cluster_labels[pred]
end

####
#### METADATA
####

metadata_pkg.(
    (KMeans, KMedoids),
    name="Clustering",
    uuid="aaaa29a8-35af-508c-8bc3-b662a17a0fe5",
    url="https://github.com/JuliaStats/Clustering.jl",
    julia=true,
    license="MIT",
    is_wrapper=false
)

metadata_model(
    KMeans,
    input = MMI.Table(Continuous),
    output = MMI.Table(Continuous),
    weights = false,
    descr = KMeansDescription,
    path = "$(PKG).KMeans"
)

metadata_model(
    KMedoids,
    input = MMI.Table(Continuous),
    output = MMI.Table(Continuous),
    weights = false,
    descr = KMedoidsDescription,
    path = "$(PKG).KMedoids"
)

end # module