text:
module Lasso
module Util
    # Extract fields from object into function locals
    # See https://github.com/JuliaLang/julia/issues/9755
    macro extractfields(from, fields...)
        esc(Expr(:block, [:($(fields[i]) = $(from).$(fields[i])) for i = 1:length(fields)]...))
    end
    export @extractfields
end

include("FusedLasso.jl")
include("TrendFiltering.jl")

using Reexport, LinearAlgebra, SparseArrays, Random, .Util, MLBase
import Random: Sampler
@reexport using GLM, Distributions, .FusedLassoMod, .TrendFiltering
using GLM: FPVector, LinPred, Link, LmResp, GlmResp, cholpred, updateμ!,
    linkfun, linkinv, LinPredModel
export RegularizationPath, LassoPath, GammaLassoPath, NaiveCoordinateDescent,
       CovarianceCoordinateDescent, fit, fit!, coef, predict,
       minAICc, hasintercept, dof, aicc, distfun, linkfun, cross_validate_path,
       SegSelect, segselect, selectmodel,
       AllSeg, MinAIC, MinAICc, MinBIC, CVSegSelect, MinCVmse, MinCV1se,
       RegularizedModel, LassoModel, GammaLassoModel

## HELPERS FOR SPARSE COEFFICIENTS

struct SparseCoefficients{T} <: AbstractVector{T}
    coef::Vector{T}              # Individual coefficient values
    coef2predictor::Vector{Int}  # Mapping from indices in coef to indices in original X
    predictor2coef::Vector{Int}  # Mapping from indices in original X to indices in coef

    SparseCoefficients{T}(n::Int) where {T} = new(T[], Int[], zeros(Int, n))
end

function LinearAlgebra.mul!(out::Vector, X::Matrix, coef::SparseCoefficients{T}) where T
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = 1:size(X, 1)
            out[i] += c*X[i, ipred]
        end
    end
    out
end

function LinearAlgebra.mul!(out::Vector, X::SparseMatrixCSC, coef::SparseCoefficients{T}) where T
    @extractfields X colptr rowval nzval
    fill!(out, zero(eltype(out)))
    @inbounds for icoef = 1:nnz(coef)
        ipred = coef.coef2predictor[icoef]
        c = coef.coef[icoef]
        @simd for i = colptr[ipred]:colptr[ipred+1]-1
            out[rowval[i]] += c*nzval[i]
        end
    end
    out
end

function LinearAlgebra.dot(x::Vector{T}, coef::SparseCoefficients{T}) where T
    v = 0.0
    @inbounds @simd for icoef = 1:nnz(coef)
        v += x[coef.coef2predictor[icoef]]*coef.coef[icoef]
    end
    v
end

Base.size(x::SparseCoefficients) = (length(x.predictor2coef),)
SparseArrays.nnz(x::SparseCoefficients) = length(x.coef)
Base.getindex(x::SparseCoefficients{T}, ipred::Int) where {T} =
    x.predictor2coef[ipred] == 0 ? zero(T) : x.coef[x.predictor2coef[ipred]]

function Base.setindex!(A::Matrix{T}, coef::SparseCoefficients, rg::UnitRange{Int}, i::Int) where T
    A[:, i] = zero(T)
    for icoef = 1:nnz(coef)
        A[rg[coef.coef2predictor[icoef]], i] = coef.coef[icoef]
    end
    A
end

function Base.copyto!(x::SparseCoefficients, y::SparseCoefficients)
    length(x) == length(y) || throw(DimensionMismatch())
    n = length(y.coef)
    resize!(x.coef, n)
    resize!(x.coef2predictor, n)
    copyto!(x.coef, y.coef)
    copyto!(x.coef2predictor, y.coef2predictor)
    copyto!(x.predictor2coef, y.predictor2coef)
    x
end

# Add a new coefficient to x, returning its index in x.coef
function addcoef!(x::SparseCoefficients{T}, ipred::Int) where T
    push!(x.coef, zero(T))
    push!(x.coef2predictor, ipred)
    coefindex = nnz(x)
    x.predictor2coef[ipred] = coefindex
end

# Add newcoef to column i of coefs
function addcoefs!(coefs::SparseMatrixCSC, newcoef::SparseCoefficients, i::Int)
    n = nnz(coefs)
    nzval = coefs.nzval
    rowval = coefs.rowval
    resize!(nzval, n+nnz(newcoef))
    resize!(rowval, n+nnz(newcoef))
    @inbounds for ipred = 1:length(newcoef.predictor2coef)
        icoef = newcoef.predictor2coef[ipred]
        if icoef != 0
            cval = newcoef.coef[icoef]
            if cval != 0
                n += 1
                nzval[n] = cval
                rowval[n] = ipred
            end
        end
    end
    resize!(nzval, n)
    resize!(rowval, n)
    coefs.colptr[i+1:end] .= n+1
end

## COEFFICIENT ITERATION IN SEQUENTIAL OR RANDOM ORDER
struct RandomCoefficientIterator
    rng::MersenneTwister
    rg::Sampler
    coeforder::Vector{Int}
end
const RANDOMIZE_DEFAULT = true

function RandomCoefficientIterator()
    rng = MersenneTwister(1337)
    RandomCoefficientIterator(rng, Sampler(rng, 1:2), Int[])
end

const CoefficientIterator = Union{UnitRange{Int},RandomCoefficientIterator}

# Iterate over coefficients in random order
function Base.iterate(x::RandomCoefficientIterator)
    if !isempty(x.coeforder)
        @inbounds for i = length(x.coeforder):-1:2
            j = rand(x.rng, x.rg)
            x.coeforder[i], x.coeforder[j] = x.coeforder[j], x.coeforder[i]
        end
        x.coeforder[1], 2
    else
        nothing
    end
end

Base.iterate(x::RandomCoefficientIterator, i) = (i > length(x.coeforder)) ? nothing : (x.coeforder[i], i += 1)

# Add an additional coefficient and return a new CoefficientIterator
function addcoef(x::RandomCoefficientIterator, icoef::Int)
    push!(x.coeforder, icoef)
    RandomCoefficientIterator(x.rng, Sampler(x.rng, 1:length(x.coeforder)), x.coeforder)
end
addcoef(x::UnitRange{Int}, icoef::Int) = 1:length(x)+1

abstract type RegularizationPath{S<:Union{LinearModel,GeneralizedLinearModel},T} <: RegressionModel end

# don't add an intercept when using a @formula because we use the intercept keyword arg to add an intercept
StatsModels.drop_intercept(::Type{R}) where R<:RegularizationPath = true

## LASSO PATH

mutable struct LassoPath{S<:Union{LinearModel,GeneralizedLinearModel},T} <: RegularizationPath{S,T}
    m::S
    nulldev::T                    # null deviance
    nullb0::T                     # intercept of null model, if one was fit
    λ::Vector{T}                  # shrinkage parameters
    autoλ::Bool                   # whether λ is automatically determined
    Xnorm::Vector{T}              # original squared norms of columns of X before standardization
    pct_dev::Vector{T}            # percent deviance explained by each model
    coefs::SparseMatrixCSC{T,Int} # model coefficients
    b0::Vector{T}                 # model intercepts
    niter::Int                    # number of coordinate descent iterations

    LassoPath{S,T}(m, nulldev::T, nullb0::T, λ::Vector{T}, autoλ::Bool, Xnorm::Vector{T}) where {S,T} =
        new(m, nulldev, nullb0, λ, autoλ, Xnorm)
end

function Base.show(io::IO, path::RegularizationPath)
    prefix = isa(path.m, GeneralizedLinearModel) ? string(typeof(distfun(path)).name.name, " ") : ""
    pathsize = size(path)
    println(io, "$(prefix)$(typeof(path).name.name) ($(pathsize[2])) solutions for $(pathsize[1]) predictors in $(path.niter) iterations):")

    if isdefined(path, :coefs)
        coefs = path.coefs
        ncoefs = zeros(Int, size(coefs, 2))
        for i = 1:size(coefs, 2)-1
            ncoefs[i] = coefs.colptr[i+1] - coefs.colptr[i]
        end
        ncoefs[end] = nnz(coefs) - coefs.colptr[size(coefs, 2)] + 1
        show(io, CoefTable(Union{Vector{Int},Vector{Float64}}[path.λ, path.pct_dev, ncoefs], ["λ", "pct_dev", "ncoefs"], []))
    else
        print(io, "    (not fit)")
    end
end

## MODEL CONSTRUCTION

# Controls early stopping criteria with automatic λ
const MIN_DEV_FRAC_DIFF = 1e-5
const MAX_DEV_FRAC = 0.999

# Compute automatic λ values based on λmax and λminratio
function computeλ(λmax, λminratio, α, nλ)
    λmax /= α
    logλmax = log(λmax)
    exp.(range(logλmax, stop=logλmax + log(λminratio), length=nλ))
end

# Compute automatic λ values based on X'y and λminratio without weights ω
function computeλ(Xy::Vector{T}, λminratio, α, nλ, ω::Nothing) where T
    λmax = zero(T)
    for i = eachindex(Xy)
        x = abs(Xy[i])
        if x > λmax
            λmax = x
        end
    end
    computeλ(λmax, λminratio, α, nλ)
end

# Compute automatic λ values based on X'y and λminratio with specific weights ω
function computeλ(Xy::Vector{T}, λminratio, α, nλ, ω::Vector) where T
    λmax = zero(T)
    for i = eachindex(Xy)
        if ω[i] > 0
            x = abs(Xy[i]) / ω[i]
            if x > λmax
                λmax = x
            end
        end
    end
    computeλ(λmax, λminratio, α, nλ)
end

function weightedresiduals(model::LinearModel)
    r = model.rr
    y = r.y
    mu = r.mu
    if isempty(r.wts)
        y - mu
    else
        wts = r.wts
        resid = similar(y)
        @simd for i = eachindex(resid,y,mu,wts)
            @inbounds resid[i] = (y[i] - mu[i]) * wts[i]
        end
        resid
    end
end

