text:
"""
Converts a Base.Dates.XX format into a string in HH:MM:SS format.
"""
function duration_to_string(duration::Base.Dates.Hour)
    duration = Dates.value(duration)
    return string(duration)
end

function duration_to_string(duration::Base.Dates.Minute)
    duration = Dates.value(duration)
    minutes = duration % 60
    hours = (duration - minutes) / 60
    minutes = Dates.value(minutes)
    hours = Dates.value(hours)
    return string(hours, ":", minutes)
end

function duration_to_string(duration::Base.Dates.Second)
    duration = Dates.value(duration)
    seconds = duration % 60
    duration = (duration - seconds) / 60
    minutes = duration % 60
    hours = (duration - minutes) / 60
    seconds = convert(Int64, seconds)
    minutes = convert(Int64, minutes)
    hours = convert(Int64, hours)
    return string(hours, ":", minutes, ":", seconds)
end

function duration_to_string(duration::Base.Dates.Millisecond)
    duration = Dates.value(duration)
    duration = floor(duration / 1000)
    seconds = duration % 60
    duration = (duration - seconds) / 60
    minutes = duration % 60
    hours = (duration - minutes) / 60
    seconds = convert(Int64, seconds)
    minutes = convert(Int64, minutes)
    hours = convert(Int64, hours)
    return string(hours, ":", minutes, ":", seconds)
end

"""
Returns only the rows (2D) or entries corresponding to the supplied indicies
"""
function subset(A::Array{T, 1}, idxs::Array{G}) where {T <: Any, G <: Integer}
    return A[idxs]
end

function subset!(A::Array{T, 1}, idxs::Array{G}) where {T <: Any, G <: Integer}
    A = A[idxs]
    nothing
end

function subset(A::Array{T, 2}, idxs::Array{G}) where {T <: Any, G <: Integer}
    return A[idxs, :]
end

function subset!(A::Array{T, 2}, idxs::Array{G}) where {T <: Any, G <: Integer}
    A = A[idxs, :]
    nothing
end


text:
include("../src/ccsa.jl")
using LinearAlgebra
using CatViews
using BenchmarkTools
using Profile
using ProfileView
using PProf
a1=2;b1=0;a2=-1;b2=1;
function f(x)
    return [sqrt(x[2]),(a1*x[1]+b1)^3-x[2],(a2*x[1]+b2)^3-x[2]]
end
function ∇f(x) #(m+1)*n
    return [0 1/2/x[2]; 3*a1*(a1*x[1]+b1)^2 -1; 3*a2*(a2*x[1]+b2)^2 -1]
end
function f_and_∇f(x)
    return f(x),∇f(x)
end
n=2
m=2
σ=[30.0, 30.0] #n
ρ=[100.0, 100.0, 100.0]
x=[1.234, 5.678]
lb=[-Inf, 0.]
opt=CCSAState(n,m,f_and_∇f,ρ,σ,x)
opt.lb=lb
opt.xtol_rel = 1e-4
@time iters=optimize(opt)
ProfileView.@profview optimize(opt)

ext: jl


text:
abstract type DiffPrecSupport end
struct ANTSupport <: DiffPrecSupport end
struct BootStdSupport <: DiffPrecSupport end
struct BootMaxSupport <: DiffPrecSupport end   # not implemented yet
struct DTraceValidationSupport <: DiffPrecSupport end

####################################
#
#   support estimator
#
####################################

function __initSupport(
    Sx::Symmetric, Sy::Symmetric, X, Y)

    p = size(X)
    np = div((p + 1)*p, 2)

    estimSupp = Vector{Vector{Int64}}(undef, np)

    # first stage
    λ = 1.01 * quantile(Normal(), 1. - 0.1 / (p * (p+1)))
    H = Symmetric( skron(Sx, Sy) )
    b = svec(Sx) - svec(Sy)
    x1 = reducedDiffEstimation(H, b, Sx, Sy, X, Y, λ)
    S1 = getReducedSupport(x1)

    # second stage
    ind = 0
    for ci=1:p
      for ri=ci:p
        indElem = sub2indLowerTriangular(p, ri, ci)
        # x2 = invQSymHessian(Sx, Sy, j, X, Y, λ)
        # S2 = getSupport(x2, p)
        x2 = invHessianReduced(H, indElem, Sx, Sy, X, Y, λ)
        S2 = getReducedSupport(x2)
        ind = ind + 1
        estimSupp[ind] = (1:np)[S1 .| S2]
      end
    end

    estimSupp
end

function __threshold(Δpoint::Vector{DiffPrecResultNormal}, p::Int, τ0::Float64)
    Δ = zeros(Float64, p, p)
    it = 0
    for col=1:p
        for row=col:p
            it = it + 1
            τ = τ0 * Δpoint[it].std
            v = abs(Δpoint[it].p) > τ ? Δpoint[it].p : 0.
            if row == col
                Δ[row, row] = v
            else
                Δ[row, col] = v
                Δ[col, row] = v
            end
        end
    end
    sparse(Δ)
end

function __threshold(Δpoint::Vector{DiffPrecResultBoot}, p::Int, τ0::Float64)
    Δ = zeros(Float64, p, p)
    it = 0
    for col=1:p
        for row=col:p
            it = it + 1
            τ = τ0 * std( Δpoint[it].boot_p )
            v = abs(Δpoint[it].p) > τ ? Δpoint[it].p : 0.
            if row == col
                Δ[row, row] = v
            else
                Δ[row, col] = v
                Δ[col, row] = v
            end
        end
    end
    sparse(Δ)
end

function __getΔ(::ANTSupport, X, Y; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    nx, p = size(X)
    ny    = size(Y, 1)

    Sx = Symmetric( cov(X) )
    Sy = Symmetric( cov(Y) )

    eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport
    Δpoint = Array{DiffPrecResultNormal}(undef, div((p+1)*p, 2))

    it = 0
    for col=1:p
        for row=col:p
            it = it + 1
            indS = getLinearSupport(row, col, eS[it])

            Δpoint[it] = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS)
        end
    end

    Δpoint, eS
end

