text:
mutable struct CRHCCache{M<:AbstractMatrix, V<:AbstractVector, C, IN}
    momentmatrix::M       # nxm
    "For regression type"
    modelmatrix::M
    residuals::V
    "Scratch matrix"
    matscr::M
    "Factorization of X'X"
    crossx::C         # Factorization of X'*X
    "Meat"
    Shat::M
    "The cluster indices (*sorted*)"
    clustersindices::IN
    "The *sorted* cluster indicator"
    f::CategoricalArray
end

"""
    bysort(x, f)

Sort each element of `x` according to f (a categorical).

# Arguments
- `x` an iterable whose elements are arrays.
- `f::CategoricalArray` a categorical array defining the sorting order

# Returns
- `Tuple`: a tuple (xs, fs) containng the sorted element of x and f
"""
function bysort(x, f)
    issorted(f) && return x, f
    @assert all(map(x->isa(x,AbstractArray), x))
    @assert all(map(x->size(x,1)==length(f), x))
    idx = sortperm(f)
    rt = map(x->similar(x), x)
    for ix in eachindex(rt)
        xin = x[ix]
        xout= rt[ix]
        for j in 1:size(xin,2)
            for i in eachindex(idx)
                xout[i,j] = xin[idx[i],j]
            end
        end
    end
    return rt, f[idx]
end

categorize(a::AbstractArray) = categorical(a)
categorize(f::CategoricalArray) = f

function clustersintervals(f::CategoricalArray)
    s = CategoricalArrays.order(f.pool)[f.refs]
    ci = collect(searchsorted(s, j) for j in unique(s))
    return ci
end

function installcache(k::CRHC, X::AbstractMatrix{T}) where T
    f = categorize(k.cl)
    (X, ), sf = bysort((X,), f)
    ci = clustersintervals(sf)
    n, p = size(X)
    chol = cholesky(diagm(0=>repeat([one(T)], inner=[p])))
    em = similar(X)
    ev = T[]
    Shat= Matrix{T}(undef,p,p)
    CRHCCache(X, em, ev, Matrix{T}(undef, 0,0), chol, Shat, ci, sf)
end

"""
    clustersindices(c::CRHCCache)

Return an array whose element `i` is a `Range{Int}` with indeces of the i-th cluster. Since
the data is sorted when cached, the indices are contigous.
"""
clustersindices(c::CRHCCache) = c.clustersindices
StatsModels.modelmatrix(c::CRHCCache) = c.modelmatrix
momentmatrix(c::CRHCCache) = c.momentmatrix
resid(c::CRHCCache) = c.residuals
invcrossx(c::CRHCCache) = inv(crossx(c))
crossx(c::CRHCCache) = c.crossx
#clusters_categorical(c::CRHCCache) = c.clustersindices ## TODO: remove it (??)

"""
    dofadjustment(k::CRHC, ::CRHCCache)

Calculate the default degrees-of-freedom adjsutment for `CRHC`

# Arguments
- `k::CRHC`: cluster robust variance type
- `c::CRHCCache`: the `CRHCCache` from which to extract the information
# Return
- `Float`: the degrees-of-fredom adjustment
# Note: the adjustment is a multyplicative factor.
"""
function dofadjustment(k::CRHC0, c::CRHCCache)
    g = length(clustersindices(c))
    return g/(g-1)
end

function dofadjustment(k::CRHC1, c::CRHCCache)
    g, (n, p) = length(clustersindices(c)), size(modelmatrix(c))
    return (n-1)/(n-p) * g/(g-1)
end

dofadjustment(k::CRHC2, c::CRHCCache) = 1

function dofadjustment(k::CRHC3, c::CRHCCache)
     g, (n, p) = length(clustersindices(c)), size(modelmatrix(c))
    return g/(g-1)
end

adjustresid!(k::CRHC0, c::CRHCCache) = resid(c)
adjustresid!(k::CRHC1, c::CRHCCache) = resid(c)

function adjustresid!(v::CRHC2, c::CRHCCache)
    n, p = size(momentmatrix(c))
    X, u = modelmatrix(c), resid(c)
    invxx, indices = invcrossx(c), clustersindices(c)
    Threads.@threads for index in indices
        Xv = view(X, index, :)
        uv = view(u, index, :)
        xAx = Xv*invxx*Xv'
        ldiv!(cholesky!(I - xAx; check=false).L, uv)
    end
    return u
end

function adjustresid!(k::CRHC3, c::CRHCCache)
    X, u = modelmatrix(c), resid(c)
    n, p = size(X)
    invxx, indices = invcrossx(c), clustersindices(c)
    Threads.@threads for index in indices
        Xv = view(X, index, :)
        uv = view(u, index, :)
        xAx = Xv*invxx*Xv'
        ldiv!(cholesky!(I - xAx; check=false), uv)
    end
    return rmul!(u, 1/sqrt(dofadjustment(k, c)))
end

Base.@propagate_inbounds function clusterize!(c::CRHCCache)
    U = momentmatrix(c)
    zeroM = zero(eltype(U))
    M = fill!(c.Shat, zeroM)
    p = size(M, 1)
    s = Array{eltype(U)}(undef, p)
    for m in clustersindices(c)
        fill!(s, zeroM)
        for j in 1:p
            for i in m
                s[j] += U[i, j]
            end
        end
        for j in 1:p
            for i in 1:j
                M[i, j] += s[i]*s[j]
            end
        end
    end
    return LinearAlgebra.copytri!(M, 'U')
end

renew(::CRHC0, id) = CRHC0(id, nothing)
renew(::CRHC1, id) = CRHC1(id, nothing)
renew(::CRHC2, id) = CRHC2(id, nothing)
renew(::CRHC3, id) = CRHC3(id, nothing)

# -----------------------------------------------------------------------------
# Very Very private api (subject to continuous change - DO NOT USE!!!)
# -----------------------------------------------------------------------------
Base.@propagate_inbounds function __vcov(k::CRHC, cache, df)
    B = inv(cache.crossx)
    res = adjustresid!(k, cache)
    cache.momentmatrix .= cache.modelmatrix.*res
    Shat = clusterize!(cache)
    return Symmetric((B*Shat*B).*df)
end


text:
DummyRegressor_ = SKDU.DummyRegressor
@sk_reg mutable struct DummyRegressor <: MMI.Deterministic
    strategy::String = "mean"::(_ in ("mean", "median", "quantile", "constant"))
    constant::Any     = nothing
    quantile::Float64 = 0.5::(0 ≤ _ ≤ 1)
end
MMI.fitted_params(m::DummyRegressor, (f, _, _)) = (
    constant  = f.constant_,
    n_outputs = f.n_outputs_
    )
metadata_model(DummyRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    weights = false,
    descr   = "DummyRegressor is a regressor that makes predictions using simple rules."
    )

# ----------------------------------------------------------------------------
DummyClassifier_ = SKDU.DummyClassifier
@sk_clf mutable struct DummyClassifier <: MMI.Probabilistic
    strategy::String  = "stratified"::(_ in ("stratified", "most_frequent", "prior", "uniform", "constant"))
    constant::Any     = nothing
    random_state::Any = nothing
end
MMI.fitted_params(m::DummyClassifier, (f, _, _)) = (
    classes   = f.classes_,
    n_classes = f.n_classes_,
    n_outputs = f.n_outputs_
    )
metadata_model(DummyClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "DummyClassifier is a classifier that makes predictions using simple rules."
    )

# ============================================================================
GaussianNBClassifier_ = SKNB.GaussianNB
@sk_clf mutable struct GaussianNBClassifier <: MMI.Probabilistic
    priors::Option{AbstractVector{Float64}} = nothing::(_ === nothing || all(_ .≥ 0))
    var_smoothing::Float64                  = 1e-9::(_ > 0)
end
MMI.fitted_params(m::GaussianNBClassifier, (f, _, _)) = (
    class_prior = f.class_prior_,
    class_count = f.class_count_,
    theta       = f.theta_,
    sigma       = f.sigma_,
    epsilon     = f.epsilon_,
    )
metadata_model(GaussianNBClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Gaussian naive bayes model."
    )

# ============================================================================
BernoulliNBClassifier_ = SKNB.BernoulliNB
@sk_clf mutable struct BernoulliNBClassifier <: MMI.Probabilistic
    alpha::Float64            = 1.0::(_ ≥ 0)
    binarize::Option{Float64} = 0.0
    fit_prior::Bool           = true
    class_prior::Option{AbstractVector} = nothing::(_ === nothing || all(_ .≥ 0))
end
MMI.fitted_params(m::BernoulliNBClassifier, (f, _, _)) = (
    class_log_prior  = f.class_log_prior_,
    feature_log_prob = f.feature_log_prob_,
    class_count      = f.class_count_,
    feature_count    = f.feature_count_
    )
metadata_model(BernoulliNBClassifier,
    input   = Table(Count),      # it expects binary but binarize takes care of that
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Binomial naive bayes classifier. It is suitable for classification with binary features; features will be binarized based on the `binarize` keyword (unless it's `nothing` in which case the features are assumed to be binary)."
    )

# ============================================================================
MultinomialNBClassifier_ = SKNB.MultinomialNB
@sk_clf mutable struct MultinomialNBClassifier <: MMI.Probabilistic
    alpha::Float64  = 1.0::(_ ≥ 0)
    fit_prior::Bool = true
    class_prior::Option{AbstractVector} = nothing::(_ === nothing || all(_ .≥ 0))
end
MMI.fitted_params(m::MultinomialNBClassifier, (f, _, _)) = (
    class_log_prior  = f.class_log_prior_,
    intercept        = f.intercept_,
    feature_log_prob = f.feature_log_prob_,
    coef             = f.coef_,
    class_count      = f.class_count_,
    feature_count    = f.feature_count_
    )
metadata_model(MultinomialNBClassifier,
    input   = Table(Count),        # NOTE: sklearn may also accept continuous (tf-idf)
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Multinomial naive bayes classifier. It is suitable for classification with discrete features (e.g. word counts for text classification)."
    )

# ============================================================================
ComplementNBClassifier_ = SKNB.ComplementNB
@sk_clf mutable struct ComplementNBClassifier <: MMI.Probabilistic
    alpha::Float64  = 1.0::(_ ≥ 0)
    fit_prior::Bool = true
    class_prior::Option{AbstractVector} = nothing::(_ === nothing || all(_ .≥ 0))
    norm::Bool      = false
end
MMI.fitted_params(m::ComplementNBClassifier, (f, _, _)) = (
    class_log_prior  = f.class_log_prior_,
    feature_log_prob = f.feature_log_prob_,
    class_count      = f.class_count_,
    feature_count    = f.feature_count_,
    feature_all      = f.feature_all_
    )
metadata_model(ComplementNBClassifier,
    input   = Table(Count),        # NOTE: sklearn may also accept continuous (tf-idf)
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Similar to Multinomial NB classifier but with more robust assumptions. Suited for imbalanced datasets."
    )

# ============================================================================
KNeighborsRegressor_ = SKNE.KNeighborsRegressor
@sk_reg mutable struct KNeighborsRegressor <: MMI.Deterministic
    n_neighbors::Int    = 5::(_ > 0)
    weights::Union{String,Function} = "uniform"::((_ isa Function) || _ in ("uniform", "distance"))
    algorithm::String   = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Int      = 30::(_ > 0)
    p::Int              = 2::(_ ≥ 0)
    metric::Any         = "minkowski"
    metric_params::Any  = nothing
    n_jobs::Option{Int} = nothing
end
MMI.fitted_params(m::KNeighborsRegressor, (f, _, _)) = (
    effective_metric        = f.effective_metric_,
    effective_metric_params = f.effective_metric_params_
    )
metadata_model(KNeighborsRegressor,
    input=Table(Continuous),
    target=AbstractVector{Continuous},
    weights=false,
    descr="K-Nearest Neighbors regressor: predicts the response associated with a new point by taking an average of the response of the K-nearest points."
    )

# ----------------------------------------------------------------------------
KNeighborsClassifier_ = SKNE.KNeighborsClassifier
@sk_clf mutable struct KNeighborsClassifier <: MMI.Probabilistic
    n_neighbors::Int    = 5::(_ > 0)
    weights::Union{String,Function} = "uniform"::((_ isa Function) || _ in ("uniform", "distance"))
    algorithm::String   = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Int      = 30::(_ > 0)
    p::Int              = 2::(_ ≥ 0)
    metric::Any         = "minkowski"
    metric_params::Any  = nothing
    n_jobs::Option{Int} = nothing
end
MMI.fitted_params(m::KNeighborsClassifier, (f, _, _)) = (
    classes                 = f.classes_,
    effective_metric        = f.effective_metric_,
    effective_metric_params = f.effective_metric_params_,
    outputs_2d              = f.outputs_2d_
    )
metadata_model(KNeighborsClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "K-Nearest Neighbors classifier: predicts the class associated with a new point by taking a vote over the classes of the K-nearest points."
    )

ext: jl


text:
doc"""
Multipole function

$f_{lm}(r,\theta,\phi) = u(r) Y_{lm}(\theta,\phi)$,

where $u(r)$ is a user-defined radial function.

If $r_{cut}<r<r_{max}$, then the long-range analytical
form $q Y_{lm}(\theta,\phi)/r^{l+1}$ is used, where
`q` is the multipole moment.

The function vanishes for $r>r_{max}$.

f = MultipoleFunc(u; l=0, m=0, R=(0.,0.,0.), rcut=Inf, rmax=Inf, q=0.)

`u` - a user-defined radial function

e.g.
```
u = r -> exp(-r^2)
u = Spline1D(r, exp(-r.^2))
```

`l`,`m` - angular quantum numbers

`R` - origin of spherical coordinate system
"""
type MultipoleFunc{T} <: ScalarFunc3d{Float64}
    u::T
    l::Int
    m::Int
    R::Point3{Float64}
    rcut::Float64
    rmax::Float64
    q::Float64
end

MultipoleFunc(u; l=0, m=0, R=(0.,0.,0.), rcut=Inf, rmax=Inf, q=0.) = MultipoleFunc{typeof(u)}(u,l,m,R,rcut,rmax,q)

MultipoleFunc(r::AbstractVector, u::AbstractVector; l=0, m=0, R=(0.,0.,0.), rcut=Inf,
        rmax=Inf, q=0., npts::Int=100) =
    MultipoleFunc(Spline1D(r,u,npts),l,m,R,rcut,rmax,q)

function Base.call(f::MultipoleFunc, x::Float64, y::Float64, z::Float64)
    r,θ,ϕ = cart2sph(x,y,z,f.R)
    r > f.rmax ? 0.0 : ( (r > f.rcut ? f.q/r^(f.l+1) : f.u(r)) * realsphharm(f.l,f.m,θ,ϕ))
end

Base.call(f::MultipoleFunc, p::Point3{Float64}) = f(p[1],p[2],p[3])

@vectorize2 Base.call MultipoleFunc Point3{Float64} Float64

doc"""
`enumeratelm(lmax)` lists all pairs `(l,m)` for `l=0,...,lmax` and `m=-l,...,l`.
"""
function enumeratelm(lmax)
    q = Point2{Int}[]
    for l in 0:lmax
        for m in -l:l
            push!(q, (l,m))
        end
    end
    q
end

typealias MultipoleFuncSpl MultipoleFunc{Spline1D{Float64,Vector{Float64}}}
typealias MultipoleBasis Vector{MultipoleFuncSpl}

function eval_basis_sum(b::MultipoleBasis, g::AbstractGrid3d, sphharmtab::AbstractArray, disttab::AbstractArray)
    n,m = length(g),length(b)
    φ = zeros(n)

    i = 1; R = b[1].R
    ind = Int[(f.R != R ? begin R = f.R; i += 1 end : i) for f in b]  # map basis function index -> center index

    for j in 1:m
        f = b[j]
        k = iq(f.l,f.m)
        R = f.R
        for i in 1:n
            r = disttab[i,ind[j]]
            if r <= f.rmax
                φ[i] += (r > f.rcut ? f.q/r^(f.l+1) : f.u(r)) * sphharmtab[i,ind[j],k]
            end
        end
    end

    φ
end

max_stars_repo_name: johncsnyder/Element


text:
### Binning
# A bin is a range of values between which values are grouped. Bins range from
# 1 to however many bins are selected. This allows for continuous data to be used
# for classification or histograms or other analysis.

SampleSpaceTooSmall = AssertionError("equalfrequencylimits(dv, bins) -> bins is not less than the length of dv. Reduce bins.")
BinsTooSmall = AssertionError("equalfrequencylimits(dv, bins) -> bins must be at least 2.")
BinWidthTooSmall = AssertionError("equalfrequencylimits(dv, bins) -> The calculated width of the bins is too small for sample space. Add more data or reduce bins.")

"equaldistancelimits(dv::Vector{<:AbstractFloat}, bins::Signed)

Produces a StepRangeLen of floats used to create bins, based on a spread
from the minimum to the maximum value of elements given."
function equaldistancelimits(dv::Vector{<:AbstractFloat}, bins::Signed)
    min = minimum(dv)
    max = maximum(dv)
    binwidth = (max-min)/(bins-1)
    return(min:binwidth:max)
end

"equalfrequencylimits(dv::Vector{<:AbstractFloat}, bins::Signed)

Produces a vector of floats used to create bins, based upon an equal frequency
division of the number of given elements after the range has been sorted.
If the number of elements is not evenly divided by the number of bins, the extra
elements are added to the last bin.

For example: 100 elements, 9 bins, produces 8 bins of 11 elements and 1 bin of 12."
function equalfrequencylimits(dv::Vector{<:AbstractFloat}, bins::Signed)
    if bins > length(dv)
        throw(SampleSpaceTooSmall)
    end
    if bins < 2
        throw(BinsTooSmall)
    end
    dv = sort(dv)
    elements = length(dv)
    binwidth = elements/bins

    if binwidth < 1.0
        throw(BinWidthTooSmall)
    end

    binrange = zeros(Float64, bins)
    for i in 1:bins
        binrange[i] = dv[Int(floor(i*binwidth))]
    end
    binrange[end] = dv[end]
    return sort(binrange)
end

"find_bin(x::AbstractFloat, limits::Vector{<:AbstractFloat})

This function finds which of the given bins the value belongs to."
function find_bin(x::AbstractFloat, limits::Vector{<:AbstractFloat})
    bin = Union{Missing, Int16}(missing)
    for i in 1:length(limits)
        if x <= limits[i]
            bin = i
            break
        end
    end
    bin
end

"find_bin(x::AbstractFloat, limits::StepRangeLen{<:AbstractFloat}=equaldistancelimits(x, 10))

This function finds which of the given bins the value belongs to."
function find_bin(x::AbstractFloat, limits::StepRangeLen{<:AbstractFloat}=equaldistancelimits(x, 10))
    bin = Union{Missing, Int16}(missing)
    for i in 1:length(limits)
        if x <= limits[i]
            bin = i
            break
        end
    end
    bin
end

"binindex(dv::Vector{<:AbstractFloat}, limits::Vector{<:AbstractFloat})

This produces a vector, the same length as the data, containing the bin values."
function binindex(dv::Vector{<:AbstractFloat}, limits::Vector{<:AbstractFloat})
    vbins = zeros(Int16, length(dv))
    for i in 1:length(dv)
        vbins[i] = find_bin(dv[i], limits)
    end
    vbins
end

"binindex(dv::Vector{<:AbstractFloat}, limits::StepRangeLen{<:AbstractFloat}=equaldistancelimits(dv, 10))

This produces a vector, the same length as the data, containing the bin values."
function binindex(dv::Vector{<:AbstractFloat}, limits::StepRangeLen{<:AbstractFloat}=equaldistancelimits(dv, 10))
    vbins = zeros(Int16, length(dv))
    for i in 1:length(dv)
        vbins[i] = find_bin(dv[i], limits)
    end
    vbins
end

"binsupport(dv::Vector{Int16}, limits::Vector{<:Signed})

This produces a vector of tuples that contains the bin and the support, or number of samples seen, for that bin."
function binsupport(dv::Vector{Int16}, limits::Vector{<:Signed})
    map(λ -> (λ, sum(dv .== λ)), limits)
end

"binsupport(dv::Vector{Int16}, limits::UnitRange{<:Signed})

This produces a vector of tuples that contains the bin and the support, or number of samples seen, for that bin."
function binsupport(dv::Vector{Int16}, limits::UnitRange{<:Signed})
    map(λ -> (λ, sum(dv .== λ)), limits)
end

# TODO: Missings support
# TODO: Group binning


text:
# Testing from REPL

using LSystem

function generate(model::LModel, n::Integer)
    state = LState(model)
    println(state)
    for i in 1:n
        state = next(state)
        println(state)
    end
end

algae_model = @lsys begin
    axiom : A
    rule  : A → AB
    rule  : B → A
end

generate(algae_model, 4)


text:
module LegacyEphemeris

using JLD

const STANDARD_EPHEMERIS = "430"
const PATH = "$(Pkg.dir())/JPLEphemeris/deps"

type Ephemeris
    id::Int
    fid::JLD.JldFile
    startepoch::Float64
    finalepoch::Float64
    dtable::Vector{Vector{Float64}}
    intervall::Vector{Float64}
    cache::Dict{AbstractString, Matrix{Float64}}
    constants::Dict{AbstractString, Float64}

    function Ephemeris(file::AbstractString)
        fid = jldopen(file, "r")
        startepoch = read(fid, "startepoch")
        finalepoch = read(fid, "finalepoch")
        dtable = read(fid, "dtable")
        intervall = read(fid, "intervall")
        cache = Dict{AbstractString, Matrix{Float64}}()
        sizehint!(cache, mapreduce(length, +, dtable))
        constants = read(fid, "constants")
        id = constants["DENUM"]
        merge!(constants, Dict("earthshare"=>1.0/(1.0+constants["EMRAT"]),
        "moonshare"=>constants["EMRAT"]/(1.0+constants["EMRAT"])))
        return new(id, fid, startepoch, finalepoch, dtable, intervall, cache, constants)
    end
end

function Ephemeris()
    return Ephemeris("$PATH/de$STANDARD_EPHEMERIS.jld")
end

function Ephemeris(denum::Integer)
    if isfile("$PATH/de$denum.jld")
        return Ephemeris("$PATH/de$denum.jld")
    else
        error("Ephemeris file '$PATH/de$denum.jld' not found.")
    end
end

const planets = Dict(
    "mercury"=>1,
    "venus"=>2,
    "earthmoon"=>3,
    "mars"=>4,
    "jupiter"=>5,
    "saturn"=>6,
    "uranus"=>7,
    "neptune"=>8,
    "pluto"=>9,
    "moon"=>10,
    "sun"=>11,
    "nutations"=>12,
    "librations"=>13,
)

function position(ephem::Ephemeris, body::AbstractString, date::Float64)
    i = planets[body]
    checkdate(ephem, date)
    c = coefficients(ephem, i, date)
    return position(c...)
end

function position(ephem::Ephemeris, body::AbstractString, date::AbstractVector{Float64})
    n = body == "nutations" ? 2 : 3
    p = zeros(n, length(date))
    for (i,d) in enumerate(date)
        p[:,i] = position(ephem, body, d)
    end
    return p
end

function velocity(ephem::Ephemeris, body::AbstractString, date::Float64)
    i = planets[body]
    checkdate(ephem, date)
    c = coefficients(ephem, i, date)
    return velocity(c...)
end

function velocity(ephem::Ephemeris, body::AbstractString, date::AbstractVector{Float64})
    n = body == "nutations" ? 2 : 3
    v = zeros(n, length(date))
    for (i,d) in enumerate(date)
        v[:,i] = velocity(ephem, body, d)
    end
    return v
end

function state(ephem::Ephemeris, body::AbstractString, date::Float64)
    nbody = planets[body]
    checkdate(ephem, date)
    c = coefficients(ephem, nbody, date)
    return [position(c...); velocity(c...)]
end

function state(ephem::Ephemeris, body::AbstractString, date::AbstractVector{Float64})
    n = body == "nutations" ? 4 : 6
    s = zeros(n, length(date))
    for (i,d) in enumerate(date)
        s[:,i] = state(ephem, body, d)
    end
    return s
end

function checkdate(ephem::Ephemeris, date::Float64)
    if (date < ephem.startepoch) | (date > ephem.finalepoch)
        error("The date must be between $(ephem.startepoch) and
        $(ephem.finalepoch).")
    end
end

function coefficients(ephem::Ephemeris, nbody::Int64, date::Float64)
    dt = ephem.intervall[nbody]
    if date == ephem.finalepoch
        frac = dt
        d = ephem.dtable[nbody][end]
    else
        index, frac = divrem(date-ephem.startepoch, dt)
        index = round(Int,(index) + 1)
        d = ephem.dtable[nbody][index]
    end

    # Check if the requested coeffcients are already cached,
    # if not retrieve them from the JLD file.
    key = "$nbody-$d"
    if ~haskey(ephem.cache, key)
        merge!(ephem.cache, Dict(key=>read(ephem.fid, key)))
    end
    c = ephem.cache[key]
    x = zeros(size(c)[1])

    # Normalized Chebyshev time
    tc = 2.0 * frac/dt - 1.0
    x[1] = 1
    x[2] = tc
    twotc = tc + tc
    for i = 3:length(x)
        x[i] = twotc*x[i-1] - x[i-2]
    end
    return c, x, dt, twotc
end

function position(c::Matrix{Float64}, x::Vector{Float64},
    dt::Float64, twotc::Float64)
    return c'*x
end

function velocity(c::Matrix{Float64}, x::Vector{Float64},
    dt::Float64, twotc::Float64)
    t = zeros(length(x))
    t[2] = 1.0
    t[3] = twotc + twotc
    for i = 4:length(t)
        t[i] = twotc*t[i-1] - t[i-2] + x[i-1] + x[i-1]
    end
    t *= 2.0
    t /= dt
    return c'*t
end

function close(ephem::Ephemeris)
    close(ephem.fid)
end

include("util.jl")

end # module


text:
using Plots, Distributions, LaTeXStrings

# Evaluates target distribution π(q)
function π_q(q)
    return 0.3*1/sqrt(2*pi)*exp(-(q+2)^2/2) + 0.7*1/sqrt(pi)*exp(-(q-3)^2)
end

# Evaluates conditional probability π(p|q)
function π_pq(p,q)
    return pdf(Normal(0, 1), p)
end

# Evaluates Kinetic enery K(p,q)
function K(p,q)
    return -log(π_pq(p,q))
end

# Evaluates Potential energy V(q)
function V(q)
    return -log(π_q(q))
end

# Evaluates dV/dq at q
function dV_dq(q)
    return 1/π_q(q) * (0.3*1/sqrt(2*pi)*(q+2)*exp(-(q+2)^2/2) +
    0.7*2/sqrt(pi)*(q-3)*exp(-(q-3)^2))
end

# Compute Hamiltonian at q,p
function H(p,q)
    return K(p,q) + V(q)
end

# Gets a trajectory of length L with n steps
function get_trajectory(q,p,L,n)
    # Flip sign of momentum
    # p *= -1

    # Step size
    ϵ = L/n

    q_trajectory = zeros(n)
    p_trajectory = zeros(n)

    q_trajectory[1] = q
    p_trajectory[1] = p

    for j=2:n
        ψ = p_trajectory[j-1] - ϵ/2 * dV_dq(q_trajectory[j-1])
        q_trajectory[j] = q_trajectory[j-1] + ϵ*ψ
        p_trajectory[j] = ψ - ϵ/2*dV_dq(q_trajectory[j])
    end

    return q_trajectory, p_trajectory
end

function hmc_alg(q₀, p₀)
    # Number of samples to return
    n = 20000

    q_samples = zeros(n)
    p_samples = zeros(n)
    q_samples[1] = q₀
    p_samples[1] = p₀

    trajectory_length = 1
    trajectory_partitions = 3

    for j=2:n

        q = q_samples[j-1]
        p = rand(Normal(0,1))
        q_proposal_trajectory, p_proposal_trajectory =
        get_trajectory(q, p, trajectory_length, trajectory_partitions)

        # Proposal
        q_L = last(q_proposal_trajectory)
        p_L = last(p_proposal_trajectory)

        # Compute acceptance probabilty
        α = min(1, exp(H(p,q) - H(-p_L,q_L)))

        if rand() < α
            q_samples[j] = q_L
            p_samples[j] = -p_L
        else
            q_samples[j] = q
            p_samples[j] = p
        end

    end

    return q_samples
end

# Generate HMC samples
samples = hmc_alg(0,1)

# Plot histogram of samples
x_vals = -5:0.01:5
y_vals = [π_q(q) for q in x_vals]
Plots.plot(x_vals, y_vals,
 bg = RGB(247/255, 236/255, 226/255),
 grid = false,
 legend = false,
 xlabel = L"q",
 ylabel = L"\pi(q)")
histogram!(samples, color = "blue", opacity = 0.15, normalize = true)

savefig("../images/hmc_ex_1_histogram.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
function mark_innerrectangle_perimetr!(r::Robot, side::HorizonSide)
    if side == Nord

        direction_of_movement=(Nord,Ost,Sud, West)
        direction_to_border=(Ost,Sud,West,Nord)
    else

        direction_of_movement=(Sud,Ost,Nord,West)
        direction_to_border=(Ost,Nord,West,Sud)
    end

    for i ∈ 1:4

        putmarkers!(r,  direction_of_movement[i], direction_to_border[i])
    end
end

function mark_innerrectangle_perimetr!(r::Robot, side::HorizonSide)
    direction_of_movement, direction_to_border = get_directions(side)
    for i ∈ 1:4
        putmarkers!(r, direction_of_movement[i], direction_to_border[i])
    end
end

get_directions(side::HorizonSide) = if (side == Nord)

        return (Nord,Ost,Sud, West), (Ost,Sud,West,Nord)
    else
        return (Sud,Ost,Nord,West), (Ost,Nord,West,Sud)
    end

function putmarkers!(r::Robot, direction_of_movement::HorizonSide, direction_to_border::HorizonSide)
    while isborder(r,direction_to_border)==true
        move!(r,direction_of_movement)
    end
end

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


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


text:
using OAuth

BOT_C_KEY = ENV["BOT_C_KEY"]
BOT_C_SEC = ENV["BOT_C_SEC"]
BOT_A_TOK = ENV["BOT_A_TOK"]
BOT_A_SEC = ENV["BOT_A_SEC"]

get_oauth(endpoint::String,
          options::Dict) = oauth_request_resource(endpoint,
                                                  "GET",
                                                  options,
                                                  BOT_C_KEY,
                                                  BOT_C_SEC,
                                                  BOT_A_TOK,
                                                  BOT_A_SEC)

post_oauth(endpoint::String,
           options::Dict) = oauth_request_resource(endpoint,
                                                   "POST",
                                                   options,
                                                   BOT_C_KEY,
                                                   BOT_C_SEC,
                                                   BOT_A_TOK,
                                                   BOT_A_SEC)

max_stars_repo_name: abhishalya/twitter-bot

max_issues_repo_name: abhishalya/twitter-bot

max_forks_repo_name: abhishalya/twitter-bot


text:
## The code in this file *does* depend on SparseMatrixCSC
## The algorithms could be factored out.

"""
    c_to_julia_index!(colptr, rowval, nzval)::SparseArrays.SparseMatrixCSC

Convert a sparse matrix with zero-based indices to a `SparseMatrixCSC`.
`colptr` and `rowval` are altered in place. `nzval` is neither copied
nor altered.
"""
function c_to_julia_index!(colptr, rowval, nzval)::SparseArrays.SparseMatrixCSC
    colptr .= colptr .+ 1
    rowval .= rowval .+ 1
    n = length(colptr)-1
    m = maximum(rowval)
    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

"""
    c_to_julia_index(colptr, rowval, nzval)::SparseArrays.SparseMatrixCSC

Convert a sparse matrix with zero-based indices to a `SparseMatrixCSC`.
`colptr` and `rowval` are not altered, but rather copied. `nzval` is neither copied, nor altered.
"""
function c_to_julia_index(colptr, rowval, nzval)::SparseArrays.SparseMatrixCSC
    c_to_julia_index!(deepcopy(colptr), deepcopy(rowval), nzval)
end

# We should use a doc like this.
# """
#     transpose(M::SparseMatrixCSC)

# Return the lazy transpose of `M`, which is a thin wrapper
# around `M`. To get the material
# transpose, which is another matrix of type `SparseMatrixCSC`,
# use `materialize(transpose(M))`.
# """

## FIXME: probably should remove this.
material_transpose(M::SparseMatrixCSC) = SparseArrays.ftranspose(M, identity)

"""
    nnz(sp::SparseMatrixCSC, colnum::Integer)

Return the number of structurally non-zero elements in the `colnum`th
column of `sp`.
"""
function SparseArrays.nnz(sp::SparseMatrixCSC, colnum::Integer)
    return length(SparseArrays.nzrange(sp, colnum))
end

"""
    nnzcounts(sp::SparseMatrixCSC; rev=true, byvalue=true)

Return a sorted countmap of the number of structural non-zeros
in the columns of `sp`. The keys controlling the sorting are
passed to `StatsBase.sort`.

# Examples
In this example, there are 132781 columns with one non-zero element,
32096 columns with two, etc.
```julia
julia> conns = nnzcounts(sp);

julia> collect(Iterators.take(conns,3))
3-element Array{Pair{Int64,Int64},1}:
 1 => 132781
 2 => 32096
 3 => 17478
```
"""
function nnzcounts(sp::SparseMatrixCSC; rev=true, byvalue=true)
    spcountmap = SparseUtils.countmap(nnzcols(sp), datatype=Int)
    sorted_spcountmap = StatsBase.sort(spcountmap, rev=rev, byvalue=rev)
    return sorted_spcountmap
end

# These may have been in SparseArrays, but were removed.
# They are very convenient at the command line.
nrows(sp::SparseMatrixCSC) = size(sp)[1]
ncols(sp::SparseMatrixCSC) = size(sp)[2]

"""
    nnzcols(sp::SparseMatrixCSC)

Return an iterator over the number of non-zero elements
in each column of `sp`.
"""
nnzcols(sp::SparseMatrixCSC) = (nnz(sp, i) for i in 1:ncols(sp))

"""
    prunecols!(sp::SparseMatrixCSC, min_connections)

For each column, set all elements of the column equal to zero if the number of non-zero
elements is less than `min_connections`. `sp` is modified in place. Here "zero" means
structural zeros.
"""
function prunecols!(sp::SparseMatrixCSC, min_connections)
    connection_flags = [nnz(sp, j) >= min_connections for j in 1:ncols(sp)] # This must not change during fkeep!
    SparseArrays.fkeep!(sp, (i, j, v) -> connection_flags[j], true)
    return sp
end

"""
    hasemptycols(sp::SparseMatrixCSC)

Return true if `sp` contains columns for which all elements are structural
zeros.

`prunecols!` will create such columns.
"""
hasemptycols(sp::SparseMatrixCSC) = length(sp.colptr) != length(unique(sp.colptr))

"""
    renumbercols(sp::SparseMatrixCSC)

If `sp` has "empty" columns as determined by `hasemptycols`, then remove these
columns and renumber the remaining columns. A new `SparseMatrixCSC` is returned.
Modifying elements of `.rowval` and `.nzval` in the output matrix will change
these values in `sp`.
"""
function renumbercols(sp::SparseMatrixCSC)
    (newcolptr, n) = _renumbercols(sp)
    SparseMatrixCSC(sp.m, n, newcolptr, sp.rowval, sp.nzval)
end

function _renumbercols(sp::SparseMatrixCSC)
    newcolptr = unique(sp.colptr)
    n = length(newcolptr) - 1
    return (newcolptr, n)
end

"""
    hasemptyrows(sp::SparseMatrixCSC)

Return true if `sp` contains rows for which all elements are structural
zeros.

`prunecols!` may create such rows.
"""
function hasemptyrows(sp::SparseMatrixCSC)
    rowlist = unique(sp.rowval)
    return length(rowlist) != maximum(rowlist)
end

"""
    renumberrows(sp::SparseMatrixCSC)

If `sp` has "empty" rows as determined by `hasemptycols`, then remove these
columns and renumber the remaining columns. A new `SparseMatrixCSC` is returned.
Modifying elements of `.rowval` and `.nzval` in the output matrix will change
these values in `sp`.
"""
function renumberrows(sp::SparseMatrixCSC)
    (renumbered_rowval, m) = _renumberrows(sp)
    SparseMatrixCSC(m, sp.n, sp.colptr, renumbered_rowval, sp.nzval)
end

function _renumberrows(sp::SparseMatrixCSC)
    sorted_rowval = sort!(unique(sp.rowval))
    rowval_type = eltype(sp.rowval)
    renumbered_rowval = rowval_type[searchsortedfirst(sorted_rowval, val) for val in sp.rowval]
    m = length(sorted_rowval)
    return (renumbered_rowval, m)
end

"""
    renumberrowscols(sp::SparseMatrixCSC)

Perform both `renumberrows` and `renumbercols` on `sp`
"""
function renumberrowscols(sp::SparseMatrixCSC)
    (newcolptr, n) = _renumbercols(sp)
    (renumbered_rowval, m) = _renumberrows(sp)
    return SparseMatrixCSC(m, n, newcolptr, renumbered_rowval, sp.nzval)
end


text:
using AdvancedMH
using Distributions
using StructArrays
using MCMCChains

using Random
using Test
using DiffResults
using ForwardDiff

include("util.jl")

@testset "AdvancedMH" begin
    # Set a seed
    Random.seed!(1234)

    # Generate a set of data from the posterior we want to estimate.
    data = rand(Normal(0, 1), 300)

    # Define the components of a basic model.
    insupport(θ) = θ[2] >= 0
    dist(θ) = Normal(θ[1], θ[2])
    density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf

    # Construct a DensityModel.
    model = DensityModel(density)

    @testset "StaticMH" begin
        # Set up our sampler with initial parameters.
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])
        spl2 = StaticMH(MvNormal([0.0, 0.0], 1))

        # Sample from the posterior.
        chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"])
        chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"])

        # chn_mean ≈ dist_mean atol=atol_v
        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
        @test mean(chain2.μ) ≈ 0.0 atol=0.1
        @test mean(chain2.σ) ≈ 1.0 atol=0.1
    end

    @testset "RandomWalk" begin
        # Set up our sampler with initial parameters.
        spl1 = RWMH([Normal(0,1), Normal(0, 1)])
        spl2 = RWMH(MvNormal([0.0, 0.0], 1))

        # Sample from the posterior.
        chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"])
        chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"])

        # chn_mean ≈ dist_mean atol=atol_v
        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
        @test mean(chain2.μ) ≈ 0.0 atol=0.1
        @test mean(chain2.σ) ≈ 1.0 atol=0.1
    end

    @testset "parallel sampling" begin
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])

        chain1 = sample(model, spl1, MCMCDistributed(), 10000, 4;
                        param_names=["μ", "σ"], chain_type=Chains)
        @test mean(chain1["μ"]) ≈ 0.0 atol=0.1
        @test mean(chain1["σ"]) ≈ 1.0 atol=0.1

        if VERSION >= v"1.3"
            chain2 = sample(model, spl1, MCMCThreads(), 10000, 4;
                            param_names=["μ", "σ"], chain_type=Chains)
            @test mean(chain2["μ"]) ≈ 0.0 atol=0.1
            @test mean(chain2["σ"]) ≈ 1.0 atol=0.1
        end
    end

    @testset "Proposal styles" begin
        m1 = DensityModel(x -> logpdf(Normal(x,1), 1.0))
        m2 = DensityModel(x -> logpdf(Normal(x[1], x[2]), 1.0))
        m3 = DensityModel(x -> logpdf(Normal(x.a, x.b), 1.0))
        m4 = DensityModel(x -> logpdf(Normal(x,1), 1.0))

        p1 = StaticProposal(Normal(0,1))
        p2 = StaticProposal([Normal(0,1), InverseGamma(2,3)])
        p3 = (a=StaticProposal(Normal(0,1)), b=StaticProposal(InverseGamma(2,3)))
        p4 = StaticProposal((x=1.0) -> Normal(x, 1))

        c1 = sample(m1, MetropolisHastings(p1), 100; chain_type=Vector{NamedTuple})
        c2 = sample(m2, MetropolisHastings(p2), 100; chain_type=Vector{NamedTuple})
        c3 = sample(m3, MetropolisHastings(p3), 100; chain_type=Vector{NamedTuple})
        c4 = sample(m4, MetropolisHastings(p4), 100; chain_type=Vector{NamedTuple})

        @test keys(c1[1]) == (:param_1, :lp)
        @test keys(c2[1]) == (:param_1, :param_2, :lp)
        @test keys(c3[1]) == (:a, :b, :lp)
        @test keys(c4[1]) == (:param_1, :lp)
    end

    @testset "Initial parameters" begin
        # Set up our sampler with initial parameters.
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])

        val = [0.4, 1.2]

        # Sample from the posterior.
        chain1 = sample(model, spl1, 10, init_params = val)

        @test chain1[1].params == val
    end

    @testset "is_symmetric_proposal" begin
        # True distributions
        d1 = Normal(5, .7)

        # Model definition.
        m1 = DensityModel(x -> logpdf(d1, x))

        # Set up the proposal (StandardNormal is a custom distribution in "util.jl").
        p1 = RandomWalkProposal(StandardNormal())

        # Implement `is_symmetric_proposal` for StandardNormal random walk proposal.
        AdvancedMH.is_symmetric_proposal(::RandomWalkProposal{<:StandardNormal}) = true

        # Make sure `is_symmetric_proposal` behaves correctly.
        @test AdvancedMH.is_symmetric_proposal(p1)

        # Sample from the posterior with initial parameters.
        chain1 = sample(m1, MetropolisHastings(p1), 100000;
                        chain_type=StructArray, param_names=["x"])

        @test mean(chain1.x) ≈ mean(d1) atol=0.05
        @test std(chain1.x) ≈ std(d1) atol=0.05
    end

    @testset "MALA" begin

        # Set up the sampler.
        sigma = 1e-1
        spl1 = MALA(x -> MvNormal((sigma^2 / 2) .* x, sigma))

        # Sample from the posterior with initial parameters.
        chain1 = sample(model, spl1, 100000; init_params=ones(2), chain_type=StructArray, param_names=["μ", "σ"])

        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
    end

    @testset "EMCEE" begin include("emcee.jl") end

end


text:
using LanguageServer
import LanguageServer.JSONRPC:Request, parse
import LanguageServer: process, lint
server = LanguageServerInstance(IOBuffer(), IOBuffer(), false)
init_request = """
{
    "jsonrpc":"2.0",
    "id":0,
    "method":"initialize",
    "params":{"processId":9902,
              "rootPath":null,
              "rootUri":"$(Pkg.dir("LanguageServer"))",
              "capabilities":{"workspace":{"applyEdit":true,"workspaceEdit":{"documentChanges":true},"didChangeConfiguration":{"dynamicRegistration":false},"didChangeWatchedFiles":{"dynamicRegistration":false},"symbol":{"dynamicRegistration":true},"executeCommand":{"dynamicRegistration":true}},"textDocument":{"synchronization":{"dynamicRegistration":true,"willSave":true,"willSaveWaitUntil":true,"didSave":true},"completion":{"dynamicRegistration":true,"completionItem":{"snippetSupport":true}},"hover":{"dynamicRegistration":true},"signatureHelp":{"dynamicRegistration":true},"references":{"dynamicRegistration":true},"documentHighlight":{"dynamicRegistration":true},"documentSymbol":{"dynamicRegistration":true},"formatting":{"dynamicRegistration":true},"rangeFormatting":{"dynamicRegistration":true},"onTypeFormatting":{"dynamicRegistration":true},"definition":{"dynamicRegistration":true},"codeAction":{"dynamicRegistration":true},"codeLens":{"dynamicRegistration":true},"documentLink":{"dynamicRegistration":true},"rename":{"dynamicRegistration":true}}},
              "trace":"off"}
}"""

process(LanguageServer.parse(Request,init_request), server)
process(LanguageServer.parse(Request, """{"jsonrpc":"2.0","method":"initialized","params":{}}"""), server)

# Workspace Symbols
r = parse(Request, """{"jsonrpc":"2.0","id":59,"method":"workspace/symbol","params":{"query":""}}""")
process(r, server);

# Document Symbols
for doc in values(server.documents)
    uri = doc._uri
    r = parse(Request, """{"jsonrpc":"2.0","id":1,"method":"textDocument/documentSymbol","params":{"textDocument":{"uri":"$(uri)"}}}""")
    process(r, server)
end

# Hovers

for doc in values(server.documents)
    uri = doc._uri
    print("Hovers: $uri ")
    for loc in 1:sizeof(doc._content)-1
        line, character = LanguageServer.get_position_at(doc, loc)
        r = parse(Request, """{"jsonrpc":"2.0","id":2,"method":"textDocument/hover","params":{"textDocument":{"uri":"$uri"},"position":{"line":$line,"character":$character}}}""")
        process(r, server)
        while !isempty(server.user_modules.data)
            take!(server.user_modules)
        end
    end
    println("($(sizeof(doc._content)))")
end

# Completions
for doc in values(server.documents)
    uri = doc._uri
    print("Completions: $uri ")
    for loc in 1:sizeof(doc._content)-1
        line, character = LanguageServer.get_position_at(doc, loc)
        r = parse(Request, """{"jsonrpc":"2.0","id":2,"method":"textDocument/completion","params":{"textDocument":{"uri":"$uri"},"position":{"line":$line,"character":$character}}}""")
        process(r, server)
        while !isempty(server.user_modules.data)
            take!(server.user_modules)
        end
    end
    println("($(sizeof(doc._content)))")
end

# Definitions
for doc in values(server.documents)
    uri = doc._uri
    print("Definitions: $uri ")
    for loc in 1:sizeof(doc._content)-1
        line, character = LanguageServer.get_position_at(doc, loc)
        r = parse(Request, """{"jsonrpc":"2.0","id":2,"method":"textDocument/definition","params":{"textDocument":{"uri":"$uri"},"position":{"line":$line,"character":$character}}}""")
        process(r, server)
        while !isempty(server.user_modules.data)
            take!(server.user_modules)
        end
    end
    println("($(sizeof(doc._content)))")
end

# Signatures
for doc in values(server.documents)
    uri = doc._uri
    print("Signatures: $uri ")
    for loc in 1:sizeof(doc._content)-1
        line, character = LanguageServer.get_position_at(doc, loc)
        r = parse(Request, """{"jsonrpc":"2.0","id":2,"method":"textDocument/signatureHelp","params":{"textDocument":{"uri":"$uri"},"position":{"line":$line,"character":$character}}}""")
        process(r, server)
        while !isempty(server.user_modules.data)
            take!(server.user_modules)
        end
    end
    println("($(sizeof(doc._content)))")
end


text:
using SpecialFunctions

mean_beta(a, b) = a / (a + b)
meanlog_beta(a, b) = digamma(a) - digamma(a + b)
mean_mirror(a, b) = mean_beta(b, a)
meanlog_mirror(a, b) = meanlog_beta(b, a)

function expected_logpi!(pi::Vector{Float64}, vm::VectorModel, w::Tw,
	min_prob::Float64=1e-3) where {Tw <: Integer}
	r = 0.
	x = 1.
	senses = 0
	pi[:] = vm.counts[:, w] #view(vm.counts, :, w)
	ts = sum(pi)
	for k in 1:T(vm)-1
		ts = max(ts - pi[k], 0.)
		a, b = 1. + pi[k] - vm.d, vm.alpha + k*vm.d + ts
		pi[k] = meanlog_beta(a, b) + r
		r += meanlog_mirror(a, b)

		pi_k = mean_beta(a, b) * x
		x = max(x - pi_k, 0.)
		if pi_k >= min_prob
			senses += 1
		end
	end
	pi[T(vm)] = r
	if x >= min_prob senses += 1 end
	return senses
end

function expected_pi!(pi::Vector{Float64}, vm::VectorModel,
		w::Tw, min_prob=1e-3) where {Tw <: Integer}
	r = 1.
	senses = 0
	ts = sum(view(vm.counts, :, w))
	for k in 1:T(vm)-1
		ts = max(ts - vm.counts[k, w], 0.)
		a, b = 1. + vm.counts[k, w] - vm.d, vm.alpha + k*vm.d + ts
		pi[k] = mean_beta(a, b) * r
		if pi[k] >= min_prob senses += 1 end
		r = max(r - pi[k], 0.)
	end
	pi[T(vm)] = r
	if r >= min_prob senses += 1 end
	return senses
end

function expected_pi(vm::VectorModel, w::Tw, min_prob=1e-3) where
	{Tw <: Integer}
	z = zeros(T(vm))
	expected_pi!(z, vm, w, min_prob)
	return z
end


text:
# SPDX-License-Identifier: MIT

struct NormalModel{T} <: AbstractModel
    pθ::T
end

#function loglikelihood(::NormalModel, x::AbstractVector, uTθ)
#    Λμ, μᵀΛμ, Λ, logdetΛ = uTθ
#    D = length(x)
#    xᵀΛx = dot(x, Λ * x)
#    μᵀΛx = dot(x, Λμ)
#    -(1/2)*(xᵀΛx + (μᵀΛμ - logdetΛ + D*log(2π))) + μᵀΛx
#end

function loglikelihood(::NormalModel, X::AbstractMatrix, uTθ)
    Λμ, μᵀΛμ, Λ, logdetΛ = uTθ
    D = size(X, 1)
    XᵀΛX = sum(X .* (Λ * X), dims=1)
    μᵀΛX = Λμ' * X
    -(1/2)*(XᵀΛX .+ (μᵀΛμ - logdetΛ + D*log(2π))) .+ μᵀΛX
end

#######################################################################
# GSM prior/posterior

#getstats(model::NormalModel, pθ::NormalWishart) = μ(pθ)
#unpack(model::NormalModel, pθ::NormalWishart, Tθ) = unpack(pθ, Tθ)
#
#function newposterior(::NormalModel, qθ::NormalWishart, ∇μ::AbstractVector;
#                      lrate=1)
#    T = typeof(qθ)
#    ηq = η(qθ)
#    ∇̃ξ = ForwardDiff.derivative(t -> ξ(qθ, ηq + t*∇μ), 0)
#    ξᵗ⁺¹ = ξ(qθ, ηq) + lrate*∇̃ξ
#    NormalWishart(unpack(qθ, ξᵗ⁺¹)...)
#end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
function correspondence(z,m₀,μ,Ω,graph)
    k = size(z,2)
    if k == 0
        return zeros(0)
    end
    n = (length(μ)-3)÷2
    qnoise = diagm( 0 => [5.0,0.02])
    co = zeros(Int,k)

    if n > 0
        t = μ[4:end]
        t₁ = t[1:2:end]
        t₂ = t[2:2:end]
        lst = findall((abs.(t₁.-μ[1]) .≤ 75) .*
                (abs.(t₂.-μ[2]) .≤ 75) .*
                ((t₁.-μ[1])*cos(μ[3]) .+ (t₂.-μ[2])*sin(μ[3]) .≥ -5))
        a = length(lst)
        jcb = jacobian_measurement(μ,lst)
        ẑ = reshape(equation_measurement(μ,lst),(2,1,a))

        q = SharedArray{Float64,2}((2a,2))

        @sync begin
            for p in procs(q)
                @async remotecall_wait(compute_cov_shared!,p,
                                        Ω,lst,m₀,graph,q,jcb,qnoise)
            end
        end

        f = reshape(permutedims(ẑ.-z,[1,3,2]),(2*a,k))
        f[2:2:end,:] = measure.(f[2:2:end,:])
        d₁ = zeros(2a,k)
        cc, ia = intersect_ind(lst,m₀,n)

        @inbounds for j ∈ ia
            d₁[2j-1:2j,:] = q[2j-1:2j,:]*f[2j-1:2j,:]
        end

        ind = [[2j-1,2j] for j ∈ ia]
        d₂ = d₁[vcat(ind...),:]
        d₂ .= d₂.^2
        d₃ = [d₂[i,j]+d₂[i+1,j] for i ∈ 1:2:2length(ia),j ∈ 1:k]

        co = argmin(d₃,dims=1)
        co = co[:]
        co = getindex.(co,1)
        m₁ = minimum(d₃,dims=1)[:]
        c₁ = cquantile(Chisq(2), .50)
        c₂ = cquantile(Chisq(2), .05)
        newl = findall(m₁ .> c₁)
        new₁ = findall(m₁[newl] .< c₂)
        ia = ia[setdiff(1:length(ia),co[newl[new₁]])]
        co = cc[co]

        if (length(newl) > 0 && length(cc) < a) || length(cc) == 0
                if length(cc) == 0
                        newl = 1:k
                end
                ib = setdiff(1:a,ia)
                lst = lst[ib]
                @inbounds for j ∈ ib
                        d₁[2j-1:2j,newl] = q[2j-1:2j,:]*f[2j-1:2j,newl]
                end
                if length(ib) > 0
                        ind = [[2j-1,2j] for j ∈ ib]
                        d₂ = d₁[vcat(ind...),newl]
                        d₂ = d₂.^2
                        d₃ = [d₂[i,j]+d₂[i+1,j] for i ∈ 1:2:2*length(ib),j ∈ 1:length(newl)]

                        co₁ = argmin(d₃,dims=1)
                        co₁ = co₁[:]
                        co[newl] = lst[getindex.(co₁,1)]
                        m₁ = minimum(d₃,dims=1)[:]
                        new₁ = []
                        newl = newl[m₁ .> c₂]
                end
        end
    else
        newl = 1:k
        new₁ = []
    end

    if length(new₁) > 0
            newl = newl[[i for i ∈ 1:length(newl) if i ∉ new₁]]
    end
    if length(newl) > 0
            co[newl] = n .+ [i for i ∈ 1:length(newl)]
    end
    return co
end

@everywhere include("markov_blanket.jl")

@everywhere function myrange(q)
        idx = indexpids(q)
        if idx == 0
            return 1:0
        end
        nchunks = length(procs(q))
        splits = [round(Int,s) for s ∈ range(0,stop=size(q,1)÷2,length=nchunks+1)]
        splits[idx]+1:splits[idx+1]
end

@everywhere function compute_cov!(Ω,lst,m₀,graph,q,jcb,qnoise,irange)
        @inbounds for i ∈ irange
            ind = markov_blanket(lst[i]+1,m₀.+1,graph)
            ind = ind[2:end] .- 1
            ind = [[2j+2,2j+3] for j ∈ ind ]
            ind = vcat(ind...)
            ind1 = findfirst(ind.==2*lst[i]+2)
            f = fxm(ind1,3+length(ind))
            j = jcb[:,:,i]
            #drastic improvements with symmetric matrices specific solver.
            s = j*f / Symmetric(Ω[[1:3;ind],[1:3;ind]]) *f'*j'
            vals,vec = eigen(Symmetric(s+qnoise))
            q[2i-1:2i,:] = diagm( 0 => sqrt.(1 ./ vals) )*vec'
        end
end

@everywhere function compute_cov_shared!(Ω,lst,m₀,graph,q,jcb,qnoise)
        compute_cov!(Ω,lst,m₀,graph,q,jcb,qnoise,myrange(q))
end

@everywhere function fxm(p,n)
        sparse(1:5,[1,2,3,3+p,4+p],1,5,n)
end

function measure(θ)
    tmp = θ%2π
    tmp-min(tmp÷π,1)*2π
end

function intersect_ind(a,b,n)
    #assume a is sorted vectors of positive integers
    u = zeros(Int,2,n)
    v = zeros(Int,n)
    u[1,a] .= 1
    u[2,a] = 1:length(a)
    v[b] .= 1
    ps = u[1,:] .* v
    res =  findall(ps .≠ 0)
    res, u[2,(ps .* u[2,:]) .> 0]
end


text:
"""
Phiefec(args...)

Calculate the Efective porosity by estimating the Geometric mean to Density and neutron Porosity

`` Phie=\\sqrt { \\frac { { RhoPhi }^{ 2 }+{ NeutronPhi }^{ 2 } }{ 2 }  } \\times \\left( 1-Vshale \\right)  ``
"""
function Phiefec(RhoPorosity,NeutronPorosity,Vshale)
    phi=(((RhoPorosity.^2 .+ NeutronPorosity.^2)./2)^0.5).*(1 .-Vshale)
end


text:
function dbeta(p::Real, a::Real, b::Real)
    any((a, b) .<= 0) ? (return NaN) : nothing
    !(0 <= p <= 1) ? (return 0.0) : nothing
    return gamma(a + b)/gamma(a)/gamma(b)*(1 - p)^(b - 1)*p^(a - 1)
end
function pbeta(p::Real, a::Real, b::Real)
    any((a, b) .<= 0) ? (return NaN) : nothing
    p <= 0  ? (return 0.0) : nothing
    p >= 1  ? (return 1.0) : nothing
    return beta_inc(a, b, p, 1 - p)[1]
end

gauss_legendre = function(low, high, order::Integer)
    pivots, weights = gauss(order)
    a, b = (high - low)/2, (high + low)/2
    scaled_pivots  = a .* pivots .+ b
    scaled_weights = a .* weights
    return(scaled_pivots, scaled_weights)
end

# precompute
gl_25_pivots, gl_25_weights = gauss_legendre(-1, 1, 25)
gauss_legendre_25 = function(low, high)
    pivots, weights = gl_25_pivots, gl_25_weights
    a, b = (high - low)/2, (high + low)/2
    scaled_pivots  = a .* pivots .+ b
    scaled_weights = a .* weights
    return(scaled_pivots, scaled_weights)
end

function one_stage_sample_size(p0::Real, α::Real, p1::Real, β::Real)
    z_1_α   = Distributions.quantile(Distributions.Normal(), 1 - α)
    z_1_β   = Distributions.quantile(Distributions.Normal(), 1 - β)
    napprox = p1*(1 - p1)*( (z_1_α + z_1_β) / (p1 - p0) )^2
    return Int(ceil(napprox))
end


text:
using Documenter, Stheno

DocMeta.setdocmeta!(
    Stheno,
    :DocTestSetup,
    :(using AbstractGPs, Stheno, Random, LinearAlgebra);
    recursive=true,
)

makedocs(
	modules = [Stheno],
    format = Documenter.HTML(),
    sitename = "Stheno.jl",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Input Types" => "input_types.md",
        "Kernel Design" => "kernel_design.md",
        "Internals" => "internals.md",
        "API" => "api.md",
    ],
)

deploydocs(repo="github.com/JuliaGaussianProcesses/Stheno.jl.git")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

"""
    get_beam_mass_matrix_2d(X1, X2, A, ro)

Function integrates mass matrix for 6 DOF Euler-Bernoulli beam element in 2d.

X1 = beams left node coordinates
X2 = beams right node coordinates
A = Cross section area
ro = Density
"""
function get_beam_mass_matrix_2d(X1,X2,A,ro)
    le=norm(X2-X1)                      # Lenght of element
    a=atan((X2[2]-X1[2])/(X2[1]-X1[1])) # Rotation angle of the element
    nn=2                            # Number of nodes
    nd=3*1+3                           # Number of DOFs
    Gp=[-1/3*sqrt(5+2*sqrt(10/7)) -1/3*sqrt(5-2*sqrt(10/7)) 0 1/3*sqrt(5-2*sqrt(10/7)) 1/3*sqrt(5+2*sqrt(10/7))] # Five Gauss integration points
    w=[(322-13*sqrt(70))/900 (322+13*sqrt(70))/900 128/225 (322+13*sqrt(70))/900 (322-13*sqrt(70))/900] # Five Gauss integration weights
    # Rotation matrix
    B=[cos(a) sin(a) 0  0     0       0;
      -sin(a) cos(a) 0  0     0       0;
       0      0      1  0     0       0;
       0      0      0  cos(a) sin(a) 0;
       0      0      0 -sin(a) cos(a) 0;
       0      0      0  0      0      1]
    # Integration of the truss mass matrix tm
    detJ=le/2
    function tmint(wi::Float64, xi::Float64)
        N1s=1+(-1/2)*(1+xi)
        N2s=(1/2)*(1 + xi)
        Ns=[N1s N2s]
        return wi*Ns'*Ns
    end
    tm=zeros(2,2)
    for i = 1:size(Gp,2)
        tm +=tmint(w[i],Gp[i])
    end
    tm=ro*A*tm*detJ
    # Integration of the beam mass matrix bm
    detJ=(le/2)
    function bmint(wi::Float64, xi::Float64)
        N1=1/4*(1-xi)^2*(2+xi)
        N2=le/8*(1-xi)^2*(xi+1)
        N3=1/4*(1+xi)^2*(2-xi)
        N4=le/8*(1+xi)^2*(xi-1)
        N=[N1 N2 N3 N4]
        return wi*N'*N
    end
    bm=zeros(4,4)
    for i = 1:size(Gp,2)
        bm +=bmint(w[i],Gp[i])
    end
    bm=ro*A*bm*detJ
    # Assembly of the 6 DOF truss-beam mass matrix m
    m=zeros(6,6)
    m[1,1]=tm[1,1];m[1,4]=tm[1,2];m[4,1]=tm[2,1];m[4,4]=tm[2,2]
    m[2,2]=bm[1,1];m[2,3]=bm[1,2];m[2,5]=bm[1,3];m[2,6]=bm[1,4]
    m[3,2]=bm[2,1];m[3,3]=bm[2,2];m[3,5]=bm[2,3];m[3,6]=bm[2,4]
    m[5,2]=bm[3,1];m[5,3]=bm[3,2];m[5,5]=bm[3,3];m[5,6]=bm[3,4]
    m[6,2]=bm[4,1];m[6,3]=bm[4,2];m[6,5]=bm[4,3];m[6,6]=bm[4,4]
    # Rotation
    m=B'*m*B
    return m
end


text:
module JLD2Writers

using JLD2

using ClimateMachine.Ocean.Domains: DiscontinuousSpectralElementGrid
using ClimateMachine.Ocean: current_step, current_time
using ..Fields: SpectralElementField

struct JLD2Writer{A, F, M, O}
    filepath::F
    model::M
    outputs::O
    array_type::A
end

function Base.show(io::IO, writer::JLD2Writer{A}) where {A}

    header = "JLD2Writer{$A}"
    filepath = "    ├── filepath: $(writer.filepath)"
    outputs = "    └── $(length(writer.outputs)) outputs: $(keys(ow.outputs))"

    print(io, header, '\n', filepath, '\n', outputs)

    return nothing
end

"""
    JLD2Writer(model, outputs=model.fields; filepath, array_type=Array, overwrite_existing=true)

Returns a utility for writing field output to JLD2 files, `overwrite_existing` file at `filepath`.

`write!(jld2_writer::JLD2Writer)` writes `outputs` to `filepath`, where `outputs` is either a `NamedTuple` or
`Dict`ionary of `fields` or functions of the form `output(model)`.

Field data is converted to `array_type` before outputting.
"""
function JLD2Writer(
    model,
    outputs = model.fields;
    filepath,
    array_type = Array,
    overwrite_existing = true,
)

    # Convert grid to CPU
    cpu_grid =
        DiscontinuousSpectralElementGrid(model.domain, array_type = array_type)

    # Initialize output
    overwrite_existing && isfile(filepath) && rm(filepath; force = true)

    file = jldopen(filepath, "a+")

    file["domain"] = model.domain
    file["grid"] = cpu_grid

    close(file)

    writer = JLD2Writer(filepath, model, outputs, array_type)

    write!(writer, first = true)

    return writer
end

initialize_output!(file, args...) = nothing

initialize_output!(file, field::SpectralElementField, name) =
    file["$name/meta/realelems"] = field.realelems

function write!(writer::JLD2Writer; first = false)

    model = writer.model
    filepath = writer.filepath
    outputs = writer.outputs

    # Add new data to file
    file = jldopen(filepath, "a+")

    N_output = first ? 0 : length(keys(file["times"]))

    step = current_step(model)
    time = current_time(model)

    file["times/$N_output"] = time
    file["steps/$N_output"] = step

    for (name, output) in zip(keys(outputs), values(outputs))
        first && initialize_output!(file, output, name)
        write_single_output!(file, output, name, N_output, writer)
    end

    close(file)

    return nothing
end

function write_field!(file, field, name, N_output, writer)
    data = convert(writer.array_type, field.data)
    file["$name/$N_output"] = data
    return nothing
end

write_single_output!(file, field::SpectralElementField, args...) =
    write_field!(file, field, args...)

function write_single_output!(file, output, name, N_output, writer)
    data = output(writer.model)
    file["$name/$N_output"] = data
    return nothing
end

struct OutputTimeSeries{F, N, D, G, T, S}
    filepath::F
    name::N
    domain::D
    grid::G
    times::T
    steps::S
end

function OutputTimeSeries(name, filepath)
    file = jldopen(filepath)

    domain, grid, times, steps = [nothing for i in 1:4]

    try
        domain = file["domain"]
        grid = file["grid"]

        output_indices = keys(file["times"])

        times = [file["times/$i"] for i in output_indices]
        steps = [file["steps/$i"] for i in output_indices]

    catch err
        @warn "Could not build time series of $name from $filepath because $(sprint(showerror, err))"

    finally
        close(file)

    end

    return OutputTimeSeries(filepath, name, domain, grid, times, steps)
end

function Base.length(timeseries::OutputTimeSeries)
    file = jldopen(timeseries.filepath)

    timeseries_length = length(keys(file["times"]))

    close(file)

    return timeseries_length
end

function Base.getindex(timeseries::OutputTimeSeries, i)

    name = timeseries.name
    domain = timeseries.domain
    grid = timeseries.grid

    file = jldopen(timeseries.filepath)

    data = file["$name/$(i-1)"]
    realelems = file["$name/meta/realelems"]

    close(file)

    realdata = view(data, :, realelems)

    return SpectralElementField(domain, grid, realdata, data, realelems)
end

end # module

max_forks_repo_name: jm-c/CLIMA


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Helper structs to dispatch known decision variable construction
struct NoSpecifiedConstraint <: MOI.AbstractSet end
Base.copy(set::NoSpecifiedConstraint) = set

struct DecisionSet <: MOI.AbstractScalarSet
    stage::Int
    constraint::Union{MOI.AbstractSet, JuMP.AbstractVectorSet}
    is_recourse::Bool

    function DecisionSet(stage::Integer; constraint::Union{MOI.AbstractSet, JuMP.AbstractVectorSet} = NoSpecifiedConstraint(), is_recourse::Bool = false)
        return new(stage, constraint, is_recourse)
    end
end

struct KnownSet <: MOI.AbstractScalarSet
    stage::Int
end

struct SingleDecisionSet{T} <: MOI.AbstractScalarSet
    stage::Int
    decision::Decision{T}
    constraint::MOI.AbstractSet
    is_recourse::Bool
end
Base.copy(set::SingleDecisionSet) = reuse(set, set.decision)

set_constraint(set::SingleDecisionSet) = set.constraint

struct MultipleDecisionSet{T} <: MOI.AbstractVectorSet
    stage::Int
    decisions::Vector{Decision{T}}
    constraint::MOI.AbstractSet
    is_recourse::Bool
end
MOI.dimension(set::MultipleDecisionSet) = length(set.decisions)
Base.copy(set::MultipleDecisionSet) = reuse(set, set.decisions)

MOIU.variable_function_type(::Type{<:SingleDecisionSet}) = SingleDecision
MOIU.variable_function_type(::Type{<:MultipleDecisionSet}) = VectorOfDecisions

is_decision_type(::Type{SingleDecisionSet}) = true
is_decision_type(::Type{MultipleDecisionSet}) = false

function JuMP.in_set_string(print_mode, set::SingleDecisionSet)
    decision_str(set) = set.is_recourse ? "RecourseDecisions" : "Decisions"
    if state(set.decision) == Taken
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))(value = $(set.decision.value))")
    end
    if state(set.decision) == Known
        return string(JuMP._math_symbol(print_mode, :in), " Known(value = $(set.decision.value))")
    end
    if set.constraint isa NoSpecifiedConstraint
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))")
    else
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function JuMP.in_set_string(print_mode, set::MultipleDecisionSet)
    decision_str(set) = set.is_recourse ? "RecourseDecisions" : "Decisions"
    if all(d -> state(d) == Taken, set.decisions)
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))(values = $([d.value for d in set.decisions]))")
    end
    if all(d -> state(d) == Known, set.decisions)
        return string(JuMP._math_symbol(print_mode, :in), " Known(values = $([d.value for d in set.decisions]))")
    end
    if set.constraint isa NoSpecifiedConstraint
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))")
    else
        return string(JuMP._math_symbol(print_mode, :in), " $(decision_str(set))($(JuMP.in_set_string(print_mode, set.constraint)))")
    end
end

function reuse(set::SingleDecisionSet, decision::Decision)
    return SingleDecisionSet(set.stage, decision, copy(set.constraint), set.is_recourse)
end

function reuse(set::MultipleDecisionSet, decisions::Vector{<:Decision})
    return MultipleDecisionSet(set.stage, decisions, copy(set.constraint), set.is_recourse)
end

function decision_set(variable::JuMP.ScalarVariable, set::DecisionSet)
    return SingleDecisionSet(set.stage, Decision(variable.info, Float64), set.constraint, set.is_recourse)
end

function decision_set(variables::Vector{<:JuMP.ScalarVariable}, set::DecisionSet)
    decisions = map(variables) do variable
        Decision(variable.info, Float64)
    end
    if set.constraint isa JuMP.AbstractVectorSet
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   JuMP.moi_set(set.constraint, length(variables)),
                                   set.is_recourse)
    else
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   set.constraint,
                                   set.is_recourse)
    end
end

function decision_set(variables::Matrix{<:JuMP.ScalarVariable}, set::DecisionSet)
    decisions = map(variables) do variable
        Decision(variable.info, Float64)
    end
    if set.constraint isa JuMP.AbstractVectorSet
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   JuMP.moi_set(set.constraint, length(variables)),
                                   set.is_recourse)
    else
        return MultipleDecisionSet(set.stage,
                                   decisions,
                                   set.constraint,
                                   set.is_recourse)
    end
end

function decision_set(variable::JuMP.ScalarVariable, set::KnownSet)
    return SingleDecisionSet(set.stage, KnownDecision(variable.info, Float64), NoSpecifiedConstraint(), false)
end

function decision_set(variables, set::KnownSet)
    decisions = map(variables) do variable
        KnownDecision(variable.info, Float64)
    end
    return MultipleDecisionSet(set.stage,
                               decisions,
                               NoSpecifiedConstraint(),
                               false)
end


text:
### solving a large system of linear equations demo.

using SigmoidNumbers

P32 = Posit{32,2}

#demonstrating what quires are.

q = Quire(P32)

add!(q, P32(1e20))

add!(q, P32(1e-20))

add!(q, P32(-1e20))

Float64(ans)

#demonstrating the power of quires.

A = random_exact_matrix(P32, 100)

y = exact_rowsum(A)

println((sum(Float64.(A),2)[1], Float64(y[1])))

println((Float64(sum(A,2)[1]), Float64(y[1])))

x = A \ y

sum(abs.(Float64.(x) - 1.0)) / 100.0

x = Float32.(A) \ Float32.(y)

sum(abs.(Float64.(x) - 1.0)) / 100.0

A[:, idx] = A[:, idx] * P32(2.0);

x1 = solve(A, y)

x2 = refine(x1, A, y)

#demonstrate that this is impossible with Float64

A = Float64.(A)
y = Float64.(y)

x1 = A \ y

x1[1]

x2 = x1 + A \ (y - A * x1)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using Tmp: f
using Test

@testset "Tmp.jl" begin
    @test f(2) == 4
    @test f(5) == 10
end


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_energy_market_at_fo" begin

    #=
        ECH = 10h
        TS: [11h, 11h15]
        S: [S1,S2]
                            bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=100    |                      |
        Csta=0, Cprop=1     |                      |
        S1: 20    S1: 15    |----------------------|
        S2: 30    S2: 30    |         35           |
                            |                      |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=10, Pmax=100   |                      | Pmin=10, Pmax=100
        Csta=45k, Cprop=10  |                      | Csta=80k, Cprop=15
    INIT: ON                |                      |INIT: ON
                            |                      | DP=DMO=2h
                            |                      |
            load(bus_1)     |                      |load(bus_2)
        S1: 10     S1: 17   |                      | S1: 40  S1: 48
        S2: 10     S2: 13   |                      | S2: 45  S2: 52

    prod_2_1's level was already decided for TS1 and TS2 : 10 and 15MW

    =#

    TS = [DateTime("2015-01-01T11:00:00"), DateTime("2015-01-01T11:15:00")]
    ech = DateTime("2015-01-01T10:00:00")
    context = PSCOPFFixtures.context_2buses_2TS_2S(TS, ech)
    context.tso_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T09:00:00"), SortedDict(
                                            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                                SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>20., "S2"=>30.)),
                                                            TS[2] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>15., "S2"=>30.))),
                                                ),
                                            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                                                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing,
                                                                                SortedDict("S1"=>PSCOPF.OFF, "S2"=>PSCOPF.OFF)),
                                                            TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(missing,
                                                                                SortedDict("S1"=>PSCOPF.OFF, "S2"=>PSCOPF.OFF))
                                                            ),
                                                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>20., "S2"=>15.)),
                                                            TS[2] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>35., "S2"=>20.))
                                                            )
                                                ),
                                            "prod_2_1" => PSCOPF.GeneratorSchedule("prod_2_1",
                                                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON,
                                                                                SortedDict("S1"=>PSCOPF.ON, "S2"=>PSCOPF.ON)),
                                                            TS[2] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON,
                                                                                SortedDict("S1"=>PSCOPF.ON, "S2"=>PSCOPF.ON))
                                                            ),
                                                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(10.,
                                                                                SortedDict("S1"=>10., "S2"=>10.)),
                                                            TS[2] => PSCOPF.UncertainValue{PSCOPF.Float64}(15.,
                                                                                SortedDict("S1"=>15., "S2"=>15.))
                                                            )
                                                )
                                            )
                                    )

    market = PSCOPF.EnergyMarketAtFO()
    next_ech = DateTime("2015-01-01T10:01:00") # does not matter : all decisions will be DECIDED or TODECIDE
    result, firmness = PSCOPF.run_step!(context, market, ech, next_ech)

    @testset "energy_market_at_fo_is_launched_at_fo" begin
        @test ech == (TS[1] - PSCOPF.get_fo_length(PSCOPF.get_management_mode(context)))
    end

    @testset "energy_market_required_decisions_are_all_firm" begin
        println(firmness)
        expected_firmness = PSCOPF.Firmness(
                        SortedDict("prod_1_1" => SortedDict(TS[1] => PSCOPF.TO_DECIDE,
                                                            TS[2] => PSCOPF.TO_DECIDE,
                                                            ),
                                    "prod_2_1" => SortedDict(TS[1] => PSCOPF.DECIDED,
                                                            TS[2] => PSCOPF.DECIDED,
                                                            ),
                                    ),
                        SortedDict("wind_1_1" => SortedDict(TS[1] => PSCOPF.TO_DECIDE,
                                                            TS[2] => PSCOPF.TO_DECIDE,
                                                            ),
                                    "prod_1_1" => SortedDict(TS[1] => PSCOPF.TO_DECIDE,
                                                            TS[2] => PSCOPF.TO_DECIDE,
                                                            ),
                                    "prod_2_1" => SortedDict(TS[1] => PSCOPF.DECIDED,
                                                            TS[2] => PSCOPF.DECIDED,
                                                            ),
                                    )
                        )
        @test expected_firmness == firmness
    end

    @testset "energy_market_at_fo_successful_launch" begin
        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
    end

    #=
    EnergyMarketAtFO returns firm decisions
    =#
    @testset "energy_market_at_fo_all_decisions_are_firm" begin
        @test context.market_schedule.decision_time == ech
        for (gen_id,gen_schedule) in context.market_schedule.generator_schedules
            if PSCOPF.Networks.get_type(PSCOPF.Networks.get_generator(PSCOPF.get_network(context), gen_id)) == PSCOPF.Networks.PILOTABLE
                for (ts, uncertain) in gen_schedule.production
                    @test PSCOPF.is_definitive(uncertain)
                end
            end
        end
        for (gen_id,gen_schedule) in context.market_schedule.generator_schedules
            for (ts, uncertain) in gen_schedule.commitment
                @test PSCOPF.is_definitive(uncertain)
            end
        end
    end

    #=
    EnergyMarketAtFO operates on an aggregated scenario
    =#
    @testset "energy_market_at_fo_aggregated_scenario" begin
        @test PSCOPF.SCENARIOS_DELIMITER == "_+_"

        expected_s_agg = "S1_+_S2"
        s_agg = PSCOPF.aggregate_scenario_name(context, ech)
        @test s_agg == expected_s_agg
    end

    #=
    EnergyMarketAtFO operates on an aggregated scenario and respects preceding decisions
    preceding decision : prod_2_1 started at 10MW for TS[1]

    The aggregate scenario :
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=100    |                      |
    Csta=0, Cprop=1     |                      |
        25      22.5    |----------------------|
                        |         35           |
                        |                      |
                        |                      |
    (pilotable) prod_1_1|                      |(pilotable) prod_2_1
    Pmin=10, Pmax=100   |                      | Pmin=10, Pmax=100
    Csta=45k, Cprop=10  |                      | Csta=80k, Cprop=15
    ON    ON       ON   |                      | ON  ON      ON
                        |                      |
           load(bus_1)  |                      |load(bus_2)
        10      15      |                      |    42.5  50

    =#
    @testset "energy_market_at_fo_production_levels" begin
        s_agg = PSCOPF.aggregate_scenario_name(context, ech)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        @test 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S2")
        @test 17.5 ≈ value(result.pilotable_model.p_injected["prod_1_1", TS[1], s_agg])
        @test 10. ≈ value(result.pilotable_model.p_injected["prod_2_1", TS[1], s_agg])
        @test value(result.limitable_model.p_global_capping[TS[1], s_agg]) < 1e-09
        @test value(result.lol_model.p_global_loss_of_load[TS[1], s_agg]) < 1e-09

        @test 22.5 ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[2], "S1")
        @test 22.5 ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[2], "S2")
        @test 27.5 ≈ value(result.pilotable_model.p_injected["prod_1_1", TS[2], s_agg])
        @test 15. ≈ value(result.pilotable_model.p_injected["prod_2_1", TS[2], s_agg])
        @test value(result.limitable_model.p_global_capping[TS[2], s_agg]) < 1e-09
        @test value(result.lol_model.p_global_loss_of_load[TS[2], s_agg]) < 1e-09

        @testset "energy_market_at_fo_respect_preceding_decisions" begin
            #prod_2_1 was already set to 10 for TS1
            @test PSCOPF.ON == PSCOPF.safeget_commitment_value(context.market_schedule, "prod_2_1", TS[1])
            @test 10. ≈ PSCOPF.safeget_prod_value(context.market_schedule, "prod_2_1", TS[1])
            #prod_2_1 was already set to 15 for TS2
            @test PSCOPF.ON == PSCOPF.safeget_commitment_value(context.market_schedule, "prod_2_1", TS[2])
            @test 15. ≈ PSCOPF.safeget_prod_value(context.market_schedule, "prod_2_1", TS[2])
        end
    end

    #=
    EnergyMarketAtFO considers aggregated scenario uncertainties,
    => It does not necessarily respect the EOD of the initial individual scenarios
    If we look at the initial scenarios we can have : overloads, overproductions or balanced scenarios
    =#
    @testset "energy_market_at_fo_does_not_necessarily_respect_EOD_on_initial_scenarios" begin
        s_agg = PSCOPF.aggregate_scenario_name(context, ech)

        #TS[1], "S1" : aggregated scenarios induces an over prodction
        @test PSCOPF.is_different(0. ,PSCOPF.compute_eod(PSCOPF.get_uncertainties(context),
                                                        PSCOPF.get_market_schedule(context),
                                                        PSCOPF.get_network(context),
                                                        ech, TS[1], "S1"))
        #TS[1], "S2" : aggregated scenarios induces an over load
        @test PSCOPF.is_different(0., PSCOPF.compute_eod(PSCOPF.get_uncertainties(context),
                                                        PSCOPF.get_market_schedule(context),
                                                        PSCOPF.get_network(context),
                                                        ech, TS[1], "S2") )
        #TS[2] : initial scenarios are balanced
        @test PSCOPF.compute_eod(PSCOPF.get_uncertainties(context),
                                PSCOPF.get_market_schedule(context),
                                PSCOPF.get_network(context),
                                ech, TS[2], "S1") < 1e-09
        @test PSCOPF.compute_eod(PSCOPF.get_uncertainties(context),
                                PSCOPF.get_market_schedule(context),
                                PSCOPF.get_network(context),
                                ech, TS[2], "S2") < 1e-09

        # capping and LoL are defined for the aggregated scenario
        @test value( result.lol_model.p_global_loss_of_load[TS[1], s_agg] ) < 1e-09
        @test value( result.limitable_model.p_global_capping[TS[1], s_agg] ) < 1e-09
        @test value( result.lol_model.p_global_loss_of_load[TS[2], s_agg] ) < 1e-09
        @test value( result.limitable_model.p_global_capping[TS[2], s_agg] ) < 1e-09
    end

    #=
    EnergyMarketAtFO satisfies the EOD of the aggregated scenario
    NOTE: If the EOD could not be satisfied, the market might use capping and LoL for the aggregated scenario
    =#
    @testset "energy_market_at_fo_respects_EOD_on_aggregated_uncertainties" begin
        s_agg, agg_uncertainties = PSCOPF.aggregate_scenarios(context, ech)
        agg_uncertainties = PSCOPF.Uncertainties(ech => agg_uncertainties)
        for ts in TS
            @test_skip 1e-09 > abs( PSCOPF.compute_eod(agg_uncertainties,
                                                PSCOPF.get_market_schedule(context),
                                                PSCOPF.get_network(context),
                                                ech, ts, s_agg) )
            @test value( result.lol_model.p_global_loss_of_load[ts, s_agg] ) < 1e-09
            @test value( result.limitable_model.p_global_capping[ts, s_agg] ) < 1e-09
        end
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
    abstract type PolynomialBasis{D,P}