# NOTE: mutates nullmodel
maxλ!(nullmodel, λ::Vector, X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T =
    convert(Vector{T}, λ)

function maxλ!(nullmodel::LinearModel, λ::Nothing,
    X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T

    # muscratch = nullmodel.rr.mu.*nullmodel.rr.wts
    Xy = X'*weightedresiduals(nullmodel)

    computeλ(Xy, λminratio, α, nλ, ω)
end

function maxλ!(nullmodel::GeneralizedLinearModel, λ::Nothing,
    X::AbstractMatrix{T}, λminratio, α, nλ, ω) where T

    Xy = X'*broadcast!(*, nullmodel.rr.wrkresid, nullmodel.rr.wrkresid, nullmodel.rr.wrkwt)
    computeλ(Xy, λminratio, α, nλ, ω)
end

# rescales A so that it sums to base
rescale(A, base) = A * (base / sum(A))

# null model's X includes unpenalized (ω_j==0) variables
nullX(X::AbstractMatrix{T}, intercept::Bool, ω::Nothing) where T =
    ones(T, size(X,1), ifelse(intercept, 1, 0))

function nullX(X::AbstractMatrix{T}, intercept::Bool, ω::Vector) where T
    ixunpenalized = findall(iszero, ω)
    [ones(T, size(X,1), ifelse(intercept, 1, 0)) view(X, :, ixunpenalized)]
end

# version for Linear models
function build_model(X::AbstractMatrix{T}, y::FPVector, d::Normal, l::IdentityLink,
                     lp::LinPred, λminratio::Real, λ::Union{Vector,Nothing},
                     wts::Union{FPVector,Nothing}, offset::Vector, α::Real, nλ::Int,
                     ω::Union{Vector,Nothing}, intercept::Bool, irls_tol::Real, dofit::Bool) where T
    mu = isempty(offset) ? copy(y) : y + offset
    nullmodel = LinearModel(LmResp{typeof(y)}(fill!(similar(y), 0), offset, wts, y),
                            GLM.cholpred(nullX(X, intercept, ω), false))
    fit!(nullmodel)
    nulldev = deviance(nullmodel)
    nullb0 = intercept ? coef(nullmodel)[1] : zero(T)

    λ = maxλ!(nullmodel, λ, X, λminratio, α, nλ, ω)

    # mu is just a placeholder here, setting it to nullmodel.mu causes problems
    model = LinearModel(LmResp{typeof(y)}(mu, offset, wts, y), lp)

    (model, nulldev, nullb0, λ)
end

# version for GLMs
function build_model(X::AbstractMatrix{T}, y::FPVector, d::UnivariateDistribution, l::Link,
                     lp::LinPred, λminratio::Real, λ::Union{Vector,Nothing},
                     wts::Union{FPVector,Nothing}, offset::Vector, α::Real, nλ::Int,
                     ω::Union{Vector, Nothing}, intercept::Bool, irls_tol::Real, dofit::Bool) where T
    # Fit to find null deviance
    # Maybe we should reuse this GlmResp object?
    nullmodel = fit(GeneralizedLinearModel, nullX(X, intercept, ω), y, d, l;
                    wts=wts, offset=offset, rtol=irls_tol, dofit=dofit)
    nulldev = deviance(nullmodel)
    nullb0 = intercept ? coef(nullmodel)[1] : zero(T)

    λ = maxλ!(nullmodel, λ, X, λminratio, α, nλ, ω)

    rr = GlmResp(y, d, l, offset, wts)
    model = GeneralizedLinearModel(rr, lp, false)

    (model, nulldev, nullb0, λ)
end

defaultalgorithm(d::Normal, l::IdentityLink, n::Int, p::Int) = p > 5n ? NaiveCoordinateDescent : CovarianceCoordinateDescent
defaultalgorithm(d::UnivariateDistribution, l::Link, n::Int, p::Int) = NaiveCoordinateDescent

# following glmnet rescale penalty factors to sum to the number of coefficients
initpenaltyfactor(penalty_factor::Nothing,p::Int,::Bool) = nothing
initpenaltyfactor(penalty_factor::Vector,p::Int,standardizeω::Bool) =
    standardizeω ? rescale(penalty_factor, p) : penalty_factor

# Standardize predictors if requested
function standardizeX(X::AbstractMatrix{T}, standardize::Bool) where T
    if standardize
        Xnorm = vec(convert(Matrix{T},std(X; dims=1, corrected=false)))
        if any(x -> x == zero(T), Xnorm)
            @warn("""One of the predicators (columns of X) is a constant, so it can not be standardized.
                  To include a constant predicator set standardize = false and intercept = false""")
        end
        for i = 1:length(Xnorm)
            @inbounds Xnorm[i] = 1/Xnorm[i]
        end
        X = X .* transpose(Xnorm)
    else
        Xnorm = T[]
    end

    X, Xnorm
end

@doc raw"""
    fit(LassoPath, X, y, d=Normal(), l=canonicallink(d); ...)

fits a linear or generalized linear Lasso path given the design
matrix `X` and response `y`:

``\underset{\beta}{\operatorname{argmin}} -\frac{1}{N} \mathcal{L}(y|X,\beta) + \lambda\left[(1-\alpha)\frac{1}{2}\|\beta\|_2^2 + \alpha\|\beta\|_1\right]``

The optional argument `d` specifies the conditional distribution of
response, while `l` specifies the link function. Lasso.jl inherits
supported distributions and link functions from GLM.jl. The default
is to fit an linear Lasso path, i.e., `d=Normal(), l=IdentityLink()`,
or ``\mathcal{L}(y|X,\beta) = -\frac{1}{2}\|y - X\beta\|_2^2 + C``

# Examples
```julia
fit(LassoPath, X, y)    # L1-regularized linear regression
fit(LassoPath, X, y, Binomial(), Logit();
    α=0.5) # Binomial logit regression with an Elastic net combination of
           # 0.5 L1 and 0.5 L2 regularization penalties
```
# Arguments
- `wts=ones(length(y))`: Weights for each observation
- `offset=zeros(length(y))`: Offset of each observation
- `λ`: can be used to specify a specific set of λ values at which models are fit.
    If λ is unspecified, Lasso.jl selects nλ logarithmically spaced λ values from
    `λmax`, the smallest λ value yielding a null model, to
    `λminratio * λmax`.
- `nλ=100` number of λ values to use
- `λminratio=1e-4` if more observations than predictors otherwise 0.001.
- `stopearly=true`: When `true`, if the proportion of deviance explained
    exceeds 0.999 or the difference between the deviance explained by successive λ
    values falls below `1e-5`, the path stops early.
- `standardize=true`: Whether to standardize predictors to unit standard deviation
    before fitting.
- `intercept=true`: Whether to fit an (unpenalized) model intercept.
- `algorithm`: Algorithm to use.
    `NaiveCoordinateDescent` iteratively computes the dot product of the
    predictors with the  residuals, as opposed to the
    `CovarianceCoordinateDescent` algorithm, which uses a precomputed Gram matrix.
    `NaiveCoordinateDescent` is typically faster when there are many
    predictors that will not enter the model or when fitting
    generalized linear models.
    By default uses `NaiveCoordinateDescent` if more than 5x as many predictors
    as observations or model is a GLM. `CovarianceCoordinateDescent` otherwise.
- `randomize=true`: Whether to randomize the order in which coefficients are
    updated by coordinate descent. This can drastically speed
    convergence if coefficients are highly correlated.
- `maxncoef=min(size(X, 2), 2*size(X, 1))`: maximum number of coefficients
    allowed in the model. If exceeded, an error will be thrown.
- `dofit=true`: Whether to fit the model upon construction. If `false`, the
    model can be fit later by calling `fit!(model)`.
- `cd_tol=1e-7`: The tolerance for coordinate descent iterations iterations in
    the inner loop.
- `irls_tol=1e-7`: The tolerance for outer iteratively reweighted least squares
    iterations. This is ignored unless the model is a generalized linear model.
- `criterion=:coef` Convergence criterion. Controls how ``cd_tol`` and ``irls_tol``
    are to be interpreted. Possible values are:
    - ``:coef``: The model is considered to have converged if the
    the maximum absolute squared difference in coefficients
    between successive iterations drops below the specified
    tolerance. This is the criterion used by glmnet.
    - ``:obj``: The model is considered to have converged if the
    the relative change in the Lasso/Elastic Net objective
    between successive iterations drops below the specified
    tolerance. This is the criterion used by GLM.jl.
- `minStepFac=0.001`: The minimum step fraction for backtracking line search.
- `penalty_factor=ones(size(X, 2))`: Separate penalty factor ``\omega_j``
    for each coefficient ``j``, i.e. instead of ``\lambda`` penalties become
    ``\lambda\omega_j``.
    Note the penalty factors are internally rescaled to sum to
    the number of variables (`glmnet.R` convention).
- `standardizeω=true`: Whether to scale penalty factors to sum to the number of
    variables (glmnet.R convention).
"""
function StatsBase.fit(::Type{LassoPath},
                       X::AbstractMatrix{T}, y::V, d::UnivariateDistribution=Normal(),
                       l::Link=canonicallink(d);
                       wts::Union{FPVector,Nothing}=ones(T, length(y)),
                       offset::V=similar(y, 0),
                       α::Number=one(eltype(y)), nλ::Int=100,
                       λminratio::Number=ifelse(size(X, 1) < size(X, 2), 0.01, 1e-4),
                       λ::Union{Vector,Nothing}=nothing, standardize::Bool=true,
                       intercept::Bool=true,
                       algorithm::Type=defaultalgorithm(d, l, size(X, 1), size(X, 2)),
                       dofit::Bool=true,
                       irls_tol::Real=1e-7, randomize::Bool=RANDOMIZE_DEFAULT,
                       maxncoef::Int=min(size(X, 2), 2*size(X, 1)),
                       penalty_factor::Union{Vector,Nothing}=nothing,
                       standardizeω::Bool=true,
                       fitargs...) where {T<:AbstractFloat,V<:FPVector}
    size(X, 1) == size(y, 1) || DimensionMismatch("number of rows in X and y must match")
    n = length(y)
    length(wts) == n || error("length(wts) = $(length(wts)) should be 0 or $n")

    X, Xnorm = standardizeX(X, standardize)

    # Lasso initialization
    α = convert(T, α)
    λminratio = convert(T, λminratio)
    coefitr = randomize ? RandomCoefficientIterator() : (1:0)

    # penalty_factor (ω) defaults to a vector of ones
    ω = initpenaltyfactor(penalty_factor, size(X, 2), standardizeω)

    cd = algorithm{T,intercept,typeof(X),typeof(coefitr),typeof(ω)}(X, α, maxncoef, 1e-7, coefitr, ω)

    # GLM response initialization
    autoλ = λ == nothing
    model, nulldev, nullb0, λ = build_model(X, y, d, l, cd, λminratio, λ, wts .* T(1/sum(wts)),
                                            Vector{T}(offset), α, nλ, ω, intercept, irls_tol, dofit)

    # Fit path
    path = LassoPath{typeof(model),T}(model, nulldev, nullb0, λ, autoλ, Xnorm)
    if dofit
        fit!(path; irls_tol=irls_tol, fitargs...)
    end
    path
end

StatsBase.nobs(path::RegularizationPath) = length(path.m.rr.y)

dispersion_parameter(path::RegularizationPath) = GLM.dispersion_parameter(distfun(path))

function StatsBase.loglikelihood(path::RegularizationPath{S,T}) where {S<:LinearModel,T}
    n = nobs(path)
    -0.5.*n.*log.(deviance(path))
end

function StatsBase.loglikelihood(path::RegularizationPath{S,T}) where {S,T}
    n = nobs(path)
    -0.5.*n.*deviance(path)
end

"""
    dof(path::RegularizationPath)

Approximates the degrees-of-freedom in each segment of the path as the number of non zero coefficients
plus a dispersion parameter when appropriate.
Note that for GammaLassoPath this may be a crude approximation, as gamlr does this differently.
"""
function StatsBase.dof(path::RegularizationPath)
    nλ = length(path.λ)
    βs = coef(path)
    df = zeros(Int,nλ)
    for s=1:nλ
        df[s] = sum(βs[:,s].!=0)
    end

    if dispersion_parameter(path)
        # add one for dispersion_parameter
        df.+=1
    end

    df
end

function infocrit(d::T,l::F,n,k) where {T,F}
    if d + one(T) > n
        floatmax(F)
    else
        -2l + (k*d + k*d*(d+1)/(n-d-1))
    end
end

function StatsBase.aicc(path::RegularizationPath;k=2)
    dfs = dof(path)
    ls = loglikelihood(path)
    n = nobs(path)
    broadcast((d,l)->infocrit(d,l,n,k), dfs, ls)
end

minAIC(path::RegularizationPath)=argmin(aic(path))
minAICc(path::RegularizationPath;k=2)=argmin(aicc(path;k=k))
minBIC(path::RegularizationPath)=argmin(bic(path))

hasintercept(path::RegularizationPath) = hasintercept(path.m.pp)

"""
size(path) returns (p,nλ) where p is the number of coefficients (including
any intercept) and nλ is the number of path segments.
If model was only initialized but not fit, returns (p,1).
"""
function Base.size(path::RegularizationPath)
  if isdefined(path,:coefs)
    p,nλ = size(path.coefs)
  else
    X = path.m.pp.X
    p = size(X,2)
    nλ = 1
  end

  if hasintercept(path)
    p += 1
  end

  p,nλ
end

"link of underlying GLM"
GLM.linkfun(path::RegularizationPath{M}) where {M<:LinearModel} = IdentityLink()
GLM.linkfun(path::RegularizationPath{GeneralizedLinearModel{GlmResp{V,D,L},L2}}) where {V<:FPVector,D<:UnivariateDistribution,L<:Link,L2<:GLM.LinPred} = L()

"""
predict(path::RegularizationPath, newX::AbstractMatrix; kwargs...)

Predicted values for a selected segment of a regularization path.

# Examples
```julia
predict(path, newX; select=MinBIC())     # predict using BIC minimizing segment
"""
function StatsBase.predict(path::RegularizationPath, newX::AbstractMatrix{T}; offset::FPVector=T[], select=AllSeg()) where {T<:AbstractFloat}
    # add an interecept to newX if the model has one
    if hasintercept(path)
        newX = [ones(eltype(newX),size(newX,1),1) newX]
    end

    # calculate etas for each obs x segment
    eta = newX * coef(path, select)

    # get model
    mm = path.m

    # adjust for any offset
    if length(mm.rr.offset) > 0
        length(offset) == size(newX, 1) ||
            throw(ArgumentError("fit with offset, so `offset` kw arg must be an offset of length `size(newX, 1)`"))
        broadcast!(+, eta, eta, offset)
    else
        length(offset) > 0 && throw(ArgumentError("fit without offset, so value of `offset` kw arg does not make sense"))
    end

    if typeof(mm) <: LinearModel
        eta
    else
        # invert all etas to mus
        μ(η) = linkinv(linkfun(path), η)
        map(μ, eta)
    end
end

"predicted values for data used to estimate `path`"
function StatsBase.predict(path::RegularizationPath; select=AllSeg())
    X = path.m.pp.X
    offset = path.m.rr.offset

    # destandardize X if needed
    if !isempty(path.Xnorm)
        X = X ./ transpose(path.Xnorm)
    end

    predict(path, X; offset=offset, select=select)
end

"distribution of underlying GLM"
distfun(path::RegularizationPath{M}) where {M<:LinearModel} = Normal()
distfun(path::RegularizationPath) = path.m.rr.d

"deviance at each segment of the path for the fitted model and data"
StatsBase.deviance(path::RegularizationPath) = (1 .- path.pct_dev) .* path.nulldev # * nobs(path)

"""
deviance at each segement of the path for (potentially new) data X and y
select=AllSeg() or MinAICc() like in coef()
"""
function StatsBase.deviance(path::RegularizationPath, X::AbstractMatrix{T}, y::V;
                    offset::FPVector=T[], select=AllSeg(),
                    wts::FPVector=ones(T, length(y))) where {T<:AbstractFloat,V<:FPVector}
    μ = predict(path, X; offset=offset, select=select)
    deviance(path, y, μ, wts)
end

"""
deviance at each segment of the path for (potentially new) y and predicted values μ
"""
function StatsBase.deviance(path::RegularizationPath, y::V, μ::AbstractArray{T},
                wts::FPVector=ones(T, length(y))) where {T<:AbstractFloat,V<:FPVector}
    # get model specs from path
    d = distfun(path)

    # rescale weights
    wts = wts .* convert(T, 1/sum(wts))

    # closure for deviance of a single observation
    dev(ys,μs,ws) = ws * devresid(d, ys, μs)

    # deviances of all obs x segment
    devresidv = broadcast(dev,y,μ,wts)

    # deviance is just their sum
    if size(μ,2) > 1
        vec(sum(devresidv,dims=1))
    else
        sum(devresidv)
    end
end

include("coordinate_descent.jl")
include("gammalasso.jl")
include("segselect.jl")
include("cross_validation.jl")
include("deprecated.jl")

end


text:
function MonotoneChain()
    return MonotoneChain(())
end

function MonotoneChain(arg0::jboolean)
    return MonotoneChain((jboolean,), arg0)
end

function MonotoneChain(arg0::jboolean, arg1::jdouble)
    return MonotoneChain((jboolean, jdouble), arg0, arg1)
end

function find_hull_vertices(obj::MonotoneChain, arg0::Collection)
    return jcall(obj, "findHullVertices", Collection, (Collection,), arg0)
end

ext: jl


text:
using JuMP,Cbc,Gadfly

function Ecualidian_distanse(x1::Float64,x2::Float64,y1::Float64,y2::Float64)
    dist=((x1-x2)^2+(y1-y2)^2)^(1/2)
    return dist
end
################ Index ###################
I=20; # set of customer (i<I)
J=3; # set of service facilities (j<J)
K=9; # set of possible locations for establishing defence
################## Parameters ########################
#costumer location
costumer_location=[ 162.312   286.268
                    197.187   842.411
                    868.178   814.67
                    462.966   680.692
                    598.665   684.046
                    27.3373  440.171
                    96.5991  718.215
                    142.16    745.368
                    454.988   731.582
                    608.22    288.027
                    117.664   791.96
                    806.964   101.707
                    220.986   400.692
                    682.917    18.6734
                    999.23    949.455
                    209.858   932.908
                    727.525   181.241
                    117.496   314.118
                    613.169   682.684
                    804.046   134.233];
#potential defence facilities location
pdef_location=[ 882.754  846.182
                843.319  215.314
                610.862  412.394
                460.733  110.689
                177.43   845.614
                819.911  452.415
                813.314  496.604
                336.915  339.528
                96.554  682.244];
#service_location
service_location= [180.963   761.906;353.318   510.514;15.5935  708.127];

# customer demand
a=[35.0 95.0 85.0 30.0 10.0 75.0 10.0 95.0 40.0 50.0 35.0 100.0 10.0 65.0 15.0 55.0 35.0 55.0 60.0 55.0];

#distance between customer and service facility
d=[476.003 294.567 446.644;
   82.1233 366.787 225.85;
   689.238 597.989 859.216;
   293.464 202.443 448.212;
   424.897 300.513 583.568;
   356.531 333.484 268.214;
   95.0059 330.219 81.6312;
   42.1797 315.823 131.932;
   275.698 243.327 440.02;
   638.052 338.342 726.423;
   70.071 367.076 132.085;
   909.802 610.669 997.001;
   363.425 171.967 369.734;
   896.857 592.066 959.514;
   839.485 780.942 1012.81;
   173.426 446.092 297.094;
   797.435 498.448 885.695;
   452.263 306.893 406.973;
   439.407 311.713 598.116;
   884.425 587.149 975.199];
# the cost of stablishing a defence facilitiy in location k
c=[18000.0,19000.0,19000.0,1000.0,5000.0,17000.0,5000.0,2000.0,6000.0]

r=500; #coverage range of each defence facility

#distance between facility j and defence k
f=[706.834 858.766 554.051 708.77 83.7825 709.958 685.75 450.249 116.064;
   626.877 572.052 275.602 414.002 378.456 470.196 460.205 171.771 308.9;
   878.081 963.325 664.683 745.038 212.353 843.988 825.287 488.992 84.9972]

#probability of service facility j defenced buy defence facility located in k th position
p=[0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15
 0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15
 0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15]
# service facilities capacity
q=[220 700 700];

#unit shipment cost of demand per unit distance
CS=2

 # unit outsourcing cost of demand (independent of distance)
CO=500

#Attacekr Budget
B=2

###################### Model ###################################

m1=Model(solver=CbcSolver())

@variable(m1,X[1:K],Bin) #if decided to establish a defence facility in site k X[k]=1...
@variable(m1,Y[1:J,1:K],Bin) # if jth facility covered by kth defence

@objective(m1,Min,sum(c[k]*X[k] for k=1:K))

for j=1:J
    @constraint(m1,sum(Y[j,k=1:K])>=1)  ##2
end

for k=1:K
    @constraint(m1,sum(Y[j=1:J,k])<=1000*X[k]) ##3
end

for j=1:J
    for k=1:K
        @constraint(m1,Y[j,k]*f[j,k]<=r) ##4
    end
end

for j=1:J
   for k=1:K
      if Ecualidian_distanse(service_location[j,1],pdef_location[k,1],service_location[j,2],pdef_location[k,2])<=r
         @constraint(m1,Y[j,k]==1)
      end
   end
end

solve(m1)
getobjectivevalue(m1)

########## level 2 model ##########
m2=Model(solver=CbcSolver());

@variable(m2,S[1:J],Bin)
@constraint(m2,sum(S[j] for j=1:J)<=B)

########## level 3 model ##########
m3=Model(solver=CbcSolver())

@variable(m3,U[1:I,1:J],Bin) # if customer i allocated to service j
#@variable(m3,T[1:J],Bin) ## this variable is equal to S[j] variable in M2 model
#@variable(m3,V[1:J,1:K],Bin) ## this vaiable is equal to Y[i,j] variable in M1 model
T=[1,0,1];
V=[1 0 0 1 0 1 0 1 1;0 0 0 0 0 1 0 1 1;1 0 0 1 1 0 0 0 0;]

@objective(m3,Min,CS*(sum(sum(a[i]*d[i,j]*U[i,j] for i=1:I) for j=1:J))+(sum(a[i] for i=1:I)*(1-(sum(sum(U[i,j] for i=1:I) for j=1:J)))))

for i=1:I
    @constraint(m3,sum(U[i,j] for j=1:J)<=1)
end

for j=1:J
    c=prod(1-p[j,:].*V[j,:])
    @constraint(m3,sum(a[i]*U[i,j] for i=1:I)<=q[j]*(1-T[j])*prod(1-p[1,:].*V[1,:])+q[j]*(1-prod(1-p[1,:].*V[1,:])))
end
@time solve(m3)
getobjectivevalue(m3)

plot(
      layer(x=costumer_location[:,1], y=costumer_location[:,2],Geom.point,color=[colorant"black"],shape=[Shape.diamond]),
      layer(x=service_location[:,1],y=service_location[:,2],Geom.point,color=[colorant"orange"],shape=[Shape.circle]),
      layer(x=pdef_location[:,1],y=pdef_location[:,2],Geom.point,color=[colorant"purple"],shape=[Shape.xcross]),
      Guide.manual_color_key("locations", ["costumer locations", "service locations","defence faciity locations"],
                           ["black", "orange","purple"]),
      Theme(background_color="white", point_size=6pt)
       )

max_stars_repo_name: imanpanahi/MilProject

max_issues_repo_name: imanpanahi/MilProject

max_forks_repo_name: imanpanahi/MilProject


text:
logspace(start, stop, n = 50; base = 10.0) = base.^range(start, stop, length = n)

"""
    Returns a new dict with symbol keys.
    It will share the source data.
"""
function to_symbol_dict(src_dict::Dict)
    dict = Dict()
    for (k, dat) in src_dict
        dict[Symbol(k)] = dat
    end
    return dict
end

"""
    Returns dict with symbol keys.
    It will share the source object data.
"""
function struct_to_dict(obj::T) where T
    dict = Dict()
    for f in fieldnames(T)
        dict[f] = getproperty(obj, f)
    end
    return dict
end

"""
    give the error text as string
"""
function err_str(err; max_len = 10000)
    s = sprint(showerror, err, catch_backtrace())
    return length(s) > max_len ? s[1:max_len] * "\n[...]" : s
end

function get!push!(d; kwargs...)
    for (k, val) in kwargs
        push!(get!(d, k, []), val)
    end
    d
end

"""
    Give the percent [0-1] of zero elements
"""
sparsity(col::AbstractArray{<:Number}) = float(count(iszero, col) / length(col))

"""
    Returns a copy of the given data. If possible the object is
    stored in a more compacted type, e.i: array -> sparse
"""
compressed_copy(dat; sparsity_th = 0.66) =
    dat isa AbstractArray{<:Number} ?
        sparsity(dat) > sparsity_th ? sparse(dat) : deepcopy(dat) :
        deepcopy(dat)

function compressed_copy(dict::Dict; sparsity_th = 0.66)
    new_dict = Dict()
    for (k, dat) in dict
        new_dict[k] = compressed_copy(dat; sparsity_th = sparsity_th)
    end
    return new_dict
end

"""
    Returns a copy of the given data in a less
    compressed format, e.i: sparse -> array
"""
uncompressed_copy(dat) = deepcopy(dat)
uncompressed_copy(dat::AbstractSparseArray) = dat |> collect

function uncompressed_copy(dict::Dict)
    new_dict = Dict()
    for (k, dat) in dict
        new_dict[k] = uncompressed_copy(dat)
    end
    return new_dict
end


text:
export Multiplication

abstract Multiplication{T} <:Operator{T}

immutable ConcreteMultiplication{D<:Space,S<:Space,V,T} <: Multiplication{T}
    f::Fun{D,V}
    space::S

    ConcreteMultiplication(f::Fun{D,V},sp::S)=new(f,sp)
end

function ConcreteMultiplication{D,T}(f::Fun{D,T},sp::Space)
    @assert domainscompatible(space(f),sp)
    ConcreteMultiplication{D,typeof(sp),T,promote_type(T,eltype(sp))}(chop(f,maxabs(f.coefficients)*40*eps(eltype(f))),sp)
end

# We do this in two stages to support Modifier spaces
# without ambiguity errors
function defaultMultiplication(f::Fun,sp::Space)
    csp=space(f)
    if csp==sp
        error("Implement Multiplication(::Fun{$(typeof(space(f)))},::$(typeof(sp)))")
    end
    MultiplicationWrapper(f,Multiplication(f,csp)*Conversion(sp,csp))
end

Multiplication(f::Fun,sp::Space)=defaultMultiplication(f,sp)

Multiplication(f::Fun,sp::UnsetSpace)=ConcreteMultiplication(f,sp)
Multiplication(f::Fun)=Multiplication(f,UnsetSpace())
Multiplication(c::Number)=Multiplication(Fun(c) )

# This covers right multiplication unless otherwise specified.
Multiplication{D,T}(S::Space,f::Fun{D,T}) = Multiplication(f,S)

function Base.convert{S,V,TT,T}(::Type{Operator{T}},C::ConcreteMultiplication{S,V,TT})
    if T==eltype(C)
        C
    else
        ConcreteMultiplication{S,V,TT,T}(C.f,C.space)
    end
end

domainspace{D,S,T,V}(M::ConcreteMultiplication{D,S,T,V})=M.space
domain(T::ConcreteMultiplication)=domain(T.f)

## Default implementation: try converting to space of M.f

# avoid ambiguity
rangespace{F,T}(D::ConcreteMultiplication{F,UnsetSpace,T}) = UnsetSpace()
getindex{F,T}(D::ConcreteMultiplication{F,UnsetSpace,T},k::Integer,j::Integer) =
    error("No range space attached to Multiplication")

##multiplication can always be promoted, range space is allowed to change
promotedomainspace(D::Multiplication,sp::UnsetSpace) = D
promotedomainspace(D::Multiplication,sp::Space) = Multiplication(D.f,sp)

choosedomainspace{D}(M::ConcreteMultiplication{D,UnsetSpace},::UnsetSpace) = space(M.f)
# we assume multiplication maps spaces to themselves
choosedomainspace{D}(M::ConcreteMultiplication{D,UnsetSpace},sp::Space) = sp

Base.diagm(a::Fun) = Multiplication(a)

immutable MultiplicationWrapper{D<:Space,O<:Operator,V,T} <: Multiplication{T}
    f::Fun{D,V}
    op::O
end

MultiplicationWrapper{D<:Space,V}(T::Type,f::Fun{D,V},op::Operator)=MultiplicationWrapper{D,typeof(op),V,T}(f,op)
MultiplicationWrapper{D<:Space,V}(f::Fun{D,V},op::Operator)=MultiplicationWrapper(eltype(op),f,op)

@wrapper MultiplicationWrapper

function Base.convert{TT,S,V,VV,T}(::Type{Operator{TT}},C::MultiplicationWrapper{S,V,VV,T})
    if TT==eltype(C)
        C
    else
        MultiplicationWrapper{S,V,VV,TT}(C.f,C.op)
    end
end

## Multiplication operators allow us to multiply two spaces

hasfasttransform(::)=false
hasfasttransform(f::Fun)=hasfasttransform(space(f))
hasfasttransformtimes(f,g)=spacescompatible(f,g) && hasfasttransform(f) && hasfasttransform(g)

# This should be overriden whenever the multiplication space is different
function .*{T,N,S,V}(f::Fun{S,T},g::Fun{V,N})
    # When the spaces differ we promote and multiply
    if domainscompatible(space(f),space(g))
        m,n = ncoefficients(f),ncoefficients(g)
        # Heuristic division of parameter space between value-space and coefficient-space multiplication.
        if hasfasttransformtimes(f,g) && log10(m)*log10(n)>4
            transformtimes(f,g)
        elseif m≤n
            coefficienttimes(f,g)
        else
            coefficienttimes(g,f)
        end
    else
        sp=union(space(f),space(g))
        Fun(f,sp).*Fun(g,sp)
    end
end

coefficienttimes(f::Fun,g::Fun) = Multiplication(f,space(g))*g

function transformtimes(f::Fun,g::Fun,n)
    @assert spacescompatible(space(f),space(g))
    f2,g2,sp = pad(f,n),pad(g,n),space(f)
    hc = transform(sp,values(f2).*values(g2))
    chop!(Fun(hc,sp),10norm(hc,Inf)*eps(eltype(hc)))
end
transformtimes(f::Fun,g::Fun)=transformtimes(f,g,ncoefficients(f) + ncoefficients(g) - 1)


text:
@recipe function plot(s::DSP.Periodograms.Spectrogram)
    seriestype := :heatmap
    xguide --> "Time"
    yguide --> "Freq"
    s.time, s.freq, log.(s.power)
end

@recipe function plot(s::DSP.Periodograms.Periodogram; normalize=true)
    yguide --> "Power"
    xguide --> "Freq"
    xscale --> :log10
    yscale --> :log10
    if normalize
        y =  s.power[2:end] ./ std(s.power)
    else
        y = s.power[2:end]
    end
    s.freq[2:end], y
end

function Base.show(io::IO, m::MIME"image/png", anim::Plots.Animation)
    path = mktemp()[1]
    path = path*".gif"
    g = gif(anim, path, fps=15)
    show(io, MIME"image/gif"(), g)
end
Base.showable(::MIME"image/gif", agif::Plots.Animation) = true
Base.showable(::MIME"text/html", agif::Plots.Animation) = true

@recipe function plot(v::Vector{<:Complex}, d = determine_domain(v); circle=true)
    @series begin
        primary := true
        seriestype := :scatter
        real.(v), imag.(v)
    end
    xguide --> "Re"
    yguide --> "Im"
    if circle && d isa Discrete
        t = LinRange(0, 2pi, 100)
        @series begin
            seriestype := :path
            marker := false
            primary := false
            linestyle := :dash
            linecolor := :black
            cos.(t), sin.(t)
        end
    end
end

@recipe function plot(r::AbstractRoots)
    title --> "Root map"
    xguide --> "Re"
    yguide --> "Im"
    r.r, domain(r)
end

"""
    assignmentplot(m1, m2)

Plots the poles of two models and the optimal assignment between them
"""
assignmentplot
@userplot AssignmentPlot

@recipe function assignmentplot(h::AssignmentPlot; d=Discrete(),p=2)
    m1,m2 = h.args[1:2]
    rf(a) = roots(d, a)
    r1 = (rf(m1))
    r2 = hungariansort(r1, (rf(m2)), p)
    @series begin
        markerstrokealpha --> 0.1
        label --> "m1"
        seriestype := :scatter
        real.(r1), imag.(r1)
    end
    @series begin
        markerstrokealpha --> 0.1
        label --> "m2"
        seriestype := :scatter
        real.(r2), imag.(r2)
    end
    @series begin
        primary := false
        seriestype := :path
        xc = [real.(r1) real.(r2) fill(Inf, length(m1.p))]'[:]
        yc = [imag.(r1) imag.(r2) fill(Inf, length(m1.p))]'[:]
        xc,yc
    end
    d isa Discrete && @series begin
        line := (:dash, :black)
        primary := false
        cos.(0:0.1:2pi),sin.(0:0.1:2pi)
    end
end

"""
    flowplot(dist,m1, m2)

Plots the poles of two models and the optimal mass flow between them
"""
flowplot
@userplot FlowPlot

@recipe function flowplot(h::FlowPlot)
    dist,m1,m2 = h.args[1:3]

    Γ = transport_plan(dist,m1,m2;d=d)
    e1,e2 = roots.(domain(dist), (m1,m2))
    # display(Γ)

    @series begin
        markerstrokealpha --> 0.1
        label --> "m1"
        seriestype := :scatter
        real.(e1), imag.(e1)
    end
    @series begin
        markerstrokealpha --> 0.1
        label --> "m2"
        seriestype := :scatter
        real.(e2), imag.(e2)
    end
    for i in eachindex(e1), j in eachindex(e2)
        @series begin
            primary := false
            seriestype := :path
            coords = [e1[i], e2[j]]
            linewidth --> 10Γ[i,j]
            linealpha --> sqrt(10Γ[i,j])
            real(coords), imag(coords)
        end
    end
    domain(dist) isa Discrete && @series begin
        l := (:dash, :black)
        primary := false
        cos.(0:0.1:2pi),sin.(0:0.1:2pi)
    end
end

@recipe function plot(m::AbstractModel, w=exp10.(LinRange(-2, log10(pi), 200)); rad=true)
    mag = evalfr.(Continuous(), Identity(), w, m, m.b)
    yscale --> :log10
    xscale --> :log10
    rad || (w = w .* 2π)
    w, sqrt.(mag)
end

@recipe function plot(
    m::TimeVaryingAR,
    w = exp10.(LinRange(-2, log10(pi), 200));
    rad = true,
)
    rad || (w = w .* 2π)
    mag =
        evalfr.(
            Continuous(),
            Identity(),
            w,
            reshape(m.models, 1, :),
            getfield.(m.models, :b)',
        )
    seriestype --> :heatmap
    log.(mag)
end


text:
using Plots, LaTeXStrings, Statistics

function Deposition(;len, tot_time, time_steps, rate)
    Time = ceil.(Int, exp.(0:(tot_time)/(time_steps):tot_time))
    surf = [0 for i=1:len]
    VarList = [0.0 for i=1:time_steps]
    for n in 1:time_steps
        randsurf = rand(1:len,(Time[n+1]-Time[n])*rate)
        for i in randsurf
            surf[i] += 1
        end
        VarList[n] = std(surf)
    end
    return VarList
end

function Linear_fit(;len, tot_time, time_steps, rate)
    A = [hcat(log.(Time)) reshape(ones(time_steps), time_steps, 1)]
    b = reshape(log.(meanVar), time_steps, 1)
    line = (A \ b)
    x = 0:tot_time
    y = x .* line[1] .+ line[2]
    return x, y, line
end

iternum = 1000
Parameters = Dict(
                :len => 200,
                :tot_time => 10,
                :time_steps => 20,
                :rate => 10)
allVar = [[0.0 for i in 1:Parameters[:time_steps]] for j = 1:iternum]
meanVar = [0.0 for i in 1:Parameters[:time_steps]]
vars = [0.0 for i in 1:Parameters[:time_steps]]
for i in 1:iternum
    VarList = Deposition(;Parameters...)
    allVar[i] = VarList
    meanVar += VarList
    print("\r$i")
end
meanVar /= iternum
for i in 1:Parameters[:time_steps]
    vars[i] = std(log.(hcat(allVar...))[i,:])
end

Time = ceil.(Int, exp.(0:(Parameters[:tot_time])/(Parameters[:time_steps]):Parameters[:tot_time]))[1:end-1]
X, Y, Line = Linear_fit(;Parameters...)

scatter(log.(Time),log.(meanVar),
    xlabel= L"Log\ Time",
    ylabel= L"Log\ W_{(t)}",
    title= L"Log-Log\ Plot\ of\ ~W_{(t)}-Time~, Particles \ number \approx e^{%$(Parameters[:tot_time])}",
    yerror = vars,
    label = L"Data\ point")
plot!(X,Y,c= :black,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))")
savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet2\\Figs\\Q2\\RBDW(t).png")


text:
test_doc(obj, doc) = @test Documented(obj).p[:docstring] == doc

@testset "Generated documentation" begin
    @testset "Low-level structs" begin
        test_doc(StructDefinition{false}(struct_by_name(:VkExtent2D)),
            """
            Intermediate wrapper for VkExtent2D.

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkExtent2D.html)
            """
        )

        test_doc(
            Constructor(StructDefinition{false}(struct_by_name(:VkExtent2D))),
            """
            Arguments:
            - `width::UInt32`
            - `height::UInt32`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkExtent2D.html)
            """,
        )

        test_doc(
            Constructor(StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo))),
            """
            Arguments:
            - `enabled_layer_names::Vector{String}`
            - `enabled_extension_names::Vector{String}`
            - `next::Ptr{Cvoid}`: defaults to `C_NULL`
            - `flags::UInt32`: defaults to `0`
            - `application_info::_ApplicationInfo`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkInstanceCreateInfo.html)
            """,
        )

        test_doc(
            Constructor(StructDefinition{false}(struct_by_name(:VkSubmitInfo2))),
            """
            Arguments:
            - `wait_semaphore_infos::Vector{_SemaphoreSubmitInfo}`
            - `command_buffer_infos::Vector{_CommandBufferSubmitInfo}`
            - `signal_semaphore_infos::Vector{_SemaphoreSubmitInfo}`
            - `next::Ptr{Cvoid}`: defaults to `C_NULL`
            - `flags::SubmitFlag`: defaults to `0`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkSubmitInfo2.html)
            """,
        )
    end

    @testset "High-level structs" begin
        test_doc(StructDefinition{true}(struct_by_name(:VkInstanceCreateInfo)),
            """
            High-level wrapper for VkInstanceCreateInfo.

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkInstanceCreateInfo.html)
            """,
        )

        test_doc(
            Constructor(StructDefinition{true}(struct_by_name(:VkInstanceCreateInfo))),
            """
            Arguments:
            - `enabled_layer_names::Vector{String}`
            - `enabled_extension_names::Vector{String}`
            - `next::Any`: defaults to `C_NULL`
            - `flags::UInt32`: defaults to `0`
            - `application_info::ApplicationInfo`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/VkInstanceCreateInfo.html)
            """,
        )

        test_doc(
            APIFunction(create_func(:vkCreateInstance), false),
            """
            Return codes:
            - `ERROR_OUT_OF_HOST_MEMORY`
            - `ERROR_OUT_OF_DEVICE_MEMORY`
            - `ERROR_INITIALIZATION_FAILED`
            - `ERROR_LAYER_NOT_PRESENT`
            - `ERROR_EXTENSION_NOT_PRESENT`
            - `ERROR_INCOMPATIBLE_DRIVER`

            Arguments:
            - `enabled_layer_names::Vector{String}`
            - `enabled_extension_names::Vector{String}`
            - `allocator::_AllocationCallbacks`: defaults to `C_NULL`
            - `next::Ptr{Cvoid}`: defaults to `C_NULL`
            - `flags::UInt32`: defaults to `0`
            - `application_info::_ApplicationInfo`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkCreateInstance.html)
            """,
        )
    end

    @testset "API functions" begin
        @testset "Intermediate functions" begin
            test_doc(
                APIFunction(func_by_name(:vkEnumerateInstanceExtensionProperties), false),
                """
                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`
                - `ERROR_LAYER_NOT_PRESENT`

                Arguments:
                - `layer_name::String`: defaults to `C_NULL`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkEnumerateInstanceExtensionProperties.html)
                """,
            )

            test_doc(
                APIFunction(func_by_name(:vkDestroyDevice), false),
                """
                Arguments:
                - `device::Device` (externsync)
                - `allocator::_AllocationCallbacks`: defaults to `C_NULL`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkDestroyDevice.html)
                """,
            )

            test_doc(
                APIFunction(func_by_name(:vkGetPipelineCacheData), false),
                """
                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`

                Arguments:
                - `device::Device`
                - `pipeline_cache::PipelineCache`

                !!! warning
                    The pointer returned by this function holds memory owned by Julia. It is therefore **your** responsibility to free it after use (e.g. with `Libc.free`).

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkGetPipelineCacheData.html)
                """,
            )

            test_doc(
                APIFunction(func_by_name(:vkWriteAccelerationStructuresPropertiesKHR), false),
                """
                Extension: VK\\\\_KHR\\\\_acceleration\\\\_structure

                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`

                Arguments:
                - `device::Device`
                - `acceleration_structures::Vector{AccelerationStructureKHR}`
                - `query_type::QueryType`
                - `data_size::UInt`
                - `data::Ptr{Cvoid}` (must be a valid pointer with `data_size` bytes)
                - `stride::UInt`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkWriteAccelerationStructuresPropertiesKHR.html)
                """,
            )

            test_doc(
                APIFunction(func_by_name(:vkEnumerateInstanceLayerProperties), false),
                """
                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkEnumerateInstanceLayerProperties.html)
                """,
            )
        end

        @testset "High-level functions" begin
            test_doc(
                VulkanGen.promote_hl(APIFunction(func_by_name(:vkWriteAccelerationStructuresPropertiesKHR), false)),
                """
                Extension: VK\\\\_KHR\\\\_acceleration\\\\_structure

                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`

                Arguments:
                - `device::Device`
                - `acceleration_structures::Vector{AccelerationStructureKHR}`
                - `query_type::QueryType`
                - `data_size::UInt`
                - `data::Ptr{Cvoid}` (must be a valid pointer with `data_size` bytes)
                - `stride::UInt`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkWriteAccelerationStructuresPropertiesKHR.html)
                """,
            )

            test_doc(
                VulkanGen.promote_hl(APIFunction(func_by_name(:vkEnumerateInstanceLayerProperties), false)),
                """
                Return codes:
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkEnumerateInstanceLayerProperties.html)
                """,
            )

            test_doc(
                VulkanGen.promote_hl(APIFunction(func_by_name(:vkAcquireNextImageKHR), false)),
                """
                Extension: VK\\\\_KHR\\\\_swapchain

                Return codes:
                - `SUCCESS`
                - `TIMEOUT`
                - `NOT_READY`
                - `SUBOPTIMAL_KHR`
                - `ERROR_OUT_OF_HOST_MEMORY`
                - `ERROR_OUT_OF_DEVICE_MEMORY`
                - `ERROR_DEVICE_LOST`
                - `ERROR_OUT_OF_DATE_KHR`
                - `ERROR_SURFACE_LOST_KHR`
                - `ERROR_FULL_SCREEN_EXCLUSIVE_MODE_LOST_EXT`

                Arguments:
                - `device::Device`
                - `swapchain::SwapchainKHR` (externsync)
                - `timeout::UInt64`
                - `semaphore::Semaphore`: defaults to `C_NULL` (externsync)
                - `fence::Fence`: defaults to `C_NULL` (externsync)

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkAcquireNextImageKHR.html)
                """,
            )

            test_doc(
                VulkanGen.promote_hl(APIFunction(func_by_name(:vkGetDescriptorSetLayoutSupport), false)),
                """
                Arguments:
                - `device::Device`
                - `create_info::DescriptorSetLayoutCreateInfo`
                - `next_types::Type...`: types of members to initialize and include as part of the `next` chain

                [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkGetDescriptorSetLayoutSupport.html)
                """
            )
        end
    end
    @testset "Handles" begin
        test_doc(
            Constructor(HandleDefinition(handle_by_name(:VkInstance)), VulkanGen.promote_hl(APIFunction(create_func(:vkCreateInstance), false))),
            """
            Arguments:
            - `enabled_layer_names::Vector{String}`
            - `enabled_extension_names::Vector{String}`
            - `allocator::AllocationCallbacks`: defaults to `C_NULL`
            - `next::Any`: defaults to `C_NULL`
            - `flags::UInt32`: defaults to `0`
            - `application_info::ApplicationInfo`: defaults to `C_NULL`

            [API documentation](https://www.khronos.org/registry/vulkan/specs/1.3-extensions/man/html/vkCreateInstance.html)
            """
        )
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
function markovmatrix(points::Array{Float64, 2},
                      image_points::Array{Float64, 2},
                      simplex_inds::Array{Int, 2})

    const n_simplices = size(simplex_inds, 1)
    const dim = size(points, 2)

    # Mismatch for the markovity of the Markov matrix is at most delta if volume_tolerance tolerance
    # is defined as follows
    const delta::Float64 = 1/10^5
    const volume_tolerance::Float64 = delta/n_simplices

    # Tolerance for similary of convex expansion coefficients of simplex vertices in simplexintersection function.
    const tolerance::Float64 = 1/10^12

    P = zeros(Float64, n_simplices, n_simplices)
    #println("There are ", n_simplices, " simplices in the triangulation.")
    for i = 1:n_simplices
        image = image_points[simplex_inds[i, :], :].'
        # Volume of the image of the simplex
        image_volume = abs(det([ones(1, dim + 1); image]))
        t1 = time_ns()
        for j = 1:n_simplices
            simplex = points[simplex_inds[j, :], :].'
            simplex_volume = abs(det([ones(1, dim + 1); simplex]))

            if (simplex_volume * image_volume > 0 && simplex_volume/image_volume > volume_tolerance)
                P[i, j] = SimplexIntersection.simplexintersection(simplex, image)/image_volume
            elseif image_volume < tolerance
                P[i, j] = 0
            end
        end
        println("Time for image simplex #", i, "\n", (time_ns() - t1)/10^9, " s")

    end

    return P
end

function markovmatrix(t::Triangulation)

    const n_simplices = size(t.simplex_inds, 1)
    const dim = size(points, 2)

    # Mismatch for the markovity of the Markov matrix is at most delta if volume_tolerance tolerance
    # is defined as follows
    const delta::Float64 = 1/10^5
    const volume_tolerance::Float64 = delta/n_simplices

    # Tolerance for similary of convex expansion coefficients of simplex vertices in simplexintersection function.
    const tolerance::Float64 = 1/10^12

    P = zeros(Float64, n_simplices, n_simplices)
    #println("There are ", n_simplices, " simplices in the triangulation.")
    for i = 1:n_simplices
        # Volume of the image of the simplex
        image_volume = t.volumes_im[i]
        for j = 1:n_simplices
            simplex_volume = t.volumes[j]
            if (simplex_volume * image_volume > 0 && simplex_volume/image_volume > volume_tolerance)
                P[i, j] = Simplices.simplexintersection(view(t.points[t.simplex_inds[j, :], :].', :, :), view(t.image_points[t.simplex_inds[i, :], :].', :, :))/image_volume
            elseif image_volume < tolerance
                P[i, j] = 0
            end
        end
        #println("Time for image simplex #", i, "\n", (time_ns() - t1)/10^9, " s")
    end

    return P