function supportEstimate(::ANTSupport, X, Y, τ0::Float64; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    p = size(X, 2)

    Δpoint, eS = __getΔ(ANTSupport(), X, Y; estimSupport=estimSupport)
    Δsupp = __threshold(Δpoint, p, τ0)

    Δsupp, Δpoint, eS
end

function supportEstimate(::ANTSupport, X, Y, τArr::Vector{Float64}; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    p = size(X, 2)

    Δpoint, eS = __getΔ(ANTSupport(), X, Y; estimSupport=estimSupport)

    Δsupp = Array{SparseMatrixCSC{Float64,Int64}}(undef, length(τArr))
    for i=1:length(τArr)
        Δsupp[i] = __threshold(Δpoint, p, τArr[i])
    end

    Δsupp, Δpoint, eS
end

function __getΔ(::BootStdSupport, X, Y; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    nx, p = size(X)
    ny    = size(Y, 1)

    Sx = Symmetric( cov(X) )
    Sy = Symmetric( cov(Y) )

    eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport
    Δpoint = Array{DiffPrecResultBoot}(undef, div((p+1)*p, 2))

    it = 0
    for col=1:p
        for row=col:p
            it = it + 1
            indS = getLinearSupport(row, col, eS[it])

            Δpoint[it] = estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, row, col, indS)
        end
    end

    Δpoint, eS
end

function supportEstimate(::BootStdSupport, X, Y, τ0::Float64; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    p = size(X, 2)

    Δpoint, eS = __getΔ(BootStdSupport(), X, Y; estimSupport=estimSupport)
    Δsupp = __threshold(Δpoint, p, τ0)

    Δsupp, Δpoint, eS
end

function supportEstimate(::BootStdSupport, X, Y, τArr::Vector{Float64}; estimSupport::Union{Array{BitArray},Nothing}=nothing)
    p = size(X, 2)

    Δpoint, eS = __getΔ(BootStdSupport(), X, Y; estimSupport=estimSupport)

    Δsupp = Array{SparseMatrixCSC{Float64,Int64}}(undef, length(τArr))
    for i=1:length(τArr)
        Δsupp[i] = __threshold(Δpoint, p, τArr[i])
    end

    Δsupp, Δpoint, eS
end

####################################
#
#   support estimator  -- DTrace
#
####################################

# Λarr should be in descending order
function supportEstimate(::DTraceValidationSupport,
    Sx::Symmetric, Sy::Symmetric,
    SxValid::Symmetric, SyValid::Symmetric,
    Λarr::Vector{Float64},
    options=CDOptions())

    numΛ = length(Λarr)
    eΔarr = Vector{SparseMatrixCSC{Float64,Int64}}(undef, numΛ)
    loss2arr = Vector{Float64}(undef, numΛ)
    lossInfarr = Vector{Float64}(undef, numΛ)

    f = CDDirectDifferenceLoss(Sx, Sy)
    x = SymmetricSparseIterate(f.p)
    g = ProxL1(Λarr[1])
    coordinateDescent!(x, f, g, options)
    eΔarr[1] = sparse(Matrix(x))
    loss2arr[1] = CovSel.diffLoss(SxValid, x, SyValid, 2)
    lossInfarr[1] = CovSel.diffLoss(SxValid, x, SyValid, Inf)

    opt = CDOptions(;
        maxIter=options.maxIter,
        optTol=options.optTol,
        randomize=options.randomize,
        warmStart=true,
        numSteps=options.numSteps)

    for i=2:numΛ
        g = ProxL1(Λarr[i])
        coordinateDescent!(x, f, g, opt)
        eΔarr[i] = sparse(Matrix(x))
        loss2arr[i] = CovSel.diffLoss(SxValid, x, SyValid, 2)
        lossInfarr[i] = CovSel.diffLoss(SxValid, x, SyValid, Inf)
    end

    # find min loss
    i2 = argmin(loss2arr)
    iInf = argmin(loss2arr)

    eΔarr, i2, iInf, loss2arr, lossInfarr
end


text:
# Linear Discriminant Analysis

#### Type to represent a linear discriminant functional

"""
A linear discriminant functional can be written as

```math
    f(\\mathbf{x}) = \\mathbf{w}^T \\mathbf{x} + b
```

Here, ``w`` is the coefficient vector, and ``b`` is the bias constant.
"""
struct LinearDiscriminant{T<:Real} <: RegressionModel
    w::Vector{T}
    b::T
end

#### function to solve linear discriminant
"""
    ldacov(C, μp, μn)

Performs LDA given a covariance matrix `C` and both mean vectors `μp` & `μn`.  Returns a linear discriminant functional of type [`LinearDiscriminant`](@ref).

*Parameters*
- `C`: The pooled covariane matrix (*i.e* ``(Cp + Cn)/2``)
- `μp`: The mean vector of the positive class.
- `μn`: The mean vector of the negative class.
"""
function ldacov(C::DenseMatrix{T},
                μp::DenseVector{T},
                μn::DenseVector{T}) where T<:Real

    w = cholesky(C) \ (μp - μn)
    ap = w ⋅ μp
    an = w ⋅ μn
    c = 2 / (ap - an)
    LinearDiscriminant(rmul!(w, c), 1 - c * ap)
end

"""
    ldacov(Cp, Cn, μp, μn)

Performs LDA given covariances and mean vectors. Returns a linear discriminant functional of type [`LinearDiscriminant`](@ref).

*Parameters*
- `Cp`: The covariance matrix of the positive class.
- `Cn`: The covariance matrix of the negative class.
- `μp`: The mean vector of the positive class.
- `μn`: The mean vector of the negative class.

**Note:** The coefficient vector is scaled such that ``w'μp + b = 1`` and ``w'μn + b = -1``.
"""
ldacov(Cp::DenseMatrix{T},
       Cn::DenseMatrix{T},
       μp::DenseVector{T},
       μn::DenseVector{T}) where T<:Real = ldacov(Cp + Cn, μp, μn)

"""
    evaluate(f, x::AbstractVector)

Evaluate the linear discriminant value, *i.e* ``w'x + b``, it returns a real value.
"""
evaluate(f::LinearDiscriminant, x::AbstractVector) = dot(f.w, x) + f.b

"""
    evaluate(f, X::AbstractMatrix)

Evaluate the linear discriminant value, *i.e* ``w'x + b``, for each sample in columns of `X`. The function returns a vector of length `size(X, 2)`.
"""
function evaluate(f::LinearDiscriminant, X::AbstractMatrix)
    R = transpose(X) * f.w
    if f.b != 0
        broadcast!(+, R, R, f.b)
    end
    return R
end

# RegressionModel interface

"""
    predict(f, x::AbstractVector)

Make prediction for the vector `x`. It returns `true` iff `evaluate(f, x)` is positive.
"""
predict(f::LinearDiscriminant, x::AbstractVector) = evaluate(f, x) > 0

"""
    predict(f, X::AbstractMatrix)

Make predictions for the matrix `X`.
"""
predict(f::LinearDiscriminant, X::AbstractMatrix) = Bool[y > 0 for y in evaluate(f, X)]

"""
    coef(f::LinearDiscriminant)

Return the coefficients of the linear discriminant model.
"""
coef(f::LinearDiscriminant) = (f.b, f.w)

"""
    coef(f::LinearDiscriminant)

Return the coefficients' names of the linear discriminant model.
"""
coefnames(f::LinearDiscriminant) = ["Bias", "Weights"]

"""
    dof(f::LinearDiscriminant)

Return the number of degrees of freedom in the linear discriminant model.
"""
dof(f::LinearDiscriminant) = length(f.w)+1

"""
    weights(f::LinearDiscriminant)

Return the linear discriminant model coefficient vector.
"""
weights(f::LinearDiscriminant) = f.w

"""
Get the length of the coefficient vector.
"""
length(f::LinearDiscriminant) = length(f.w)

"""
    fit(LinearDiscriminant, Xp, Xn; covestimator = SimpleCovariance())

Performs LDA given both positive and negative samples. The function accepts follwing parameters:

**Parameters**
- `Xp`: The sample matrix of the positive class.
- `Xn`: The sample matrix of the negative class.

**Keyword arguments:**
- `covestimator`: Custom covariance estimator for between-class covariance. The covariance matrix will be calculated as `cov(covestimator_between, #=data=#; dims=2, mean=zeros(#=...=#)`. Custom covariance estimators, available in other packages, may result in more robust discriminants for data with more features than observations.
"""
function fit(::Type{LinearDiscriminant}, Xp::DenseMatrix{T}, Xn::DenseMatrix{T};
             covestimator::CovarianceEstimator = SimpleCovariance()) where T<:Real
    μp = vec(mean(Xp, dims=2))
    μn = vec(mean(Xn, dims=2))
    Zp = Xp .- μp
    Zn = Xn .- μn
    Cp = calcscattermat(covestimator, Zp)
    Cn = calcscattermat(covestimator, Zn)
    ldacov(Cp, Cn, μp, μn)
end

#==============================================================================#

#### Multiclass LDA Stats

"""
Resulting statistics of the multi-class LDA evaluation.
"""
mutable struct MulticlassLDAStats{T<:Real, M<:AbstractMatrix{T}, N<:AbstractMatrix{T}}
    dim::Int              # sample dimensions
    nclasses::Int         # number of classes
    cweights::Vector{T}   # class weights
    tweight::T            # total sample weight
    mean::Vector{T}       # overall sample mean
    cmeans::Matrix{T}     # class-specific means
    Sw::M                 # within-class scatter matrix
    Sb::N                 # between-class scatter matrix
end

mean(S::MulticlassLDAStats) = S.mean
classweights(S::MulticlassLDAStats) = S.cweights
classmeans(S::MulticlassLDAStats) = S.cmeans

withclass_scatter(S::MulticlassLDAStats) = S.Sw
betweenclass_scatter(S::MulticlassLDAStats) = S.Sb

function MulticlassLDAStats(cweights::Vector{T},
                            mean::Vector{T},
                            cmeans::Matrix{T},
                            Sw::AbstractMatrix{T},
                            Sb::AbstractMatrix{T}) where T<:Real
    d, nc = size(cmeans)
    length(mean) == d || throw(DimensionMismatch("Incorrect length of mean"))
    length(cweights) == nc || throw(DimensionMismatch("Incorrect length of cweights"))
    tw = sum(cweights)
    size(Sw) == (d, d) || throw(DimensionMismatch("Incorrect size of Sw"))
    size(Sb) == (d, d) || throw(DimensionMismatch("Incorrect size of Sb"))
    MulticlassLDAStats(d, nc, cweights, tw, mean, cmeans, Sw, Sb)
end

function multiclass_lda_stats(nc::Int, X::AbstractMatrix{T}, y::AbstractVector{Int};
                              covestimator_within::CovarianceEstimator=SimpleCovariance(),
                              covestimator_between::CovarianceEstimator=SimpleCovariance()) where T<:Real
    # check sizes
    d = size(X, 1)
    n = size(X, 2)
    n ≥ nc || throw(ArgumentError("The number of samples is less than the number of classes"))
    length(y) == n || throw(DimensionMismatch("Inconsistent array sizes."))

    # compute class-specific weights and means
    cmeans, cweights, Z = center(X, y, nc)

    Sw = calcscattermat(covestimator_within, Z)

    # compute between-class scattering
    mean = cmeans * Vector{T}(cweights ./ n)
    U = rmul!(cmeans .- mean, Diagonal(sqrt.(cweights)))
    Sb = calcscattermat(covestimator_between, U)

    return MulticlassLDAStats(Vector{T}(cweights), mean, cmeans, Sw, Sb)
end

#### Multiclass LDA

"""
A multi-class linear discriminant model type has following fields:
- `proj` is the projection matrix
- `pmeans` is the projected means of all classes
- `stats` is an instance of [`MulticlassLDAStats`](@ref) type that captures all statistics computed to train the model (which we will discuss later).
"""
mutable struct MulticlassLDA{T<:Real} <: RegressionModel
    proj::Matrix{T}
    pmeans::Matrix{T}
    stats::MulticlassLDAStats{T}
end

"""
    size(M::MulticlassLDA)

Get the input (*i.e* the dimension of the observation space) and output (*i.e* the dimension of the transformed features) dimensions of the model `M`.
"""
size(M::MulticlassLDA) = size(M.proj)
"""
    length(M::MulticlassLDA)

Get the sample dimensions.
"""
length(M::MulticlassLDA) = M.stats.dim
"""
    projection(M::MulticlassLDA)

Get the projection matrix (of size *d x p*).
"""
projection(M::MulticlassLDA) = M.proj
"""
    mean(M::MulticlassLDA)

Get the overall sample mean vector (of length *d*).
"""
mean(M::MulticlassLDA) = mean(M.stats)

"""
    classmeans(M)

Get the matrix comprised of class-specific means as columns (of size ``(d, m)``).
"""
classmeans(M::MulticlassLDA) = classmeans(M.stats)
"""
    classweights(M)

Get the weights of individual classes (a vector of length ``m``). If the samples are not weighted,
the weight equals the number of samples of each class.
"""
classweights(M::MulticlassLDA) = classweights(M.stats)
"""
withinclass_scatter(M)

    Get the within-class scatter matrix (of size ``(d, d)``).
"""
withclass_scatter(M::MulticlassLDA) = withclass_scatter(M.stats)
"""
betweenclass_scatter(M)

    Get the between-class scatter matrix (of size ``(d, d)``).
"""
betweenclass_scatter(M::MulticlassLDA) = betweenclass_scatter(M.stats)

"""
    predict(M::MulticlassLDA, x)

Transform input sample(s) in `x` to the output space of MC-LDA model `M`. Here, `x` can be either a sample vector or a matrix comprised of samples in columns.
"""
predict(M::MulticlassLDA, x::AbstractVecOrMat{T}) where {T<:Real} = M.proj'x

"""
    fit(MulticlassLDA, nc, X, y; ...)

Perform multi-class LDA over a given data set `X` and collecttion of labels `y`.

This function returns the resultant multi-class LDA model as an instance of [`MulticlassLDA`](@ref).

*Parameters*

- `nc`:  the number of classes
- `X`:   the matrix of input samples, of size `(d, n)`. Each column in `X` is an observation.
- `y`:   the vector of class labels, of length `n`. Each element of `y` must be an integer between `1` and `nc`.

**Keyword arguments**

- `method`: The choice of methods:
    - `:gevd`: based on generalized eigenvalue decomposition (*default*).
    - `:whiten`: first derive a whitening transform from `Sw` and then solve the problem based on eigenvalue
    decomposition of the whiten `Sb`.
- `outdim`: The output dimension, i.e. dimension of the transformed space `min(d, nc-1)`
- `regcoef`: The regularization coefficient (*default:* `1.0e-6`). A positive value `regcoef * eigmax(Sw)`
    is added to the diagonal of `Sw` to improve numerical stability.
- `covestimator_between`: Custom covariance estimator for between-class covariance (*default:* `SimpleCovariance()`).
    The covariance matrix will be calculated as `cov(covestimator_between, #=data=#; dims=2, mean=zeros(#=...=#))`.
    Custom covariance estimators, available in other packages, may result in more robust discriminants for data
    with more features than observations.
- `covestimator_within`:  Custom covariance estimator for within-class covariance (*default:* `SimpleCovariance()`).
    The covariance matrix will be calculated as `cov(covestimator_within, #=data=#; dims=2, mean=zeros(nc))`.
    Custom covariance estimators, available in other packages, may result in more robust discriminants for data
    with more features than observations.

**Notes:**

The resultant projection matrix ``P`` satisfies:
```math
\\mathbf{P}^T (\\mathbf{S}_w + \\kappa \\mathbf{I}) \\mathbf{P} = \\mathbf{I}
```
Here, ``\\kappa`` equals `regcoef * eigmax(Sw)`. The columns of ``P`` are arranged in descending order of
the corresponding generalized eigenvalues.

Note that [`MulticlassLDA`](@ref) does not currently support the normalized version using ``\\mathbf{S}_w^*`` and
``\\mathbf{S}_b^*`` (see [`SubspaceLDA`](@ref)).
"""
function fit(::Type{MulticlassLDA}, nc::Int, X::DenseMatrix{T}, y::AbstractVector{Int};
             method::Symbol=:gevd,
             outdim::Int=min(size(X,1), nc-1),
             regcoef::T=T(1.0e-6),
             covestimator_within::CovarianceEstimator=SimpleCovariance(),
             covestimator_between::CovarianceEstimator=SimpleCovariance()) where T<:Real

    multiclass_lda(multiclass_lda_stats(nc, X, y;
                                        covestimator_within=covestimator_within,
                                        covestimator_between=covestimator_between);
                   method=method,
                   regcoef=regcoef,
                   outdim=outdim)
end

function multiclass_lda(S::MulticlassLDAStats{T};
                        method::Symbol=:gevd,
                        outdim::Int=min(S.dim, S.nclasses-1),
                        regcoef::T=T(1.0e-6)) where T<:Real

    P = mclda_solve(S.Sb, S.Sw, method, outdim, regcoef)
    MulticlassLDA(P, P'S.cmeans, S)
end

mclda_solve(Sb::AbstractMatrix{T}, Sw::AbstractMatrix{T}, method::Symbol, p::Int, regcoef::T) where T<:Real =
    mclda_solve!(copy(Sb), copy(Sw), method, p, regcoef)

function mclda_solve!(Sb::AbstractMatrix{T},
                      Sw::AbstractMatrix{T},
                      method::Symbol, p::Int, regcoef::T) where T<:Real

    p <= size(Sb, 1) || throw(ArgumentError("p cannot exceed sample dimension."))

    if method == :gevd
        regularize_symmat!(Sw, regcoef)
        E = eigen!(Symmetric(Sb), Symmetric(Sw))
        ord = sortperm(E.values; rev=true)
        P = E.vectors[:, ord[1:p]]

    elseif method == :whiten
        W = _lda_whitening!(Sw, regcoef)
        wSb = transpose(W) * (Sb * W)
        Eb = eigen!(Symmetric(wSb))
        ord = sortperm(Eb.values; rev=true)
        P = W * Eb.vectors[:, ord[1:p]]

    else
        throw(ArgumentError("Invalid method name $(method)"))
    end
    return P::Matrix{T}
end

function _lda_whitening!(C::AbstractMatrix{T}, regcoef::T) where T<:Real
    n = size(C,1)
    E = eigen!(Symmetric(C))
    v = E.values
    a = regcoef * maximum(v)
    for i = 1:n
        @inbounds v[i] = 1.0 / sqrt(v[i] + a)
    end
    return rmul!(E.vectors,  Diagonal(v))
end

#### SubspaceLDA

"""Subspace LDA model type has following fields:
- `projw`: the projection matrix of the subspace spanned by the between-class scatter
- `projLDA`: the projection matrix of the subspace spanned by the within-class scatter
- `λ`: the projection eigenvalues
- `cmeans`: the class centroids
- `cweights`: the class weights
"""
struct SubspaceLDA{T<:Real} <: RegressionModel
    projw::Matrix{T}
    projLDA::Matrix{T}
    λ::Vector{T}
    cmeans::Matrix{T}
    cweights::Vector{Int}
end

"""
    size(M)

Get the input (*i.e* the dimension of the observation space) and output (*i.e* the dimension of the subspace projection)
dimensions of the model `M`.
"""
size(M::SubspaceLDA) = (size(M.projw,1), size(M.projLDA, 2))
"""
    length(M)

Get dimension of the LDA model.
"""
length(M::SubspaceLDA) = size(M.projLDA, 2)
"""
    predict(M::SubspaceLDA, x)

Transform input sample(s) in `x` to the output space of LDA model `M`.
Here, `x` can be either a sample vector or a matrix comprised of samples in columns.
"""
predict(M::SubspaceLDA, x::AbstractVecOrMat{T}) where {T<:Real} = M.projLDA' * (M.projw' * x)
"""
    projection(M)

Get the projection matrix.
"""
projection(M::SubspaceLDA) = M.projw * M.projLDA

"""
    mean(M::SubspaceLDA)

Returns the mean vector of the subspace LDA model `M`.
"""
mean(M::SubspaceLDA) = vec(sum(M.cmeans * Diagonal(M.cweights / sum(M.cweights)), dims=2))

"""
    eigvals(M::SubspaceLDA)

Get the eigenvalues of the subspace LDA model `M`.
"""
eigvals(M::SubspaceLDA) = M.λ

classmeans(M::SubspaceLDA) = M.cmeans
classweights(M::SubspaceLDA) = M.cweights
fit(::Type{F}, X::AbstractMatrix{T}, nc::Int, label::AbstractVector{Int}) where {T<:Real, F<:SubspaceLDA} =
    fit(F, X, label, nc)

"""
    fit(SubspaceLDA, X, labels; normalize=true)

Fit an subspace projection of LDA model using the equivalent of ``\\mathbf{S}_w^*`` and ``\\mathbf{S}_b^*```.

Note: Subspace LDA also supports the normalized version of LDA via the `normalize` keyword.
"""
function fit(::Type{F}, X::AbstractMatrix{T},
             label::AbstractVector{Int},
             nc=maximum(label);
             normalize::Bool=false) where {T<:Real, F<:SubspaceLDA}
    d, n = size(X, 1), size(X, 2)
    n ≥ nc || throw(ArgumentError("The number of samples is less than the number of classes"))
    length(label) == n || throw(DimensionMismatch("Inconsistent array sizes."))
    # Compute centroids, class weights, and deviation from centroids
    # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = center(X, label, nc)
    dmeans = cmeans .- (normalize ? mean(cmeans, dims=2) : cmeans * (cweights / T(n)))
    Hb = normalize ? dmeans : dmeans * Diagonal(convert(Vector{T}, sqrt.(cweights)))
    if normalize
        Hw /= T(sqrt(n))
    end
    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full=false)
    keep = Σw .> sqrt(eps(T)) * maximum(Σw)
    projw = Uw[:,keep]
    pHb = projw' * Hb
    pHw = projw' * Hw
    λ, G = lda_gsvd(pHb, pHw, cweights)
    SubspaceLDA(projw, G, λ, cmeans, cweights)
end

# Reference: Howland & Park (2006), "Generalizing discriminant analysis
# using the generalized singular value decomposition", IEEE
# Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
function lda_gsvd(Hb::AbstractMatrix{T}, Hw::AbstractMatrix{T}, cweights::AbstractVector{Int}) where T<:Real
    nc = length(cweights)
    K = vcat(Hb', Hw')
    P, R, Q = svd(K, full=false)
    keep = R .> sqrt(eps(T))*maximum(R)
    R = R[keep]
    Pk = P[1:nc, keep]
    U, ΣA, W = svd(Pk)
    ncnz = sum(cweights .> 0)
    G = Q[:,keep]*(Diagonal(1 ./ R) * W[:,1:ncnz-1])
    # Normalize
    Gw = G' * Hw
    nrm = Gw * Gw'
    G = G ./ reshape(sqrt.(diag(nrm)), 1, ncnz-1)
    # Also get the eigenvalues
    Gw = G' * Hw
    Gb = G' * Hb
    λ = diag(Gb * Gb')./diag(Gw * Gw')
    λ, G
end

function center(X::AbstractMatrix{T}, label::AbstractVector{Int}, nc=maximum(label)) where T<:Real
    d, n = size(X,1), size(X,2)
    # Calculate the class weights and means
    cmeans = zeros(T, d, nc)
    cweights = zeros(Int, nc)
    for j = 1:n
        k = label[j]
        for i = 1:d
            cmeans[i,k] += X[i,j]
        end
        cweights[k] += 1
    end
    for j = 1:nc
        cw = cweights[j]
        cw == 0 && continue
        for i = 1:d
            cmeans[i,j] /= cw
        end
    end
    # Compute differences from the means
    dX = Matrix{T}(undef, d, n)
    for j = 1:n
        k = label[j]
        for i = 1:d
            dX[i,j] = X[i,j] - cmeans[i,k]
        end
    end
    cmeans, cweights, dX
end


text:
module SLFN

using Distances
import Optim

# stdlib
using LinearAlgebra, Statistics, InteractiveUtils, Printf, Random

include("stable_regression.jl")

export AbstractLinReg, AbstractLSMethod, AbstractLADMethod,
    OLS, LSSVD, LSLdiv, RLSTikhonov, RLST, LADPP, LADDP,
    RLADPP, RALDDP, RLSSVD, regress

abstract type AbstractSLFN  end

export AbstractActivation, Sigmoid, SoftPlus, Tanh, Relu, AbstractSLFN,
    fit!, isexact, input_to_node, hidden_out, n_coefs, activate, activate!,
    AlgebraicNetwork, ELM, TSELM, EIELM, OSELM, ROSELM, LSIELM,
    Linear, RBF, Gaussian

const WARNINGS = [false]
dont_warn_me() = (WARNINGS[1] = false; return)
warn_me() = (WARNINGS[1] = true; return)

isexact(elm::AbstractSLFN) = false

function input_to_node(elm::AbstractSLFN, x, Wt=elm.Wt, d=elm.d)
    input_to_node(elm.neuron_type, x, Wt, d)
end

function hidden_out(elm::AbstractSLFN, x, Wt=elm.Wt, d=elm.d)
    hidden_out(elm.neuron_type, x, Wt, d)
end

n_coefs(elm::AbstractSLFN) = length(elm.Wt) + length(elm.d) + length(elm.v)

include("activations.jl")
include("node_inputs.jl")
include("elm.jl")
include("oselm.jl")
include("tselm.jl")
include("eielm.jl")
include("roselm.jl")
include("algebraic.jl")
include("lsielm.jl")
include("optimelm.jl")

standardize(x::AbstractMatrix, μ::AbstractVector, σ::AbstractVector) = (x .- μ') ./ σ'

standardize(x::AbstractVector, μ::AbstractVector, σ::AbstractVector) =
    length(μ) == 1 ? (x .- μ') ./ σ' : (x .- μ) ./ σ

function standardize(x::Number, μ, σ)
    length(μ) == length(σ) == 1 || error("x shouldn't be scalar")
    (x - μ[1]) / σ[1]
end

for T in InteractiveUtils.subtypes(AbstractSLFN)
    @eval function (elm::$(T))(_x′) where T
        @assert size(_x′, 2) == elm.q "wrong input dimension"
        x′ = isdefined(elm, :μx) ? standardize(_x′, elm.μx, elm.σx) : _x′
        if length(elm.v) == size(elm.Wt, 2) + 1 # has intercept
            hidden_out(elm, x′) * elm.v[2:end] .+ elm.v[1]
        else
            hidden_out(elm, x′) * elm.v
        end
    end

    @eval function (elm::$(T))(_x′::Number) where T
        x′ = isdefined(elm, :μx) ? standardize(_x′, elm.μx, elm.σx) : _x′
        if length(elm.v) == size(elm.Wt, 2) + 1 # has intercept
            out_vec = hidden_out(elm, x′) * elm.v[2:end] + elm.v[1]
        else
            out_vec = hidden_out(elm, x′) * elm.v
        end

        @assert length(out_vec) == 1 "There's a bug -- file an issue"
        out_vec[1]
    end
end

Base.getindex(elm::AbstractSLFN, x) = elm(x)

end # module


text:
# Julia wrapper for header: /usr/include/libavutil/camellia.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_camellia_alloc,
    av_camellia_init,
    av_camellia_crypt

function av_camellia_alloc()
    ccall((:av_camellia_alloc, libavutil), Ptr{AVCAMELLIA}, ())
end

function av_camellia_init(ctx, key, key_bits::Integer)
    ccall((:av_camellia_init, libavutil), Cint, (Ptr{AVCAMELLIA}, Ptr{UInt8}, Cint), ctx, key, key_bits)
end

function av_camellia_crypt(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_camellia_crypt, libavutil), Cvoid, (Ptr{AVCAMELLIA}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end


text:
#------------------------------------------------------------------------------
#  Example.jl
#
#
#
#
#
#
#
#
#
#
#
#
#
#  Paul.Soderlind@unisg.ch   May 2013, to Julia Jan 2016
#------------------------------------------------------------------------------

using Optim, SpecialFunctions, Statistics
include("HistogramFitfn.jl")         #functions used

CaseQ = 3               #try 1,2, or 3 (see Prob below)
#------------------------------------------------------------------------------

                                                #DATA
#catogories, -Inf<x<=Bound(1),Bond(1)<x<=Bound(2),..,Bound(n)<x<=Inf
CatBounds = [-2.00;-1.00; 0.00; 1.00; 2.00; 3.00; 4.00; 5.00; 6.00]
#mid points, first and last interval are artificially closed
CatMid    = [-2.50;-1.50;-0.50; 0.50; 1.50; 2.50; 3.50; 4.50; 5.50; 6.50]

if CaseQ == 1          #only one interval has non-zero prob
  Prob = [ 0.00; 0.00; 0.00; 0.00; 1.00; 0.00; 0.00; 0.00; 0.00; 0.00]
elseif CaseQ == 2      #two intervals have non-zero prob
  Prob = [ 0.00; 0.00; 0.00; 0.00; 0.30; 0.70; 0.00; 0.00; 0.00; 0.00]
elseif CaseQ == 3      #3 or more intervals have non-zero prob
  Prob = [ 0.00; 0.00; 0.00; 0.10; 0.30; 0.40; 0.15; 0.05; 0.00; 0.00]   #probabilities
end
#------------------------------------------------------------------------------

if sum(Prob) != 1
  error("probabilities do not sum to one")
end

MeanCrude = sum(Prob .* CatMid)                    #crude mean and variance
VarCrude  = sum(Prob .* (CatMid.-MeanCrude).^2)

BinWidth    = mean(diff(CatBounds))
VarSheppard = VarCrude - BinWidth^2/12       #variance, Sheppard's correction
#------------------------------------------------------------------------------

NActiveCat = sum(Prob .> 0)      #no. intervals with non-zero probabilities

if NActiveCat == 1        #if only one active intervals: triangular distribution
   parM = Any[MeanCrude;sqrt(BinWidth^2/24);NActiveCat]
elseif NActiveCat == 2        #if only two active intervals: N(MeanCrude,VarSheppard)
   parM = Any[MeanCrude;sqrt(VarSheppard);NActiveCat]
elseif NActiveCat > 2         #if three or more active intervals: N(estimate,estimate)
  par0 = [MeanCrude;sqrt(VarSheppard)]
  Sol  = optimize(par->NormalHistLoss(par,Prob,CatBounds),par0)
  par1 = Optim.minimizer(Sol)
  parM = Any[par1;NActiveCat]
end
println("\n[mean,std]=$(round.(parM[1:2]',digits=3)), no. active intervals: $(parM[3])")
#------------------------------------------------------------------------------

#Comment out this if you do not have PyPlot installed.

using PyPlot
#PyPlot.svg(true)           #for ipynb notebooks
close("all")

if NActiveCat == 1                           #assumed triangular dist
  whichBin = findfirst(Prob .== 1)
  (a,b,c) = (CatBounds[whichBin-1],CatBounds[whichBin],CatMid[whichBin]) #for distribution
  y    = range(a,stop=b,length=101)
  pdfy = TriangularPdfPs.(y,a,b,c)
  figure()
    bar(CatMid,Prob/BinWidth,width=BinWidth,color="lightgray",edgecolor="black",align="center")
    xlim(-2,6)
    title("Histogram and assumed triangular distribution")
    plot(y,pdfy)
    #display(gcf())          #uncomment in Atom/Juno
elseif NActiveCat >= 2                           #fitted N(mu,s^2)
  y    = range(-2,stop=6,length=101)
  pdfy = NormPdfPs.(y,parM[1],parM[2]^2)
  figure()
    bar(CatMid,Prob/BinWidth,width=BinWidth,color="lightgray",edgecolor="black",align="center")
    xlim(-2,6)
    title("Histogram and fitted \$N(\\mu,\\sigma^2)\$")
    plot(y,pdfy)
    #display(gcf())          #uncomment in Atom/Juno
end

#-------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
using BenchmarkTools
using JLD

include("stinger.jl")
include("parallelbfslg.jl")

function createsuite(bench::GraphBenchmarkSpec, scale::Int, edgefactor::Int)
    suite = BenchmarkGroup([])
    suite["lg"] = bfsbenchsuite(DiGraph, bench, scale, edgefactor)
    suite["stinger"] = bfsbenchsuite(StingerGraph, bench, scale, edgefactor)
    suite
end

function bench(bench::GraphBenchmarkSpec, scaleRange::Range{Int}, edgefactor::Int)
    benchname = split("$(typeof(bench))", '.')[end]
    suite = BenchmarkGroup([benchname])
    for scale in scaleRange
        suite["$(scale)_$(edgefactor)"] = createsuite(bench, scale, edgefactor)
    end
    threads = Threads.nthreads()
    info("Running the $benchname benchmarks for threads=$threads")
    results = run(suite)
    jldopen(joinpath(dirname(@__FILE__), "bfs_$threads.jld"), "w") do f
        write(f, "bfssuite_$(benchname)", results)
    end
end

bench(Graph500(), 10:12, 16)


text:
"""
    mutable struct Dense <: Block

Applies y = fn(w * x .+ b) transformation to the incoming data x
"""
mutable struct Dense <: Block
    w::VarOrNil
    b::VarOrNil
    f::Function
    # type specilized may be CuArray/AFArray/ClArray/Array etc
    function Dense(inputSize::Int, hiddenSize::Int, fn::Function=relu; type::Type=Array{Float32})
        T = eltype(type)
        a = sqrt(T(2/inputSize))
        w = randn(T, hiddenSize, inputSize) .* a
        b = randn(T, hiddenSize,         1) .* a
        new(Variable{type}(w,true,true,true), Variable{type}(b,true,true,true), fn)
    end
    function Dense(fn::Function)
        new(nothing, nothing, fn)
    end
end

function clone(this::Dense; type::Type=Array{Float32})
    cloned = Dense(this.f)
    cloned.w = clone(this.w, type=type)
    cloned.b = clone(this.b, type=type)
    return cloned
end

# pretty show
function Base.show(io::IO, m::Dense)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    print(io, "Dense($(SIZE[2]), $(SIZE[1]), $(m.f); type=$TYPE)")
end

mutable struct MLP <: Block
    layers::Vector{Dense}
    function MLP(topology::Vector{Int}; type::Type=Array{Float32})
        n = length(topology) - 1
        layers = Vector{Dense}(undef, n)
        for i = 1:n
            layers[i] = Dense(topology[i], topology[i+1], relu; type=type)
        end
        new(layers)
    end

    function MLP(topology::Vector{Int}, fn::Vector{F}; type::Type=Array{Float32}) where F
        n = length(topology) - 1
        layers = Vector{Dense}(undef, n)
        for i = 1:n
            layers[i] = Dense(topology[i], topology[i+1], fn[i]; type=type)
        end
        new(layers)
    end
end

# pretty show
function Base.show(io::IO, m::MLP)
    print(io, "MLP\n      (\n          ")
    join(io,    m.layers, "\n          ")
    print(io,                "\n      )")
end

Base.length(m::MLP)             = length(m.layers)
Base.lastindex(m::MLP)          = length(m.layers)
Base.getindex(m::MLP, k...)     =        m.layers[k...]
Base.setindex!(m::MLP, v, k...) =       (m.layers[k...] = v)
Base.firstindex(m::MLP) = 1
Base.iterate(m::MLP, i=firstindex(m)) = i>length(m) ? nothing : (m[i], i+1)

function forward(m::Dense, x::Variable)
    f = m.f
    w = m.w
    b = m.b
    return f( matAddVec(w*x, b) )
end

function forward(m::MLP, x::Variable)
    for i = 1:length(m)
        x = forward(m[i], x)
    end
    return x
end

function predict(m::Dense, x)
    f = m.f
    w = m.w.value
    b = m.b.value
    return f(w * x .+ b)
end

function predict(m::MLP, x)
    for i = 1:length(m)
        x = predict(m[i], x)
    end
    return x
end

"""
    unbiasedof(m::Dense)

unbiased weights of Dense block
"""
function unbiasedof(m::Dense)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function weightsof(m::Dense)
    weights = Vector(undef, 2)
    weights[1] = m.w.value
    weights[2] = m.b.value
    return weights
end

"""
    unbiasedof(m::MLP)

unbiased weights of MLP block
"""
function unbiasedof(m::MLP)
    weights = Vector(undef,0)
    for i = 1:length(m)
        append!(weights, unbiasedof(m[i]))
    end
    return weights
end

function weightsof(m::MLP)
    weights = Vector(undef,0)
    for i = 1:length(m)
        append!(weights, weightsof(m[i]))
    end
    return weights
end

function gradsof(m::Dense)
    grads = Vector(undef, 2)
    grads[1] = m.w.delta
    grads[2] = m.b.delta
    return grads
end

function gradsof(m::MLP)
    grads = Vector(undef,0)
    for i = 1:length(m)
        append!(grads, gradsof(m[i]))
    end
    return grads
end

function zerograds!(m::Dense)
    for v in gradsof(m)
        v .= 0.0
    end
end

function zerograds!(m::MLP)
    for v in gradsof(m)
        v .= 0.0
    end
end

function paramsof(m::Dense)
    params = Vector{Variable}(undef,2)
    params[1] = m.w
    params[2] = m.b
    return params
end

function xparamsof(m::Dense)
    xparams = Vector{XVariable}(undef,2)
    xparams[1] = ('w', m.w)
    xparams[2] = ('b', m.b)
    return xparams
end

function paramsof(m::MLP)
    params = Vector{Variable}(undef,0)
    for i = 1:length(m)
        append!(params, paramsof(m[i]))
    end
    return params
end

function xparamsof(m::MLP)
    params = Vector{XVariable}(undef,0)
    for i = 1:length(m)
        append!(params, xparamsof(m[i]))
    end
    return params
end

function nparamsof(m::Dense)
    lw = length(m.w)
    lb = length(m.b)
    return (lw + lb)
end

function bytesof(model::Dense, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function nparamsof(m::MLP)
    num = 0
    for i = 1:length(m)
        num += nparamsof(m[i])
    end
    return num
end

function bytesof(model::MLP, unit::String="MB")
    n = nparamsof(model) * elsizeof(model[1].w)
    return blocksize(n, uppercase(unit))
end

function to(type::Type, m::Dense)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return m
end

function to!(type::Type, m::Dense)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return nothing
end

function to(type::Type, m::MLP)
    for layer in m
        layer = to(type, layer)
    end
    return m
end

function to!(type::Type, m::MLP)
    for layer in m
        to!(type, layer)
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
module Optimizers
using Statistics
export Optimizer, OptimizerContext, accuracy, logcrossentropy, fit!, σ, mse, logcrossentropy!, sigmoid!, stopearly!, addvalidationset!, postprocess!

"""
    Optimizer

Abstract type for all optimizers
"""
abstract type Optimizer end
abstract type OptimizerContext end

function fit!(optim::Optimizer, β, X, y)
    @error "Not implemented"
end

function addvalidationset!(optim::Optimizer, X, y, earlystopping=false)
    optim.validationset = (X, y)
    optim.validationloss = zeros(Float32, optim.context.maxiterations)
    optim.earlystopping = earlystopping
    optim.hasvalidationset = true
end

function stopearly!(optim::Optimizer, β, step::Int)::Bool
    if !optim.hasvalidationset
        return false
    end

    history = 3
    yhat = σ.(optim.validationset[1]*β)
    optim.validationloss[step] = logcrossentropy(optim.validationset[2],
                                                 yhat)
    if step ≤ history || !optim.earlystopping
        return false
    end
    limit = mean(optim.validationloss[(step - history):(step - 1)])
    if optim.validationloss[step] > limit
        return true
    else
        return false
    end
end

function postprocess!(optim::Optimizer)::Nothing
    resize!(optim.loss, optim.iterations)
    if optim.hasvalidationset
        if optim.converged
            resize!(optim.validationloss, optim.iterations-1)
        else
            resize!(optim.validationloss, optim.iterations)
        end
    end
    return
end

function accuracy(observations, predictions)
    truepositive = observations .== Int.(predictions .> 0.5)
    sum(truepositive)/length(observations)
end

function lognull(x::T)::T where T<:Real
    x == zero(x) && return zero(x)
    log(x)
end

function logcrossentropy(observations, predictions)
    -mean(@. observations * lognull(predictions) + (1-observations) * lognull(1-predictions))
end

function logcrossentropy!(output, tmp, observations, predictions)
    @. tmp = - observations * lognull(predictions) - (1-observations) * lognull(1-predictions)
    mean!(output, tmp)
end

function mse(observations, predictions)
    mean((observations .- predictions).^2)
end

sigmoid(x) = 1 /(1 + exp(-x))
σ = sigmoid
function sigmoid!(y, x)
    for i in eachindex(y)
        y[i] = 1/(1 + exp(-x[i]))
    end
end

end


text:
using AxisArrays
using HTTP
using JSON
using JuMP
using MathOptInterface
using MutableArithmetics
using Cbc


text:
"""
      SeisPatch(in::String,out::String;<keyword arguments>)

Creates overlapping 5d patches from a 5d volume

# Arguments
* `in::String`: input filename (data should have grid information in headers)
* `out::String`: prefix for output filenames

# Keyword arguments
* `style="sxsygxgy"`: bin style. Options: "mxmyhxhy","mxmyhaz","sxsyhxhy","gxgyhxhy","sxsyhaz","gxgyhaz"
* `min_isx=0`,`max_isx=0`,`min_isy=0`,`max_isy=0`: grid extreme values for sources
* `min_igx=0`,`max_igx=0`,`min_igy=0`,`max_igy=0`: grid extreme values for receivers
* `min_imx=0`,`max_imx=0`,`min_imy=0`,`max_imy=0`: grid extreme values for midpoints
* `min_ihx=0`,`max_ihx=0`,`min_ihy=0`,`max_ihy=0`: grid extreme values for offsets
* `min_ih=0`,`max_ih=0`,`min_iaz=0`,`max_iaz=0`: grid extreme values for azimuth and offset
* `it_WL=9e9`,`it_WO=0` : length and overlapping samples in time patches
* `ix1_WL=9e9`,`ix1_WO=0`:length and overlapping samples in first space dimension
* `ix2_WL=9e9`,`ix2_WO=0`,`ix3_WL=9e9`,`ix3_WO=0`,`ix4_WL=9e9`,`ix4_WO=0`

# Output
`filename,npatch`: String Array with the file name of the data patches, number of patches created

*Credits: A. Stanton, F. Carozzi, 2017*
"""
function SeisPatch(in::String, out::String; style="sxsygxgy", min_isx=0,
                   max_isx=0, min_isy=0, max_isy=0, min_igx=0, max_igx=0,
                   min_igy=0, max_igy=0, min_imx=0, max_imx=0, min_imy=0,
                   max_imy=0, min_ihx=0, max_ihx=0, min_ihy=0, max_ihy=0,
                   min_ih=0, max_ih=0, min_iaz=0, max_iaz=0, it_WL=9e9, it_WO=0,
                   ix1_WL=9e9, ix1_WO=0, ix2_WL=9e9, ix2_WO=0, ix3_WL=9e9,
                   ix3_WO=0, ix4_WL=9e9, ix4_WO=0)

if (style == "sxsygxgy")
  key = ["t","isx","isy","igx","igy"]
	min_ix1 = min_isx
	max_ix1 = max_isx
	min_ix2 = min_isy
	max_ix2 = max_isy
	min_ix3 = min_igx
	max_ix3 = max_igx
	min_ix4 = min_igy
	max_ix4 = max_igy

elseif (style=="mxmyhxhy")
	key = ["t","imx","imy","ihx","ihy"]
	min_ix1 = min_imx
	max_ix1 = max_imx
	min_ix2 = min_imy
	max_ix2 = max_imy
	min_ix3 = min_ihx
	max_ix3 = max_ihx
	min_ix4 = min_ihy
	max_ix4 = max_ihy

elseif (style=="mxmyhaz")
	key = ["t","imx","imy","ih","iaz"]
	min_ix1 = min_imx
	max_ix1 = max_imx
	min_ix2 = min_imy
	max_ix2 = max_imy
	min_ix3 = min_ih
	max_ix3 = max_ih
	min_ix4 = min_iaz
	max_ix4 = max_iaz

elseif (style=="sxsyhxhy")
	key = ["t","isx","isy","ihx","ihy"]
	min_ix1 = min_isx
	max_ix1 = max_isx
	min_ix2 = min_isy
	max_ix2 = max_isy
	min_ix3 = min_ihx
	max_ix3 = max_ihx
	min_ix4 = min_ihy
	max_ix4 = max_ihy

elseif (style=="gxgyhxhy")
	key = ["t","igx","igy","ihx","ihy"]
	min_ix1 = min_igx
	max_ix1 = max_igx
	min_ix2 = min_igy
	max_ix2 = max_igy
	min_ix3 = min_ihx
	max_ix3 = max_ihx
	min_ix4 = min_ihy
	max_ix4 = max_ihy

elseif (style=="sxsyhaz")
	key = ["t","isx","isy","ih","iaz"]
	min_ix1 = min_isx
	max_ix1 = max_isx
	min_ix2 = min_isy
	max_ix2 = max_isy
	min_ix3 = min_ih
	max_ix3 = max_ih
	min_ix4 = min_iaz
	max_ix4 = max_iaz

elseif (style=="gxgyhaz")
	key = ["t","igx","igy","ih","iaz"]
	min_ix1 = min_igx
	max_ix1 = max_igx
	min_ix2 = min_igy
	max_ix2 = max_igy
	min_ix3 = min_ih
	max_ix3 = max_ih
	min_ix4 = min_iaz
	max_ix4 = max_iaz

else
	error("style not defined.")
end

  nx1 = max_ix1 - min_ix1 + 1
  nx2 = max_ix2 - min_ix2 + 1
  nx3 = max_ix3 - min_ix3 + 1
  nx4 = max_ix4 - min_ix4 + 1

  filename_data = ParseDataName(in)
  filename_headers = ParseHeaderName(in)
  extent = ReadTextHeader(in)

  nt = extent.n1
  dt = extent.d1
  ot = extent.o1

  it_WL  = it_WL  > nt  ? nt  : it_WL
  ix1_WL = ix1_WL > nx1 ? nx1 : ix1_WL
  ix2_WL = ix2_WL > nx2 ? nx2 : ix2_WL
  ix3_WL = ix3_WL > nx3 ? nx3 : ix3_WL
  ix4_WL = ix4_WL > nx4 ? nx4 : ix4_WL

  tmax = ot + dt*nt
  it_NW = Int(floor(nt/(it_WL-it_WO)))
  ix1_NW = Int(floor(nx1/(ix1_WL-ix1_WO)))
  ix2_NW = Int(floor(nx2/(ix2_WL-ix2_WO)))
  ix3_NW = Int(floor(nx3/(ix3_WL-ix3_WO)))
  ix4_NW = Int(floor(nx4/(ix4_WL-ix4_WO)))

if (ot + dt*(it_NW-1)*(it_WL-it_WO) + dt*it_WL < tmax)
	it_NW += 1
end
if (min_ix1 + (ix1_NW-1)*(ix1_WL-ix1_WO) + ix1_WL < max_ix1)
	ix1_NW += 1
end
if (min_ix2 + (ix2_NW-1)*(ix2_WL-ix2_WO) + ix2_WL < max_ix2)
	ix2_NW += 1
end
if (min_ix3 + (ix3_NW-1)*(ix3_WL-ix3_WO) + ix3_WL < max_ix3)
	ix3_NW += 1
end
if (min_ix4 + (ix4_NW-1)*(ix4_WL-ix4_WO) + ix4_WL < max_ix4)
	ix4_NW += 1
end

println("Patches in each dimension ",it_NW," ",ix1_NW," ",ix2_NW," ",ix3_NW," ",ix4_NW)

#split cube into multipatches
npatch=0

NW=it_NW*ix1_NW*ix2_NW*ix3_NW*ix4_NW
println("Number of patches=",NW)

patch_list = Patch[]
patch_names = String[]
nmbr = 1
for it_W = 1 : it_NW
	mint = ot + dt*(it_W-1)*(it_WL-it_WO)

	maxt = mint + dt*(it_WL - 1)

 	if (maxt >= tmax)
	    maxt = tmax
	end

	for ix1_W = 1 : ix1_NW
	    minx1 = min_ix1 + (ix1_W-1)*(ix1_WL-ix1_WO)
	    maxx1 = minx1 + ix1_WL - 1
	    if (maxx1 >= max_ix1)
		      maxx1 =  max_ix1
	    end
	    for ix2_W = 1 : ix2_NW
		      minx2 = min_ix2 + (ix2_W-1)*(ix2_WL-ix2_WO)
		      maxx2 = minx2 + ix2_WL - 1
		      if (maxx2 >= max_ix2)
		          maxx2 = max_ix2
		      end
		      for ix3_W = 1 : ix3_NW
		          minx3 = min_ix3 + (ix3_W-1)*(ix3_WL-ix3_WO)
		          maxx3 = minx3 + ix3_WL - 1
		          if (maxx3 >= max_ix3)
			             maxx3 = max_ix3
		          end
		          for ix4_W = 1 : ix4_NW
			             minx4 = min_ix4 + (ix4_W-1)*(ix4_WL-ix4_WO)
			             maxx4 = minx4 + ix4_WL - 1
			             if (maxx4 >= max_ix4)
			                  maxx4 = max_ix4
			             end
                         patch_name = join([out "_" nmbr])
			             minval=[mint minx1 minx2 minx3 minx4]
			             maxval=[maxt maxx1 maxx2 maxx3 maxx4]

			             patch_names = push!(patch_names,patch_name)
				         npatch += 1
                   		 nmbr +=1
						 it_nt=it_WL

				         push!(patch_list, Patch(in, patch_name, key, mint, maxt,minx1, maxx1, minx2, maxx2,minx3, maxx3, minx4, maxx4,it_nt))

		          end
		      end
	      end
	   end
   end

   pmap(grab_patch,patch_list)
   return patch_names,npatch

end

#*******************************************************************************
mutable struct Patch
     in
     name
     key
     mint
     maxt
     minx1
     maxx1
     minx2
     maxx2
     minx3
     maxx3
     minx4
     maxx4
     it_nt

end

#*******************************************************************************
function grab_patch(patch)
	 minval=[patch.mint patch.minx1 patch.minx2 patch.minx3 patch.minx4]
	 maxval=[patch.maxt patch.maxx1 patch.maxx2 patch.maxx3 patch.maxx4]
	 SeisWindowPatch(patch.in,patch.name,key=patch.key,minval=minval,maxval=maxval,it_nt=patch.it_nt)
end


text:
"""
    HalfInteger <: Real

Abstract supertype for half-integers.
"""
abstract type HalfInteger <: Real end

(T::Type{<:Rational})(x::HalfInteger)      = (t=twice(x); T(t, oftype(t,2)))
(T::Type{<:AbstractFloat})(x::HalfInteger) = T(float(x))
(T::Type{<:Integer})(x::HalfInteger)       = T(Integer(x))

Base.Bool(x::HalfInteger) = Bool(Integer(x))

function Base.Integer(x::HalfInteger)
    if isinteger(x)
        t = twice(x)
        t ÷ oftype(t,2)
    else
        throw(InexactError(:Integer, Integer, x))
    end
end

Base.:<(x::HalfInteger, y::HalfInteger)    = twice(x) < twice(y)
Base.:≤(x::HalfInteger, y::HalfInteger)    = twice(x) ≤ twice(y)
Base.:(==)(x::HalfInteger, y::HalfInteger) = twice(x) == twice(y)

# For correct hashing
Base.decompose(x::HalfInteger) = twice(x), 0, 2

Base.float(x::HalfInteger) = twice(x)/2

Base.isinteger(x::HalfInteger) = iseven(twice(x))

"""
    HalfInt <: HalfInteger

Number type for representing values `n/2` where `n` is an `Int`.

!!! warning

    Since the implementation of `HalfInt` is based on standard `Int` arithmetic, it is
    subject to integer overflow:

    ```jldoctest
    julia> x = typemax(HalfInt)
    $(typemax(Int))/2

    julia> x + 1
    -$(typemax(Int))/2

    julia> y = typemax(Int64)
    9223372036854775807

    julia> HalfInt(y) # 2*y == -2
    -1
    ```
"""
struct HalfInt <: HalfInteger
    twice::Int # the number 2x, where x is the value represented by the HalfInt object
    @doc """
        HalfInt(num, den)

    Create a `HalfInt` with the numerator `num` and denominator `den`. The denominator
    must be equal to `1` or `2`.

    # Examples

    ```jldoctest
    julia> HalfInt(5, 2)
    5/2

    julia> HalfInt(-3, 1)
    -3

    julia> HalfInt(1, 3)
    ERROR: DomainError with 3:
    denominator must be 1 or 2.
    Stacktrace:
    [...]
    ```
    """
    function HalfInt(num, den)
        den == 1 && return new(twice(num))
        den == 2 && return new(num)
        throw(DomainError(den, "denominator must be 1 or 2."))
    end
end

"""
    HalfInt(x)

Create a `HalfInt` equal to the number `x`.

# Examples

```jldoctest
julia> HalfInt(2)
2

julia> HalfInt(1.5)
3/2

julia> HalfInt(7//2)
7/2
```
"""
HalfInt(x)           = HalfInt(x, 1)
HalfInt(x::Rational) = HalfInt(numerator(x), denominator(x))
HalfInt(x::HalfInt)  = x

Base.:+(x::HalfInt, y::HalfInt) = HalfInt(twice(x)+twice(y), 2)
Base.:-(x::HalfInt, y::HalfInt) = HalfInt(twice(x)-twice(y), 2)
Base.:-(x::HalfInt)             = HalfInt(-twice(x), 2)

# Multiplication of two `HalfInt`s yields a Float64,
Base.:*(x::HalfInt, y::HalfInt) = twice(x)*twice(y)/4
# Multiplication of `HalfInt` and `Integer` yields a `HalfInt`
Base.:*(x::HalfInt, y::Integer) = HalfInt(twice(x)*y, 2)
Base.:*(x::Integer, y::HalfInt) = HalfInt(x*twice(y), 2)

Base.:/(x::HalfInt, y::HalfInt) = twice(x)/twice(y)

Base.://(x::HalfInt, y::HalfInt) = twice(x)//twice(y)
# Promotion doesn’t work for `Base.://`
Base.://(a::HalfInteger, b) = twice(a)//twice(b)
Base.://(a, b::HalfInteger) = twice(a)//twice(b)

Base.:^(x::Real, y::HalfInt) = x^float(y)

Base.gcd(x::HalfInt, y::HalfInt) = halfint_gcd(x,y)
Base.gcd(x::HalfInt, y::Integer) = halfint_gcd(x,y)
Base.gcd(x::Integer, y::HalfInt) = halfint_gcd(x,y)

Base.gcdx(x::HalfInt, y::HalfInt) = halfint_gcdx(x,y)
Base.gcdx(x::HalfInt, y::Integer) = halfint_gcdx(x,y)
Base.gcdx(x::Integer, y::HalfInt) = halfint_gcdx(x,y)

Base.lcm(x::HalfInt, y::HalfInt) = halfint_lcm(x,y)
Base.lcm(x::HalfInt, y::Integer) = halfint_lcm(x,y)
Base.lcm(x::Integer, y::HalfInt) = halfint_lcm(x,y)

# gcd where at least one of the arguments is a HalfInt
halfint_gcd(x, y) = HalfInt(gcd(twice(x), twice(y)), 2)
# lcm where at least one of the arguments is a HalfInt
halfint_lcm(x, y) = HalfInt(lcm(twice(x), twice(y)), 2)
# gcdx where at least one of the arguments is a HalfInt
function halfint_gcdx(x, y)
    d, u, v = gcdx(twice(x), twice(y))
    HalfInt(d, 2), u, v
end

Base.promote_rule(::Type{HalfInt}, T::Type)           = promote_type(Rational{Int}, T)
Base.promote_rule(::Type{HalfInt}, ::Type{<:Integer}) = HalfInt

Base.rem(x::HalfInt, y::HalfInt) = HalfInt(rem(twice(x),twice(y)), 2)

Base.round(x::HalfInt, ::RoundingMode{:Down}) = ifelse(isinteger(x), x, HalfInt(twice(x)-1, 2))

Base.show(io::IO, x::HalfInt) = isinteger(x) ? print(io, Integer(x)) : print(io, twice(x), "/2")

Base.typemax(::Type{HalfInt}) = HalfInt(typemax(Int), 2)
Base.typemin(::Type{HalfInt}) = HalfInt(typemin(Int), 2)

Base.intersect(r::StepRange{HalfInt}, s::StepRange{HalfInt})   = halfint_intersect(r, s)
Base.intersect(r::StepRange{HalfInt}, s::StepRange{<:Integer}) = halfint_intersect(r, s)
Base.intersect(r::StepRange{<:Integer}, s::StepRange{HalfInt}) = halfint_intersect(r, s)

Base.intersect(r::OrdinalRange{HalfInt}, s::OrdinalRange{HalfInt})   = halfint_intersect(r, s)
Base.intersect(r::OrdinalRange{HalfInt}, s::OrdinalRange{<:Integer}) = halfint_intersect(r, s)
Base.intersect(r::OrdinalRange{<:Integer}, s::OrdinalRange{HalfInt}) = halfint_intersect(r, s)

# Intersection of two ranges where the eltype of at least on of the ranges is HalfInt
function halfint_intersect(r::AbstractUnitRange, s::AbstractUnitRange)
    ifelse(isinteger(first(r)) ⊻ isinteger(first(s)),
           one(HalfInt):zero(HalfInt),
           max(first(r),first(s)):min(last(r),last(s)))
end

function halfint_intersect(r::OrdinalRange, s::OrdinalRange)
    twice_r = twice(first(r)):twice(step(r)):twice(last(r))
    twice_s = twice(first(s)):twice(step(s)):twice(last(s))
    twice_range = intersect(twice_r, twice_s)
    HalfInt(first(twice_range),2):HalfInt(step(twice_range),2):HalfInt(last(twice_range),2)
end

"""
    twice(x)

Return a value equal to `2*x`. For `HalfInteger`s, this returns an integer type. (not
exported)

# Examples

```jldoctest
julia> twice(3)
6

julia> twice(1.5)
3.0

julia> twice(HalfInt(1.5))
3
```
"""
twice(x)          = x + x
twice(x::HalfInt) = x.twice


text:
using WordCloud
wc = wordcloud(
    processtext(open(pkgdir(WordCloud) * "/res/alice.txt"), stopwords=WordCloud.stopwords_en ∪ ["said"]),
    mask=loadmask(pkgdir(WordCloud) * "/res/alice_mask.png", color="#faeef8"),
    colors=:seaborn_dark,
    angles=(0, 90),
    density=0.5) |> generate!
println("results are saved to alice.png")
paint(wc, "alice.png", background=outline(wc.mask, color="purple", linewidth=4))
wc
#eval# runexample(:alice)
#md# ![alice](alice.png)


text:
using StatsBase
using Test

X = Float64[1 0; 2 1; 3 0; 4 1; 5 10]

x1 = X[:,1]
x2 = X[:,2]
y = [5, 3, 4, 2, 5]

# corspearman

@test corspearman(x1, y) ≈ -0.102597835208515
@test corspearman(x2, y) ≈ -0.081110710565381

@test corspearman(X, y) ≈ [-0.102597835208515, -0.081110710565381]
@test corspearman(y, X) ≈ [-0.102597835208515 -0.081110710565381]

c11 = corspearman(x1, x1)
c12 = corspearman(x1, x2)
c22 = corspearman(x2, x2)
@test c11 ≈ 1.0
@test c22 ≈ 1.0
@test corspearman(X, X) ≈ [c11 c12; c12 c22]
@test corspearman(X)    ≈ [c11 c12; c12 c22]

# corkendall

@test corkendall(x1, y) ≈ -0.105409255338946
@test corkendall(x2, y) ≈ -0.117851130197758

@test corkendall(X, y) ≈ [-0.105409255338946, -0.117851130197758]
@test corkendall(y, X) ≈ [-0.105409255338946 -0.117851130197758]

c11 = corkendall(x1, x1)
c12 = corkendall(x1, x2)
c22 = corkendall(x2, x2)

@test c11 ≈ 1.0
@test c22 ≈ 1.0
@test corkendall(X, X) ≈ [c11 c12; c12 c22]
@test corkendall(X)    ≈ [c11 c12; c12 c22]


text:
using RecyclePermanent
using BenchmarkTools
using JLD2
include("../models/bpModel.jl")

@load "data/bpData.jld" ys θ
model = makeBPModel(θ, ys)

N = 40 * model.n

rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@btime simpleEstimate($model, $rpio)
@btime biasedEstimate($model, $rpio)
@btime biasedEstimateDiscrete($model, $rpio)
@btime rpEstimate($model, $rpio)
@btime rpEstimateDiscrete($model, $rpio)

rpio = RPIO{model.Tx}(N, model.n, 1)
@btime simpleEstimate($model, $rpio)
@btime biasedEstimate($model, $rpio)
@btime biasedEstimateDiscrete($model, $rpio)
@btime rpEstimate($model, $rpio)
@btime rpEstimateDiscrete($model, $rpio)


text:
module TestGapLas

using Test
import TreeLas: GapLas
import GraphIdx

@testset "GapLas  : demo 3x7             " begin
    include("demo3x7.jl")
    graph = GraphIdx.Grid.GridGraph(size(y)...)
    λ = GraphIdx.Ones{Float64}()
    μ = GraphIdx.Ones{Float64}()
    mem = GapLas.GapMem(y, graph, λ)
    GapLas.gaplas!(mem, y, graph, λ, μ)
    @testset "tree edges" begin
        parent = copy(mem.tree.parent)
        tree_edges = Dict{NTuple{2,Int},Int}()
        GraphIdx.enumerate_edges(graph) do e::Int, u::Int, v::Int
            if parent[v] == u || parent[u] == v
                tree_edges[(u, v)] = e
            end
        end
        tree_sel = Dict(mem.egraph.edges[s] => s for s in @view mem.mst.selected[2:end])
        @test tree_sel == tree_edges
        # be able to compare to c++ implementation
        @test join(
            (
                "{{$(e[1] -1), $(e[2] -1)}, $a}"
                for (e, a) in zip(mem.egraph.edges, round.(mem.sol.α, digits=3))
            ),
            ",\n"
        ) == """
        {{0, 1}, 0.0},
        {{3, 4}, 0.248},
        {{6, 7}, -0.273},
        {{9, 10}, 0.0},
        {{12, 13}, 0.0},
        {{15, 16}, -0.171},
        {{18, 19}, -0.238},
        {{1, 2}, 0.0},
        {{4, 5}, -0.505},
        {{7, 8}, -1.0},
        {{10, 11}, 0.0},
        {{13, 14}, -0.758},
        {{16, 17}, -0.266},
        {{19, 20}, 0.0},
        {{0, 3}, 0.179},
        {{3, 6}, 0.0},
        {{6, 9}, 0.405},
        {{9, 12}, -0.254},
        {{12, 15}, -0.582},
        {{15, 18}, 0.0},
        {{1, 4}, 0.179},
        {{4, 7}, 0.0},
        {{7, 10}, 0.618},
        {{10, 13}, 0.0},
        {{13, 16}, 0.0},
        {{16, 19}, -0.223},
        {{2, 5}, -0.101},
        {{5, 8}, -0.127},
        {{8, 11}, 0.151},
        {{11, 14}, 0.0},
        {{14, 17}, -0.997},
        {{17, 20}, -0.442}
        """ |> chomp
    end
    @test mem.dp_mem.proc_order == [
        12, 11, 19, 20, 21, 14, 15, 18, 17, 16, 13, 10, 7, 8, 9, 3, 6, 2, 5, 4, 1
    ]
    @test round.(mem.α_tree, digits=3) ≈ [
        -0.0, 0.179, -0.101, -0.179, -0.248, 0.505, -0.273,
        -1.0, 0.127, -0.405, -0.618, -0.151, 0.254, -0.758,
        -0.997, 0.582, 0.171, 0.266, -0.238, 0.223, 0.442
    ]
    @test round.(mem.sol.α, digits=3) ≈ [
        +0.0,   +0.248, -0.273, +0.0,   +0.0,   -0.171, -0.238, +0.0,
        -0.505, -1.0,   +0.0,   -0.758, -0.266, +0.0,   +0.179, +0.0,
        +0.405, -0.254, -0.582, +0.0,   +0.179, +0.0,   +0.618, +0.0,
        +0.0,   -0.223, -0.101, -0.127, +0.151, +0.0,   -0.997, -0.442
    ]
end

@testset "GapLas  : demo 3x7 iter 5      " begin
    include("demo3x7.jl")
    graph = GraphIdx.Grid.GridGraph(size(y)...)
    lam = GraphIdx.Const(0.1);
    GapLas.gaplas(y, graph, lam; max_iter=5, verbose=false)
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# Note: tests are sensitive to the line numbers of statements below
function foo()
    x = 5
    @showln x
    x = 7
    @showln x
end

function recurses(n)
    @showlnt n
    n += 1
    @showlnt n
    if n < 10
        n = recurses(n+1)
    end
    return n
end


text:
# ------------ Diff system ------------ #

struct ScalarDiff{K} <: Diff
    diff::K
end

struct BoolDiff <: Diff
    new::Bool
end

struct Change <: Diff end
struct Consumed <: Diff end

@inline tupletype(dfs::Diffed...) = Tuple{map(d -> valtype(d), dfs)...}

valtype(d::Diffed{V, DV}) where {V, DV} = V

# Diff propagation is basically a form of linear typing - Change types get "consumed" when they enter into a trace statement.

function propagate(args...)
    unwrapped = map(args) do a
        _lift(unwrap(a)) <: Change
    end
    any(unwrapped) && return Change
    return NoChange
end

function consume(args...)
    if any(args) do a
            a isa Change || _lift(a) <: Change
        end
        Consumed
    else
        NoChange
    end
end

struct DiffInterpreter <: InterpretationContext end

include("lib/numeric.jl")
include("lib/distributions.jl")
include("lib/base.jl")


text:
module PlaneStrainSolver

function bulk_modulus(l, m)
    return l + 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct CylindricalSolver
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function CylindricalSolver(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::CylindricalSolver, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::CylindricalSolver)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if r ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function shell_radial_stress(A::CylindricalSolver, r)
    return (A.ls + 2A.ms) * (A.A1s - A.A2s / r^2) +
           A.ls * (A.A1s + A.A2s / r^2) - (A.ls + 2A.ms / 3) * A.theta0
end

function shell_circumferential_stress(A::CylindricalSolver, r)
    return A.ls * (A.A1s - A.A2s / r^2) +
           (A.ls + 2A.ms) * (A.A1s + A.A2s / r^2) -
           (A.ls + 2A.ms / 3) * A.theta0
end

function shell_out_of_plane_stress(A::CylindricalSolver)
    return 2 * A.ls * A.A1s - (A.ls + 2A.ms / 3) * A.theta0
end

function shell_stress(A::CylindricalSolver, r)

    srr = shell_radial_stress(A, r)
    stt = shell_circumferential_stress(A, r)
    s33 = shell_out_of_plane_stress(A)

    return [srr, stt, s33]
end

function shell_pressure(A::CylindricalSolver, r)
    stress = shell_stress(A, r)
    return pressure(stress)
end

function shell_radial_strain(A::CylindricalSolver, r)
    return A.A1s - A.A2s / r^2 - A.theta0 / 3
end

function shell_circumferential_strain(A::CylindricalSolver, r)
    return A.A1s + A.A2s / r^2 - A.theta0 / 3
end

function shell_out_of_plane_strain(A::CylindricalSolver)
    return -A.theta0 / 3
end

function shell_strain(A::CylindricalSolver, r)
    err = shell_radial_strain(A, r)
    ett = shell_circumferential_strain(A, r)
    ezz = shell_out_of_plane_strain(A)
    return [err, ett, ezz]
end

function shell_strain_energy(A::CylindricalSolver, r, V0)
    stress = shell_stress(A, r)
    strain = shell_strain(A, r)

    return V0 * 0.5 * (sum(stress .* strain))
end

function shell_dilatation(A::CylindricalSolver)
    return 2 * A.A1s - A.theta0
end

function shell_compression_work(A::CylindricalSolver, r, V0)
    ekk = shell_dilatation(A)
    srr = shell_radial_stress(A, r)
    return V0 * (1 + ekk) * srr
end

function shell_potential(A::CylindricalSolver, r, V0)
    strainenergy = shell_strain_energy(A, r, V0)
    compwork = shell_compression_work(A, r, V0)
    return strainenergy - compwork
end

###########################################################################
function core_in_plane_stress(A::CylindricalSolver)
    return (A.lc + 2A.mc) * A.A1c + A.lc * A.A1c
end

function core_in_plane_strain(A::CylindricalSolver)
    return A.A1c
end

function core_strain(A::CylindricalSolver)
    err = core_in_plane_strain(A)
    return [err, err, 0.0]
end

function core_out_of_plane_stress(A::CylindricalSolver)
    return 2 * A.lc * A.A1c
end

function core_stress(A::CylindricalSolver)
    s11 = core_in_plane_stress(A)
    s33 = core_out_of_plane_stress(A)
    return [s11, s11, s33]
end

function core_pressure(A::CylindricalSolver)
    stress = core_stress(A)
    return pressure(stress)
end

function core_strain_energy(A::CylindricalSolver, V0)
    srr = core_in_plane_stress(A)
    err = core_in_plane_strain(A)
    stt = srr
    ett = err
    return V0 * 0.5 * (srr * err + stt * ett)
end

function core_dilatation(A::CylindricalSolver)
    return 2 * A.A1c
end

function core_compression_work(A::CylindricalSolver, V0)
    srr = core_in_plane_stress(A)
    ekk = core_dilatation(A)
    return V0 * (1 + ekk) * srr
end

function core_potential(A::CylindricalSolver, V0)
    strainenergy = core_strain_energy(A, V0s)
    compwork = core_compression_work(A, V0)
    return strainenergy - compwork
end

###########################################################################

function interface_potential_difference(solver,V01,V02)
    cse = core_strain_energy(solver,V02)
    ccw = core_compression_work(solver,V02)

    sse = shell_strain_energy(solver,solver.inradius,V01)
    scw = shell_compression_work(solver,solver.inradius,V01)

    return (sse - scw) - (cse - ccw)
end

function pressure(stress)
    return -1.0 / 3.0 * sum(stress)
end

end


text:
using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")
SnoopCompile.@snoop ["--project"] abspath(@__DIR__, "compiles.csv") begin
    include(abspath(@__DIR__, "..", "run.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts
blacklist = []

data = SnoopCompile.read(abspath(@__DIR__, "compiles.csv"))
pc = SnoopCompile.parcel(reverse!(data[2]), blacklist=blacklist)
SnoopCompile.write(abspath(@__DIR__), pc, always=true)


text:
"""

Construct command line for chain id.

$(SIGNATURES)

### Required arguments
```julia
* `m::SampleModel`                     : SampleModel
* `id::Int`                            : Chain id
```
Not exported
"""
function cmdline(m::SampleModel, id)

    #= cmdline with default parameter values
        `./bernoulli3 num_threads=4
        sample num_chains=4
        num_samples=1000 num_warmup=1000 save_warmup=0
        thin=1
        adapt engaged=1 gamma=0.05 delta=0.8 kappa=0.75
        t0=10.0 init_buffer=75 term_buffer=50 window=25
        algorithm=hmc engine=nuts
        max_depth=10 metric=diag_e stepsize=1.0 stepsize_jitter=1.0
        random seed=-1
        init=bernoulli3_1.init.R
        id=1
        data file=bernoulli3_1.data.R
        output file=bernoulli3_samples_1.csv
        refresh=100`
    =#

    cmd = ``
    # Handle the model name field for unix and windows
    cmd = `$(m.exec_path)`

    if m.use_cpp_chains
        cmd = `$cmd num_threads=$(m.num_threads)`
        cmd = `$cmd sample num_chains=$(m.num_cpp_chains)`
    else
        #cmd = `$cmd sample num_chains=1`
        cmd = `$cmd sample`
    end

    cmd = `$cmd num_samples=$(m.num_samples) num_warmup=$(m.num_warmups)`

    if m.save_warmup
        cmd = `$cmd save_warmup=1`
    else
        cmd = `$cmd save_warmup=0`
    end

    # Common to all models
    cmd = `$cmd thin=$(m.thin)`

    # Adapt section
    if m.engaged
        cmd = `$cmd adapt engaged=1`
    else
        cmd = `$cmd adapt engaged=0`
    end
    cmd = `$cmd gamma=$(m.gamma) delta=$(m.delta) kappa=$(m.kappa)`
    cmd = `$cmd t0=$(m.t0) init_buffer=$(m.init_buffer)`
    cmd = `$cmd term_buffer=$(m.init_buffer) window=$(m.window)`

    # Algorithm section
    cmd = `$cmd algorithm=$(string(m.algorithm)) engine=$(string(m.engine))`
    if m.engine == :nuts
        cmd = `$cmd max_depth=$(m.max_depth)`
    elseif m.engine == :static
        cmd = `$cmd int_time=$(m.int_time)`
    end
    cmd = `$cmd metric=$(string(m.metric)) stepsize=$(m.stepsize)`
    cmd = `$cmd stepsize_jitter=$(m.stepsize_jitter)`

    cmd = `$cmd random seed=$(m.seed)`

    # Init file required?
    if length(m.init_file) > 0 && isfile(m.init_file[id])
      cmd = `$cmd init=$(m.init_file[id])`
    else
      cmd = `$cmd init=$(m.init_bound)`
    end

    cmd = `$cmd id=$(id)`

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file[id])
      cmd = `$cmd data file=$(m.data_file[id])`
    end

    # Output files
    cmd = `$cmd output`
    if length(m.sample_file[id]) > 0
      cmd = `$cmd file=$(m.sample_file[id])`
    end
    if length(m.diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(m.diagnostic_file[id])`
    end

    # Refresh rate
    cmd = `$cmd refresh=$(m.refresh)`

    cmd
  end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
# helpers from http://nbviewer.jupyter.org/url/www.maths.usyd.edu.au/u/olver/teaching/MATH3976/notes/04.ipynb
# by Sheehan Olver
using Printf
printred(x)=print("\x1b[1m\x1b[31m"*x*"\x1b[0m")
printgreen(x)=print("\x1b[7m\x1b[32m"*x*"\x1b[0m")
printblue(x)=print("\x1b[4m\x1b[34m"*x*"\x1b[0m")

function decode_Float64(x::Float64; name="")
    s = bitstring(x)
    sgn = parse(Int,s[1],base=2)
    e = parse(Int,s[2:12],base=2)
    f = parse(Int,s[13:64],base=2)
    ffrac = f/2^52

    ebias = 1023

    if length(name) > 0
        @printf("The computer representation of %s is:\n", name)
        print("  ")
        printred(s[1:1])
        printgreen(s[2:12])
        printblue(s[13:end])
        println()
        println("  which decodes to ")
    end

    println("  | sign |   exponent  |    mantissa    ")
    print("     ");printred(s[1:1]);print("     ")
    printgreen(s[2:12]);print("   ")
    printblue(s[13:end])
    println()

    fstr = (@sprintf( "%.17f", ffrac) )[3:end]

    if e == 0
        println( @sprintf("= (-1)^(%i) x 2^(%+4i)   x 0.%s = %s ", sgn, 1-ebias, fstr, string(x)) )
    elseif e == 2^11 - 1
        println( @sprintf("= (-1)^(%i) x 2^( Inf)   x 1.%s = %s ", sgn, fstr, string(x)) )
    else
        println( @sprintf("= (-1)^(%i) x 2^(%+5i)   x 1.%s = %s ", sgn, e-ebias, fstr, string(x)) )
    end
end

macro show_float(var)
    varname = string(var)
    :(decode_Float64($var;name=$varname))
end

##
@show_float(Float64(pi))
@show_float(Float64(1.0))


text:
using Distributed
using Test

@testset "Noiseless" begin
# In this test we check that the noiseless case gets to probability 1 in the
# expected time (up to tolerance depending on maxdt)
    N = 20
    maxdt = 0.01
    psi0 = 1. / sqrt(N) * ones(ComplexF64, N)
    @testset "Complete graph" begin
        res = spatialsearch(psi0, complete_graph_Ad;
                    maxdt=maxdt, # maximum allowed time for dt
                    gamma=1. / N,
                    posW=1,
                    time=2*sqrt(N),
                    mu=1.,
                    coupling=1.,
                    noiseStrength=0.,
                    noiseRealizations=1,
                    dysonOrder=4);

        @test res.p >= zero(res.p) # Positive probability
        @test res.pmax ≈ 1 atol = 1e-5 # Maximum is 1
        @test res.tmax ≈ pi*sqrt(N)/2  atol = maxdt # For the correct time

    end

    @testset "Star graph" begin
        @testset "Central target" begin
            res = spatialsearch(psi0, star_graph_Ad;
                        maxdt=maxdt, # maximum allowed time for dt
                        gamma=1. / N,
                        posW=1,
                        time=2*sqrt(N),
                        mu=1.,
                        coupling=1.,
                        noiseStrength=0.,
                        noiseRealizations=1,
                        dysonOrder=4);

            @test res.p >= zero(res.p) # Positive probability
            @test res.pmax ≈ 1 atol = 1e-5
            @test res.tmax ≈ pi*sqrt(N)/2 atol = maxdt
        end

        @testset "External target" begin
            res = spatialsearch(psi0, star_graph_Ad;
                        maxdt=maxdt, # maximum allowed time for dt
                        gamma=1.,
                        posW=2,
                        time=2*sqrt(N),
                        mu=1.,
                        coupling=1.,
                        noiseStrength=0.,
                        noiseRealizations=1,
                        dysonOrder=4);

            @test res.p >= zero(res.p) # Positive probability
            @test res.pmax ≈ 1  atol = 1 / N^2
            @test res.tmax ≈ pi*sqrt(N)/2  atol = maxdt
        end
    end
end

@testset "Noisy" begin
    N = 10
    maxdt = 0.01
    psi0 = 1/sqrt(N) * ones(ComplexF64, N)
    @testset "Star graph" begin
        @testset "Central target" begin
            res = spatialsearch(psi0, star_graph_Ad;
                        maxdt=maxdt, # maximum allowed time for dt
                        gamma=1. / N,
                        posW=1,
                        time=2*sqrt(N),
                        mu=0.01,
                        coupling=1.,
                        noiseStrength=1.,
                        noiseRealizations=1000,
                        dysonOrder=4);

            @test res.p >= zero(res.p) # Positive probability
            @test res.pmax ≈ 0.5 atol = 5e-2
        end
    end
end

@testset "Distributed" begin
    addprocs(2)
    @everywhere using QuantumSpatialSearch
    N = 10
    maxdt = 0.01
    psi0 = 1/sqrt(N) * ones(ComplexF64, N)
    @testset "Star graph" begin
        @testset "Central target" begin
            res = spatialsearch(psi0, star_graph_Ad;
                        maxdt=maxdt, # maximum allowed time for dt
                        gamma=1. / N,
                        posW=1,
                        time=2*sqrt(N),
                        mu=0.01,
                        coupling=1.,
                        noiseStrength=1.,
                        noiseRealizations=1000,
                        dysonOrder=4);

            @test res.p >= zero(res.p) # Positive probability
            @test res.pmax ≈ 0.5 atol = 5e-2
        end
    end
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
# ------------ Choice sites ------------ #

@inline function (ctx::SimulateContext)(call::typeof(trace),
                                        addr::T,
                                        d::Distribution{K}) where {T <: Address, K}
    visit!(ctx.visited, addr)
    s = rand(d)
    add_choice!(ctx, addr, logpdf(d, s), s)
    return s
end

# ------------ Learnable ------------ #

@inline function (ctx::SimulateContext)(model::typeof(learnable), addr::T) where T <: Address
    visit!(ctx, addr)
    haskey(ctx.params, addr) && return getindex(ctx.params, addr)
    error("Parameter not provided at address $addr.")
end

# ------------ Black box call sites ------------ #

@inline function (ctx::SimulateContext)(c::typeof(trace),
                                        addr::T,
                                        call::Function,
                                        args...) where T <: Address
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ret, cl = simulate(ps, call, args...)
    add_call!(ctx, addr, cl)
    return ret
end

@inline function (ctx::SimulateContext)(c::typeof(trace),
                                        addr::T,
                                        call::G,
                                        args...) where {G <: GenerativeFunction,
                                                        T <: Address}
    visit!(ctx, addr)
    tr = simulate(call, args)
    ret = get_retval(tr)
    add_call!(ctx, addr, DynamicCallSite(get_choices(tr), get_score(tr), get_gen_fn(tr), get_args(tr), ret))
    ret
end

# ------------ Convenience ------------ #

function simulate(opt::J, params::P, model::Function, args...) where {J <: CompilationOptions, P <: AddressMap}
    ctx = Simulate(opt, DynamicTrace(), params)
    ret = ctx(model, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, model, args, ret)
end

function simulate(params::P, model::Function, args...) where P <: AddressMap
    ctx = Simulate(DefaultPipeline(), DynamicTrace(), params)
    ret = ctx(model, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, model, args, ret)
end

function simulate(model::Function, args...)
    ctx = Simulate(DefaultPipeline(), DynamicTrace(), Empty())
    ret = ctx(model, args...)
    return ret, DynamicCallSite(ctx.tr, ctx.score, model, args, ret)
end

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct SimulateContext{T <: AddressMap, P <: AddressMap} <: ExecutionContext
    tr::T
    visited::Visitor
    params::P
    SimulateContext(params) where T <: AddressMap = new{T}(AddressMap(), Visitor(), params)
end
```

`SimulateContext` is used to simulate traces without recording likelihood weights. `SimulateContext` can be instantiated with custom `AddressMap` instances, which is useful when used for gradient-based learning.

Inner constructors:
```julia
SimulateContext(params) = new{DynamicAddressMap}(AddressMap(), Visitor(), params)
```
""", SimulateContext)

@doc(
"""
```julia
ret, cl = simulate(fn::Function, args...; params = LearnableByAddress())
ret, cl = simulate(fn::typeof(trace), d::Distribution{T}; params = LearnableByAddress()) where T
```

`simulate` function provides an API to the `SimulateContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, and a `RecordSite` instance specialized to the call. `simulate` is used to express unconstrained generation of a probabilistic program trace, without likelihood weight recording.
""", simulate)


text:
using Test
using Nystrom
using StaticArrays
using Random
Random.seed!(1)

@testset "Interior representation" begin
    # test the interior representation formula u(r) = 𝒮[γ₁u](r) - 𝒟[γ₀u](r)
    @testset "2d" begin
        Geometry.clear_entities!()
        Ω    = ParametricSurfaces.Circle() |> Geometry.Domain
        Γ    = boundary(Ω)
        M    = ParametricSurfaces.meshgen(Γ,(10,))
        mesh = NystromMesh(view(M,Γ),order=5)
        x₀   = SVector(0.1,-0.2)
        xout = SVector(3,3)
        ops = (
            Laplace(;dim=2),
            Helmholtz(;dim=2,k=1.2),
            Elastostatic(;dim=2,μ=2,λ=3)
        )
        for pde in ops
            S    = SingleLayerPotential(pde,mesh)
            D    = DoubleLayerPotential(pde,mesh)
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xout,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xout,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            uₐ(x) = S[γ₁u](x) - D[γ₀u](x)
            @test isapprox(u(x₀),uₐ(x₀),rtol=1e-3)
        end
    end
    @testset "3d" begin
        Geometry.clear_entities!()
        Ω  = ParametricSurfaces.Sphere() |> Geometry.Domain
        Γ    = boundary(Ω)
        M    = ParametricSurfaces.meshgen(Γ,(4,4))
        mesh = NystromMesh(view(M,Γ),order=5)
        x₀   = SVector(0.1,-0.2,0.1)
        xout = SVector(3,3,3)
        ops = (
            Laplace(;dim=3),
            Helmholtz(;dim=3,k=1.2),
            Elastostatic(;dim=3,μ=2,λ=3),
            Maxwell(;k=1)
        )
        for pde in ops
            S    = SingleLayerPotential(pde,mesh)
            D    = DoubleLayerPotential(pde,mesh)
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xout,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xout,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            uₐ(x) = S[γ₁u](x) - D[γ₀u](x)
            @test isapprox(u(x₀),uₐ(x₀),rtol=1e-3)
        end
    end
end

@testset "Exterior representation" begin
    # test the exterior representation formula -u(r) = 𝒮[γ₁u](r) - 𝒟[γ₀u](r)
    @testset "2d" begin
        Geometry.clear_entities!()
        Ω  = ParametricSurfaces.Circle() |> Geometry.Domain
        Γ    = boundary(Ω)
        M    = ParametricSurfaces.meshgen(Γ,(10,))
        mesh = NystromMesh(view(M,Γ),order=5)
        x₀   = SVector(3,3)
        xin = SVector(0.1,0.2)
        ops = (
            Laplace(;dim=2),
            Helmholtz(;dim=2,k=1.2),
            Elastostatic(;dim=2,μ=2,λ=3)
        )
        for pde in ops
            S    = SingleLayerPotential(pde,mesh)
            D    = DoubleLayerPotential(pde,mesh)
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xin,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xin,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            uₐ(x) = -S[γ₁u](x) + D[γ₀u](x)
            @test isapprox(u(x₀),uₐ(x₀),rtol=1e-3)
        end
    end
    @testset "3d" begin
        Geometry.clear_entities!()
        Ω  = ParametricSurfaces.Sphere() |> Geometry.Domain
        Γ    = boundary(Ω)
        M    = ParametricSurfaces.meshgen(Γ,(4,4))
        mesh = NystromMesh(view(M,Γ),order=5)
        x₀   = SVector(3,3,3)
        xin = SVector(0.1,-0.2,0)
        ops = (
            Laplace(;dim=3),
            Helmholtz(;dim=3,k=1.2),
            Elastostatic(;dim=3,μ=2,λ=3),
            Maxwell(;k=1)
        )
        for pde in ops
            S    = SingleLayerPotential(pde,mesh)
            D    = DoubleLayerPotential(pde,mesh)
            T    = Nystrom.default_density_eltype(pde)
            c    = rand(T)
            u    = (qnode) -> SingleLayerKernel(pde)(xin,qnode)*c
            dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xin,qnode))*c
            γ₀u   = Density(u,mesh)
            γ₁u   = Density(dudn,mesh)
            uₐ(x) = -S[γ₁u](x) + D[γ₀u](x)
            @test isapprox(u(x₀),uₐ(x₀),rtol=1e-3)
        end
    end
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
#
# Routines for modified quadrature weights originally written by Johan Helsing
#

using FastGaussQuadrature
using Vandermonde

@inline function add_logcorr_block!(B, ztrg, zsrc, nsrc, wsrc, scale, trans, alpha)
    numpoints = length(wsrc)
    glpoints, glweights = gausslegendre(numpoints)
    if trans==0
        glpoints_scaled = glpoints
    else
        glpoints_scaled = @. trans + scale*glpoints
    end
    #WfrakL = WfrakLcomp(trans, scale, glpoints)
    WfrakL = WfrakLcomp(glpoints_scaled, glpoints)
    # Assemble block (including diagonal)
    KL = Array{Float64}(2,2)
    tdist = Array{Float64}(numpoints)
    for i=1:numpoints
        # tdist = Distance in parametrization of src panel
        @. tdist = abs(glpoints_scaled[i] - glpoints[:])
        # Very empirical limit for when neighbors need special quad due to log kernel:
        # 0.2 for 32, 0.4 for 16
        # Too gives nearly sing error
        # Too large makes special quad inaccurate, though not very sensitive
        if (trans == 0) || (minimum(tdist) < 0.2*32/numpoints)
            for j=1:numpoints
                nb_wcorr = WfrakL[i, j]/glweights[j] - log(tdist[j])
                # i==j term on self panel omitted since GL(0)=0 for modified stresslet
                zi = ztrg[1:2, i]
                zj = zsrc[1:2, j]
                nj = nsrc[1:2, j]
                wj = wsrc[j]
                rij = zj.-zi
                dblkernel_log!(rij, nj, alpha, KL)
                block_2x2_add!(B, i, j, KL.*(nb_wcorr*wj))
            end
        end
    end
end

"""
The choice of input arguments trans=0 and scale=1 gives WL
for ri and rj on the same quadrature panel γp. The choice
trans=±2 and scale=1 gives WL for ri on a neighboring panel
γp±1, assuming it is equal in parameter length. The input
argument tfrak is a column vector whose entries contain the
canonical nodes ti.
"""
function WfrakLcomp(trans::Float64,
                    scale::Float64,
                    tfrak::Array{Float64,1})
    npt = length(tfrak)
    tt=trans+scale*tfrak
    Q=zeros(npt, npt)
    p=zeros(1,npt+1)
    for j=1:npt
        p[1]=log(abs((1-tt[j])/(1+tt[j])))
        for k=1:npt
            c=(1-(-1)^k)/k
            p[k+1]=p[k]*tt[j]+c
        end
        for i=1:2:npt-1
            Q[i,j]=log(abs(1-tt[j]^2))-p[i+1] # 1658880
            Q[i+1,j]=p[1]-p[i+2]
        end
        for i=1:npt
            Q[i,j]=Q[i,j]/i
        end
    end

    # Direct solution:
    # A=ones(npt, npt)
    # for k=2:npt
    #     A[k,:]=tfrak.*A[k-1,:]
    # end
    # WfrakL =(A\Q).'

    # Bjorck-Pereyra:
    WfrakL = zeros(npt, npt);
    for i=1:npt
        WfrakL[i,:] = pvand(tfrak, Q[:,i]);
    end

    return WfrakL
end

function WfrakLcomp(targets, sources)
    tt = targets
    tfrak = sources
    nsrc = length(sources)
    ntrg = length(targets)
    Q=zeros(nsrc, ntrg)
    p=zeros(1,nsrc+1)
    for j=1:ntrg
        p[1]=log(abs((1-tt[j])/(1+tt[j])))
        for k=1:nsrc
            c=(1-(-1)^k)/k
            p[k+1]=p[k]*tt[j]+c
        end
        for i=1:2:nsrc-1
            Q[i,j]=log(abs(1-tt[j]^2))-p[i+1] # 1658880
            Q[i+1,j]=p[1]-p[i+2]
        end
        for i=1:nsrc
            Q[i,j]=Q[i,j]/i
        end
    end
    # Direct solution:
    # A=ones(npt, npt)
    # for k=2:npt
    #     A[k,:]=tfrak.*A[k-1,:]
    # end
    # WfrakL =(A\Q).'

    # Bjorck-Pereyra:
    WfrakL = zeros(ntrg, nsrc);
    for i=1:ntrg
        WfrakL[i,:] = pvand(tfrak, Q[:,i]);
    end

    return WfrakL
end

"""
    This function relies on complex arithmetic and takes, as
    input, points and vec- tors in R2 represented as points in
    C. Otherwise the notation follows Section 6: input
    parameters za and zb correspond to tau(ta) and tau(tb); z is the
    target point z ∈ E;and zj, nuj,and zpwj are column vectors
    whose entries contain the points zj , the exterior unit
    normals nu at rj , and the weighted velocity function zpwj
    , j = 1,...,npt.
    """
function modifiedweights(za::Complex{Float64},
                         zb::Complex{Float64},
                         z::Complex{Float64},
                         zj::Array{Complex{Float64}}
                         )
    npt = length(zj)
    # Transform endpoints to [-1,1]
    dz = (zb-za)/2
    ztr = (z-(zb+za)/2)/dz
    zjtr = (zj-(zb+za)/2)/dz
    p = Array{Complex{Float64}}(npt+1)
    q = Array{Complex{Float64}}(npt)
    r = Array{Complex{Float64}}(npt)
    s = Array{Complex{Float64}}(npt)

    c = (1-(-1).^(1:npt))./(1:npt)
    p[1] = log(1-ztr) - log(-1-ztr)
    p1 = log(1-ztr) + log(-1-ztr)
    # Signs here depend on how we have defined
    # direction of curve and normals. I think.
    if imag(ztr)<0 && abs(real(ztr))<1
        p[1] += 2im*pi
        p1 -= 2im*pi
    end
    for k=1:npt
	p[k+1] = ztr*p[k] + c[k]
    end

    r[1] = -1/(1+ztr)-1/(1-ztr)
    s[1] = -1/2*(1/(1-ztr)^2 - 1/(-1-ztr)^2)
    for k=1:npt-1
	r[k+1] = ztr*r[k] + p[k]
        s[k+1] = ztr*s[k] + r[k]
    end

    q[1:2:npt-1] = p1 - p[2:2:npt]
    q[2:2:npt] = p[1] - p[3:2:npt+1]
    q = q./(1:npt)
    # Bjorck-Pereyra solve
    wL = pvand(zjtr, q)*dz
    wC = pvand(zjtr, p[1:npt])
    wQ = pvand(zjtr, r) / dz
    wT = pvand(zjtr, s) / dz^2

    # Direct solve:
    # A = Array{Complex{Float64}}(npt, npt)
    # A[1,:] = 1.0
    # for k=2:npt
    #     for j=1:npt
    #         A[k,j]=zjtr[j]*A[k-1,j] # vandermonde transposed
    #     end
    # end
    # wL = (A\q)*dz
    # wC = A\p[1:npt]
    # wQ = (A\r) / dz
    # Only do one solve
    # sol = A\[q p[1:npt] r];
    # wL = sol[:,1]*dz
    # wC = sol[:,2]
    # wQ = sol[:,3] / dz

    return (wL, wC, wQ, wT)
end

# Compute near eval weights for log, Cauchy and stresslet singularities
function near_eval_weights(za, zb, z, zj, nuj, zpwj)
    # Compute modified weights
    wL, wC, wQ = modifiedweights(za, zb, z, zj)
    # Create corrections
    return near_eval_weights(za, zb, z, zj, nuj, zpwj, wL, wC, wQ)
end

function near_eval_weights(za, zb, z, zj, nuj, zpwj, wL, wC, wQ)
    # Create corrections
    dz = (zb-za)/2
    nuj = -nuj # Helsing has n = -i*t
    ## Log
    wcorrL = @. imag(wL*conj(nuj)) + log(abs(dz))*abs(zpwj)
    ## Cauchy
    wcorrC = @. -imag(wC)
    ## Unmodified stresslet
    wcorrQ1 = -1/2*imag(wC)
    wcorrQ2 = @. 1/4im*conj(conj(nuj)^2*wC + conj(zj-z)*wQ)
    return wcorrL, wcorrC, wcorrQ1, wcorrQ2
end


text:
"""
    DiscreteDomain{T <: Number} <: AbstractDomain
An abstract supertype for discrete domains (set, range).
"""
abstract type DiscreteDomain{T <: Number} <: AbstractDomain end

"""
    SetDomain{T <: Number} <: DiscreteDomain{T}
Domain that stores discrete values as a set of (unordered) points.
"""
struct SetDomain{T <: Number} <: DiscreteDomain{T}
    domain::Set{T}
end
SetDomain(values) = SetDomain(Set(values))

"""
    RangeDomain
A discrete domain defined by a `range <: AbstractRange{Real}`. As ranges are immutable in Julia, changes in `RangeDomain` must use [`set_domain!`](@ref).
"""
struct RangeDomain{T <: Real, R <: AbstractRange{T}} <: DiscreteDomain{T}
    domain::R
end

"""
    domain(values)
    domain(range::R) where {T <: Real, R <: AbstractRange{T}}
Construct either a [`SetDomain`](@ref) or a [RangeDomain](@ref).
```julia
d1 = domain(1:5)
d2 = domain([53.69, 89.2, 0.12])
d3 = domain([2//3, 89//123])
d4 = domain(4.3)
```
"""
domain(values) = SetDomain(values)
domain(range::R) where {T <: Real, R <: AbstractRange{T}} = RangeDomain{T, R}(range)

"""
    Base.length(d::D) where D <: DiscreteDomain
Return the number of points in `d`.
"""
Base.length(d::D) where D <: DiscreteDomain = length(get_domain(d))

"""
    Base.rand(d::D) where D <: DiscreteDomain
Draw randomly a point in `d`.
"""
Base.rand(d::D) where D <: DiscreteDomain = rand(get_domain(d))

"""
    Base.in(value, d::D) where D <: DiscreteDomain
Return `true` if `value` is a point of `d`.
"""
Base.in(value, d::D) where D <: DiscreteDomain = value ∈ get_domain(d)

"""
    domain_size(d::D) where D <: DiscreteDomain
Return the maximum distance between two points in `d`.
"""
domain_size(d::D) where D <: DiscreteDomain = δ_extrema(get_domain(d))

"""
    add!(d::SetDomain, value)
Add `value` to the list of points in `d`.
"""
add!(d::SetDomain, value) = !(value ∈ d) && push!(get_domain(d), value)

"""
    Base.delete!(d::SetDomain, value)(d::SetDomain, value)
Delete `value` from the list of points in `d`.
"""
Base.delete!(d::SetDomain, value) = pop!(get_domain(d), value)

function merge_domains(rd₁::RangeDomain, rd₂::RangeDomain)
    d₁ = get_domain(rd₁)
    d₂ = get_domain(rd₂)
    if step(d₁) == step(d₂) == 1
        a₁, b₁ = extrema(d₁)
        a₂, b₂ = extrema(d₂)
        a₂ < a₁ && return merge_domains(rd₂, rd₁)
        b₂ ≤ b₁ && return d₁
        a₂ ≤ b₁ + 1 && return a₁:b₂
    end
    return merge_domains(domain(Set(d₁)), rd₂)
end
merge_domains(rd::RangeDomain, d::D) where {D <: DiscreteDomain} = merge_domains(d, rd)
function merge_domains(d₁::D1, d₂::D2) where {D1 <: DiscreteDomain, D2 <: DiscreteDomain}
    return union(get_domain(d₁), get_domain(d₂))
end

function intersect_domains(d₁::D1,d₂::D2
) where {D1 <: DiscreteDomain, D2 <: DiscreteDomain}
    return intersect(get_domain(d₁), get_domain(d₂))
end


text:
# Autogenerated wrapper script for normaliz_jll for powerpc64le-linux-gnu-cxx11
export libnormaliz, normaliz

using GMP_jll
using CompilerSupportLibraries_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libnormaliz`
const libnormaliz_splitpath = ["lib", "libnormaliz.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libnormaliz_path = ""

# libnormaliz-specific global declaration
# This will be filled out by __init__()
libnormaliz_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libnormaliz = "libnormaliz.so.3"

# Relative path to `normaliz`
const normaliz_splitpath = ["bin", "normaliz"]

# This will be filled out by __init__() for all products, as it must be done at runtime
normaliz_path = ""

# normaliz-specific global declaration
function normaliz(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(normaliz_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"normaliz")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (GMP_jll.PATH_list, CompilerSupportLibraries_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (GMP_jll.LIBPATH_list, CompilerSupportLibraries_jll.LIBPATH_list,))

    global libnormaliz_path = normpath(joinpath(artifact_dir, libnormaliz_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libnormaliz_handle = dlopen(libnormaliz_path)
    push!(LIBPATH_list, dirname(libnormaliz_path))

    global normaliz_path = normpath(joinpath(artifact_dir, normaliz_splitpath...))

    push!(PATH_list, dirname(normaliz_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
julia> sum(k -> 1/k^2, 1:1000)
1.643934566681559

julia> pi^2/6
1.6449340668482264

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode