text:
mutex = Threads.ReentrantLock()
#iomutex = Threads.ReentrantLock()

function setindex_taskthreads_encode_worker( channel::Channel{Tuple}, buf::Array{T,N},
                                        ba::BigArray{D,T}, iter::Tuple) where {D,T,N}
    println("encode worker at ", Threads.threadid())
    C = get_encoding(ba)
    (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) = iter
    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba,
                                    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
    chk = buf[rangeInBuffer]
    mipLevelName = get_mip_level_name(ba)
    key = joinpath(mipLevelName, cartesian_range2string( chunkGlobalRange ) )
    data = encode( chk, C)

    println("encode worker at ", Threads.threadid(), " start ingesting data to channel")
    Threads.lock(mutex)
    put!(channel, (data, iter))
    Threads.unlock(mutex)
    println("encode worker at ", Threads.threadid(), " finished work.")
end

function upload_worker(channel::Channel{Tuple}, ba::BigArray)
    println("upload worker at ", Threads.threadid())
    #Threads.lock(mutex)
    (data, iter) = take!(channel)
    #Threads.unlock(mutex)
    println("upload worker at ", Threads.threadid(), " have get the data from channel.")
    (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) = iter
    ba.kvStore[ cartesian_range2string(chunkGlobalRange) ] = data
    println("upload worker at ", Threads.threadid(), " finished work.")
end