end


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

mutable struct IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    fields::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Fields } # spec name: fields
    manager::Any # spec type: Union{ Nothing, String } # spec name: manager
    operation::Any # spec type: Union{ Nothing, String } # spec name: operation
    time::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: time

    function IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry(;apiVersion=nothing, fields=nothing, manager=nothing, operation=nothing, time=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("fields"), fields)
        setfield!(o, Symbol("fields"), fields)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("manager"), manager)
        setfield!(o, Symbol("manager"), manager)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("operation"), operation)
        setfield!(o, Symbol("operation"), operation)
        validate_property(IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry

const _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("fields")=>Symbol("fields"), Symbol("manager")=>Symbol("manager"), Symbol("operation")=>Symbol("operation"), Symbol("time")=>Symbol("time"))
const _property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("fields")=>"IoK8sApimachineryPkgApisMetaV1Fields", Symbol("manager")=>"String", Symbol("operation")=>"String", Symbol("time")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1ManagedFieldsEntry)
    true
end

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


text:
α = 0.025

@testset "gaussian copula" begin
  @testset "exceptions" begin
    @test_throws DomainError Gaussian_cop([1. 0.; .1 1.])
    @test_throws DomainError Gaussian_cop([1. 2.; .2 1.])
    @test_throws DomainError Gaussian_cop([2. 1.; 1. 1.])
  end
  @testset "examples on data" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Gaussian_cop([1. 0.; 0. 1.]); rng = rng) ≈ [0.44953 0.74757] atol=1.0e-5

    Random.seed!(43)
    x = simulate_copula(350000, Gaussian_cop([1. 0.5; 0.5 1.]))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0
    @test tail(x[:,1], x[:,2], "r", 0.00001) ≈ 0
  end