A polynomial basis over `D` of order `≤P`.
"""
abstract type PolynomialBasis{D,P} end

"""
        struct MonomialBasis{D,P}

A basis of monomials spanning all polynomials of *order* `P` over `D`.

For multi-dimensional domains `D`, the precise definion of *order* depends in
fact on `D`. When `D` is a unit hypercube, the order of ``xᵅ`` where `α` is a
multi-index (e.g. an `NTuple`) is infinity norm. When `D` is a reference
simplex, the norm is the *L¹-norm*.
"""
struct MonomialBasis{D,P} <: PolynomialBasis{D,P} end
MonomialBasis(d::AbstractReferenceShape,p::Int) = MonomialBasis{typeof(d),p}()

Base.getindex(b::MonomialBasis,::Colon) = SVector{length(b)}(collect(b))

# monomial basis over reference segment
Base.length(::MonomialBasis{ReferenceLine,P}) where {P} = P+1
Base.eltype(::MonomialBasis{ReferenceLine}) = Monomial{1}
function Base.iterate(basis::MonomialBasis{ReferenceLine,P},state=0) where {P}
    if state > P
        return nothing
    else
        return Monomial(state),state+1
    end
end

# monomial basis over reference triangle
Base.length(::MonomialBasis{ReferenceTriangle,P}) where {P} = Int((P+1)*(P+2)/2)
Base.eltype(::MonomialBasis{ReferenceTriangle}) = Monomial{2}
function Base.iterate(basis::MonomialBasis{ReferenceTriangle,P},state=(0,0)) where {P}
    i,j = state
    if j > P
        return nothing
    else
        if i == P - j
            next = (0,j+1)
        else
            next = (i+1,j)
        end
        return Monomial(state),next
    end
end

# monomial basis over reference square
Base.length(::MonomialBasis{ReferenceSquare,P}) where {P} = (P)^2 + 1
Base.eltype(::MonomialBasis{ReferenceSquare}) = Monomial{2}
function Base.iterate(basis::MonomialBasis{ReferenceSquare,P},state=(0,0)) where {P}
    if state[1] > P || state[2] > P
        return nothing
    else
        if state[1]==state[2]
            return Monomial(state),state .+ (1,0)
        elseif state[1]>state[2]
            return Monomial(state),state .+ (-1,1)
        elseif state[2]>state[1]
            return Monomial(state),state .+ (1,0)
        end
    end
end

# vandermond matrix
function vandermond(xi,basis::PolynomialBasis)
    return [p(x) for x in xi, p in basis]
end

function lagrange_basis(xi,m::MonomialBasis)
    L = vandermond(xi,m)
    C = L\I
    return C
end

"""
    struct LagrangeBasis{D,P}

A Lagrange basis of order `P` over the reference shape
`D<:AbstractReferenceShape`. Spans the same space [`MonomialBasis{D,P}`](@ref)

Can be used as a function through `(b::LagrangeBasis)(x)`, where it Evaluate all
base elements in `b` at the point `x`. Returns a `StaticVector` of length
`length(b)`.
"""
struct LagrangeBasis{D,P} <: PolynomialBasis{D,P}
end

Base.getindex(b::LagrangeBasis,I) = b[:][I]

Base.length(::LagrangeBasis{ReferenceLine,P}) where {P} = P+1 |> Int
Base.length(::LagrangeBasis{ReferenceTriangle,P}) where {P} = (P+1)*(P+2)/2 |> Int
Base.length(::LagrangeBasis{ReferenceTetrahedron,P}) where {P} = (P+1)*(P+2)*(P+3)/6 |> Int

function (b::LagrangeBasis{D,0})(x::SVector{<:Any,<:Number}) where {D}
    SVector(1)
end
Base.getindex(b::LagrangeBasis{D,0},::Colon) where {D} = (x->SVector{1},)

function (b::LagrangeBasis{ReferenceLine,1})(x::SVector{<:Any,<:Number})
    map(f->f(x),b[:])
end
Base.getindex(b::LagrangeBasis{ReferenceLine,1},::Colon) = (x->1-x[1],x->x[1])

function (b::LagrangeBasis{ReferenceTriangle,1})(x::SVector{<:Any,<:Number})
    SVector(1 - x[1] - x[2], x[1], x[2])
end
(b::LagrangeBasis{ReferenceTriangle,1})()              = (x -> 1 - x[1] - x[2], x -> x[1], x -> x[2])
Base.getindex(b::LagrangeBasis{ReferenceTriangle,1},I) = b()[I]

function (b::LagrangeBasis{ReferenceTetrahedron,1})(x::SVector{<:Any,<:Number})
    SVector(1 - x[1] - x[2] - x[3], x[1], x[2], x[3])
end
Base.getindex(b::LagrangeBasis{ReferenceTetrahedron,1},::Colon) = (x -> 1 - x[1] - x[2] - x[3], x -> x[1], x -> x[2], x -> x[3])

function (b::LagrangeBasis{D,p})(x̂::SVector{<:Any,<:SVector}) where {D,p}
    mapreduce(x->b(x),hcat,x̂)
end

trace(::LagrangeBasis{ReferenceTriangle,p}) where {p} = LagrangeBasis{ReferenceLine,p}()
trace(::LagrangeBasis{ReferenceTetrahedron,p}) where {p} = LagrangeBasis{ReferenceTriangle,p}()

struct gradLagrangeBasis{D,P} <: PolynomialBasis{D,P}
end

Base.length(::gradLagrangeBasis{ReferenceLine,P}) where {P} = P+1 |> Int
Base.length(::gradLagrangeBasis{ReferenceTriangle,P}) where {P} = (P+1)*(P+2)/2 |> Int
Base.length(::gradLagrangeBasis{ReferenceTetrahedron,P}) where {P} = (P+1)*(P+2)*(P+3)/6 |> Int

function (b::gradLagrangeBasis{ReferenceLine,1})(x::SVector{<:Any,<:Number})
    SVector(SVector(-1),
            SVector(1))
end

function (b::gradLagrangeBasis{ReferenceTriangle,1})(x::SVector{<:Any,<:Number})
    SVector(SVector(-1, -1),
            SVector(1, 0),
            SVector(0, 1))
end
(b::gradLagrangeBasis{ReferenceTriangle,1})()               = ( x -> SVector(-1, -1),
                                                                x -> SVector(1, 0),
                                                                x-> SVector(0, 1)
                                                                )
Base.getindex(b::gradLagrangeBasis{ReferenceTriangle,1},I)  = b()[I]

function (b::gradLagrangeBasis{ReferenceTetrahedron,1})(x::SVector{<:Any,<:Number})
    SVector(SVector(-1, -1, -1),
            SVector(1, 0, 0),
            SVector(0, 1, 0),
            SVector(0, 0, 1))
end

function (b::gradLagrangeBasis{D,p})(x̂::SVector{<:Any,<:SVector}) where {D,p}
    mapreduce(x->b(x),hcat,x̂)
end

grad(::LagrangeBasis{D,p}) where {D,p} = gradLagrangeBasis{D,p}()

# FIXME: need a better way to related grad to functions in a basis
blist = (LagrangeBasis{ReferenceTriangle,1}(),)
for b in blist
    for i in 1:length(b)
        f = b[i]
        @eval grad(::typeof($f)) = grad($b)[$i]
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using OpenTelemetry
using Test

@testset "OpenTelemetry.jl" begin
    # Write your tests here.
end


text:
# ---
# title: 526. Beautiful Arrangement
# id: problem526
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/beautiful-arrangement/description/>
# hidden: true
# ---
#
# Suppose you have **N** integers from 1 to N. We define a beautiful arrangement
# as an array that is constructed by these **N** numbers successfully if one of
# the following is true for the ith position (1 <= i <= N) in this array:
#
#   1. The number at the ith position is divisible by **i**.
#   2. **i** is divisible by the number at the ith position.
#
#
#
# Now given N, how many beautiful arrangements can you construct?
#
# **Example 1:**
#
#
#
#     Input: 2
#     Output: 2
#     Explanation:
#
#     The first beautiful arrangement is [1, 2]:
#
#     Number at the 1st position (i=1) is 1, and 1 is divisible by i (i=1).
#
#     Number at the 2nd position (i=2) is 2, and 2 is divisible by i (i=2).
#
#     The second beautiful arrangement is [2, 1]:
#
#     Number at the 1st position (i=1) is 2, and 2 is divisible by i (i=1).
#
#     Number at the 2nd position (i=2) is 1, and i (i=2) is divisible by 1.
#
#
#
#
# **Note:**
#
#   1. **N** is a positive integer and will not exceed 15.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end