"""
    put array in RAM to a BigArray backend
this version uses channel to control the number of asynchronized request
"""
function setindex_taskthreads!( ba::BigArray{D,T}, buf::Array{T,N},
                       idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    idxes = colon2unit_range(buf, idxes)
    offset = get_offset(ba)
    chunkSize = BigArrays.get_chunk_size(ba)
    # check alignment
    @assert all(map((x,y,z)->mod(first(x) - 1 - y, z),
                    idxes, offset.I, chunkSize).==0)
                    "the start of index should align with BigArray chunk size"
    t1 = time()
    baIter = ChunkIterator(idxes, chunkSize; offset=offset)
    channel = Channel{Tuple}( CHUNK_CHANNEL_SIZE )
    @sync begin
        @async begin
            for iter in baIter
                Threads.@spawn setindex_taskthreads_encode_worker(channel, buf, ba, iter)
            end
        end
        for iter in baIter
            @async upload_worker(channel, ba)
        end
        close(channel)
    end
    elapsed = time() - t1 # sec
    println("saving speed: $(sizeof(buf)/1024/1024/elapsed) MB/s")
end

function getindex_taskthreads_download_worker(channel::Channel{Tuple},
                                        ba::BigArray, iter::Tuple)
    baRange = CartesianIndices(ba)
    blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = iter
    chunkSize = (last(chunkGlobalRange) - first(chunkGlobalRange) + one(CartesianIndex{N})).I
    if any(map((x,y)->x>y, first(globalRange).I, last(baRange).I)) ||
                        any(map((x,y)->x<y, last(globalRange).I, first(baRange).I))
        @warn("out of volume range, keep it as zeros")
        return
    end
    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba,
                                    chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
    data = ba.kvStore[ cartesian_range2string(chunkGlobalRange) ]

    Threads.lock(mutex)
    put!(channel, (data, iter))
    Threads.unlock(mutex)
end

function getindex_taskthreads_decode_worker!(dataChannel::Channel{Tuple}, ret::OffsetArray{T,N}) where {T,N}
    Threads.lock(mutex)
    data, iter = take!(dataChannel)
    Threads.unlock(mutex)

    blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = iter
    chk = Codings.decode(data, get_encoding(ba))
    chk = reshape(reinterpret(T, chk), chunkSize)
    chk = chk[rangeInChunk]
    block = OffsetArray(chk, cartesian_range2unit_range(globalRange)...)
    ret[axes(block)...] = parent(block)
end

function getindex_taskthreads( ba::BigArray, idxes::Union{UnitRange, Int}...)
    t1 = time()
    sz = map(length, idxes)
    ret = OffsetArray(zeros(T, sz), idxes...)

    channelSize = cld( nworkers(), 2 )
    channel = Channel{OffsetArray}( channelSize );

    baIter = ChunkIterator(idxes, ba.chunkSize; offset=get_offset(ba))

    @sync begin
        @async begin
            for iter in baIter
                @async getindex_taskthreads_download_worker(channel, ba)
            end
            close(channel)
        end
        # control the number of concurrent requests here
        for iter in baIter
            Threads.@spawn getindex_taskthreads_decode_worker(ret, channel)
        end

        @async begin
            for iter in baIter
                block = take!(results)
                ret[axes(block)...] = parent(block)
            end
        end
    end
    elapsed = time() - t1 # seconds
    println("cutout speed: $(sizeof(ret)/1024/1024/elapsed) MB/s")
    # handle single element indexing, return the single value
    ret
end


text:
import Base: convert

using FluxionArrays

# Abstract types

# Implements:
# `convert(FieldExpr, x)`
abstract type FieldExpr end

abstract type Dimension <: FieldExpr end
abstract type Field <: FieldExpr end
abstract type KnownDimension <: Dimension end

# Dimensions

struct UnknownDimension <: Dimension
    id :: Symbol

    UnknownDimension() = new(gensym("UD"))
end

struct BasicDimension <: KnownDimension
    id :: Symbol
    grid :: Grid

    BasicDimension(grid::Grid) = new(gensym("BD"), grid)
end

struct TransformedDimension <: KnownDimension
    id :: Symbol
    grid :: Grid
    base_dimension :: BasicDimension

    TransformedDimension(grid::Grid, base_dimension::BasicDimension) =
        new(gensym("TD"), grid, base_dimension)
end

# Stochastic metadata

# This object represents three possible states:
# a non-stochastic field expression, an unknown stochastic field expression,
# and a known stochastic field expression (which implies a known number of trajectories)
struct StochasticInfo
    trajectories :: Nullable{Int64}

    StochasticInfo() = new(nothing)
    StochasticInfo(trajectories::Int64) = new(trajectories)
end

struct StochasticMean <: FieldExpr
    fexpr :: FieldExpr
end

struct Sample <: FieldExpr
    fexpr :: FieldExpr
    trajectories :: Int64
end

struct StochasticDimension <: KnownDimension
end

const STOCHASTIC_DIMENSION = StochasticDimension()

# Fields

struct UnknownField <: Field
    id :: Symbol
    dimensions :: Array{Dimension, 1}
    stochastic :: Bool

    UnknownField(dimensions, stochastic=false) = new(gensym("UF"), dimensions, stochastic)
end

struct KnownField <: Field
    dimensions :: Array{KnownDimension, 1}
    data :: BackendArray
    stochastic_info :: StochasticInfo

    KnownField(dimensions, data) = new(dimensions, data, StochasticInfo())
    KnownField(dimensions, data, stochastic_info) = new(dimensions, data, stochastic_info)
end

struct FieldExprNode <: FieldExpr
    head :: Symbol
    args :: Array{FieldExpr, 1}
end

struct Fourier <: FieldExpr
    fexpr :: FieldExpr
    inverse :: Bool
    dimensions :: Dict{Dimension, Dimension}
end

struct Diff <: FieldExpr
    field :: UnknownField
    dimensions :: Dict{Dimension, Int}
end

struct RandomVariable <: FieldExpr
    distribution :: Symbol
    dimensions :: Array{Dimension, 1}
end

# Methods

grid(dim::BasicDimension) = dim.grid
grid(dim::TransformedDimension) = dim.grid

convert(::Type{FieldExpr}, val::Number) = KnownField([], to_backend(reshape([val], ())))
convert(::Type{FieldExpr}, val::FieldExpr) = val

to_subnodes(fexpr::FieldExpr) = nothing
to_subnodes(fexpr::FieldExprNode) = fexpr.args
to_subnodes(fexpr::Fourier) = [fexpr.fexpr]
to_subnodes(fexpr::Diff) = [fexpr.field]
to_subnodes(fexpr::StochasticMean) = [fexpr.fexpr]

from_subnodes(fexpr::FieldExprNode, subnodes) = FieldExprNode(fexpr.head, subnodes)
from_subnodes(fexpr::Fourier, subnodes) = Fourier(subnodes[1], fexpr.inverse, fexpr.dimensions)
from_subnodes(fexpr::Diff, subnodes) = Diff(subnodes[1], fexpr.dimensions)
from_subnodes(fexpr::StochasticMean, subnodes) = StochasticMean(subnodes[1])

known_dimensions(x) = []
known_dimensions(x::UnknownField) = [dim for dim in x.dimensions if isa(dim, KnownDimension)]
known_dimensions(x::KnownField) = x.dimensions
known_dimensions(x::KnownDimension) = [x]
known_dimensions(x::TransformedDimension) = [x]

unknown_dimensions(x::KnownField) = []
unknown_dimensions(x::UnknownField) = [dim for dim in x.dimensions if isa(dim, UnknownDimension)]
unknown_dimensions(x::KnownDimension) = []
unknown_dimensions(x::TransformedDimension) = []
unknown_dimensions(x::UnknownDimension) = [x]

function known_dimensions_expr(fexpr)
    inspect(PersistentSet{Dimension}(), fexpr) do controls, dims, x
        if typeof(fexpr) == Fourier
            skip_children!(controls)
            inner_dims = known_dimensions_expr(fexpr.fexpr)
            union(dims, [get(fexpr.dimensions, dim, dim) for dim in inner_dims])
        else
            union(dims, known_dimensions(x))
        end
    end
end

function unknown_dimensions_expr(fexpr)
    inspect(PersistentSet{Dimension}(), fexpr) do controls, dims, x
        union(dims, unknown_dimensions(x))
    end
end

function remove_unknown_dimensions(f::UnknownField)
    remove_dimensions(f, unknown_dimensions(f)...)
end

function remove_dimensions(f::UnknownField, dims...)
   to_remove = Set{Dimension}(dims)
   UnknownField([dim for dim in f.dimensions if !in(dim, to_remove)])
end

function replace_expr(fexpr::FieldExpr, replacements...)
    modify(fexpr, Dict(replacements...)) do controls, x, replacement_dict
        get(replacement_dict, x, x)
    end
end

ext: jl


text:
module ObaraSaika

function get_ijk_list(m)
    l = Array{Int, 1}[]
    for a = 1 : m + 1
        for b = 1 : a
            i = m + 1 - a
            j = a - b
            k = b - 1
            push!(l, [i, j, k])
        end
    end
    return l
end

function get_shell4(a, b, c, d)
    components = Array{Int, 1}[]
    for p in get_ijk_list(a)
        for q in get_ijk_list(b)
            for r in get_ijk_list(c)
                for s in get_ijk_list(d)
                    push!(components, vcat(p, q, r, s))
                end
            end
        end
    end
    return components
end

function get_shell2(a, b)
    components = Array{Int, 1}[]
    for p in get_ijk_list(a)
        for q in get_ijk_list(b)
            push!(components, vcat(p, q))
        end
    end
    return components
end

function list_is_flat(l)
    for x in l
        if x <: Array
            return false
        end
    end
    return true
end

function flatten_sub(l)
    return [item for sublist in l for item in sublist]
end

function flatten(l)
    l_out = deepcopy(l)
    max_tries = 50
    for _ = 1 : max_tries
        if list_is_flat(l_out)
            return l_out
        else
            l_out = flatten_sub(l_out)
        end
    end
    # thrown an error?
    # sys.stderr.write('ERROR: max depth reached in flatten\n')
    # sys.exit(1)
end

function find_fun_to_lower(q, n)

    # Determine the total angular momentum on each center.
    l = Int[]
    for i = 1 : n
        push!(l, q[i*3] + q[i*3 + 1] + q[i*3 + 2])
    end

    # find function to lower
    # start with lowest angular momentum above s
    fun = -1
    kmax = max(l) + 1
    for i = 1 : n
        k = l[i]
        # If we're larger than an s-function...
        if k > 0
            if k < kmax
                kmax = k
                fun = i
            end
        end
    end

    return fun

end

function find_component_to_lower(fun)
    for (i, c) in enumerate(fun)
        if c > 0
            return i
        end
    end
    return -1
end

function apply_os4(x4)
    # TODO
end

function apply_os2(x, kind)
    # TODO
end

function get_r12_squared(r1, r2)
    # have an assertion that they're both length 3?
    # do some array trick rather than this crap?
    return (r1[0] - r2[0])^2.0 + (r1[1] - r2[1])^2.0 + (r1[2] - r2[2])^2.0
end

function get_k(z1, z2, r1, r2)
    r12 = get_r12_squared(r1, r2)
    f0 = z1 + z2
    if r12 > 0.0
        f1 = -z1*z2*r12/f0
        f2 = exp(f1)
    else
        f2 = 1.0
    end
    return sqrt(2.0)*f2*pi^(5.0/4.0)/f0
end

function get_rho(za, zb, zc, zd)
    z = za + zb
    n = zc + zd
    return z * n / (z + n)
end

function get_bi_center(z1, z2, r1, r2)
    z = z1 + z2
    rx = (z1*r1[0] + z2*r2[0])/z
    ry = (z1*r1[1] + z2*r2[1])/z
    rz = (z1*r1[2] + z2*r2[2])/z
    return [rx, ry, rz]
    # return (z1*r1 + z2*r2) / z
end

using GSL

function boys(n, x)
    if x > 0.0
        f = 2.0*x^(n + 0.5)
        g = gamma(n + 0.5)
        gi = 1.0 - gammainc(n + 0.5, x, regularized=True)
        return g*gi/f
    else
        return 1.0/(n*2 + 1)
    end
end

end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Generate documentation

push!(LOAD_PATH, "../src/")

using Pkg
Pkg.activate(normpath(@__DIR__))
# Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using Documenter

using Raytracer
using Raytracer:
    SimpleShape, CompositeShape

makedocs(
    sitename = "Raytracer.jl",
    pages = [
        "Introduction" => "index.md",
        "Quickstart" => [
            "Basic CLI usage" => "quickstart/cli.md",
            "Basic SceneLang usage" => "quickstart/scenelang.md",
            "Basic API usage" => "quickstart/api.md"
        ],
        "SceneLang" => "scenelang.md",
        "CLI tool" => "cli.md",
        "API" => [
            "High-level API" => "api/high-level.md",
            "Low-level API" => "api/low-level.md"
        ],
        "Extendability" => "extendability.md",
        "For devs" => [
            "Collaboration instructions" => "devs/collab.md",
            "Private API" => "devs/private-api.md",
            "SceneLang interpreter API" => "devs/scenelang-api.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/Paolo97Gll/Raytracer.jl.git",
    push_preview = true
)


text:
#using LIBLINEAR, UnicodePlots
#include("utils.jl")
#include("main.jl")

# USAGE: copy and paste the following line to work with this script
# include("utils.jl"); include("main.jl"); using Revise, LIBLINEAR, UnicodePlots; includet("inspect-models.jl")
#inspect("_entropy_pan20-en.params", "PAN20/en.train.json", "PAN20/en.test.json")

function vocab(config, dataname)
    D = loadjson(dataname, textkey="text", labelkey="klass")
    textconfig = config.textconfig
    #textconfig = copy(config.textconfig, nlist=[3], qlist=[])
    config = MicroTC_Config(textconfig, config.textmodel, config.cls)
    c = MicroTC(config, D.corpus, D.labels, tok=Tokenizer(config.textconfig))
    Dict(decode(c.tok, k) => v.weight for (k, v) in c.textmodel.tokens)
end

function topk(vec, k=length(vec))
    s = sort!(collect(vec), by=p->p[end], rev=true)
    if k !== length(vec)
        s[1:k]
    else
        s
    end
end

function prune(vec, k::Integer)
    Dict(topk(vec, k))
end

function prune(vec, minw::Real)
    Dict(k => w for (k, w) in vec if w >= minw)
end

function merge(A, B)
    if length(B) < length(A)
        A, B = B, A
    end

    Dict(k => sqrt(v * B[k]) for (k, v) in A if haskey(B, k))
end

function plotweights(vec, k=0)
    L = sort!(collect(vec), by=p->p[end], rev=true)
    if k > 0
        L = L[1:k]
    end
    lineplot(last.(L))
end

function inspect(paramsfile, trainfile, testfile)
    P = loadparams(paramsfile)
    vtrain = vocab(P[1][1], trainfile)
    vtest = vocab(P[1][1], testfile)
    vtrain, vtest
end

max_stars_repo_name: sadit/text-classification-scripts

max_issues_repo_name: sadit/text-classification-scripts

max_forks_repo_name: sadit/text-classification-scripts


text:
### ADAPTIVE GAUSS-LOBATTO FROM GANDER AND GAUTSCHI (2000)

### constants from section 4 of the paper
const α = sqrt(2/3)
const β = 1/sqrt(5)

const X1 = 0.94288241569547971905635175843185720232
const X2 = 0.64185334234578130578123554132903188354
const X3 = 0.23638319966214988028222377349205292599

const A = 0.015827191973480183087169986733305510591
const B = 0.094273840218850045531282505077108171960
const C = 0.155071987336585396253963597980210298680
const D = 0.18882157396018245442000533937297167125
const E = 0.199773405226858526792068022066048840246
const F = 0.22492646533333952701601768799639508076
const G = 0.24261107190140773379964095790325635233

"""
    adaptlob(f,a,b; tol::Real=1e-10, trace::Bool=false)

Numerically integrate `f` over the closed real interval `[a,b]` to within tolerance `tol` using the adaptive Gauss-Lobatto method from Gander and Gautschi (2000).

Passing `trace=true` prints out information about where the function is being sampled.

`f` must be defined at `a` and `b`.

```
using QuadGG
adaptlob(sqrt,0,1)
```
"""
function adaptlob(f,a,b; tol::Real=1e-10, trace::Bool=false)
    m = (a+b)/2
    h = (b-a)/2

    x = (a,m-X1*h,m-α*h,m-X2*h,m-β*h,m-X3*h,m,m+X3*h,m+β*h,m+X2*h,m+α*h,m+X1*h,b)
    y = f.(x)

    fa,fb = y[1],y[end]

    _i2 = (h/6) * (y[1] + y[13] + 5*(y[5]+y[9]))
    _i1 = (h/1470) * (77*(y[1]+y[13]) + 432*(y[3]+y[11]) + 625 * (y[2]+y[12]) + 672*y[7])

    _is = h * (A*(y[1]+y[13]) + B*(y[2]+y[12]) + C*(y[3]+y[11]) + D*(y[4]+y[10]) + E*(y[5]+y[9]) + F*(y[6]+y[8]) + G*y[7])

    sgn = sign(_is)
    sgn = ifelse(iszero(sgn),1,sgn)
    erri1 = abs(_i1 - _is)
    erri2 = abs(_i2 - _is)

    R = ifelse(!iszero(erri2), erri1/erri2, 1)

    if 0 < R < 1
        tol = tol/R
    end

    _is_new = sgn * abs(_is) * tol / eps()

    if iszero(_is_new)
        _is_new = b-a
    end

    flag = WarnFlag(false)
    Q = adaptlobstp(f,a,b,fa,fb,_is_new,trace,flag)
    return Q
end

function adaptlobstp(f,a,b,fa,fb,_is,trace,flag)
    m = (a+b)/2
    h = (b-a)/2

    mll = m-α*h
    ml = m-β*h
    mr = m+β*h
    mrr = m+α*h

    fmll = f(mll)
    fml = f(ml)
    fm = f(m)
    fmr = f(mr)
    fmrr = f(mrr)

    i2=(h/6)*(fa+fb+5*(fml+fmr))
    i1=(h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr)+672*fm)

    if ≈(_is+i1-i2, _is, atol=1e-16) || mll <= a || b <= mrr || !isfinite(_is)
        if (m <= a || b <= m) && !flag.nowarn
            @warn "Interval contains no more machine numbers. Required tolerance may not be met."
            flag.nowarn = true
        end
        Q = i1
        trace && @printf "%20.15f    %20.15f    %20.15f\n" a b-a Q
    else
        Q = adaptlobstp(f,a,mll,fa,fmll,_is,trace,flag) +
            adaptlobstp(f,mll,ml,fmll,fml,_is,trace,flag) +
            adaptlobstp(f,ml,m,fml,fm,_is,trace,flag) +
            adaptlobstp(f,m,mr,fm,fmr,_is,trace,flag) +
            adaptlobstp(f,mr,mrr,fmr,fmrr,_is,trace,flag) +
            adaptlobstp(f,mrr,b,fmrr,fb,_is,trace,flag)
    end
    return Q
end


text:
using MAT
using Random
include("../src/HyperLocal.jl")

## Read in the matrix
s = time()
@time M = matread("../data/AmazonReview5core_H.mat")
NodeLabels = vec(M["NodeLabels"])
NodeNames = M["NodeNames"]
LabelNames = M["LabelNames"]
H = M["H"]
order = round.(Int64,vec(sum(H,dims=2)))

## Extract hyperedges with only one node
e = findall(x->x>1,order)
H = H[e,:]
d = vec(sum(H,dims=1))
order = round.(Int64,vec(sum(H,dims=2)))
volA = sum(d)
m,n = size(H)
Ht = sparse(H')
toc = time()-s
println("Done loading things into memory in $toc seconds.")

## Run several sets of experiments
include("Smallest_Experiments.jl")
include("Medium_Experiments.jl")
include("Largest_Experiments.jl")

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
"""
```
  lsq_constrsparsereg(
      X         :: AbstractMatrix{T},
      y         :: AbstractVector{T},
      ρ         :: Union{AbstractVector, Number} = zero(T);
      Aeq       :: AbstractMatrix = zeros(T, 0, size(X, 2)),
      beq       :: Union{AbstractVector, Number} = zeros(T, size(Aeq, 1)),
      Aineq     :: AbstractMatrix = zeros(T, 0, size(X, 2)),
      bineq     :: Union{AbstractVector, number} = zeros(T, size(Aineq, 1)),
      obswt     :: AbstractVector = ones(T, length(y)),
      penwt     :: AbstractVector = ones(T, size(X, 2)),
      warmstart :: Bool = false,
      solver = ECOSSolver(maxit=10e8, verbose=0)
    )
```
Fit constrained lasso at fixed tuning parameter value(s) by minimizing
    `0.5sumabs2(√obswt .* (y - X * β)) + ρ * sumabs(penwt .* β)`
subject to linear constraints, using `Convex.jl`.

### Arguments
- `X`       : predictor matrix.
- `y`       : response vector.
- `ρ`       : tuning parameter. Can be a number or a list of numbers. Default 0.

### Optional arguments
- `Aeq`       : equality constraint matrix.
- `beq`       : equality constraint vector.
- `Aineq`     : inequality constraint matrix.
- `bineq`     : inequality constraint vector.
- `obswt`     : observation weights. Default is `[1 1 1 ... 1]`.
- `penwt`     : predictor penalty weights. Default is `[1 1 1 ... 1]`.
- `warmstart` :
- `solver`    : a solver Convex.jl supports. Default is ECOS.
              Note that Mosek and Gurobi are more robust than ECOS. Unlike ECOS or
              SCS, both Mosek and Gurobi require a license (free for academic
              use). <http://convexjl.readthedocs.io/en/latest/solvers.html>

### Returns
- `β`       : estimated coefficents.
- `objval`  : optimal objective value.
- `problem` : Convex.jl problem.
"""

function lsq_constrsparsereg(
    X::AbstractMatrix{T},
    y::AbstractVector{T},
    ρ::Union{AbstractVector, Number} = zero(T);
    Aeq::AbstractMatrix = zeros(T, 0, size(X, 2)),
    beq::Union{AbstractVector, Number} = zeros(T, size(Aeq, 1)),
    Aineq::AbstractMatrix = zeros(T, 0, size(X, 2)),
    bineq::Union{AbstractVector, Number} = zeros(T, size(Aineq, 1)),
    obswt::AbstractVector = ones(T, length(y)),
    penwt::AbstractVector = ones(T, size(X, 2)),
    warmstart::Bool = false,
    solver = ECOSSolver(maxit=10e8, verbose=0)
    ) where T

    n, p = size(X)

    β̂ = zeros(T, p, length(ρ))
    optval_vec = zeros(length(ρ))
    prob_vec = Convex.Problem[]

    β = Variable(p)
    loss = (1//2)sumsquares(sqrt(obswt) .* (y - X * β)) # loss term
    pen  = dot(penwt, abs(β)) # penalty term
    eqconstr   = Aeq * β == beq
    ineqconstr = Aineq * β <= bineq

    for i in 1:length(ρ)

      ρi = ρ[i]
      if ρi == typemax(typeof(ρi))
          problem = minimize(pen, eqconstr, ineqconstr)
      elseif ρi ≤ 0
          problem = minimize(loss, eqconstr, ineqconstr)
      else
          problem = minimize(loss + ρi * pen, eqconstr, ineqconstr)
      end

      if warmstart
          solve!(problem, solver; warmstart = i > 1? true : false)
          push!(prob_vec, problem)
          optval_vec[i] = problem.optval
      else
          solve!(problem, solver)
          push!(prob_vec, problem)
          optval_vec[i] = problem.optval
      end

      if length(ρ) == 1
        return β.value, problem.optval, problem
      end

      β̂[:, i] = β.value

    end # end of for loop

    return β̂, optval_vec, prob_vec

end # end of function

max_issues_repo_name: Hua-Zhou/ConstrainedLasso


text:
@testset "ZGenus" begin
  A = matrix(ZZ, 2, 2, [1, 1, 1, 1])
  @test (false, 1) == Hecke._iseven(A)
  A = matrix(ZZ, 2, 2, [2, 1, 1, 2])
  @test (true, -1) == Hecke._iseven(A)

  A = matrix(ZZ, 2, 2, [1, 2, 2, 3])
  @test (false, 1) == Hecke._iseven(A)
  @test (1, ZZ[-1;]) == Hecke._split_odd(A)

  A = matrix(ZZ, 2, 2, [1, 2, 2, 3])
  @test (1, ZZ[-1;]) ==  Hecke._split_odd(A)
  @test 0== Hecke._trace_diag_mod_8(A)

  A = matrix(ZZ, 2, 2, [1, 2, 2, 5])
  @test (1, ZZ[1;]) ==  Hecke._split_odd(A)
  @test (false, 1) == Hecke._iseven(A)
  @test 2== Hecke._trace_diag_mod_8(A)

  A = 2*diagonal_matrix(map(ZZ,[1, 2, 3, 4]))
  @test Hecke._p_adic_symbol(A, 2, 2) == [[1, 2, 3, 1, 4], [2, 1, 1, 1, 1], [3, 1, 1, 1, 1]]
  @test Hecke._p_adic_symbol(A, 3, 1) == [[0, 3, 1], [1, 1, -1]]

  A = diagonal_matrix(map(ZZ, [1, 2, 3, 4]))
  @test Hecke._two_adic_symbol(A, 2) == [[0, 2, 3, 1, 4], [1, 1, 1, 1, 1], [2, 1, 1, 1, 1]]

  #equality testing
  g1 = ZpGenus(2,[[0, 2, 7, 0, 0], [3, 1, 7, 1, 7]])
  g2 = ZpGenus(2,[[0, 2, 3, 0, 0], [3, 1, 3, 1, 3]])
  @test g1 != g2

  g1 = genus(A, 3)
  g2 = genus(A, 3)
  @test g1 == g2
  @test iseven(g1)
  @test 4==dim(g1)
  @test 4 == rank(g1)
  @test 4 == dim(g1)
  @test 6 == det(g1)
  g1 = genus(A,2)
  @test det(g1) == 3*8

  A = diagonal_matrix(map(ZZ,[1,3,-3]))
  @test excess(genus(A, 2)) == 2
  @test excess(genus(A, 3)) == 0
  @test excess(genus(A, 5)) == 0
  @test excess(genus(A, 7)) == 0
  @test excess(genus(A, 11)) == 0
  A = 2*diagonal_matrix(map(ZZ,[1,3,-3]))
  @test excess(genus(A, 2)) == 2
  @test excess(genus(A, 3)) == 0
  @test excess(genus(A, 5)) == 0
  @test excess(genus(A, 7)) == 0
  @test excess(genus(A, 11)) == 0

  A = 2*diagonal_matrix(map(ZZ,[1,2,3,4]))
  @test excess(genus(A, 2)) == 2
  @test excess(genus(A, 3)) == 6
  @test excess(genus(A, 5)) == 0
  @test excess(genus(A, 7)) == 0
  @test excess(genus(A, 11)) == 0

  A = diagonal_matrix(map(ZZ,[2, 4, 18]))
  G = genus(A, 2)
  @test level(G) == 4

  A = diagonal_matrix(map(ZZ,[2, 4, 18]))
  G = genus(A, 2)
  @test norm(G) == 2
  @test scale(G) == 2
  @test scale(genus(A, 3)) == 1
  A = ZZ[0 1; 1 0]
  @test norm(genus(A, 2)) == 2

  A = diagonal_matrix(map(ZZ, [1, 2, 3]))
  g = genus(A, 2)
  @test signature(g) == 5
  g = genus(A, 3)
  @test signature(g) == 1

  A = matrix(ZZ, 2, 2, [1,1,1,2])
  G2 = genus(A, 2)
  @test Hecke._is2adic_genus(G2._symbol)
  @test Hecke._is2adic_genus(G2)
  @test sprint(show, G2) isa String
  A = matrix(ZZ, 2, 2, [1,0,0,2])
  G = genus(A)
  @test 2 == det(G)
  @test !iseven(G)

  output =[[15, 2, 3, 0, 0],
         [15, 2, 7, 0, 0],
         [15, 2, 1, 1, 2],
         [15, 2, 5, 1, 6],
         [15, 2, 1, 1, 6],
         [15, 2, 5, 1, 2],
         [15, 2, 7, 1, 0],
         [15, 2, 3, 1, 4]]
  @test output == Hecke._blocks([15, 2, 0, 0, 0])

  @test size(Hecke._local_genera(2, 3, 1, 2, false))[1]==12
  @test size(Hecke._local_genera(2, 3, 1, 2, true))[1]==4
  @test size(Hecke._local_genera(5, 2, 2, 2, true))[1]==6

  @test length(genera((2,2), 10, even=true))==0  # check that a bug in catesian_product_iterator is fixed
  @test 4 == length(genera((4,0), 125, even=true))
  G = genus(diagonal_matrix(map(ZZ,[2, 4, 18])))
  @test 36 == level(G)
  G = genus(diagonal_matrix(map(ZZ,[2, 4, 18])))
  @test 2 == scale(G)
  @test 2 == norm(G)

  G = genus(diagonal_matrix(map(ZZ,[6, 4, 18])))
  @test norm(G) == 2
  @test sprint(show, G) isa String
  G = genus(matrix(ZZ, 2, 2, [0, 1, 1, 0]))
  @test norm(G) == 2

  @test all(hasse_invariant(g) == hasse_invariant(ambient_space(representative(g)),prime(g)) for g in Hecke._local_genera(2,3,1,2,false))
  @test all(hasse_invariant(g) == hasse_invariant(ambient_space(representative(g)),prime(g)) for g in Hecke._local_genera(2,5,4,4,true))
  @test all(hasse_invariant(g) == hasse_invariant(ambient_space(representative(g)),prime(g)) for g in Hecke._local_genera(3,2,2,2,true))
  @test all(hasse_invariant(g) == hasse_invariant(ambient_space(representative(g)),prime(g)) for g in Hecke._local_genera(3,3,4,4,true))
  @test all(hasse_invariant(g) == hasse_invariant(ambient_space(representative(g)),prime(g)) for g in Hecke._local_genera(5,2,2,2,true))

  A = diagonal_matrix(fmpz[2, -4, 6, 8])
  G = genus(A)
  q1 = discriminant_group(G)
  q2 = discriminant_group(Zlattice(gram=A))

  A = matrix(ZZ, 0, 0, [])
  g2 = genus(A, 2)
  g3 = genus(A, 3)
  g = genus(A)
  @test norm(g) == 0
  @test scale(g) == 0
  @test det(g) == 1
  @test rank(g) == 0
  @test level(g) == 1
  @test iseven(g)
  @test signature(g) == 0
  @test rank(representative(g))==0
  @test norm(g2) == 0
  @test scale(g2) == 0
  @test det(g2) == 1
  @test rank(g2) == 0
  @test iseven(g2)
  @test level(g2) == 1
  @test signature(g2) == 0
  @test excess(g2) == 0
  @test norm(g3) == 0
  @test scale(g3) == 0
  @test det(g3) == 1
  @test rank(g3) == 0
  @test signature(g3) == 0
  @test excess(g3) == 0
  @test rank(representative(g2))==0

  g3 = genus(diagonal_matrix(map(ZZ,[1,3,27])), 3)
  n3 = genus(matrix(ZZ,0,0,[]),3)
  @test g3 == orthogonal_sum(n3, g3)
  @test Hecke._species_list(g3) == [1, 1, 1]
  h3 = genus(diagonal_matrix(map(ZZ,[1,3,9])), 3)
  @test Hecke._standard_mass(h3) ==  9//16
  @test orthogonal_sum(g3,h3)==direct_sum(h3,g3)

  # These examples are taken from Table 2 of [CS1988]_::

  M_p = Hecke._M_p
  @test M_p(0, 2) == 1
  @test M_p(1, 2) == 1//2
  @test M_p(-2, 2) == 1//3
  @test M_p(2, 2) == 1
  @test M_p(3, 2) == 2//3
  @test M_p(-4, 2) == 8//15
  @test M_p(4, 2) == 8//9
  @test M_p(5, 2) == 32//45
  @test M_p(0, 3) == 1
  @test M_p(1, 3) == 1//2
  @test M_p(-2, 3) == 3//8
  @test M_p(2, 3) == 3//4
  @test M_p(3, 3) == 9//16
  @test M_p(-4, 3) == 81//160
  @test M_p(4, 3) == 81//128
  @test M_p(5, 3) == 729//1280
  @test M_p(0, 5) == 1
  @test M_p(1, 5) == 1//2
  @test M_p(-2, 5) == 5//12
  @test M_p(2, 5) == 5//8
  @test M_p(3, 5) == 25//48
  @test M_p(-4, 5) == 625//1248
  @test M_p(4, 5) == 625//1152

  A = diagonal_matrix(fmpz[1, 1, 1, 1])
  G = genus(A)
  @test Hecke._standard_mass_squared(G) == (1//48)^2

  @test Hecke._quadratic_L_function_squared(1, -4) == (1//4)^2
  @test Hecke._quadratic_L_function_squared(-4, -4) == (5//2)^2
  @test Hecke._quadratic_L_function_squared(2, 1) == (1//6)^2

  @test Hecke._zeta_exact(4)==1//90
  @test Hecke._zeta_exact(-3)==1//120
  @test Hecke._zeta_exact(0)==-1//2

  @test Hecke._gamma_exact(4)==6
  @test Hecke._gamma_exact(3)==2
  @test Hecke._gamma_exact(2)==1
  @test Hecke._gamma_exact(1)==1

  @test Hecke._gamma_exact(1//2)==1
  @test Hecke._gamma_exact(3//2)==1//2
  @test Hecke._gamma_exact(5//2)==3//4
  @test Hecke._gamma_exact(7//2)==15//8

  @test Hecke._gamma_exact(-1//2)==-2
  @test Hecke._gamma_exact(-3//2)==4//3
  @test Hecke._gamma_exact(-5//2)==-8//15
  @test Hecke._gamma_exact(-7//2)==16//105

  A = matrix(ZZ, 2, 2, [1, 1, 1, 2])
  G = genus(A)
  @test Hecke._isglobal_genus(G)
  G = genus(diagonal_matrix(fmpz[2, 2, 2, 2]))
  G._symbols[1]._symbol=[[0,2,3,0,0], [1,2,5,1,0]]
  @test !Hecke._isglobal_genus(G)

  L = Zlattice(gram=matrix(ZZ, 2, 2, [0, 1, 1, 0]))
  G = genus(L)
  g2 = genus(L, 2)
  g7 = genus(L, 7)
  @test local_symbol(G, 2) == g2
  @test local_symbol(G, 7) == g7
  q = discriminant_group(G) # corner case
  @test order(q) == 1
  L2 = Zlattice(gram=2*ZZ[2 1; 1 2])
  G2 = genus(L2)
  @test genus(orthogonal_sum(L,L2)[1]) == orthogonal_sum(G, G2)
  @test length(representatives(G2)) == 1

  G = genera((8,0), 1, even=true)[1]
  @test mass(G) == 1//696729600

  G = genus(diagonal_matrix(fmpz[1, 3, 9]),3)
  @test Hecke._mass_squared(G) == (9//8)^2

  # representatives, mass and genus enumeration
  DB = lattice_database()
  for i in 1:(long_test ? 200 : 10)
    L = lattice(DB,i)
    G = genus(L)
    q1 = quadratic_space(G)
    q2 = rational_span(L)
    @test Hecke.isisometric(q1, q2)
    L2 = representative(G)
    G2 = genus(L2)
    @test G==G2
  end

  genus_orders_sage = [[1, 1], [1], [1, 1], [1, 1, 1, 1], [1, 1], [1, 1], [1, 1], [1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 2, 2], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 2], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 2], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 2, 2, 2, 2], [1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [2, 2], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2], [1, 2], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1], [1, 1, 1, 1, 2, 2], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2], [1, 2], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 3, 3, 3, 3], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 2, 2, 2], [1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 2, 2], [2, 2], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 2, 2], [1, 2], [1, 1, 1, 1, 1, 1], [1, 2, 2, 2], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 2, 2], [1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 3, 3, 4, 4], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 2], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 2, 2], [1, 1, 1, 1], [1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1], [1, 2], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2], [1, 2], [1, 1, 1, 1, 1, 1], [1, 1]]

  for d in 1:(long_test ? 199 : 50)
    L = [length(representatives(G)) for G in genera((1,1), d)]
    @test genus_orders_sage[d] == sort!(L)
  end

  # some randomly chosen unimodular matrices
  B = [ZZ[7 -4; 2 -1],
  matrix(ZZ,3,3,[1, 1, -4, 0, 1, -3, 0, -1, 4]),
  matrix(ZZ,4,4, [0, -4, -3, -1, 1, -2, -3, -1, -1, -1, 1, 0, 1, -2, -4, 1])]

  sigdet = []
  for d  in 1:(long_test ? 32 : 10)
    for sig in [(2,0), (1,1), (1,2),(4,0)]
      push!(sigdet, (sig, d))
    end
  end

  if long_test
    for sig in [(0,3), (1,2), (2,2)]
      push!(sigdet, (sig, 2^7))
    end
  end

  for (sig,d) in sigdet
    for G in genera(sig, d)
      L = representative(G)
      spL = ambient_space(L)
      b = B[rank(L)-1]
      spLt = quadratic_space(QQ, b*gram_matrix(L)*transpose(b))
      # Our isisometric_with_isometry is too slow to handle the other cases
      if rank(L) <= 2
        flag, iso = Hecke.isisometric_with_isometry(spL,spLt)
        @test flag
        @test iso*gram_matrix(spLt)*transpose(iso) == gram_matrix(spL)
      end
      if isdefinite(L)
        # compare the two algorithms used to calculate the mass
        @test mass(L) == mass(G)
      end
      @test genus(L)==G
      D = discriminant_group(L)
      q1, _ = normal_form(D)
      q1 = Hecke.gram_matrix_quadratic(q1)
      q2 = discriminant_group(G)
      q2, _ = normal_form(q2)
      q2 = Hecke.gram_matrix_quadratic(q2)
      @test q1 == q2
      G2 = genus(D, sig)
      if iseven(G) == true
        @test isgenus(D, sig) == true
      end
      @test G == G2
      # test local representations
      if rank(L) >= 2
        diag = diagonal_matrix(fmpq[1, 2])*basis_matrix(L)[1:2,1:end]
        sub = lattice(ambient_space(L), diag)
        g = genus(sub)
        @test Hecke.represents(G, genus(sub))
      end
      if rank(L) >= 3
        diag = diagonal_matrix(fmpq[1, 2, 4])*basis_matrix(L)[1:3,1:end]
        sub = lattice(ambient_space(L), diag)
        g = genus(sub)
        @test Hecke.represents(G, genus(sub))
      end
      if rank(L) >= 3
        diag = diagonal_matrix(fmpq[4, 2, 2])*basis_matrix(L)[1:3,1:end]
        sub = lattice(ambient_space(L), diag)
        g = genus(sub)
        @test Hecke.represents(G, genus(sub))
      end
    end
  end

  for d in 1:(long_test ? 50 : 10)
    for sig in [(2,0),(3,0),(4,0)]
      for G in genera(sig,d)
        m = mass(G)
        L = representative(G)
        @test genus(L)==G
        @test mass(L)==m
        rep = genus_representatives(L)
        @test sum(1//automorphism_group_order(M) for M in rep)==m
         q = ambient_space(L)
         for r in rep
           qr = ambient_space(r)
           #b, i = Hecke.isisometric_with_isometry(q,qr)
           @test isisometric(q, qr)
           #@test b
           #@test i*gram_matrix(qr)*transpose(i) == gram_matrix(q)
         end
      end
    end
  end

end


text:
## 16.jl
# ファイルをn分割

fname = "./dat/hightemp.txt"
f = open(fname)
data = readlines(f)
n = parse(Int, ARGS[1]) #分割数
sep = Int(floor(length(data)/n)) #分割幅

for i in 1:n
  o = open(replace(fname, ".txt", "_part"*string(i)*".txt"), "w")
  i != n && print(o, join(data[(i-1)*sep+1:sep*i],""))
  i == n && print(o, join(data[sep*i:end],""))
  close(o)
end
close(f)


text:
module MultivariateStats_

export RidgeRegressor, PCA, KernelPCA, ICA, LDA, BayesianLDA

import MLJBase
import MLJBase: @mlj_model, metadata_model, metadata_pkg
import StatsBase: proportions, CovarianceEstimator
using Distances, LinearAlgebra
using Tables, ScientificTypes

import MultivariateStats

const MS = MultivariateStats

struct LinearFitresult{F} <: MLJBase.MLJType
    coefficients::Vector{F}
    intercept::F
end

const RIDGE_DESCR = "Ridge regressor with regularization parameter lambda. Learns a linear regression with a penalty on the l2 norm of the coefficients."
const PCA_DESCR = "Principal component analysis. Learns a linear transformation to project the data  on a lower dimensional space while preserving most of the initial variance."
const KPCA_DESCR = "Kernel principal component analysis."
const ICA_DESCR = "Independent component analysis."
const LDA_DESCR = "Multiclass linear discriminant analysis. The algorithm learns a projection matrix `W` that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance."
const BayesianLDA_DESCR= "Bayesian Multiclass linear discriminant analysis. The algorithm learns a projection matrix `W` that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance and classifies using Bayes rule."
const BayesianSubspaceLDA_DESCR = "Bayesian Multiclass linear discriminant analysis. Suitable for high dimensional data. The algorithm learns a projection matrix `W` (without computing covariance matrices Sw ,Sb) in the subspace spanned by the within class scatter matrix, that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance and classifies using Bayes rule."
const SubspaceLDA_DESCR = "Multiclass linear discriminant analysis. Suitable for high dimensional data. The algorithm learns a projection matrix `W` (without computing covariance matrices Sw ,Sb) in the subspace spanned by the within class scatter matrix, that projects the feature matrix `Xtrain` onto a lower dimensional space of dimension `out_dim` such that the between-class variance in the transformed space is maximized relative to the within-class variance."

####
#### RIDGE
####

"""
RidgeRegressor(; lambda=1.0)

$RIDGE_DESCR

## Parameters

* `lambda=1.0`: non-negative parameter for the regularization strength.
"""
@mlj_model mutable struct RidgeRegressor <: MLJBase.Deterministic
    lambda::Real = 1.0::(_ ≥ 0)
end

function MLJBase.fit(model::RidgeRegressor, verbosity::Int, X, y)
    Xmatrix   = MLJBase.matrix(X)
    features  = Tables.schema(X).names
    θ         = MS.ridge(Xmatrix, y, model.lambda)
    coefs     = θ[1:end-1]
    intercept = θ[end]

    fitresult = LinearFitresult(coefs, intercept)
    report    = NamedTuple()
    cache     = nothing

    return fitresult, cache, report
end

MLJBase.fitted_params(::RidgeRegressor, fr) =
    (coefficients=fr.coefficients, intercept=fr.intercept)

function MLJBase.predict(::RidgeRegressor, fr, Xnew)
    Xmatrix = MLJBase.matrix(Xnew)
    return Xmatrix * fr.coefficients .+ fr.intercept
end

####
#### PCA
####

const PCAFitResultType = MS.PCA

"""
PCA(; maxoutdim=nothing, method=:auto, pratio=0.99, mean=nothing)

$PCA_DESCR

## Parameters

* `maxoutdim=nothing`: maximum number of output dimensions, unconstrained if nothing.
* `method=:auto`: method to use to solve the problem, one of `:auto`, `:cov` or `:svd`
* `pratio=0.99`: ratio of variance preserved
* `mean=nothing`: if set to nothing centering will be computed and applied, if set to `0` no centering (assumed pre-centered), if a vector is passed, the centering is done with that vector.
"""
@mlj_model mutable struct PCA <: MLJBase.Unsupervised
    maxoutdim::Union{Nothing,Int} = nothing::(_ === nothing || _ ≥ 1)
    method::Symbol  = :auto::(_ in (:auto, :cov, :svd))
    pratio::Float64 = 0.99::(0.0 < _ ≤ 1.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MLJBase.fit(model::PCA, verbosity::Int, X)
    Xarray = MLJBase.matrix(X)
    mindim = minimum(size(Xarray))

    maxoutdim = model.maxoutdim === nothing ? mindim : model.maxoutdim

    # NOTE: copy/transpose
    fitresult = MS.fit(MS.PCA, permutedims(Xarray);
                       method=model.method,
                       pratio=model.pratio,
                       maxoutdim=maxoutdim,
                       mean=model.mean)

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult),
              principalvars=MS.principalvars(fitresult),
              tprincipalvar=MS.tprincipalvar(fitresult),
              tresidualvar=MS.tresidualvar(fitresult),
              tvar=MS.tvar(fitresult))

    return fitresult, cache, report
end

MLJBase.fitted_params(::PCA, fr) = (projection=fr,)

function MLJBase.transform(::PCA, fr::PCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MLJBase.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MLJBase.table(Xnew, prototype=X)
end

####
#### KernelPCA
####

const KernelPCAFitResultType = MS.KernelPCA

default_kernel = (x, y) -> x'y

"""
KernelPCA(; maxoutdim=nothing, kernel=(x,y)->x'y, solver=:auto, pratio=0.99, mean=nothing)

$KPCA_DESCR

## Parameters

* `maxoutdim=nothing`: maximum number of output dimensions, unconstrained if nothing.
* `kernel=nothing`: kernel function of 2 vector arguments x and y, returns a scalar value, (x,y)->x'y if nothing
* `solver=:auto`: solver to use for the eigenvalues, one of `:eig` (default), `:eigs`
* `inverse=false`: perform calculation for inverse transform
* `beta=1.0`: strength of the ridge regression that learns the inverse transform when inverse is true
* `tol=0.0`: Convergence tolerance for eigs solver
* `maxiter=300`: maximum number of iterations for eigs solver
"""
@mlj_model mutable struct KernelPCA <: MLJBase.Unsupervised
    maxoutdim::Union{Nothing,Int} = nothing::(_ === nothing || _ ≥ 1)
    kernel::Function = default_kernel
    solver::Symbol   = :eig::(_ in (:eig, :eigs))
    inverse::Bool    = false
    beta::Real       = 1.0::(_ ≥ 0.0)
    tol::Real        = 1e-6::(_ ≥ 0.0)
    maxiter::Int     = 300::(_ ≥ 1)
end

function MLJBase.fit(model::KernelPCA, verbosity::Int, X)
    Xarray = MLJBase.matrix(X)
    mindim = minimum(size(Xarray))
    # default max out dim if not given
    maxoutdim = model.maxoutdim === nothing ? mindim : model.maxoutdim

    fitresult = MS.fit(MS.KernelPCA, permutedims(Xarray);
                       kernel=model.kernel,
                       maxoutdim=maxoutdim,
                       solver=model.solver,
                       inverse=model.inverse,
                       β=model.beta,
                       tol=model.tol,
                       maxiter=model.maxiter)

    cache  = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              projection=MS.projection(fitresult),
              principalvars=MS.principalvars(fitresult))

    return fitresult, cache, report
end

MLJBase.fitted_params(::KernelPCA, fr) = (projection=fr,)

function MLJBase.transform(::KernelPCA, fr::KernelPCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MLJBase.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MLJBase.table(Xnew, prototype=X)
end

####
#### ICA
####

const ICAFitResultType = MS.ICA

"""
ICA(; k=0, alg=:fastica, fun=:tanh, do_whiten=true, maxiter=100,
      tol=1e-6, mean=nothing, winit=zeros(0,0))

$ICA_DESCR

## Parameters

* `k=0`: number of independent components to recover, set automatically if `0`
* `alg=:fastica`: algorithm to use (only `:fastica` is supported at the moment)
* `fun=:tanh`: approximate neg-entropy functor, via the function `MultivariateStats.icagfun`, one of `:tanh` and `:gaus`
* `do_whiten=true`: whether to perform pre-whitening
* `maxiter=100`: maximum number of iterations
* `tol=1e-6`:  convergence tolerance for change in matrix W
* `mean=nothing`: mean to use, if nothing centering is computed and applied, if zero, no centering, a vector of means can be passed
* `winit=zeros(0,0)`: initial guess for matrix `W` either an empty matrix (random init), a matrix of size `k × k` (if `do_whiten` is true), a matrix of size `m × k` otherwise
"""
@mlj_model mutable struct ICA <: MLJBase.Unsupervised
    k::Int          = 0::(_ ≥ 0)
    alg::Symbol     = :fastica::(_ in (:fastica,))
    fun::Symbol     = :tanh::(_ in (:tanh, :gaus))
    do_whiten::Bool = true
    maxiter::Int    = 100
    tol::Real       = 1e-6::(_ ≥ 0.0)
    winit::Union{Nothing,Matrix{<:Real}} = nothing
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MLJBase.fit(model::ICA, verbosity::Int, X)
    Xarray = MLJBase.matrix(X)
    n, p   = size(Xarray)

    m = min(n, p)
    k = ifelse(model.k ≤ m, model.k, m)

    fitresult = MS.fit(MS.ICA, permutedims(Xarray), k;
                       alg=model.alg,
                       fun=MS.icagfun(model.fun),
                       do_whiten=model.do_whiten,
                       maxiter=model.maxiter,
                       tol=model.tol,
                       mean=model.mean,
                       winit=ifelse(model.winit === nothing, zeros(0,0), model.winit))

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult))

    return fitresult, cache, report
end

MLJBase.fitted_params(::ICA, fr) = (projection=fr,)

function MLJBase.transform(::ICA, fr::ICAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MLJBase.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MLJBase.table(Xnew, prototype=X)
end

####
#### MulticlassLDA
####

"""
LDA(; kwargs...)

$LDA_DESCR

## Parameters

* `method=:gev`: choice of solver, one of `:gevd` or `:whiten` methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance, by default set to the standard `MultivariateStats.CovarianceEstimator` but could be set to any robust estimator from `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class covariance.
* `out_dim`: the output dimension, i.e dimension of the transformed space, automatically set if 0 is given (default).
* `regcoef`: regularization coefficient. A positive value `regcoef * eigmax(Sw)` where `Sw` is the within-class covariance estimator, is added to the diagonal of Sw to improve numerical stability. This can be useful if using the standard covariance estimator.
* `dist=SqEuclidean`: the distance metric to use when performing classification (to compare the distance between a new point and centroids in the transformed space), an alternative choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct LDA <: MLJBase.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Real    = 0.0::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MLJBase.fit(model::LDA, ::Int, X, y)
    class_list = MLJBase.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MLJBase.matrix(X, transpose=true) # now p x n matrix
    yplain = MLJBase.int(y) # vector of n ints in {1,..., nclasses}
    p      = size(Xm_t, 1)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    core_res = MS.fit(MS.MulticlassLDA, nclasses, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef=model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    cache     = nothing
    report    = NamedTuple{}()
    fitresult = (core_res, class_list)

    return fitresult, cache, report
end

function MLJBase.fitted_params(::LDA, (core_res, class_list))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res))
end

function MLJBase.predict(m::LDA, (core_res, class_list), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    XWt = MLJBase.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `n x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)
    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return [MLJBase.UnivariateFinite(class_list, P[j, :]) for j in 1:size(P, 1)]
end

####
#### BayesianLDA
####

"""
BayesianLDA(; kwargs...)

$BayesianLDA_DESCR

## Parameters

* `method=:gev`: choice of solver, one of `:gevd` or `:whiten` methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance, by default set to the standard `MultivariateStats.CovarianceEstimator` but could be set to any robust estimator from `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class covariance.
* `out_dim`: the output dimension, i.e dimension of the transformed space, automatically set if 0 is given (default).
* `regcoef`: regularization coefficient. A positive value `regcoef * eigmax(Sw)` where `Sw` is the within-class covariance estimator, is added to the diagonal of Sw to improve numerical stability. This can be useful if using the standard covariance estimator.
* `priors=nothing`: if `priors = nothing` estimates the prior probabilities from the data else it uses the user specified `UnivariateFinite` prior probabilities.
See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct BayesianLDA <: MLJBase.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Real    = 0.0::(_ ≥ 0)
    priors::Union{Nothing, MLJBase.UnivariateFinite} = nothing
end

function MLJBase.fit(model::BayesianLDA, ::Int, X, y)
    class_list = MLJBase.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MLJBase.matrix(X, transpose=true) # now p x n matrix
    yplain = MLJBase.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    ## Estimates prior probabilities if unspecified by user.
    if model.priors == nothing
        priors = proportions(yplain)
    else
        #check if the length of priors is same as nclasses
        length(MLJBase.levels(model.priors)) == nclasses || throw(ArgumentError("Invalid size of `priors`"))
        priors = MLJBase.pdf.(model.priors, class_list)
    end

    core_res = MS.fit(MS.MulticlassLDA, nclasses, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef=model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    ## The original projection matrix satisfies Pᵀ*Sw*P=I
    ## scaled projection_matrix and core_res.proj by multiplying by sqrt(n - nclasses) this ensures Pᵀ*Σ*P=I
    ## where covariance estimate Σ = Sw / (n - nclasses)
    core_res.proj   .*= sqrt(n - nclasses)
    core_res.pmeans .*= sqrt(n - nclasses)

    cache     = nothing
    report    = NamedTuple{}()
    fitresult = (core_res, class_list, priors)

    return fitresult, cache, report
end

function MLJBase.fitted_params(::BayesianLDA, (core_res, class_list, priors))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res),
            priors            = priors)
end

function MLJBase.predict(m::BayesianLDA, (core_res, class_list, priors), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    XWt = MLJBase.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    n  = core_res.stats.tweight # n is the Number of training examples
    nclasses = size(class_list)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `P` is of dimension `n x nc`
    #  P[i,k] = -0.5*(xᵢ −  µₖ)ᵀΣ⁻¹(xᵢ −  µₖ) + log(priorsₖ) and Σ⁻¹ = I due to the nature of the projection_matrix
    P = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    P .*= -0.5
    P .+= log.(priors)'

    # apply a softmax transformation to convert P to a probability matrix
    P .-= maximum(P, dims=2)
    P  .= exp.(P)
    P ./= sum(P, dims=2)

    return [MLJBase.UnivariateFinite(class_list, P[j, :]) for j in 1:size(P, 1)]
end

####
#### BayesianSubspaceLDA
####

"""
BayesianSubspaceLDA(; kwargs...)

$BayesianSubspaceLDA_DESCR

## Parameters

* `normalize=true`: Option to normalize the between class variance by number of observations in each class, one of `true` or `false`.
* `out_dim`: the output dimension, i.e dimension of the transformed space, automatically set if 0 is given (default).
* `priors=nothing`: if `priors = nothing` estimates the prior probabilities from the data, else it uses the user specified `UnivariateFinite` prior probabilities.
See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006), "Generalizing discriminant analysis using the generalized singular value decomposition", IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct BayesianSubspaceLDA <: MLJBase.Probabilistic
    normalize::Bool=true
    out_dim::Int   = 0::(_ ≥ 0)
    priors::Union{Nothing, MLJBase.UnivariateFinite} = nothing
end

function MLJBase.fit(model::BayesianSubspaceLDA, ::Int, X, y)
    class_list = MLJBase.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MLJBase.matrix(X, transpose=true) # now p x n matrix
    yplain = MLJBase.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    # check to make sure sample size is greater or equals nclasses
    n ≥ nclasses || throw(ArgumentError("The number of samples is less than the number of classes"))

    ## Estimates prior probabilities if unspecified by user.
    if model.priors == nothing
        priors = proportions(yplain)
    else
        #check if the length of priors is same as nclasses
        length(MLJBase.levels(model.priors)) == nclasses || throw(ArgumentError("Invalid size of `priors`"))
        priors = MLJBase.pdf.(model.priors, class_list)
    end

    # # Code gotten from JuliaStats/MultivariateStats.jl
    # # Compute centroids, class weights, and deviation from centroids
    # # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = MS.center(Xm_t, Int.(yplain), nclasses) # Hw is p x n matrix , cmeans is p x nclasses matrix, cweights is nclasses x 1 vector
    dmeans = cmeans .- (model.normalize ? MS.mean(cmeans, dims=2) : cmeans * (cweights / (n))) #dmeans is p x nclasses
    Hb = model.normalize ? dmeans .* sqrt(n)  : dmeans .* sqrt.(cweights) #Hb is p x nclasses

    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full = false)
    keep = Σw .> sqrt(eps()) * maximum(Σw)
    projw = Uw[:, keep]
    pHb = projw' * Hb
    pHw = projw' * Hw

    # compute LDA projection G in the subspace spanned by the within-class variance.
    λ, G = MS.lda_gsvd(pHb, pHw, cweights) # λ is a (nclasses -1) x 1 vector containing the eigen values sorted in descending order.

    G = G[:, 1:out_dim]

    core_res = MS.SubspaceLDA(projw, G, λ, cmeans, cweights)

    ## scaled the overall projection_matrix (P = projw * G) by multiplying by sqrt(n - nclasses) this ensures Pᵀ*Σ*P=I
    ## where covariance estimate Σ = Sw / (n - nclasses)
    core_res.projLDA   .*= sqrt(n - nclasses)

    vproportions = λ ./ sum(λ) #proportions of variance
    cache     = nothing
    report    = (vproportions = vproportions,)
    fitresult = (core_res, class_list, priors)

    return fitresult, cache, report
end

function MLJBase.fitted_params(::BayesianSubspaceLDA, (core_res, class_list, priors))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res),
            priors            = priors)
end

function MLJBase.predict(m::BayesianSubspaceLDA, (core_res, class_list, priors), Xnew)
    # projection of Xnew XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MLJBase.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `P` is of dimension `n x nc`
    #  P[i,k] = -0.5*(xᵢ −  µₖ)ᵀΣ⁻¹(xᵢ −  µₖ) + log(priorsₖ) and Σ⁻¹ = I due to the nature of the projection_matrix (P = projw * G)
    P = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    P .*= -0.5
    P .+= log.(priors)'

    # apply a softmax transformation to convert P to a probability matrix
    P .-= maximum(P, dims=2)
    P  .= exp.(P)
    P ./= sum(P, dims=2)

    return [MLJBase.UnivariateFinite(class_list, P[j, :]) for j in 1:size(P, 1)]
end

####
#### SubspaceLDA
####

"""
SubspaceLDA(; kwargs...)

$SubspaceLDA_DESCR

## Parameters

* `normalize=true`: Option to normalize the between class variance by number of observations in each class, one of `true` or `false`.
* `out_dim`: the output dimension, i.e dimension of the transformed space, automatically set if 0 is given (default).
* `dist=SqEuclidean`: the distance metric to use when performing classification (to compare the distance between a new point and centroids in the transformed space), an alternative choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006), "Generalizing discriminant analysis using the generalized singular value decomposition", IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct SubspaceLDA <: MLJBase.Probabilistic
    normalize::Bool=true
    out_dim::Int   = 0::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MLJBase.fit(model::SubspaceLDA, ::Int, X, y)
    features   = MLJBase.schema(X).names
    class_list = MLJBase.classes(y[1]) #class list containing unique entries in y
    nclasses   = length(class_list)

    # NOTE: copy/transpose
    Xm_t   = MLJBase.matrix(X, transpose=true) # now p x n matrix
    yplain = MLJBase.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)

     #check to make sure we have more than one class
    nclasses >= 2 ||  throw(ArgumentError("The number of classes has to be greater than one"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nclasses - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ def_outdim || throw(ArgumentError("`out_dim` must not be larger than `min(p, nc-1)` where `p` is the dimension of `X` and `nc` is the number of classes."))

    # check to make sure sample size is greater or equals nclasses
    n ≥ nclasses || throw(ArgumentError("The number of samples is less than the number of classes"))

    # # Code gotten from JuliaStats/MultivariateStats.jl
    # # Compute centroids, class weights, and deviation from centroids
    # # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = MS.center(Xm_t, Int.(yplain), nclasses) # Hw is p x n matrix , cmeans is p x nclasses matrix, cweights is nclasses x 1 vector
    dmeans = cmeans .- (model.normalize ? MS.mean(cmeans, dims=2) : cmeans * (cweights / (n))) #dmeans is p x nclasses
    Hb = model.normalize ? dmeans .* sqrt(n)  : dmeans .* sqrt.(cweights)

    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full = false)
    keep = Σw .> sqrt(eps()) * maximum(Σw)
    projw = Uw[:, keep]
    pHb = projw' * Hb
    pHw = projw' * Hw

    # compute LDA projection G in the subspace spanned by the within-class variance.
    λ, G = MS.lda_gsvd(pHb, pHw, cweights) # λ is a (nclasses -1) x 1 vector containing the eigen values sorted in descending order.

    G = G[:, 1:out_dim]

    core_res = MS.SubspaceLDA(projw, G, λ, cmeans, cweights)
    vproportions = λ ./ sum(λ) #proportions of variance

    cache     = nothing
    report    = (vproportions = vproportions,)
    fitresult = (core_res, class_list, features)

    return fitresult, cache, report
end

function MLJBase.fitted_params(::SubspaceLDA, (core_res, class_list))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res))
end

function MLJBase.predict(m::SubspaceLDA, (core_res, class_list), Xnew)
    # projection of Xnew, XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MLJBase.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `n x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)
    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return [MLJBase.UnivariateFinite(class_list, P[j, :]) for j in 1:size(P, 1)]
end

function MLJBase.transform(m::T, (core_res,), X) where T<:Union{LDA, SubspaceLDA, BayesianLDA, BayesianSubspaceLDA}
    # projection of X, XWt is n x o  where o = number of out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MLJBase.matrix(X) * proj
    return MLJBase.table(XWt, prototype = X)
end

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

metadata_pkg.((RidgeRegressor, PCA, KernelPCA, ICA, LDA, BayesianLDA, SubspaceLDA, BayesianSubspaceLDA),
              name="MultivariateStats",
              uuid="6f286f6a-111f-5878-ab1e-185364afe411",
              url="https://github.com/JuliaStats/MultivariateStats.jl",
              license="MIT",
              julia=true,
              is_wrapper=false)

metadata_model(RidgeRegressor,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{MLJBase.Continuous},
               weights=false,
               descr=RIDGE_DESCR)

metadata_model(PCA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=MLJBase.Table(MLJBase.Continuous),
               weights=false,
               descr=PCA_DESCR)

metadata_model(KernelPCA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=MLJBase.Table(MLJBase.Continuous),
               weights=false,
               descr=KPCA_DESCR)

metadata_model(ICA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=MLJBase.Table(MLJBase.Continuous),
               weights=false,
               descr=ICA_DESCR)

metadata_model(LDA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{<:MLJBase.Finite},
               weights=false,
               output=MLJBase.Table(MLJBase.Continuous),
               descr=LDA_DESCR)
metadata_model(BayesianLDA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{<:MLJBase.Finite},
               weights=false,
               output=MLJBase.Table(MLJBase.Continuous),
               descr=BayesianLDA_DESCR)
metadata_model(SubspaceLDA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{<:MLJBase.Finite},
               weights=false,
               output=MLJBase.Table(MLJBase.Continuous),
               descr=SubspaceLDA_DESCR)
metadata_model(BayesianSubspaceLDA,
               input=MLJBase.Table(MLJBase.Continuous),
               target=AbstractVector{<:MLJBase.Finite},
               weights=false,
               output = MLJBase.Table(MLJBase.Continuous),
               descr=BayesianSubspaceLDA_DESCR)

end # of module+


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule obstack_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("obstack")
JLLWrappers.@generate_main_file("obstack", UUID("c88a4935-d25e-5644-aacc-5db6f1b8ef79"))
end  # module obstack_jll


text:
using WriteVTK
#using VTKDataIO

"""
    DataStorage

Storage for paraview input/output.
"""
mutable struct DataStorage
    file::WriteVTK.CollectionFile
    path::String
    frame::Int64
end

"""
    new_pvd_file(path::String)::DataStorage

Cteates new `DataStorage` in a given `path`.
"""
function new_pvd_file(path::String)::DataStorage
    if !ispath(path)
        mkpath(path)
        @info("created new path: "*path)
    end
    return DataStorage(paraview_collection(path*"/result"), path, 0)
end

"""
    save_pvd_file(data::DataStorage)

Saves and closes `DataStorage`.
"""
function save_pvd_file(data::DataStorage)
    vtk_save(data.file)
end

function capture_frame(sys::ParticleSystem, data::DataStorage)
    N = length(sys.particles)
    points = zeros(Float64, 3, N)
    for k in 1:N, i in 1:3
        points[i,k] = sys.particles[k].x[i]
    end
    cells = [MeshCell(PolyData.Verts(), [i]) for i in 1:N]
    vtk_file = vtk_grid(data.path*"/frame"*string(data.frame), points, cells)
    return vtk_file
end

"""
    save_frame!(data::DataStorage, sys::ParticleSystem, vars::DataField...)

Inserts one time frame into a `DataStorage` that includes all `vars...` as fields.
"""
function save_frame!(data::DataStorage, sys::ParticleSystem, vars::Symbol...)
    frame = capture_frame(sys, data)
    for var in vars
        field = ParticleField(sys, var)
        Type = attribute_type(get_particle_type(sys), var)
        N = length(sys.particles)
        if Type <: Number
            frame[string(var)] = field
        elseif Type <: RealVector
            vals = zeros(Float64, 3, N)
            for k in 1:N, i in 1:3
                vals[i,k] = field[k][i]
            end
            frame[string(var)] = vals
        elseif Type <: RealMatrix
            vals = zeros(Float64, 3, 3, N)
            for k in 1:N, i in 1:3, j in 1:3
                vals[i,j,k] = field[k][i,j]
            end
            frame[string(var)] = vals
        else
            @error("Cannot export type "*string(Type)*" to VTK.")
        end
    end
    data.file[data.frame] = frame
    data.frame += 1
end

#=
"""
    import_particles!(sys::ParticleSystem, path::String, particle_constructor::Function)

Imports particles from a vtk file in 'path' using 'constructor'.
"""
function import_particles!(sys::ParticleSystem, path::String, particle_constructor::Function)
    input = read_vtk(path)
    (_, N) = size(input.point_coords)
    resize!(sys.particles, N)
    for i in 1:N
        x = [input.point_coords[1,i] for i in 1:3]
        sys.particles[i] = particle_constructor(RealVector(x))
        for field in fieldnames(sys.particle_type)
            key = string(field)
            #skip particle data not present in file
            if !haskey(input.point_data, key)
                continue
            end
            Type = attribute_type(get_particle_type(sys), field)
            if Type <: Number
                val = input.point_data[key][i]
                setproperty!(sys.particles[i], field, val)
            elseif Type <: RealVector
                val = Type([input.point_data[key][j,i] for j in 1:3])
                setproperty!(sys.particles[i], field, val)
            elseif Type <: RealMatrix
                val = Type([input.point_data[key][j,k] for j in 1:3, k in 1:3])
                setproperty!(sys.particles[i], field, val)
            else
                throw("Cannot import data field with type:" *string(Type))
            end
        end
    end
end
=#


text:
function check_int(x)
    if x > 0
        println("$x is positive")
    elseif x < 0
        println("$x is negative")
    else
        println("$x is zero")
    end
end

check_int(10)
check_int(-10)
check_int(0)

for x in 1:5
    @show x
end

for x=1:5
    @show x
end

x = 0
while x < 5
    @show x
    x += 1
end


text:
export ViscoelasticitySolver

@doc raw"""
    ViscoelasticitySolver(globdat::GlobalData, domain::Domain,
        d0::Union{Array{Float64, 1}, PyObject},
        v0::Union{Array{Float64, 1}, PyObject},
        a0::Union{Array{Float64, 1}, PyObject},
        σ0::Union{Array{Float64, 2}, PyObject},
        ε0::Union{Array{Float64, 2}, PyObject},
        Δt::Float64, NT::Int64,
        μ::Union{Array{Float64, 1}, PyObject},
        λ::Union{Array{Float64, 1}, PyObject},
        η::Union{Array{Float64, 1}, PyObject},
        Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
        abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)

Solves the Maxwell viscoelasticity model using the generalized $\alpha$ scheme.

The constitutive relation has the form

$$\dot \sigma_{ij} + \frac{\mu}{\eta} \left( \sigma_{ij} - \frac{\sigma_{kk}}{3}\delta_{ij} \right) = 2\mu \dot \epsilon_{ij} + \lambda \dot\epsilon_{kk}\delta_{ij}$$

along with the balance of linear momentum equation

$$\mathrm{div}\ \sigma_{ij,j} + \rho f_i = \rho \ddot u_i$$

Users need to provide $\mu$, $\eta$, $\lambda$ and appropriate boundary/initial conditions.
"""
function ViscoelasticitySolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject},
    v0::Union{Array{Float64, 1}, PyObject},
    a0::Union{Array{Float64, 1}, PyObject},
    σ0::Union{Array{Float64, 2}, PyObject},
    ε0::Union{Array{Float64, 2}, PyObject},
    Δt::Float64, NT::Int64,
    μ::Union{Array{Float64, 1}, PyObject},
    λ::Union{Array{Float64, 1}, PyObject},
    η::Union{Array{Float64, 1}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)
    N = getNGauss(domain)
    @assert length(μ)==length(λ)==length(η)==N
    S, H = compute_maxwell_viscoelasticity_matrices(μ, λ, η, Δt)

    assembleMassMatrix!(globdat, domain)
    @assert 0<=ρ<=1
    init_nnfem(domain)
    αm = (2ρ-1)/(1+ρ)
    αf = ρ/(1+ρ)
    β2 = 0.5*(1 - αm + αf)^2
    γ = 0.5 - αm + αf

    Fext, ubd, abd = convert_to_tensor([Fext, ubd, abd], [Float64, Float64, Float64])
    M = constant(globdat.M)

    stiff = s_compute_stiffness_matrix(H, domain)
    A = M*(1 - αm) +  (1 - αf) * 0.5 * β2 * Δt^2 * stiff
    A = factorize(A)
    bddof = findall(domain.EBC[:] .== -2)
    fixed_bddof = findall(domain.EBC[:] .== -1)
    nbddof = findall(domain.dof_to_eq)

    function condition(i, tas...)
        i<=NT
    end
    function body(i, tas...)
        d_arr, v_arr, a_arr, σ_arr, ε_arr = tas
        u, ∂u, ∂∂u, σc, εc = read(d_arr, i), read(v_arr, i), read(a_arr, i), read(σ_arr, i), read(ε_arr, i)
        if ismissing(Fext)
            fext = zeros(domain.neqs)
        else
            fext = Fext[i]
        end
        ∂∂up = ∂∂u
        up =  (1 - αf)*(u + Δt*∂u + 0.5 * Δt * Δt * ((1 - β2)*∂∂u + β2*∂∂up)) + αf*u
        if !ismissing(abd)
            up = scatter_update(up, bddof, ubd[i])
        end
        if length(sum(fixed_bddof))>0
            up = scatter_update(up, fixed_bddof, d0[fixed_bddof])
        end

        # Here ϵc should be interpreted as the true values at (1-αf)*Δt
        ε = s_eval_strain_on_gauss_points(up, domain)
        σ = batch_matmul(H, ε-εc) + batch_matmul(S, σc)

        fint  = s_compute_internal_force_term(σ, domain)

        res = M * (∂∂up[nbddof] *(1 - αm) + αm*∂∂u[nbddof])  + fint - fext
        Δ = -(A\res)
        ∂∂up= scatter_add(∂∂up, nbddof, Δ)
        if !ismissing(abd)
            ∂∂up = scatter_update(∂∂up, bddof, abd[i])
        end

        # updaet
        u += Δt * ∂u + Δt^2/2 * ((1 - β2) * ∂∂u + β2 * ∂∂up)
        ∂u += Δt * ((1 - γ) * ∂∂u + γ * ∂∂up)

        # We use ϵ to approximate values at (1-αf)*Δt at this time step.
        ε = s_eval_strain_on_gauss_points(u, domain)
        σ = batch_matmul(H, ε-εc) + batch_matmul(S, σc)

        if length(sum(fixed_bddof))>0
            u = scatter_update(u, fixed_bddof, d0[fixed_bddof])
        end

        i+1, write(d_arr, i+1, u), write(v_arr, i+1, ∂u), write(a_arr, i+1, ∂∂up), write(σ_arr, i+1, σ), write(ε_arr, i+1, ε)
    end

    arr_d = TensorArray(NT+1); arr_d = write(arr_d, 1, d0)
    arr_v = TensorArray(NT+1); arr_v = write(arr_v, 1, v0)
    arr_a = TensorArray(NT+1); arr_a = write(arr_a, 1, a0)
    arr_σ = TensorArray(NT+1); arr_σ = write(arr_σ, 1, σ0)
    arr_ε = TensorArray(NT+1); arr_ε = write(arr_ε, 1, ε0)
    i = constant(1, dtype=Int32)
    tas = [arr_d, arr_v, arr_a, arr_σ, arr_ε]
    _, d, v, a, σ, ε = while_loop(condition, body, [i, tas...])
    d, v, a, σ, ε = stack(d), stack(v), stack(a), stack(σ), stack(ε)
    sp = (NT+1, 2domain.nnodes)
    sp2 = (NT+1, getNGauss(domain), 3)
    set_shape(d, sp), set_shape(v, sp), set_shape(a, sp), set_shape(σ, sp2), set_shape(ε, sp2)
end


text:
Base.@kwdef struct JLArgument
    name::Symbol
    type::Any = Any
    require::Bool = true
    vararg::Bool = false
    default::Maybe{String} = nothing
end

struct JLOption
    name::Symbol
    type::Any
    hint::String
end

struct JLFlag
    name::Symbol
end

struct JLMDOption
    hint::Maybe{String}
    desc::String
    short::Bool
end

struct JLMDFlag
    desc::String
    short::Bool
end

Base.@kwdef struct JLMD
    desc::String = ""
    arguments::Dict{String, String} = Dict{String, String}()
    options::Dict{String, JLMDOption} = Dict{String, JLMDOption}()
    flags::Dict{String, JLMDFlag} = Dict{String, JLMDFlag}()
end


text:
#
# mutex.jl --
#
# POSIX mutexes and condition variables for Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of IPC.jl released under the MIT "expat" license.
# Copyright (C) 2016-2019, Éric Thiébaut (https://github.com/emmt/IPC.jl).
#

"""
```julia
IPC.Mutex()
```

yields an initialized POSIX mutex.  Associated ressources are automatically
destroyed when the returned object is garbage collected.  It is however the
user's responsability to ensure that the object is eventually unlocked.

See also: [`IPC.Condition`](@ref), [`lock`](@ref), [`unlock`](@ref),
          [`trylock`](@ref).

"""
mutable struct Mutex
    handle::Ptr{Cvoid}
    function Mutex()
        buf = Libc.malloc(_sizeof_pthread_mutex_t)
        buf != C_NULL || throw(OutOfMemoryError())
        code = ccall(:pthread_mutex_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                     buf, C_NULL)
        if code != SUCCESS
            # In principle, `pthread_mutex_init` should always return 0.
            Libc.free(buf)
            throw_system_error("pthread_mutex_init", code)
        end
        return finalizer(_destroy, new(buf))
    end
end

"""
```julia
IPC.Condition()
```

yields an initialized POSIX condition variable.  Associated ressources are
automatically destroyed when the returned object is garbage collected.

See also: [`IPC.Mutex`](@ref).

"""
mutable struct Condition
    handle::Ptr{Cvoid}
    function Condition()
        buf = Libc.malloc(_sizeof_pthread_cond_t)
        buf != C_NULL || throw(OutOfMemoryError())
        code = ccall(:pthread_cond_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                     buf, C_NULL)
        if code != SUCCESS
            Libc.free(buf)
            throw_system_error("pthread_cond_init", code)
        end
        return finalizer(_destroy, new(buf))
    end
end

function _destroy(obj::Mutex)
    if (ptr = obj.handle) != C_NULL
        obj.handle = C_NULL # to not free twice
        ccall(:pthread_mutex_destroy, Cint, (Ptr{Cvoid},), ptr)
        Libc.free(ptr)
    end
end

function _destroy(obj::Condition)
    if (ptr = obj.handle) != C_NULL
        obj.handle = C_NULL # to not free twice
        ccall(:pthread_cond_destroy, Cint, (Ptr{Cvoid},), ptr)
        Libc.free(ptr)
    end
end

Base.pointer(obj::Mutex) = obj.handle
Base.pointer(obj::Condition) = obj.handle

function Base.lock(mutex::Mutex)
    code = ccall(:pthread_mutex_lock, Cint, (Mutex,), mutex)
    code == 0 || throw_system_error("pthread_mutex_lock", code)
    nothing
end

function Base.unlock(mutex::Mutex)
    code = ccall(:pthread_mutex_unlock, Cint, (Mutex,), mutex)
    code == 0 || throw_system_error("pthread_mutex_unlock", code)
    nothing
end

function Base.trylock(mutex::Mutex)
    code = ccall(:pthread_mutex_trylock, Cint, (Mutex,), mutex)
    return (code == 0 ? true :
            code == Libc.EBUSY ? false :
            throw_system_error("pthread_mutex_trylock", code))
end

function signal(cond::Condition)
    code = ccall(:pthread_cond_signal, Cint, (Condition,), cond)
    code == 0 || throw_system_error("pthread_cond_signal", code)
    nothing
end

function Base.broadcast(cond::Condition)
    code = ccall(:pthread_cond_broadcast, Cint, (Condition,), cond)
    code == 0 || throw_system_error("pthread_cond_broadcast", code)
    nothing
end

function Base.wait(cond::Condition, mutex::Mutex)
    code = ccall(:pthread_cond_wait, Cint, (Condition, Mutex), cond, mutex)
    code == 0 || throw_system_error("pthread_cond_wait", code)
    nothing
end

"""

```julia
timedwait(cond, mutex, sec)        -> bool
timedwait(cond, mutex, sec, nsec)  -> bool
timedwait(cond, mutex, abstime)    -> bool
```

waits for condition variable `cond` to be signaled using lock `mutex` for
synchronization but no longer than a time limit.

A timeout relative to the current time can be specified either by the number of
seconds `sec` (can be fractional) to wait or by the numbers of seconds and
nanoseconds `sec` and `sec` (both integers).  Otherwise, an absolute timeout
can be specified by `abstime` (of type [`IPC.TimeSpec`](@ref)).

The returned value is `true` if the condition is signalled before expiration of
the timeout, and `false` otherwise.

"""
Base.timedwait(cond::Condition, mutex::Mutex, sec::Integer, nsec::Integer = 0) =
    timedwait(cond, mutex, Int(sec), Int(nsec))

Base.timedwait(cond::Condition, mutex::Mutex, timeout::AbstractFloat) =
    timedwait(cond, mutex, Float64(timeout))

function Base.timedwait(cond::Condition, mutex::Mutex, sec::Int, nsec::Int)
    now = gettimeofday()
    sec += now.tv_sec
    nsec += 1_000*now.tv_usec
    inc, nsec = divrem(nsec, 1_000_000_000)
    sec += inc
    timedwait(cond, mutex, TimeSpec(sec, nsec))
end

function Base.timedwait(cond::Condition, mutex::Mutex, timeout::Float64)
    local sec::Int, nsec::Int
    timeout = max(timeout, 0.0)
    timedwait(cond, mutex, trunc(Int, timeout),
              round(Int, (timeout - trunc(timout))*1E9))
end

function Base.timedwait(cond::Condition, mutex::Mutex, abstime::TimeSpec)
    code = ccall(:pthread_cond_timedwait, Cint,
                 (Condition, Mutex, Ptr{TimeSpec}),
                 cond, mutex, Ref(abstime))
    code == 0 && return true
    code == Libc.ETIMEDOUT && return false
    throw_system_error("pthread_cond_timedwait", code)
end


text:
function HHop(sbasis,facequad,facescale)
    sdim = dimension(sbasis)
    vdim = sdim + 1
    NHF = number_of_basis_functions(sbasis)
    HH = zeros(vdim*NHF,vdim*NHF)
    update_mass_matrix!(HH,sbasis,facequad,facescale,vdim)
    return HH
end

function HHop!(HH,sbasis,squad,components,scale,nhdofs)

    NQ = length(squad)
    @assert length(components) == nhdofs
    projector = components*components'

    for (idx,(p,w)) in enumerate(squad)
        svals = sbasis(p)
        NP = make_row_matrix(svals,projector)
        NI = interpolation_matrix(svals,nhdofs)

        HH .+= scale*NP'*NI*w
    end
end

function HHop(sbasis,squad,components,scale)
    facedim = dimension(sbasis)
    dim = facedim + 1
    NHF = number_of_basis_functions(sbasis)
    HH = zeros(dim*NHF,dim*NHF)

    HHop!(HH,sbasis,squad,components,scale,dim)
    return HH
end

function HHop!(HH,sbasis,iquad,imap::InterpolatingPolynomial,scale,nhdofs)

    NQ = length(iquad)
    @assert length(scale) == NQ

    for (idx,(p,w)) in enumerate(iquad)
        vals = sbasis(p)
        detjac = determinant_jacobian(imap,p)
        N = interpolation_matrix(vals,nhdofs)
        HH .+= N'*N*detjac*scale[idx]*w
    end
end

function HHop(sbasis,iquad,imap::InterpolatingPolynomial,facescale)
    facedim = dimension(sbasis)
    dim = facedim + 1
    NHF = number_of_basis_functions(sbasis)
    HH = zeros(dim*NHF,dim*NHF)

    HHop!(HH,sbasis,iquad,imap,facescale,dim)
    return HH
end

function HHop(sbasis,iquad,imap,inormals,cellmap)

    facescale = scale_area(cellmap,inormals)
    return HHop(sbasis,iquad,imap,facescale)
end

function HHop!(HH,sbasis,iquad,imap,components,scale,nhdofs)

    NQ = length(iquad)
    @assert length(scale) == NQ
    @assert size(components) == (nhdofs,NQ)

    for (idx,(p,w)) in enumerate(iquad)
        t = components[:,idx]
        projector = t*t'
        vals = sbasis(p)
        NP = make_row_matrix(vals,projector)
        NI = interpolation_matrix(vals,nhdofs)
        detjac = determinant_jacobian(imap,p)
        HH .+= NP'*NI*detjac*scale[idx]*w
    end
    return HH
end

function HHop(sbasis,iquad,imap,inormals,components,cellmap)

    facedim = dimension(sbasis)
    dim = facedim + 1
    NHF = number_of_basis_functions(sbasis)
    HH = zeros(dim*NHF,dim*NHF)

    facescale = scale_area(cellmap,inormals)
    HHop!(HH,sbasis,iquad,imap,components,facescale,dim)
    return HH
end


text:
using MyInterpolations: ValsVector

@testset "Testing lin_interp.jl" begin

    @testset "MyLinInterp" begin
        grid = 0:2
        vals = [1, 3, 2.5]
        itp = MyLinInterp(grid, vals)

        x, y = 1.2, 2.9
        @test isapprox(@inferred(itp(x)), y)

        # Grid points
        for (i, x) in enumerate(grid)
            @test itp(x) == vals[i]
        end

        # Dot overload
        xys = [
            ([0.5, 1.2], [2, 2.9]),  # Vector
            (grid, vals),  # UnitRange
            (0:2:2, vals[[1, 3]]),  # StepRange
            (0:0.5:1.5, [1, 2, 3, 2.75]),  # FloatRange
            (linspace(0, 1.5, 4), [1, 2, 3, 2.75])  # LinSpace
        ]

        for (x, y) in xys
            @test isapprox(
                @eval(@inferred($(expand(:($(itp).($x)))))),
                y
            )
        end
    end

    @testset "MyLinInterp extrapolation" begin
        # Current implementation returns vals[1] if x < lower bound, and
        # vals[end] if x > upper bound
        grid = 0:2
        vals = [1, 3, 2.5]
        itp = MyLinInterp(grid, vals)

        @test itp(-1) == vals[1]
        @test itp(10) == vals[end]
    end

    @testset "ValsVector" begin
        vals = [1, 3, 2.5]
        n = length(vals)
        vv = ValsVector(vals)

        @test length(vv) == n
        @test @inferred(vv(1, 0)) == vals[1]
        @test vv(n, 0) == vals[end]
        @test vv(n, 0.5) == vals[end]
        @test vv(1, 0.5) == 2

        @test_throws DomainError vv(-1, 0.5)
        @test_throws DomainError vv(n+1, 0.5)
    end

end


text:
const _EPOUT_BETA0_CACHE_ = Dict()

"""
    Make MaxEnt-EP but applaying MaxEnt starting from a EP solution at beta zero.
    It is important to pass a epout solution computed with beta zero.
    It is faster than the fully integrated version 'maxent_ep' but less accurate.
    It is used for secundary computations as estimating a beta which lead to a
    given objective value and so use this beta for the fully integrated version of
    the algorithm
"""
function fast_maxent_ep(model::MetNet, epout_β0::EPOut,
        alpha::Real,  beta_vec::AbstractVector)

    #= Anna's code scales fluxes to [-1,1], but it
    doesn't rescale a,d back to the original range =#
    maxflux = max(maximum(abs.(model.lb)), maximum(abs.(model.ub)))
    a = epout_β0.sol.a * maxflux;
    d = epout_β0.sol.b * maxflux^2;
    #= a, d are the mean and variance of the univariate Gaussians
    used as site approximations in the EP =#

    #= v, Σ are the mean vector and covariance matrix of Q,
    the full multivariate Gaussian (in Braunstein et al paper) =#
    if haskey(_EPOUT_BETA0_CACHE_, epout_β0)
        Σ = _EPOUT_BETA0_CACHE_[epout_β0]["Σ"]
        v = _EPOUT_BETA0_CACHE_[epout_β0]["v"]
    else
        S = model.S
        b = model.b
        D = Diagonal(1.0 ./ d)
        invΣ = Matrix(alpha*S'*S + D)
        Σ = inv(invΣ)
        v = Σ*(alpha*S'*b + D*a)

        empty!(_EPOUT_BETA0_CACHE_)
        _EPOUT_BETA0_CACHE_[epout_β0] = Dict()
        _EPOUT_BETA0_CACHE_[epout_β0]["Σ"] = Σ
        _EPOUT_BETA0_CACHE_[epout_β0]["v"] = v
    end

    #= From Cossios papar (see README) =#
    w = v + Σ * beta_vec
    wn = (d .* w - diag(Σ) .* a) ./ (d .- diag(Σ)) # mean of the non-truncated marginals
    Σnn = d .* diag(Σ) ./ (d .- diag(Σ)) # variances of the non-truncated marginals

    # Fix case of blocked reactoins
    tns = Truncated.(Normal.(wn, sqrt.(Σnn)), model.lb, model.ub)

    return EPOut(wn, Σnn, mean.(tns), var.(tns), epout_β0.sol, epout_β0.status, epout_β0.iter)
end


text:
const POSITIVES = 1:∞

struct Skip{Synthesizer}
    synthesizer::Synthesizer
    time::FLOAT_SECONDS
end

function make_series(skip::Skip, sample_rate)
    make_series(skip.synthesizer, sample_rate)[
        (round(Int, skip.time * sample_rate)+1):end
    ]
end

"""
    Map(a_function, synthesizers...)

Map `a_function` over `synthesizers`.
Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> first(make_series(Map(sin, Cycles(440Hz)), 44100Hz))
0.06264832417874369
```
"""
struct Map{AFunction,Synthesizers}
    a_function::AFunction
    synthesizers::Synthesizers
    function Map(a_function::AFunction, synthesizers...) where {AFunction}
        new{AFunction,typeof(synthesizers)}(a_function, synthesizers)
    end
end
export Map

"""
    make_series(synthesizer, sample_rate)

Return an iterator that will the play the `synthesizer` at `sample_rate` (with frequency units, like `Hz`).
The iterator should yield `Float64`s between -1 and 1.
Assumes that iterators will never end while they are scheduled.
"""
function make_series(a_map::Map, sample_rate)
    broadcast(
        a_map.a_function,
        map(let sample_rate = sample_rate
            function (synthesizer)
                make_series(synthesizer, sample_rate)
            end
        end, a_map.synthesizers)...,
    )
end
export make_series

"""
    Line(start_level, duration, end_level)

A line from `start_level` to `end_level` with a `duration` in time units like `s`.
Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> first(make_series(Line(0, 1s, 1), 44100Hz))
2.2675736961451248e-5
```
"""
struct Line
    start_level::Float64
    slope::FLOAT_PER_SECOND
end
export Line

function make_series(line::Line, sample_rate)
    line.start_level .+ (line.slope / sample_rate) .* POSITIVES
end

precompile(make_series, (Line, FLOAT_HERTZ))

function Line(start_level, duration, end_level)
    Line(start_level, (end_level - start_level) / duration)
end

precompile(Line, (Float64, FLOAT_SECONDS, Float64))

const τ = 2 * π

# TODO: No need to cycle here; sin is periodic?
"""
    Cycles(frequency)

Cycles from 0 to 2π to repeat at a `frequency` (with frequency units, like `Hz`).
Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz

julia> first(make_series(Cycles(440Hz), 44100Hz))
0.06268937721449021
```
"""
struct Cycles
    frequency::FLOAT_HERTZ
end
export Cycles

function make_series(cycles::Cycles, sample_rate)
    ((cycles.frequency / sample_rate * τ) .* POSITIVES) .% τ
end

precompile(make_series, (Cycles, FLOAT_HERTZ))

"""
    Grow(start_level, duration, end_level)

Exponentially grow or decay from `start_level` to `end_level` over a `duration` in time units like `s`.
Supports [`make_series`](@ref).

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz, s

julia> first(make_series(Grow(0.1, 1s, 1), 44100Hz))
0.10000522141770128
```
"""
struct Grow
    start_level::Float64
    rate::FLOAT_PER_SECOND
end
export Grow

function make_series(grow::Grow, sample_rate)
    grow.start_level .* exp(grow.rate / sample_rate) .^ POSITIVES
end

precompile(make_series, (Grow, FLOAT_HERTZ))

function Grow(start_level, duration, end_level)
    Grow(start_level, log(end_level / start_level) / duration)
end

precompile(Grow, (Float64, FLOAT_SECONDS, Float64))


text:
using SwitchTimeOpt

# Plotting
using PyPlot, PyCall
close("all")
# Import Seaborn for nice plots
@pyimport seaborn as sns
sns.set_context("paper", font_scale=1.5)  # Set default style seaborn
sns.set_style("whitegrid")
# Use Latex Labels in Plots
plt[:rc]("text", usetex=true)
plt[:rc]("font", family="serif")   # Use Serif math

using Ipopt
# solver = IpoptSolver(
#           print_level = 0,
#           linear_solver="ma57")
solver = IpoptSolver(print_level = 0)

### Define system parameters
# Size of the state space
nx = 2;

# Number of switches
N = 5;

# Initial state
x0 = [1.0; 1.0]

# Define initial and final time
t0 = 0.0
tf = 1.0

### Define system dynamics
A = zeros(nx, nx, N+1)

A[:, :, 1] = [-1 0;
              1  2]
A[:, :, 2] = [1 1;
              1 -2]
for i = 3:N+1
  A[:, :, i] = A[:, :, mod(i+1, 2)+1]
end

### Define and solve switching time optimization problem
# m = stoproblem(x0, A, solver=solver)
#
m = stoproblem(x0, A)

# Solve problem
solve!(m)

# Obtain values
tauopt = gettau(m)
Jopt = getobjval(m)
soltime = getsoltime(m)

# Simulate linear system
xsim, xpts, Jsim, t = simulate(m)

### Print results
@printf("OBTAINED RESULTS\n")
@printf("----------------\n")
@printf("Objective Function at the Optimum:              J = %.3f\n", Jopt)
@printf("Simulated Objective Function at the Optimum: Jsim = %.3f\n", Jsim)
@printf("Elapsed Time:                                time = %.2f [ms]\n", mean(soltime)*10^3)
@printf("Optimum:                                      tau = "); show(round.(tauopt,3)); @printf("\n")

### Plot results
figure()
subplot(2,1,1)
for i = 1:m.STOev.ngrid
  axvline(x=m.STOev.tgrid[i], color="lightgrey", linestyle="--")
end
plot(t, xsim[1,:])
ylabel(L"x_1")
yticks([1; 1.5; 2; 2.5])
ylim(0.5, 3)

subplot(2,1,2)
for i = 1:m.STOev.ngrid
  axvline(x=m.STOev.tgrid[i], color="lightgrey", linestyle="--")
end
plot(t, xsim[2,:])
yticks([1; 1.5; 2; 2.5])
ylim(0.5, 3)
ylabel(L"x_2")
xlabel(L"$\mathrm{Time}\; [s]$")

# Save Plot
# tight_layout()
# savefig("linearExample.pdf")

nothing  # Don't print anything


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StorageManagementCost <: OperationalCost
        variable::VariableCost
        fixed::Float64
        start_up::Float64
        shut_down::Float64
        energy_shortage_cost::Float64
        energy_surplus_cost::Float64
    end

Data Structure for Operational Cost Data like variable cost and start - stop costs and energy storage cost.

# Arguments
- `variable::VariableCost`: variable cost
- `fixed::Float64`: fixed cost
- `start_up::Float64`: start-up cost, validation range: `(0, nothing)`, action if invalid: `warn`
- `shut_down::Float64`: shut-down cost, validation range: `(0, nothing)`, action if invalid: `warn`
- `energy_shortage_cost::Float64`: Cost incurred by the model for being short of the energy target., validation range: `(0, nothing)`, action if invalid: `warn`
- `energy_surplus_cost::Float64`: Cost incurred by the model for surplus energy stored., validation range: `(0, nothing)`, action if invalid: `warn`
"""
mutable struct StorageManagementCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "fixed cost"
    fixed::Float64
    "start-up cost"
    start_up::Float64
    "shut-down cost"
    shut_down::Float64
    "Cost incurred by the model for being short of the energy target."
    energy_shortage_cost::Float64
    "Cost incurred by the model for surplus energy stored."
    energy_surplus_cost::Float64
end

function StorageManagementCost(; variable, fixed, start_up, shut_down, energy_shortage_cost, energy_surplus_cost, )
    StorageManagementCost(variable, fixed, start_up, shut_down, energy_shortage_cost, energy_surplus_cost, )
end

# Constructor for demo purposes; non-functional.
function StorageManagementCost(::Nothing)
    StorageManagementCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
        start_up=0.0,
        shut_down=0.0,
        energy_shortage_cost=0.0,
        energy_surplus_cost=0.0,
    )
end

"""Get [`StorageManagementCost`](@ref) `variable`."""
get_variable(value::StorageManagementCost) = value.variable
"""Get [`StorageManagementCost`](@ref) `fixed`."""
get_fixed(value::StorageManagementCost) = value.fixed
"""Get [`StorageManagementCost`](@ref) `start_up`."""
get_start_up(value::StorageManagementCost) = value.start_up
"""Get [`StorageManagementCost`](@ref) `shut_down`."""
get_shut_down(value::StorageManagementCost) = value.shut_down
"""Get [`StorageManagementCost`](@ref) `energy_shortage_cost`."""
get_energy_shortage_cost(value::StorageManagementCost) = value.energy_shortage_cost
"""Get [`StorageManagementCost`](@ref) `energy_surplus_cost`."""
get_energy_surplus_cost(value::StorageManagementCost) = value.energy_surplus_cost

"""Set [`StorageManagementCost`](@ref) `variable`."""
set_variable!(value::StorageManagementCost, val) = value.variable = val
"""Set [`StorageManagementCost`](@ref) `fixed`."""
set_fixed!(value::StorageManagementCost, val) = value.fixed = val
"""Set [`StorageManagementCost`](@ref) `start_up`."""
set_start_up!(value::StorageManagementCost, val) = value.start_up = val
"""Set [`StorageManagementCost`](@ref) `shut_down`."""
set_shut_down!(value::StorageManagementCost, val) = value.shut_down = val
"""Set [`StorageManagementCost`](@ref) `energy_shortage_cost`."""
set_energy_shortage_cost!(value::StorageManagementCost, val) = value.energy_shortage_cost = val
"""Set [`StorageManagementCost`](@ref) `energy_surplus_cost`."""
set_energy_surplus_cost!(value::StorageManagementCost, val) = value.energy_surplus_cost = val