end
@testset "t-student copula" begin
  @testset "exceptions" begin
    @test_throws DomainError Student_cop([1. 0.; .1 1.], 2)
    @test_throws DomainError Student_cop([1. 2.; .2 1.], 2)
    @test_throws DomainError Student_cop([2. 1.; 1. 1.], 2)
    @test_throws DomainError Student_cop([1. 0.; 0. 1.], 0)
  end

  @testset "example on data" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Student_cop([1. 0.; 0. 1.], 1); rng = rng) ≈ [0.39949 0.832173] atol=1.0e-5

    ν = 10
    rho = 0.5
    λ = 2*pdf(TDist(ν+1), -sqrt.((ν+1)*(1-rho)/(1+rho)))
    Random.seed!(43)
    xt = simulate_copula(350000, Student_cop([1. rho; rho 1.], ν))
    @test pvalue(ExactOneSampleKSTest(xt[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(xt[:,2], Uniform(0,1))) > α
    @test tail(xt[:,1], xt[:,2], "l") ≈ λ atol=1.0e-1
    @test tail(xt[:,1], xt[:,2], "r") ≈ λ atol=1.0e-1

    convertmarg!(xt, Normal)
    @test cor(xt) ≈ [1. rho; rho 1.] atol=1.0e-2
  end
end

@testset "Frechet copula" begin
  @testset "Frechet copula exceptions" begin
    @test_throws DomainError Frechet_cop(3, -0.8)
    @test_throws DomainError Frechet_cop(1, 0.8)

    @test_throws AssertionError Frechet_cop(3, 0.8, 0.1)
    @test_throws DomainError Frechet_cop(2, 0.8, 0.5)
    @test_throws DomainError Frechet_cop(2, -0.2, 0.8)
    @test_throws DomainError Frechet_cop(2, 0.8, -0.5)

    u = zeros(3, 5)
    cop = Frechet_cop(3, 0.8)
    @test_throws AssertionError simulate_copula!(u, cop)
  end
  @testset "small example" begin
    u = [0.3, 0.5, 0.7]
    frechet_el!(u, 0.3, 0.5)
    @test u == [0.3, 0.5, 0.7]
    frechet_el!(u, 0.7, 0.5)
    @test u == [0.7, 0.7, 0.7]
    u = [0.1, 0.8]
    frechet_el2!(u, 0.2, 0.6, 0.6)
    @test u ≈ [0.2, 0.8]
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Frechet_cop(4, 1.); rng = rng) ≈ [0.042730 0.04273 0.04273 0.042730] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Frechet_cop(4, 1.); rng = rng) ≈ [0.042730 0.042730 0.04273 0.04273] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Frechet_cop(2, 0.4, 0.4); rng = rng) ≈ [0.632510 0.36748] atol=1.0e-5

    u = zeros(1,4)
    u1 = zeros(1,2)
    Random.seed!(43)
    rng = StableRNG(123)
    simulate_copula!(u, Frechet_cop(4, 1.); rng = rng)
    @test u ≈ [0.042730 0.042730 0.04273 0.04273] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    simulate_copula!(u1, Frechet_cop(2, 0.4, 0.4); rng = rng)
    @test u1 ≈ [0.63251 0.36748] atol=1.0e-5

  end
  @testset "examples on larger data" begin
    Random.seed!(43)
    x = simulate_copula(350000, Frechet_cop(3, 0.3))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l") ≈ 0.3 atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r") ≈ 0.3 atol=1.0e-1
    @test corspearman(x) ≈ [1. 0.3 0.3; 0.3 1. 0.3; 0.3 0.3 1.] atol=1.0e-2

    x = simulate_copula(350000, Frechet_cop(3, 1.))
    ret = (x[:, 1] .<  0.2).*(x[:, 2] .<  0.3).*(x[:, 3] .<  0.5)
    @test length(findall(ret))/size(x,1) ≈ 0.2 atol=1.0e-3
    ret = (x[:, 1] .<  0.8).*(x[:, 2] .<  0.3).*(x[:, 3] .<  0.5)
    @test length(findall(ret))/size(x,1) ≈ 0.3 atol=1.0e-3
    ret = (x[:, 1] .<  0.8).*(x[:, 2] .<  0.9).*(x[:, 3] .<  0.5)
    @test length(findall(ret))/size(x,1) ≈ 0.5 atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(1500000, Frechet_cop(2, .0))
    @test corspearman(x) ≈ [1. 0.; 0. 1.] atol=1.0e-2

    x = simulate_copula(350000, Frechet_cop(2, 0.8, 0.1));
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test corspearman(x) ≈ [1. 0.7; 0.7 1.] atol=1.0e-2
    @test tail(x[:,1], x[:,2], "l") ≈ 0.8 atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r") ≈ 0.8 atol=1.0e-1
  end
end

@testset "Marshall-Olkin copula" begin
  @testset "helpers and exceptions" begin
    @test_throws AssertionError Marshall_Olkin_cop([1., -2., 3.])
    @test_throws AssertionError Marshall_Olkin_cop([1., 3.])
    @test τ2λ([0.4], [4.2, 1.2]) ≈ [4.2, 1.2, 3.6]
    @test τ2λ([0.5, 0.6, 0.5], [1., 2., 3., 4.]) ≈ [1, 2, 3, 0, 0.5, 0, 4]
    @test moρ2τ(0.6) ≈ 0.5 atol=1.0e-2

    u = zeros(3,3)
    cop = Marshall_Olkin_cop([1.,2.,3.])
    @test_throws AssertionError simulate_copula!(u, cop)
  end
  @testset "small example" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Marshall_Olkin_cop([1., 2., 3.]); rng = rng) ≈ [0.585174 0.511807] atol=1.0e-5

    u = zeros(1,2)
    Random.seed!(43)
    rng = StableRNG(123)
    simulate_copula!(u, Marshall_Olkin_cop([1., 2., 3.]); rng=rng)
    @test u ≈ [0.585174 0.511807] atol=1.0e-5

    m = [0.252982 0.201189;  0.464758 0.409039; 0.585662 0.5357]
    @test mocopula([0.2 0.3 0.4; 0.3 0.4 0.6; 0.4 0.5 0.7], 2, [1., 1.5, 2.]) ≈ m atol=1.0e-4

    s = collect(combinations(1:2))
    mocopula_el([0.1, 0.2, 0.3], 2, [1., 2., 3.], s) ≈ [0.20082988502465082, 0.1344421423967149]
  end
  @testset "example on larger data" begin
    Random.seed!(43)
    x = simulate_copula(100000, Marshall_Olkin_cop([1.1, 0.2, 0.6]))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    a2 = 0.6/0.8
    a1 = 0.6/1.7
    @test corkendall(x)[1,2]≈ a1*a2/(a1+a2-a1*a2) atol=1.0e-3
    @test tail(x[:,1], x[:,2], "r") ≈ a1 atol=1.0e-1
    Random.seed!(42)
    x = simulate_copula(20000, Marshall_Olkin_cop([1.1, 0.2, 2.1, 0.6, 0.5, 3.2, 7.1, 2.1]))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  end
end

@testset "tests on Big Float" begin
  α = 0.025
  if false
    Random.seed!(1234)
    Σ = BigFloat.([1. 0.; 0. 1.])
    x = simulate_copula(1000, Gaussian_cop(Σ))

    println(typeof(x) == Array{BigFloat,2})
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

    Random.seed!(1234)
    Σ = BigFloat.([1. 0.; 0. 1.])
    x = simulate_copula(1000, Student_cop(Σ, 2))

    println(typeof(x) == Array{BigFloat,2})
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  end

  λs = BigFloat.([1., 2., 3.])
  Random.seed!(1234)
  x = simulate_copula(1000, Marshall_Olkin_cop(λs))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

  a = b = BigFloat(0.4)
  x = simulate_copula(1000, Frechet_cop(2, a, b))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/coeffs_calcs.jl
# Defines functions used to compute taylor coefficients using static
# univariate Taylor series.
#############################################################################

@generated function truncuated_STaylor1(x::STaylor1{N,T}, v::Int) where {N,T}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:N])
    syms = Symbol[Symbol("c$i") for i in 1:N]
    for i = 1:N
        sym = syms[i]
        ex_line = :($(sym) = $i <= v ? x[$(i-1)] : zero($T))
        ex_calc.args[i] = ex_line
    end
    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end

@generated function copy_recurse(dx::STaylor1{N,T}, x::STaylor1{N,T}, ord::Int, ford::Float64) where {N,T}
    ex_calc = quote end
    append!(ex_calc.args, Any[nothing for i in 1:(N+1)])
    syms = Symbol[Symbol("c$i") for i in 1:N]
    ex_line = :(nv = dx[ord-1]/ford) #/ord)
    ex_calc.args[1] = ex_line
    for i = 0:(N-1)
        sym = syms[i+1]
        ex_line = :($(sym) = $i < ord ? x[$i] : (($i == ord) ? nv : zero(T))) # nv))
        ex_calc.args[i+2] = ex_line
    end
    exout = :(($(syms[1]),))
    for i = 2:N
        push!(exout.args, syms[i])
    end
    return quote
               Base.@_inline_meta
               $ex_calc
               return STaylor1{N,T}($exout)
            end
end

"""
$(TYPEDSIGNATURES)

A variant of the jetcoeffs! function used in TaylorIntegration.jl
(https://github.com/PerezHz/TaylorIntegration.jl/blob/master/src/explicitode.jl)
which preallocates taux and updates taux coefficients to avoid further allocations.

The TaylorIntegration.jl package is licensed under the MIT "Expat" License:
Copyright (c) 2016-2020: Jorge A. Perez and Luis Benet. Permission is hereby
granted, free of charge, to any person obtaining a copy of this software and
associated documentation files (the "Software"), to deal in the Software without
restriction, including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the following
conditions: The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
"""
function jetcoeffs!(eqsdiff!, t::T, x::Vector{STaylor1{N,U}}, xaux::Vector{STaylor1{N,U}},
                    dx::Vector{STaylor1{N,U}}, order::Int, params,
                    vnxt::Vector{Int}, fnxt::Vector{Float64}) where {N, T<:Number, U<:Number}

      ttaylor = STaylor1(t, Val{N-1}())
      for ord = 1:order

          fill!(vnxt, ord)
          fill!(fnxt, Float64(ord))

          map!(truncuated_STaylor1, xaux, x, vnxt)
          eqsdiff!(dx, xaux, params, ttaylor)::Nothing
          x .= copy_recurse.(dx, x, vnxt, fnxt)
      end
      nothing
end


text:
"""
https://github.com/owensgroup/RXMesh

https://docs.google.com/document/d/1uu9dFQbsJOSDWkxnNjpp-xhIMjiTfGd12UqpSmPq8RA/edit#

"""


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const Expat_INCLUDED = 1
const Expat_External_INCLUDED = 1

# Skipping MacroDefinition: XMLPARSEAPI ( type ) XMLIMPORT type XMLCALL
# Skipping MacroDefinition: XML_TRUE ( ( XML_Bool ) 1 )
# Skipping MacroDefinition: XML_FALSE ( ( XML_Bool ) 0 )

# begin enum XML_Status
typealias XML_Status UInt32
const XML_STATUS_ERROR = (UInt32)(0)
const XML_STATUS_OK = (UInt32)(1)
const XML_STATUS_SUSPENDED = (UInt32)(2)
# end enum XML_Status

# Skipping MacroDefinition: XML_GetUserData ( parser ) ( * ( void * * ) ( parser ) )

#=
const XML_GetErrorLineNumber = XML_GetCurrentLineNumber
const XML_GetErrorColumnNumber = XML_GetCurrentColumnNumber
const XML_GetErrorByteIndex = XML_GetCurrentByteIndex
=#
const XML_MAJOR_VERSION = 2
const XML_MINOR_VERSION = 1
const XML_MICRO_VERSION = 0

typealias XML_Char UInt8
typealias XML_LChar UInt8
typealias XML_Index Clong
typealias XML_Size Culong

type XML_ParserStruct
end

typealias XML_Parser Ptr{XML_ParserStruct}
typealias XML_Bool Cuchar

# begin enum XML_Error
typealias XML_Error UInt32
const XML_ERROR_NONE = (UInt32)(0)
const XML_ERROR_NO_MEMORY = (UInt32)(1)
const XML_ERROR_SYNTAX = (UInt32)(2)
const XML_ERROR_NO_ELEMENTS = (UInt32)(3)
const XML_ERROR_INVALID_TOKEN = (UInt32)(4)
const XML_ERROR_UNCLOSED_TOKEN = (UInt32)(5)
const XML_ERROR_PARTIAL_CHAR = (UInt32)(6)
const XML_ERROR_TAG_MISMATCH = (UInt32)(7)
const XML_ERROR_DUPLICATE_ATTRIBUTE = (UInt32)(8)
const XML_ERROR_JUNK_AFTER_DOC_ELEMENT = (UInt32)(9)
const XML_ERROR_PARAM_ENTITY_REF = (UInt32)(10)
const XML_ERROR_UNDEFINED_ENTITY = (UInt32)(11)
const XML_ERROR_RECURSIVE_ENTITY_REF = (UInt32)(12)
const XML_ERROR_ASYNC_ENTITY = (UInt32)(13)
const XML_ERROR_BAD_CHAR_REF = (UInt32)(14)
const XML_ERROR_BINARY_ENTITY_REF = (UInt32)(15)
const XML_ERROR_ATTRIBUTE_EXTERNAL_ENTITY_REF = (UInt32)(16)
const XML_ERROR_MISPLACED_XML_PI = (UInt32)(17)
const XML_ERROR_UNKNOWN_ENCODING = (UInt32)(18)
const XML_ERROR_INCORRECT_ENCODING = (UInt32)(19)
const XML_ERROR_UNCLOSED_CDATA_SECTION = (UInt32)(20)
const XML_ERROR_EXTERNAL_ENTITY_HANDLING = (UInt32)(21)
const XML_ERROR_NOT_STANDALONE = (UInt32)(22)
const XML_ERROR_UNEXPECTED_STATE = (UInt32)(23)
const XML_ERROR_ENTITY_DECLARED_IN_PE = (UInt32)(24)
const XML_ERROR_FEATURE_REQUIRES_XML_DTD = (UInt32)(25)
const XML_ERROR_CANT_CHANGE_FEATURE_ONCE_PARSING = (UInt32)(26)
const XML_ERROR_UNBOUND_PREFIX = (UInt32)(27)
const XML_ERROR_UNDECLARING_PREFIX = (UInt32)(28)
const XML_ERROR_INCOMPLETE_PE = (UInt32)(29)
const XML_ERROR_XML_DECL = (UInt32)(30)
const XML_ERROR_TEXT_DECL = (UInt32)(31)
const XML_ERROR_PUBLICID = (UInt32)(32)
const XML_ERROR_SUSPENDED = (UInt32)(33)
const XML_ERROR_NOT_SUSPENDED = (UInt32)(34)
const XML_ERROR_ABORTED = (UInt32)(35)
const XML_ERROR_FINISHED = (UInt32)(36)
const XML_ERROR_SUSPEND_PE = (UInt32)(37)
const XML_ERROR_RESERVED_PREFIX_XML = (UInt32)(38)
const XML_ERROR_RESERVED_PREFIX_XMLNS = (UInt32)(39)
const XML_ERROR_RESERVED_NAMESPACE_URI = (UInt32)(40)
# end enum XML_Error

# begin enum XML_Content_Type
typealias XML_Content_Type UInt32
const XML_CTYPE_EMPTY = (UInt32)(1)
const XML_CTYPE_ANY = (UInt32)(2)
const XML_CTYPE_MIXED = (UInt32)(3)
const XML_CTYPE_NAME = (UInt32)(4)
const XML_CTYPE_CHOICE = (UInt32)(5)
const XML_CTYPE_SEQ = (UInt32)(6)
# end enum XML_Content_Type

# begin enum XML_Content_Quant
typealias XML_Content_Quant UInt32
const XML_CQUANT_NONE = (UInt32)(0)
const XML_CQUANT_OPT = (UInt32)(1)
const XML_CQUANT_REP = (UInt32)(2)
const XML_CQUANT_PLUS = (UInt32)(3)
# end enum XML_Content_Quant

#=
typealias XML_Content XML_cp
=#

type XML_cp
    _type::XML_Content_Type
    quant::XML_Content_Quant
    name::Ptr{XML_Char}
    numchildren::UInt32
    children::Ptr{XML_cp}
end

typealias XML_ElementDeclHandler Ptr{Void}
typealias XML_AttlistDeclHandler Ptr{Void}
typealias XML_XmlDeclHandler Ptr{Void}

type XML_Memory_Handling_Suite
    malloc_fcn::Ptr{Void}
    realloc_fcn::Ptr{Void}
    free_fcn::Ptr{Void}
end

typealias XML_StartElementHandler Ptr{Void}
typealias XML_EndElementHandler Ptr{Void}
typealias XML_CharacterDataHandler Ptr{Void}
typealias XML_ProcessingInstructionHandler Ptr{Void}
typealias XML_CommentHandler Ptr{Void}
typealias XML_StartCdataSectionHandler Ptr{Void}
typealias XML_EndCdataSectionHandler Ptr{Void}
typealias XML_DefaultHandler Ptr{Void}
typealias XML_StartDoctypeDeclHandler Ptr{Void}
typealias XML_EndDoctypeDeclHandler Ptr{Void}
typealias XML_EntityDeclHandler Ptr{Void}
typealias XML_UnparsedEntityDeclHandler Ptr{Void}
typealias XML_NotationDeclHandler Ptr{Void}
typealias XML_StartNamespaceDeclHandler Ptr{Void}
typealias XML_EndNamespaceDeclHandler Ptr{Void}
typealias XML_NotStandaloneHandler Ptr{Void}
typealias XML_ExternalEntityRefHandler Ptr{Void}
typealias XML_SkippedEntityHandler Ptr{Void}

type XML_Encoding
    map::NTuple{256,Cint}
    data::Ptr{Void}
    convert::Ptr{Void}
    release::Ptr{Void}
end

typealias XML_UnknownEncodingHandler Ptr{Void}

# begin enum XML_Parsing
typealias XML_Parsing UInt32
const XML_INITIALIZED = (UInt32)(0)
const XML_PARSING = (UInt32)(1)
const XML_FINISHED = (UInt32)(2)
const XML_SUSPENDED = (UInt32)(3)
# end enum XML_Parsing

type XML_ParsingStatus
    parsing::XML_Parsing
    finalBuffer::XML_Bool
end

# begin enum XML_ParamEntityParsing
typealias XML_ParamEntityParsing UInt32
const XML_PARAM_ENTITY_PARSING_NEVER = (UInt32)(0)
const XML_PARAM_ENTITY_PARSING_UNLESS_STANDALONE = (UInt32)(1)
const XML_PARAM_ENTITY_PARSING_ALWAYS = (UInt32)(2)
# end enum XML_ParamEntityParsing

type XML_Expat_Version
    major::Cint
    minor::Cint
    micro::Cint
end

# begin enum XML_FeatureEnum
typealias XML_FeatureEnum UInt32
const XML_FEATURE_END = (UInt32)(0)
const XML_FEATURE_UNICODE = (UInt32)(1)
const XML_FEATURE_UNICODE_WCHAR_T = (UInt32)(2)
const XML_FEATURE_DTD = (UInt32)(3)
const XML_FEATURE_CONTEXT_BYTES = (UInt32)(4)
const XML_FEATURE_MIN_SIZE = (UInt32)(5)
const XML_FEATURE_SIZEOF_XML_CHAR = (UInt32)(6)
const XML_FEATURE_SIZEOF_XML_LCHAR = (UInt32)(7)
const XML_FEATURE_NS = (UInt32)(8)
const XML_FEATURE_LARGE_SIZE = (UInt32)(9)
const XML_FEATURE_ATTR_INFO = (UInt32)(10)
# end enum XML_FeatureEnum

type XML_Feature
    feature::XML_FeatureEnum
    name::Ptr{XML_LChar}
    value::Clong
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
#module Utils
#using ..Types
#using LinearAlgebra

#export machNumber, density

function machNumber(v::Float64)
    ma = v / 340.0
    return ma
end

function machNumber(v::Float64,t::Float64)
    #tv = temperature(t, RH, P)
    a = speedOfSound(t)
    ma = v / a
    return ma
end

function temperature(y::Float64;y0 = 0.0, t0 = 288.15, β = -6.5)
    t = (t0+β * (y-y0)*1e-3)-273.15 #deg Celsius
    return t
end
function pressvap(t::Float64;A₀=6.1117675,A₁=0.443986062,A₂=1.43053301*1e-2,A₃=2.65027242*1e-4,A₄=3.02246994*1e-6,A₅=2.03886313*1e-8,A₆=6.38780966*1e-11)
    p = (A₀+A₁*t+A₂*t^2+A₃*t^3+A₄*t^4+A₅*t^5+A₆*t^6)*100
end
function temperature(t::Float64, RH::Float64, P::Float64;MadMav=0.622004944)
    ta = t + 273.15 #kelvin
    Ps = pressvap(t)
    t=ta/(1-(Ps*RH)/(100*P)*(1-MadMav))
end

function speedOfSound(t::Float64; γ = 1.4, R = 287.05287)
    #a = sqrt(γ*R*t)
    a=20.046796*sqrt(t+273.15)
    return a
end

function gravity(y::Float64; g0 = 9.80665, r = 6356766.0)
    g = g0*(r/(r+y))^2
    return g
end
function pressure(y::Float64;p0 = 101325.0, y0 = 0.0, t0 = 288.15, β = -6.5, R = 287.05287)
    p = p0*(1+β/t0*(y-y0)*1e-3)^(-gravity(y)/(β*R))
    return p
end
function density(y::Float64; R = 287.05287)
    t = temperature(y)
    p = pressure(y)
    #ρ = p/(R*t)
    ρ = 0.003483678761*p/(t+273.15)
    return ρ
end

function grav0(lat::Float64)
    g₀ =  9.80665*(1-0.0026*cosd(2*lat))
    return  g₀
end

#end


text:
# Helper functions for common query operators

macro modifier(name::Symbol)
  quote
    $(esc(name))(kv::Pair...) = ($(string("\$", name)) => Dict(kv))
    export $name
  end
end

@modifier query
@modifier orderby
import Base.max
import Base.min
@modifier max
@modifier min
@modifier set
@modifier inc

macro operator(name::Symbol)
  quote
    $(esc(name))(x) = Dict($(string("\$", name)) => x)
    export $name
  end
end

@operator gt
@operator lt
@operator lte
@operator gte
import Base.in
@operator in
@operator nin
@operator ne
@operator eq


text:
module Autoscale

using LightXML
using AWS.Crypto
using Requests
using AWS.AWSEnv
using AWS
using Compat

import AWS: AbstractAWSType, AWSError, AWSResponse, _ename, _parse_from_xml, _prep_params, aws_execute

const AutoScalingError = AWSError
const AutoScalingResponse = AWSResponse

parse_from_xml{T<:AbstractAWSType}(::Type{T}, types::Dict{Symbol,Vector{Pair{Symbol,DataType}}}, pd::Union{LightXML.XMLElement,Void}) = _parse_from_xml(T, types, pd; member_pfx="member.")
prep_params{T<:AbstractAWSType}(obj::T) = _prep_params(obj; member_pfx="member.")

autoscaling_execute(env::AWSEnv, action::String, params::Vector{Tuple}=Vector{Tuple}()) = aws_execute(env, action, "autoscaling", AutoScalingApiVer, params)
autoscaling_execute(env::AWSEnv, action::String, req::AbstractAWSType) = autoscaling_execute(env, action, prep_params(req))
autoscaling_execute(env::AWSEnv, action::String, req::Void) = autoscaling_execute(env, action)

include("autoscaling_types.jl")
include("autoscaling_operations.jl")
include("autoscaling_spec.jl")

export AutoScalingError, AutoScalingResponse

end # module Autoscale


text:
module TestMeanVar

using TapirBenchmarks
using Test

@testset begin
    xs = randn(1000)
    m1, v1 = meanvar_seq(xs)
    m2, v2 = meanvar_threads(xs)
    m3, v3 = meanvar_tapir(xs)
    @test m1 ≈ m2 ≈ m3
    @test v1 ≈ v2 ≈ v3
end

end