text:
using Printf,   StatsBase , Random, Distributions, LinearAlgebra,Base.Threads

function c_factor(n::Int)
    return 2.0*(log(n-1)+0.5772156649) - (2.0*(n-1.)/(n*1.0))
end

function depth_unsuccessful_search(n::Int)
    return convert(Int,ceil(log(2,n)))
end

function len(X::Array)
    return size(X)[1]
end
function dim(X::Array)
    return size(X)[2]
end
mutable struct Node

    X::Array # Data at the node.
    n::Array # Normal vector used to build the hyperplane that splits the data in the node.
    p::Array # Intercept point through which the hyperplane passes.
    e::Float64 # Depth of the tree to which the node belongs.

    size::Int # Size of the dataset present at the node.
    ntype::String #Specifies if the node is external or internal. Takes two values: 'exNode', 'inNode'.

    left::Node # Left child node.
    right::Node # Right child node.

    Node(X, n, p, e) = new(X,n,p,e,len(X),"")

    Node(X, n, p, e, node_type,left, right) = new(X,n,p,e,len(X),node_type,left,right)

    Node(X, n, p, e,node_type) = new(X,n,p,e,len(X),node_type)

end

mutable struct iTree
    X::Array # Data at the node.
    e::Int # Depth of the tree to which the node belongs.
    l::Int # Maxium depth a tree can reach before its creation is terminated.
    exlevel::Int # Extension level used in the splitting criteria.

    size::Int # Size of the dataset.
    dim::Int # Dimension of the dataset.
    Q::Array # List of ordered integers smaller than dim.

    exnodes::Int # The number of external nodes this tree has.

    n::Array # Normal vector at the root of this tree, which is used in creating hyperplanes for splitting critera
    p::Array # Intercept point at the root of this tree through which the splitting hyperplane passes.

    root::Node #  At each node create a new tree.

    iTree(X,e,l,exlevel) = new(X,e,l,exlevel,len(X),dim(X),Array(range(1, stop=len(X), step=1)),0)
end

function(t::iTree)()
    t.root=make_tree!(t,t.X,t.e,t.l)
end

function make_tree!(t::iTree,X,e,l)
    t.e = e
    if e >= l || len(X) <= 1                                              # A point is isolated in traning data, or the depth limit has been reached.
        t.exnodes += 1
        return Node(X, t.n, t.p, e,"exNode")
    else                                                                   # Building the tree continues. All these nodes are internal.
        mins =  minimum(X, dims=1)
        maxs =  maximum(X, dims=1)
        idxs = sample(1:t.dim, t.dim-t.exlevel-1,replace=false)            # Pick the indices for which the normal vector elements should be set to zero acccording to the extension level.
        t.n = rand(Normal(0, 1),t.dim)                                     # A random normal vector picked form a uniform n-sphere. Note that in order to pick uniformly from n-sphere, we need to pick a random normal for each component of this vector.
        t.n[idxs] .= 0
        t.p = [if min < max rand(Uniform(min, max)) else min end for (min, max) in zip(mins,maxs)]   # Picking a random intercept point for the hyperplane splitting data.
        w0=X.-t.p
        w=w0 * t.n .< 0                                                    # Criteria that determines if a data point should go to the left or right child node.
        return Node(X, t.n, t.p, e,"inNode",
        make_tree!(t,X[w,:],e+1,l),
        make_tree!(t,X[.!w,:],e+1,l),
        )
    end
    return t
end

mutable struct PathFactor

    x::Array # A single data point, which is represented as a list of floats.
    itree::iTree # Given Tree
    path_list::Array{String} #Path List
    e::Float64 #  The depth of a given node in the tree.
    path::Float64 # Path Length
    PathFactor(x::Array,itree::iTree) = new(x,itree,String[],0,0)
end

function (f::PathFactor)()
    f.path=find_path!(f,f.itree.root)
end

function find_path!(f::PathFactor,T::Node)

    if T.ntype == "exNode"

        if T.size <= 1
            return f.e
        else
            f.e = f.e + c_factor(T.size)
            return f.e
        end

    else
        p = T.p                                                             # Intercept for the hyperplane for splitting data at a given node.
        n = T.n                                                             # Normal vector for the hyperplane for splitting data at a given node.

        f.e += 1

        if ((transpose(f.x)-p) * n)[1] < 0
            push!(f.path_list,"L")
            return find_path!(f,T.left)
        else
            push!(f.path_list,"R")
            return find_path!(f,T.right)
        end
    end

end

struct iForest
    X::Array # Data used for training. It is a list of list of floats.
    nobjs::Int # Size of the dataset.
    sample::Int # Size of the sample to be used for tree creation.
    ntrees::Int # A list of tree objects.
    limit::Int # Maximum depth a tree can have.
    exlevel::Int # Exention level to be used in the creating splitting critera.
    c::Float64 # Multiplicative factor used in computing the anomaly scores.
    Trees::Array{iTree}  #Array of trees
    iForest(X::Array,nobjs::Int,sample::Int,ntrees::Int,limit::Int,exlevel::Int,c::Float64 ) = new(X,nobjs,sample,ntrees,limit,exlevel,c,iTree[])
end

iForest(X::Array,sample::Int,ntrees::Int,limit::Int,exlevel::Int) = iForest(X,len(X),sample,ntrees,limit,exlevel,c_factor(sample))

iForest(X::Array,sample::Int,ntrees::Int,limit::Int) = iForest(X,len(X),sample,ntrees,limit,0,c_factor(sample))

iForest(X::Array,sample::Int,ntrees::Int) = iForest(X,len(X),sample,ntrees,depth_unsuccessful_search(sample),0,c_factor(sample))

iForest(X::Array,sample::Int) = iForest(X,len(X),sample,100,depth_unsuccessful_search(sample),0,c_factor(sample))

function check_extension_level(forest::iForest)

    dim = size(forest.X)[1]
    if forest.exlevel < 0 || forest.exlevel > dim-1
        throw(ArgumentError( @printf("Extension level has to be an integer between 0 and  %f \n", dim-1)))
    end

end

function (f::iForest)()

    check_extension_level(f)

    Threads.@threads for i in 1:f.ntrees                                          # This loop builds an ensemble of iTrees (the forest).
        ix = sample(1:f.nobjs, f.sample)
        X_p = X[ix,:]
        tree=iTree(X_p, 0, f.limit, f.exlevel)
        tree()
        push!(f.Trees,tree)
    end
    return f
end

function compute_paths(f::iForest,X_in::Array)
    S = zeros(len(X_in))
    for i in  1:len(X_in)
        h_temp = 0
        for j in 1:f.ntrees
            pf=PathFactor(X_in[i,:],f.Trees[j])
            pf()
            h_temp += pf.path*1.0              # Compute path length for each point
        Eh = h_temp/f.ntrees                                             # Average of path length travelled by the point in all trees.
        S[i] = 2.0^(-Eh/f.c)                                            # Anomaly Score
        end
    end

    return S
end

function compute_paths(f::iForest)
    return compute_paths(f,f.X)
end

function compute_paths_single_tree(f::iForest,X_in::Array,tree_indx::Int)
    S = zeros(len(X_in))
    for i in  1:len(X_in)
        pf=PathFactor(X_in[i,:],f.Trees[tree_indx])
        pf()
        h_temp = pf.path*1.0              # Compute path length for each point
        S[i] = h_temp                                                 # Anomaly Score
    end

    return S
end

function compute_paths_single_tree(f::iForest,tree_indx::Int)
    return compute_paths_single_tree(f,f.X,tree_indx)
end

max_stars_repo_name: vikashranjan/isolation-forest

max_issues_repo_name: vikashranjan/isolation-forest

max_forks_repo_name: vikashranjan/isolation-forest


text:
# Extract a column as a Dagger array

function Iterators.partition(t::DDataset, n::Integer)
    PartitionIterator(t, n)
end

struct PartIteratorState{T}
    chunkno::Int
    chunk::T
    used::Int
end

function iterate(p::PartitionIterator{<:DDataset})
    state = if !isempty(p.c.chunks)
        PartIteratorState(1, collect(p.c.chunks[1]), 0)
    else
        PartIteratorState(1, collect(p.c), 0)
    end
    return iterate(p, state)
end

function iterate(t::PartitionIterator{<:DDataset}, p::PartIteratorState)
    if p.chunkno == length(t.c.chunks) && p.used >= length(p.chunk)
        return nothing
    end
    if p.used + t.n <= length(p.chunk)
        # easy
        nextpart = subtable(p.chunk, p.used+1:(p.used+t.n))
        return nextpart, PartIteratorState(p.chunkno, p.chunk, p.used + t.n)
    else
        part = subtable(p.chunk, p.used+1:length(p.chunk))
        required = t.n - length(part)
        r = required
        chunkno = p.chunkno
        used = length(p.chunk)
        nextchunk = p.chunk
        while r > 0
            chunkno += 1
            if chunkno > length(t.c.chunks)
                # we're done, last chunk
                return part, PartIteratorState(chunkno-1, nextchunk, used)
            else
                nextchunk = collect(t.c.chunks[chunkno])
                if r > length(nextchunk)
                    part = _merge(part, nextchunk)
                    r -= length(nextchunk)
                    used = length(nextchunk)
                else
                    part = _merge(part, subtable(nextchunk, 1:r))
                    used = r
                    r = 0
                end
            end
        end
        return part, PartIteratorState(chunkno, nextchunk, used)
    end
end

Base.eltype(iter::PartitionIterator{<:DIndexedTable}) = IndexedTable
Base.eltype(iter::PartitionIterator{<:DNDSparse}) = NDSparse

function DColumns(arrays::Tup)
    if length(arrays) == 0
        error("""DColumns must be constructed with at least
                 one column.""")
    end

    i = findfirst(x->isa(x, ArrayOp), arrays)
    wrap = isa(arrays, Tuple) ? tuple :
                                namedtuple(keys(arrays)...)∘tuple
    if i == 0
        error("""At least 1 array passed to
                 DColumns must be a DArray""")
    end

    darrays = asyncmap(arrays) do x
        isa(x, ArrayOp) ? compute(get_context(), x) : x
    end

    dist = domainchunks(darrays[i])
    darrays = map(darrays) do x
        if isa(x, DArray)
            domainchunks(x) == dist ?
                x : error("Distribution incompatible")
        else
            Distribute(dist, x)
        end
    end

    darrays = asyncmap(darrays) do x
        compute(get_context(), x)
    end

    if length(darrays) == 1
        cs = chunks(darrays[1])
        chunkmatrix = reshape(cs, length(cs), 1)
    else
        chunkmatrix = reduce(hcat, map(chunks, darrays))
    end
    cs = mapslices(x -> delayed((c...) -> Columns(wrap(c...)))(x...), chunkmatrix, 2)[:]
    T = isa(arrays, Tuple) ? Tuple{map(eltype, arrays)...} :
        wrap{map(eltype, arrays)...}

    DArray(T, domain(darrays[1]), domainchunks(darrays[1]), cs, dvcat)
end

function itable(keycols::DArray, valuecols::DArray)
    cs = map(delayed(itable), chunks(keycols), chunks(valuecols))
    cs1 = compute(get_context(),
                  delayed((xs...) -> [xs...]; meta=true)(cs...))
    fromchunks(cs1)
end

function extractarray(t, f)
    fromchunks(map(delayed(f), t.chunks))
end

# TODO: do this lazily after compute
# technically it's not necessary to communicate here

function Base.getindex(d::ColDict{<:DArray})
    rows(table(d.columns...; names=d.names))
end

function columns(t::Union{DDataset, DArray})
    cs = delayedmap(t.chunks) do c
        x = columns(c)
        if isa(x, AbstractArray)
            tochunk(x)
        elseif isa(x, Tup)
            map(tochunk, IndexedTables.astuple(x))
        else
            # this should never happen
            error("Columns $which could not be extracted")
        end
    end
    tuples = collect(get_context(), treereduce(delayed(vcat), cs))
    if length(cs) == 1
        tuples = [tuples]
    end

    if isa(tuples[1], Tup)
        arrays = map((xs...)->fromchunks([xs...]), tuples...)
        if t isa DDataset
            names = colnames(t)
        else
            names = fieldnames(eltype(t))
        end
        if all(x -> x isa Symbol, names)
            IndexedTables.namedtuple(names...)(arrays)
        else
            arrays
        end
    else
        fromchunks(tuples)
    end
end

function columns(t::Union{DDataset, DArray}, which::Tuple)
    columns(rows(t, which))
end

# TODO: make sure this is a DArray of Columns!!
Base.@pure IndexedTables.colnames(t::DArray{T}) where T<:Tup = fieldnames(T)

isarrayselect(x) = x isa AbstractArray || x isa Pair{<:Any, <:AbstractArray}

function dist_selector(t, f, which::Tup)
    if any(isarrayselect, which)
        refholder = []
        t1 = compute(t)
        w1 = map(which) do x
            isarrayselect(x) ? distfor(t1, x, refholder) : x
        end
        # this repeats the non-chunks to all other chunks,
        # then queries with the corresponding chunks
        broadcast(t1.chunks, w1...) do x...
            delayed((inp...)->f(inp[1], inp[2:end]))(x...)
        end |> fromchunks
    else
        extractarray(t, x->f(x, which))
    end
end

function dist_selector(t, f, which::AbstractArray)
    which
end

function dist_selector(t, f, which)
    extractarray(t, x->f(x,which))
end

function distfor(t, x::AbstractArray, refholder)
    y = rows(t)
    if length(y) != length(x)
        error("Input column is not the same length as the table")
    end
    d = distribute(x, domainchunks(y))
    push!(refholder, d)
    d.chunks
end

function distfor(t, x::Pair{<:Any, <:AbstractArray}, refholder)
    cs = distfor(t, x[2], refholder)
    [delayed(c->x[1]=>c)(c) for c in cs]
end

function rows(t::Union{DDataset, DArray})
    extractarray(t, rows)
end

function rows(t::Union{DDataset, DArray}, which)
    dist_selector(t, rows, which)
end

function pkeys(t::DIndexedTable)
    if isempty(t.pkey)
        Columns((Base.OneTo(length(compute(t))),))
    else
        extractarray(t, pkeys)
    end
end
pkeys(t::DNDSparse) = keys(t)

for f in [:keys, :values]
    @eval function $f(t::DNDSparse)
        extractarray(t, x -> $f(x))
    end

    @eval function $f(t::DNDSparse, which)
        dist_selector(t, $f, which)
    end
end

function column(t::DDataset, name)
    extractarray(t, x -> column(x, name))
end
function column(t::DDataset, xs::AbstractArray)
    # distribute(xs, rows(t).subdomains)
    xs
end

function pairs(t::DNDSparse)
    extractarray(t, x -> map(Pair, x.index, x.data))
end

ext: jl


text:
module TestGaussianProcesses

# using Revise
using MLJBase
using Test
using Random:seed!
import CategoricalArrays

seed!(113355)

task = load_crabs()

X, y = X_and_y(task)

# load code to be tested:
import MLJModels
import GaussianProcesses # MLJModels.GaussianProcesses_ now available for loading
using MLJModels.GaussianProcesses_

baregp = GPClassifier(target_type=String)

# split the rows:
allrows = eachindex(y)
train, test = partition(allrows, 0.7, shuffle=true)
@test sort(vcat(train, test)) == allrows

fitresult, cache, report = MLJBase.fit(baregp, 1, MLJBase.selectrows(X, train), y[train])
@test fitresult isa MLJBase.fitresult_type(baregp)
yhat = predict(baregp, fitresult, MLJBase.selectrows(X, test));

@test sum(yhat .== y[test]) / length(y[test]) >= 0.7 # around 0.7

fitresult, cache, report = MLJBase.fit(baregp, 1, X, y)
yhat2 = predict(baregp, fitresult, MLJBase.selectrows(X, test));

# gp = machine(baregp, X, y)
# fit!(gp)
# yhat2 = predict(gp, MLJBase.selectrows(X, test))

@test sum(yhat2 .== y[test]) / length(y[test]) >= 0.7

info(baregp)

end # module
true

max_stars_repo_name: ablaom/MLJModels

max_issues_repo_name: ablaom/MLJModels

max_forks_repo_name: ablaom/MLJModels


text:
"""
A library of functions used in the Julia Companion to
the book **Vectors, Matrices, and Least Squares**.
"""
module VMLS

using LinearAlgebra
using SparseArrays
using DelimitedFiles
export avg, rms, stdev, ang, correl_coef, eye, speye, kmeans, gram_schmidt
export diagonal, spdiagonal, vandermonde, toeplitz, linspace
export confusion_matrix, row_argmax, one_hot, mols_solve, cls_solve
export levenberg_marquardt, aug_lag_method
export house_sales_data, population_data, petroleum_consumption_data
export vehicle_miles_data, temperature_data, iris_data, ozone_data
export regularized_fit_data, portfolio_data, lq_estimation_data
export orth_dist_reg_data

"""
    avg(x)

Returns the average of the elements of the vector or matrix `x`.
"""
avg(x) = sum(x) / length(x)

"""
    rms(x)

Returns the RMS value of the elements of the vector or matrix `x`.
"""
rms(x) = norm(x) / sqrt(length(x))

"""
    stdev(x)

Returns the standard deviation of the elements of the vector or matrix `x`.
"""
stdev(x) = rms(x .- avg(x))

"""
    ang(x,y)

Returns the angle in radians between non-zero vectors `x` and `y`.
"""
ang(x,y) = acos(x'*y  / (norm(x) * norm(y)))

"""
    correl_coef(x,y)

Returns the correlation coefficient between non-constant vectors `x` and `y`.
"""
function correl_coef(x,y)
    xdm = x .- avg(x)
    ydm = y .- avg(y)
    return xdm'*ydm / (norm(xdm) * norm(ydm))
end

"""
    eye(n)

Returns an `n` times `n` identity matrix.
"""
eye(n) = Matrix(1.0I, n, n)

"""
    speye(n)

Returns an `n` times `n` sparse identity matrix.
"""
speye(n) = sparse(1.0I, n, n)

"""
    kmeans(X, k; maxiters = 100, tol = 1e-5)

Applies the k-means algorithm for `k` clusters to the vectors stored in `X`.

The argument `X` is a one-dimensional array of N n-vectors or an `n` times `N` matrix.

`kmeans` returns a tuple with two elements: `assignment` is an array of N
integers between 1 and k with the cluster assignment for the N data points.
`reps` is an array of k n-vectors with the k cluster representatives.
"""
function kmeans(X, k; maxiters = 100, tol = 1e-5)
    if ndims(X) == 2
        X = [X[:,i] for i in 1:size(X,2)]
    end;
    N = length(X)
    n = length(X[1])
    distances = zeros(N)
    reps = [zeros(n) for j=1:k]
    assignment = [ rand(1:k) for i in 1:N ]
    Jprevious = Inf
    for iter = 1:maxiters
        for j = 1:k
            group = [i for i=1:N if assignment[i] == j]
            reps[j] = sum(X[group]) / length(group);
        end;
        for i = 1:N
            (distances[i], assignment[i]) =
                findmin([norm(X[i] - reps[j]) for j = 1:k])
        end;
        J = norm(distances)^2 /  N
        println("Iteration ", iter, ": Jclust = ", J, ".")
        if iter > 1 && abs(J - Jprevious) < tol * J
            return assignment, reps
        end
        Jprevious = J
    end
end

"""
    gram_schmidt(a; tol = 1e-10)

Applies the Gram-Schmidt algorithm to the vectors stored in the array `a`
and returns the result as an array of (orthonormal) vectors.  Issues a
warning if the original vectors are linearly dependent.
"""
function gram_schmidt(a; tol = 1e-10)
    q = []
    for i = 1:length(a)
        qtilde = a[i]
        for j = 1:i-1
            qtilde -= (q[j]'*a[i]) * q[j]
        end
        if norm(qtilde) < tol
            println("Vectors are linearly dependent.")
            return q
        end
        push!(q, qtilde/norm(qtilde))
    end;
    return q
end

"""
    diagonal(x)

Returns a diagonal matrix with the entries of the vector `x` on its diagonal.
"""
diagonal(x) = diagm(0 => x)

"""
    spdiagonal(x)

Returns a sparse diagonal matrix with the entries of the vector `x` on its diagonal.
"""
spdiagonal(x) = spdiagm(0 => x)

"""
    vandermonde(t, n)

Returns the Vandermonde matrix with `n` columns and i-th column `t`^(i-1).
"""
vandermonde(t,n) = hcat( [t.^i for i = 0:n-1]... )

"""
    toeplitz(a, n)

Returns the Toeplitz matrix with `n` columns and the vector `a` in the
leading positions of the first column.
"""
function toeplitz(a,n)
    m = length(a)
    T = zeros(n+m-1,n)
    for i=1:m
        T[i : n+m : end] .= a[i]
    end
    return T
end

"""
    linspace(a,b,n)

Returns a vector with `n` equally spaced numbers between `a` and `b`.
"""
linspace(a,b,n) = Vector(range(a, stop = b, length = n))

"""
    confusion_matrix(y, yhat, K=2)

Returns the confusion matrix for a data vector `y` and the vector of predictions
`yhat`.  If `K` is 2, the vectors are Boolean.  If `K` is greater than 2, they
contain integers from 1,...,K.
"""
function confusion_matrix(y, yhat, K=2)
    if K==2
        tp = sum( (y .== true) .& (yhat .== true) )
        fn = sum( (y .== true) .& (yhat .== false) )
        fp = sum( (y .== false) .& (yhat .== true) )
        tn = sum( (y .== false) .& (yhat .== false) )
        return [tp fn; fp tn]
    else
        C = zeros(K,K)
        for i = 1:K  for j=1:K
            C[i,j] = sum( (y .== i) .& (yhat .== j) )
        end end
        return C
    end;
end;

"""
    row_argmax(A)

Returns a `size(A,1)` vector with as i-th element the column index of the
largest element in row i of `A`.
"""
row_argmax(A) = [ argmax(A[i,:]) for i = 1:size(A,1) ]

"""
    one_hot(x, K)

Returns the one-hot encoding of the vector `x`.
"""
function one_hot(y, K)
    n = length(y)
    Y = zeros(n, K)
    for j in 1:K
       Y[findall(y .== j), j] .= 1.0
    end
    return Y
end

"""
    mols_solve(As, bs, lambdas)

Returns the solution of the multi-objective least squares problem
with coefficient matrices in the array `As`, right-hand side vectors in the
array `bs`, and weights in the array `lambdas`.
"""
function mols_solve(As,bs,lambdas)
    k = length(lambdas);
    Atil = vcat([sqrt(lambdas[i])*As[i] for i=1:k]...)
    btil = vcat([sqrt(lambdas[i])*bs[i] for i=1:k]...)
    return Atil \ btil
end

"""
   cls_solve(A, b, C, d)

Returns the solution of the constrained least squares problem with coefficient
matrices `A` and `C`, and right-hand side vectors or matrices `b` and `d`.
"""
function cls_solve(A,b,C,d)
    m, n = size(A)
    p, n = size(C)
    Q, R = qr([A; C])
    Q = Matrix(Q)
    Q1 = Q[1:m,:]
    Q2 = Q[m+1:m+p,:]
    Qtil, Rtil = qr(Q2')
    Qtil = Matrix(Qtil)
    w = Rtil \ (2*Qtil'*Q1'*b - 2*(Rtil'\d))
    return R \ (Q1'*b - Q2'*w/2)
end

"""
    levenberg_marquardt(f, Df, x1, lambda1; kmax=500, tol=1e-6)

Applies the Levenberg-Marquardt algorithm to the function defined in `f` and
`Df`, with starting point `x1` and initial regularization parameter `lambda1`.
The function returns the final iterate `x` and a dictionary with the convergence
history.
"""
function levenberg_marquardt(f, Df, x1, lambda1; kmax=100, tol=1e-6)
    n = length(x1)
    x = x1
    lambda = lambda1
    objectives = zeros(0,1)
    residuals = zeros(0,1)
    for k = 1:kmax
         fk = f(x)
         Dfk = Df(x)
         objectives = [objectives; norm(fk)^2]
         residuals = [residuals; norm(2*Dfk'*fk)]
         if norm(2*Dfk'*fk) < tol
             break
         end;
         xt = x - [ Dfk; sqrt(lambda)*eye(n) ] \ [ fk; zeros(n) ]
         if norm(f(xt)) < norm(fk)
             lambda = 0.8*lambda
             x = xt
         else
             lambda = 2.0*lambda
         end
    end
    return x, Dict([ ("objectives", objectives), ("residuals", residuals) ])
end

"""
    aug_lag_method(f, Df, g, Dg, x1, lambda1; kmax=100, feas_tol=1e-4, oc_tol=1e-4)

Applies the augmented Lagrangian method to the constrained nonlinear least
squares problem defined by `f`, `Df`, `g`, `Dg`, with starting point `x1`.
The subproblems are solved using the Levenberg-Marquardt method with
initial regularization parameter `lambda1`.  Returns the final iterate `x`,
multiplier `z`, and a dictionary with the convergence history.
"""
function aug_lag_method(f, Df, g, Dg, x1, lambda1;  kmax = 100,
    feas_tol = 1e-4, oc_tol = 1e-4)
    x = x1
    z = zeros(length(g(x)))
    mu = 1.0
    feas_res = [norm(g(x))]
    oc_res = [norm(2*Df(x)'*f(x) + 2*mu*Dg(x)'*z)]
    lm_iters = zeros(Int64,0,1);
    mus = [mu]
    for k=1:kmax
        F(x) = [f(x); sqrt(mu)*(g(x) + z/(2*mu))]
        DF(x) = [Df(x); sqrt(mu)*Dg(x)]
        x, hist = levenberg_marquardt(F, DF, x, lambda1, tol = oc_tol)
        z = z + 2*mu*g(x)
        feas_res = [feas_res; norm(g(x))]
        oc_res = [oc_res; hist["residuals"][end]]
        lm_iters = [lm_iters; length(hist["residuals"])]
        if norm(g(x)) < feas_tol
            break
        end
        mu = (norm(g(x)) < 0.25*feas_res[end-1]) ? mu : 2*mu
        mus = [mus; mu]
    end
    return x, z, Dict([ ("lm_iterations", lm_iters),
         ("feas_res", feas_res), ("oc_res", oc_res), ("mus", mus) ])
end

"""
    petroleum_consumption_data()

Returns a 34-vector with the world annual petroleum consumption between
1980 and 2013, in thousand barrels/day (discussed on page 252).
"""
petroleum_consumption_data() = [
    63122, 60953, 59551, 58785, 59795, 60083, 61819, 63107, 64978, 66090,
    66541, 67186, 67396, 67619, 69006, 70258, 71880, 73597, 74274, 75975,
    76928, 77732, 78457, 80089, 83063, 84558, 85566, 86724, 86046, 84972,
    88157, 89105, 90340, 91195 ]

include("house_sales_data.jl")
include("population_data.jl")
include("vehicle_miles_data.jl")
include("temperature_data.jl")
include("iris_data.jl")
include("ozone_data.jl")
include("regularized_fit_data.jl")
include("portfolio_data.jl")
include("lq_estimation_data.jl")
include("orth_dist_reg_data.jl")

end # module

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
# ArgParse can't be a dependency here,
# so users need to install it.
using ArgParse, PackageCompiler

Base.@ccallable function julia_main(args::Vector{String})::Cint

    s = ArgParseSettings("Static Julia Compiler",
                         version = "$(basename(@__FILE__)) version 0.7",
                         add_version = true)

    @add_arg_table s begin
        "juliaprog"
            arg_type = String
            required = true
            help = "Julia program to compile"
        "cprog"
            arg_type = String
            help = "C program to compile (required only when building an executable, if not provided a minimal driver program is used)"
        "--verbose", "-v"
            action = :store_true
            help = "increase verbosity"
        "--quiet", "-q"
            action = :store_true
            help = "suppress non-error messages"
        "--builddir", "-d"
            arg_type = String
            metavar = "<dir>"
            help = "build directory"
        "--outname", "-n"
            arg_type = String
            metavar = "<name>"
            help = "output files basename"
        "--snoopfile", "-p"
            arg_type = String
            metavar = "<file>"
            help = "specify script calling functions to precompile"
        "--clean", "-c"
            action = :store_true
            help = "remove build directory"
        "--autodeps", "-a"
            action = :store_true
            help = "automatically build required dependencies"
        "--object", "-o"
            action = :store_true
            help = "build object file"
        "--shared", "-s"
            action = :store_true
            help = "build shared library"
        "--init-shared", "-i"
            action = :store_true
            help = "add `init_jl_runtime` and `exit_jl_runtime` to shared library for runtime initialization"
        "--executable", "-e"
            action = :store_true
            help = "build executable file"
        "--rmtemp", "-t"
            action = :store_true
            help = "remove temporary build files"
        "--copy-julialibs", "-j"
            action = :store_true
            help = "copy Julia libraries to build directory"
        "--copy-file", "-f"
            arg_type = String
            action = :append_arg
            dest_name = "copy-files"
            metavar = "<file>"
            help = "copy file to build directory, can be repeated for multiple files"
        "--release", "-r"
            action = :store_true
            help = "build in release mode, implies `-O3 -g0` unless otherwise specified"
        "--Release", "-R"
            action = :store_true
            help = "perform a fully automated release build, equivalent to `-atjr`"
        "--sysimage", "-J"
            arg_type = String
            metavar = "<file>"
            help = "start up with the given system image file"
        "--home", "-H"
            arg_type = String
            metavar = "<dir>"
            help = "set location of `julia` executable"
        "--startup-file"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "load `~/.julia/config/startup.jl`"
        "--handle-signals"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "enable or disable Julia's default signal handlers"
        "--sysimage-native-code"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "use native code from system image if available"
        "--compiled-modules"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "enable or disable incremental precompilation of modules"
        "--depwarn"
            arg_type = String
            metavar = "{yes|no|error}"
            range_tester = (x -> x ∈ ("yes", "no", "error"))
            help = "enable or disable syntax and method deprecation warnings"
        "--warn-overwrite"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "enable or disable method overwrite warnings"
        "--compile"
            arg_type = String
            metavar = "{yes|no|all|min}"
            range_tester = (x -> x ∈ ("yes", "no", "all", "min"))
            help = "enable or disable JIT compiler, or request exhaustive compilation"
        "--cpu-target", "-C"
            arg_type = String
            metavar = "<target>"
            help = "limit usage of CPU features up to <target> (implies default `--sysimage-native-code=no`)"
        "--optimize", "-O"
            arg_type = Int
            metavar = "{0,1,2,3}"
            range_tester = (x -> x ∈ (0, 1, 2, 3))
            help = "set the optimization level"
        "--debug", "-g"
            arg_type = Int
            metavar = "<level>"
            range_tester = (x -> x ∈ (0, 1, 2))
            help = "enable / set the level of debug info generation"
        "--inline"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "control whether inlining is permitted"
        "--check-bounds"
            arg_type = String
            metavar = "{yes|no}"
            range_tester = (x -> x ∈ ("yes", "no"))
            help = "emit bounds checks always or never"
        "--math-mode"
            arg_type = String
            metavar = "{ieee,fast}"
            range_tester = (x -> x ∈ ("ieee", "fast"))
            help = "disallow or enable unsafe floating point optimizations"
        "--cc"
            arg_type = String
            metavar = "<cc>"
            help = "system C compiler"
        "--cc-flag"
            arg_type = String
            action = :append_arg
            dest_name = "cc-flags"
            metavar = "<flag>"
            help = "pass custom flag to the system C compiler when building a shared library or executable, can be repeated for multiple flags"
    end

    s.epilog = """
        examples:\n
        \ua0\ua0juliac.jl -vae hello.jl        # verbose, build executable and deps\n
        \ua0\ua0juliac.jl -vae hello.jl prog.c # embed into user defined C program\n
        \ua0\ua0juliac.jl -qo hello.jl         # quiet, build object file only\n
        \ua0\ua0juliac.jl -vosej hello.jl      # build all and copy Julia libs\n
        \ua0\ua0juliac.jl -vRe hello.jl        # fully automated release build
        """

    parsed_args = parse_args(args, s)

    parsed_args["copy-files"] == String[] && (parsed_args["copy-files"] = nothing)

    # TODO: in future it may be possible to broadcast dictionary indexing, see: https://discourse.julialang.org/t/accessing-multiple-values-of-a-dictionary/8648
    if getindex.(Ref(parsed_args), ["clean", "object", "shared", "executable", "rmtemp", "copy-julialibs", "copy-files"]) == [false, false, false, false, false, false, nothing]
        parsed_args["quiet"] || println("nothing to do, exiting\ntry \"$(basename(@__FILE__)) -h\" for more information")
        exit(0)
    end

    juliaprog = pop!(parsed_args, "juliaprog")
    filter!(kv -> kv.second !== nothing && kv.second !== false, parsed_args)
    kw_args = [Symbol(replace(kv.first, "-" => "_")) => kv.second for kv in parsed_args]

    static_julia(juliaprog; kw_args...)

    return 0
end

if get(ENV, "COMPILE_STATIC", "false") == "false"
    julia_main(ARGS)
end


text:
include("functions.jl")

"""
    c4_traving(line::AtomicLine)
Calculate the \$C_4\$ interaction constant (quadratic Stark effect) using the
recipe of Traving (1960), "Uber die Theorie der Druckverbreiterung von Spektrallinien",
p 93.
"""
function c4_traving(line::HydrogenicLine)
    n_eff_u = Transparency.n_eff(line.χ∞, line.χj, line.Z)
    n_eff_l = Transparency.n_eff(line.χ∞, line.χi, line.Z)
    C4 = (e^2 * inv_4πε0 * a_0^3 * 2 * π / (h * 18 * line.Z^4) *
        ((n_eff_u * (5 * n_eff_u^2 + 1))^2 - (n_eff_l * (5 * n_eff_l^2 + 1))^2))
    return C4 |> u"m^4 / s"
end

"""
    function const_unsold(line::HydrogenicLine; H_scaling::Real=1, He_scaling::Real=1)
Compute atmosphere-independent constant for γ_unsold, to be used in function `γ_unsold`.
Based on expressions from RH broad.c, which uses formula in Mihalas (1978),
pp 282, 286-287, eq. (9-50) for v_rel, table 9-1 and eq. (9-76) for the interaction
coefficient C6.  Input is an `AtomicLine`, which contains the necessary energies and
element weight. The van der Waals broadening can be scaled for both H and He perturbers
using `H_scaling` and `He_scaling`.
"""
function const_unsold(line::HydrogenicLine;
                      H_scaling::Real=1,
                      He_scaling::Real=1)

    Δr = (Ry^2 * (1 / (line.χ∞ - line.χj)^2 - 1 / (line.χ∞ - line.χi)^2)) |> u"J/J"
    C6 = ustrip((2.5 * e^2 * αp * inv_4πε0^2 * 2 * π *
                 (line.Z * a_0)^2 / h * Δr) |> u"C^2 * m^6 / (F * J * s)")
    v_rel_const = ustrip(8 * k_B / (π * line.atom_weight) |> u"J/(K * kg)")
    v_rel_H = v_rel_const * (1 + line.atom_weight / mass_H)
    v_rel_He = v_rel_const * (1 + line.atom_weight / mass_He)
    return 8.08 * (H_scaling * v_rel_H^0.3 + He_scaling * abund_He * v_rel_He^0.3) * C6^0.4
end

# TODO
# Verify mean_atomic_weight::Unitful.Mass=28 * m_u
"""
    const_quadratic_stark(line::AtomicLine;
                          mean_atomic_weight::Unitful.Mass=28 * m_u,
                          scaling::Real=1)
Calculate height-independent constant to use in `γ_quadratic_stark`, using the recipe
from RH, which is based on the following estimate:
\$\$
\\gamma = 11.37 \\cdot vrel^{1/3} * C_4^{2/3} * (n_e + n_{ion}),
\$\$
Using the estimate for \$C_4\$ from Traving (1960), "Uber die Theorie der
Druckverbreiterung von Spektrallinien", p 93., and \$n_{ion}\\approx n_e\$
(following Gray).
"""
function const_quadratic_stark(line::HydrogenicLine;
                               mean_atomic_weight::Unitful.Mass=28 * m_u,
                               scaling::Real=1)
    C = ustrip(8 * k_B / (π * line.atom_weight) |> u"J/(K * kg)")
    Cm = ((1 + line.atom_weight / m_e)^(1/6) +
          (1 + line.atom_weight / mean_atomic_weight)^(1/6))
    C4 = ustrip(c4_traving(line) |> u"m^4 / s")
    cStark23 = 11.37u"m^3 / s" * (scaling * C4)^(2/3)
    return C^(1/6) * cStark23 * Cm
end

function γ_constant(line::HydrogenicLine,
                    temperature::Array{<:Unitful.Temperature},
                    neutral_hydrogen_density::Array{<:NumberDensity},
                    electron_density::Array{<:NumberDensity})

    u=2 #line.j
    l=1 #line.i

    unsold_const = const_unsold(line)
    quad_stark_const = const_quadratic_stark(line)

    γ = γ_unsold.(unsold_const, temperature, neutral_hydrogen_density)
    # γ .+= line.Aji
    γ .+= 4.702e8u"s^-1" # Almost the same as Aji
    γ .+= γ_linear_stark.(electron_density, u, l)
    γ .+= γ_quadratic_stark.(electron_density, temperature, stark_constant=quad_stark_const)

    return γ

end

"""
Compute damping parameter.
"""
function damping(γ::Unitful.Frequency, λ::Unitful.Length, ΔλD::Unitful.Length)
    return (γ * λ^2 / (4 * π * c_0 * ΔλD)) |> Unitful.NoUnits
end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
function corrp_test()
	x = [1.,2,3,5,10,22,0,1,-3,4,0,8,10,3,4,5,5,5,-10,0];
	y = [2.,3,4,5,7,10,10,11,-3,4,1,-1,0,1,2,3,4,4,4,-5];
	r1,p1 = corrp(x,y);
	# compare to Matlab results ([r,p] = corr(x',y'))
	@test r1 ≈ 0.2923 atol = .0001
	@test p1 ≈ 0.2111 atol = .0001

	x[2] = NaN;
	y[[2,4]] .= NaN;
	r2,p2 = corrp(x,y);
	@test r2 ≈ 0.2892 atol = .0001
	@test p2 ≈ 0.2444 atol = .0001
end

corrp_test();

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
using Nystrom
using StaticArrays
using LinearAlgebra
using IterativeSolvers
using Plots

##
k = 3.3
pde = Maxwell(;dim=3,k)

# exact solution
G    = SingleLayerKernel(pde)
xs   = SVector(0.1,0.2,0.3)
x0   = SVector(5,5,5)
pol    = SVector(1+im,-2.,3.)
E    = (dof) -> G(dof,xs)*pol
exa  = E(x0)

# geometry and parameters
sph_radius = 2
geo = ParametricSurfaces.Sphere(;radius=sph_radius)
Ω = Geometry.Domain(geo)
Γ = boundary(Ω)

## CFIE
α = 1
β = im*k/sph_radius
order = 5         # quadrature order 1D
ndofs = Float64[]
errs = Float64[]
iterative = false;
##
for n in [4]
    M     = meshgen(Γ,(n,n))
    nmesh  = NystromMesh(view(M,Γ);order)
    N     = Nystrom.diagonal_ncross_matrix(nmesh)
    J     = Nystrom.diagonal_jacobian_matrix(nmesh)
    dualJ = Nystrom.diagonal_dualjacobian_matrix(nmesh)
    γ₀E   = γ₀(E,nmesh)
    S,D   = Nystrom.single_doublelayer_dim(pde,nmesh)
    Sm = Nystrom.materialize(S)
    Dm = Nystrom.materialize(D)
    rhs   = dualJ * ncross(γ₀E)

    @info "Assembling matrix..."
    L     = Nystrom.assemble_dim_nystrom_matrix_Luiz(nmesh, α, β, Dm, Sm)
    @info "Solving..."
    if iterative
        ϕ_coeff = gmres(L,rhs; verbose=true,maxiter=600,restart=600,abstol=1e-6)
    else
        ϕ_coeff = L \ rhs
    end
    ϕ     = J*ϕ_coeff
    Spot  = Nystrom.MaxwellEFIEPotential(pde,nmesh)
    Dpot  = Nystrom.MaxwellMFIEPotential(pde,nmesh)
    Eₐ    = (x) -> α*Dpot[ϕ](x) + β*Spot[ncross(ϕ)](x)
    er    = (Eₐ(x0) - exa)/norm(exa)

    ndof = length(Nystrom.dofs(nmesh))
    err = norm(er,Inf)

    push!(ndofs, ndof)
    push!(errs, err)
    @info "" ndof err
end

## Plot
sqrt_ndofs = sqrt.(ndofs)
fig = plot(sqrt_ndofs,errs,xscale=:log10,yscale=:log10,m=:o,label="error",lc=:black)
title = "k=$k, α=$α, β=$β, qorder=$order"
plot!(xlabel="√ndofs",ylabel="error",title=title)
for p in 1:5
    cc = errs[end]*sqrt_ndofs[end]^p
    plot!(fig,sqrt_ndofs,cc./sqrt_ndofs.^p,xscale=:log10,yscale=:log10,label="h^$p",ls=:dash)
end
display(fig)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using jInv.Mesh
using Multigrid
using SparseArrays
using LinearAlgebra
println("************************************************* GMG for Elasticity 2D using Vanka ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [128,128];
Mr = getRegularMesh(domain,n)
mu = 1.0*ones(prod(Mr.n));
lambda = 1.0.*mu;

mutable struct ElasticityParam
	Mesh   			:: RegularMesh;
	lambda          :: Array{Float64};
	mu				:: Array{Float64}
	MixedFormulation:: Bool
end

Hparam = ElasticityParam(Mr,lambda,mu,true);

function restrictParam(mesh_fine,mesh_coarse,param_fine,level)
	mu_c 	 	= restrictCellCenteredVariables(param_fine.mu,mesh_fine.n);
	lambda_c 	= restrictCellCenteredVariables(param_fine.lambda,mesh_fine.n);
	return ElasticityParam(mesh_coarse,lambda_c,mu_c,true);
end

# function GetVectorLapOperatorMixedFormulation(M::RegularMesh, mu::Vector,lambda::Vector)
	# vecGRAD,~,Div,nf, = getDifferentialOperators(M,2);
	# Mu     			  = getTensorMassMatrix(M,mu[:],saveAvMat=false)[1];
	# Rho	   			  = getFaceMassMatrix(M,1e-5*mu, saveMat = false, avN2C = avN2C_Nearest);
	# A 				  = vecGRAD'*Mu*vecGRAD + Rho;
	# A 				  = spdiagm(0=>ones(size(Rho,1)));
	# C 				  = -(spdiagm(0=>0.0./(lambda[:])));
	# Div.nzval .*= -1.0;
	# H = [A  Div' ; -Div  -C;];
	# return H;
# end

# function GetVectorLapOperator(M::RegularMesh, mu::Vector,lambda::Vector)
	# vecGRAD,~,Div,nf, = getDifferentialOperators(M,2);
	# Mu     			  = getTensorMassMatrix(M,mu[:],saveAvMat=false)[1];
	# Rho	   			  = getFaceMassMatrix(M,1e-5*mu, saveMat = false, avN2C = avN2C_Nearest);
	# A 				  = vecGRAD'*Mu*vecGRAD + Rho;
	# return A;
# end

function getOperator(Mr,Hparam)
	# Ar = GetVectorLapOperator(Mr, Hparam.mu,Hparam.lambda);
	# Ar = GetVectorLapOperatorMixedFormulation(Mr, Hparam.mu,Hparam.lambda);
	Ar = GetLinearElasticityOperatorMixedFormulation(Mr, Hparam.mu,Hparam.lambda);
	Ar = Ar + 1e-3*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
	return Ar;
end

levels      = 5;
numCores 	= 2;
maxIter     = 5;
relativeTol = 1e-10;
relaxType   = "VankaFaces";
# relaxType   = "Jac";
relaxParam  = 0.5;
relaxPre 	= 1;
relaxPost   = 1;
cycleType   ='V';
coarseSolveType = "Julia";
transferOperatorType = "SystemsFacesMixedLinear";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,transferOperatorType);

MGsetup(getMultilevelOperatorConstructor(Hparam,getOperator,restrictParam),Mr,MG,1,true);
Ar = MG.As[1]; # It is actually transposed, but matrix is symmetric.
# MGsetup(Ar,Mr,MG,1,true);

N = size(Ar,2);
sol = randn(N);
sol = sol/norm(sol);
sol = sol .- sum(sol)./N;
b = Ar*sol; b = b/norm(b);
x = zeros(N);

println("****************************** Stand-alone GMG: ******************************")
x, = solveMG(MG,copy(b),x,true);
@test norm(Ar*x - b) < 0.05;

println("****************************** GMG + CG: ******************************")
x[:].=0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.01;

println("TODO: see why the code below doesn't work.")
clear!(MG)

println("************************************************* GMG for Elasticity 3D  using Vanka ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0,0.0,1.0];
n = [8,8,8];
Mr = getRegularMesh(domain,n)
mu = 1.0*ones(prod(Mr.n));
lambda = 0.5.*mu;

Hparam = ElasticityParam(Mr,lambda,mu,true);

MGsetup(getMultilevelOperatorConstructor(Hparam,getOperator,restrictParam),Mr,MG,size(b,1),true);
Ar = MG.As[1]; # It is actually transposed, but matrix is symmetric.
N = size(Ar,2);
b = Ar*rand(N); b = b/norm(b);
x = zeros(N);

println("****************************** Stand-alone 3D GMG: ******************************")
solveMG(MG,b,x,true);
@test norm(Ar*x - b) < 0.1;
println("****************************** 3D GMG + CG: ******************************")
x[:].=0.0;
solveCG_MG(Ar,MG,b,x,true)
@test norm(Ar*x - b) < 0.1;


text:
#i64o-helpers.jl
#various utility functions that operate on ArrayNums as well as metaprogramming
#macros to help duplicate functions for both standard 64-bit integers and
#arraynums.

doc"""
`Unums.max_fsize(::Int64)` retrieves the maximum possible fsize value based on
the FSS value.
"""
max_fsize(FSS::Int64) = to16((1 << FSS) - 1)

#two helper functions to enable fracproc.
function __ffcall(s::Symbol, inner_block, tname, p...)
  fracname = Symbol("frac_", s, "!")
  #outer_curly is the curly variadic notation for the function name.
  outer_curly = Expr(:curly)
  append!(outer_curly.args, [fracname, :ESS, :FSS])
  #build up the outer call, which is the representation for the full function definition
  outer_call = Expr(:call)
  append!(outer_call.args, [outer_curly, :(x::$tname{ESS,FSS}), p...])
  #build up the outer block, which is the actual code.
  outer_block = Expr(:block)
  push!(outer_block.args, inner_block)
  #build up the full outer function, which merges the call and the block.
  outer_function = Expr(:function)
  push!(outer_function.args, outer_call)
  push!(outer_function.args, outer_block)
  outer_function
end

function __smallexpr(s::Symbol, p...)
  #build up the inner call mechanism
  inner_call = Expr(:call)
  inner_call.args = [s, :(x.fraction), p...]
  :(x.fraction = $inner_call; x)
end

function __largeexpr(s::Symbol, p...)
  bangsym = Symbol(s, "!")
  inner_call = Expr(:call)
  inner_call.args = [bangsym, :(x.fraction), p...]
  :($inner_call; x)
end

doc"""
`Unums.@fracproc(fn)` generates two functions that takes Unum{ESS,FSS} and executes
fn on the fraction part of the unum, if FSS < 7.  if FSS >= 7 then it executes
fn! on the fraction part of the unum.
"""
macro fracproc(name, params...)
  fracname = Symbol("frac_", name, "!")
  bangname = Symbol(name, "!")
  smallexpr = __smallexpr(name, params...)
  largeexpr = __largeexpr(name, params...)
  smallcall = __ffcall(name, smallexpr, :UnumSmall, params...)
  largecall = __ffcall(name, largeexpr, :UnumLarge, params...)
  esc(quote
    Base.@__doc__ $smallcall
    $largecall
  end)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Writer
# ======

"""
    BigBedWriter(output::IO)

Create a data writer of the BigBed file format.

# Arguments
* `output`: data sink
"""
type BigBedWriter{T<:IO} <: Bio.IO.AbstractWriter
    output::T
    block_size::Int
    items_per_slot::Int
    compress::Bool
end

function Bio.IO.stream(writer::BigBedWriter)
    return writer.output
end

function BigBedWriter(output::IO)
    # These parameters are constants for now since valid values are uncertain.
    block_size = 256
    items_per_slot = 512
    compress = true
    return BigBedWriter(output, block_size, items_per_slot, compress)
end

function Base.write(writer::BigBedWriter, intervals::IntervalCollection{BEDMetadata})
    return write_bigbed_bigwig(
        writer.output,
        :bigbed,
        intervals,
        writer.block_size,
        writer.items_per_slot,
        writer.compress)
end

immutable BigBedChromInfo
    name::StringField   # chromosome name
    id::UInt32          # unique id for chromosome
    size::UInt32        # size of chromosome
    item_count::UInt64  # number of items
end

immutable BigBedBounds
    offset::UInt64
    chrom_ix::UInt32
    start::UInt32
    stop::UInt32
end

function bigbed_chrom_info(intervals::IntervalCollection, chrom_sizes::Dict)
    info = Array(BigBedChromInfo, length(intervals.trees))
    seqnames = collect(StringField, keys(intervals.trees))
    sort!(seqnames)
    for (i, seqname) in enumerate(seqnames)
        info[i] = BigBedChromInfo(
            seqname,
            i - 1,
            chrom_sizes[seqname],
            length(intervals.trees[seqname]))
    end
    return info
end

function bigbed_btree_count_levels(max_block_size, item_count)
    levels = 1
    while item_count > max_block_size
        item_count = div(item_count + max_block_size - 1, max_block_size)
        levels += 1
    end
    return levels
end

function bigbed_btree_write_index_level(out::IO, block_size,
                                        items::Vector{BigBedChromInfo},
                                        key_size, index_offset, level)
    # See writeIndexLevel in bPlusTree.c

    # calculate number of nodes to write at this level
    slot_size_per = block_size ^ level
    node_size_per = slot_size_per * block_size
    node_count = div(length(items) + node_size_per - 1, node_size_per)

    # calculate sizes and offsets
    block_header_size = 4
    value_size = 8
    bytes_in_index_block = block_header_size + block_size * (key_size + sizeof(UInt64))
    bytes_in_leaf_block = block_header_size + block_size * (key_size + value_size)
    bytes_in_next_level_block = level == 1 ? bytes_in_leaf_block : bytes_in_index_block
    level_size = node_count * bytes_in_index_block
    end_level = index_offset + level_size
    next_child = end_level

    for i in 1:node_size_per:length(items)
        # calculate size of this block
        count_one = min(block_size, div(length(items) - (i-1) + slot_size_per - 1, slot_size_per))

        # write block header
        write(out, BigBedBTreeNode(0, 0, count_one))

        # write used slots
        slots_used = 0
        end_ix = min(length(items), i + node_size_per)
        for j in i:slot_size_per:end_ix
            item = items[j]
            write(out, item.name)
            # pad name
            for k in length(item.name)+1:key_size
                write(out, 0x00)
            end
            write(out, convert(UInt64, next_child))
            next_child += bytes_in_next_level_block
            slots_used += 1
        end
        @assert slots_used == count_one "Incorrect number of items written to index"

        # write out empty slots as all zero
        slot_size = key_size + 8
        for j in count_one:(block_size-1)
            for k in 1:slot_size
                write(out, 0x00)
            end
        end
    end

    return end_level
end

function bigbed_btree_write_leaf_level(out::IO, block_size,
                                       items::Vector{BigBedChromInfo},
                                       key_size, val_size)
    # See writeLeafLevel in bPlusTree.c

    count_one = 0
    count_left = length(items)
    i = 0
    while i < length(items)
        count_one = min(count_left, block_size)

        # write block header
        write(out, BigBedBTreeNode(1, 0, count_one))

        # write out position in genome and in file for each item
        for j in 1:count_one
            @assert i + j <= length(items) "Incorrect item count when writing B-tree leaves"
            item = items[i + j]
            write(out, item.name)
            # pad name
            for k in length(item.name)+1:key_size
                write(out, 0x00)
            end
            write(out, convert(UInt32, item.id))
            write(out, convert(UInt32, item.size))
        end

        # pad out any unused bits of last block with zeroes
        slot_size = key_size + val_size
        for j in count_one:(block_size-1)
            for k in 1:slot_size
                write(out, 0x00)
            end
        end

        count_left -= count_one
        i += count_one
    end
end

function bigbed_write_chrom_tree(out::IO, intervals::IntervalCollection,
                                 block_size,
                                 chrom_sizes::Dict=Dict{AbstractString, Int64}())
    # See bbiWriteChromInfo in bbiWrite.c
    length(intervals.trees)
    for seqname in keys(intervals.trees)
        if !haskey(chrom_sizes, seqname)
            chrom_sizes[seqname] =
                length(intervals.trees[seqname]) > 0 ?
                    intervals.trees[seqname].root.maxend : 0
        end
    end

    max_chrom_name_size = 0
    for seqname in keys(intervals.trees)
        max_chrom_name_size = max(max_chrom_name_size, length(seqname))
    end

    info = bigbed_chrom_info(intervals, chrom_sizes)
    item_count = length(info)
    chrom_block_size = min(block_size, item_count)
    value_size = 8; # sizeof(chrom_id) + sizeof(chrom_size)

    # write header. See bptFileBulkIndex in bPlusTree.c
    header = BigBedBTreeHeader(
        BIGBED_BTREE_MAGIC,
        chrom_block_size,
        max_chrom_name_size,
        value_size,
        item_count,
        0)
    write(out, header)
    index_offset = position(out)

    # write non-leaf nodes.
    levels = bigbed_btree_count_levels(chrom_block_size, item_count)
    for i in (levels-1):-1:1
        end_level_offset = bigbed_btree_write_index_level(out, chrom_block_size,
                                                          info, max_chrom_name_size,
                                                          index_offset, i)
        index_offset = position(out)
        @assert index_offset == end_level_offset "Incorrect file offsets reported."
    end

    # write leaf nodes
    bigbed_btree_write_leaf_level(out, chrom_block_size, info,
                                  max_chrom_name_size, value_size)
    return info
end

function bigbed_count_sections_needed(items::Vector{BigBedChromInfo}, items_per_slot)
    count = 0
    for item in items
        count_one = div(item.item_count + items_per_slot - 1, items_per_slot)
        count += count_one
    end
    return count
end

function bigbed_write_blocks(out::IO, intervals::IntervalCollection,
                             chrom_info::Vector{BigBedChromInfo},
                             items_per_slot, bounds::Vector{BigBedBounds},
                             section_count, compressed::Bool,
                             res_try_count, res_scales, res_sizes)
    # See writeBlocks in bedToBigBed.c

    section_ix = 0
    buf = IOBuffer()
    res_ends = zeros(Int, res_try_count)
    max_block_size = 0

    for info in chrom_info
        tree = intervals.trees[info.name]
        it_state = start(tree)
        fill!(res_ends, 0)

        while !done(tree, it_state)
            section_ix += 1
            item_ix = 1
            block_offset = position(out)
            first_interval = first(intervals.trees[info.name])
            start_pos = -1
            end_pos  = first_interval.last

            while item_ix <= items_per_slot && !done(tree, it_state)
                interval, it_state = next(tree, it_state)

                interval_start = interval.first - 1 # bed is 0-based
                interval_end   = interval.last
                if start_pos < 0
                    start_pos = interval_start
                end
                end_pos = max(end_pos, interval_end)

                write(buf, convert(UInt32, info.id))
                write(buf, convert(UInt32, interval_start))
                write(buf, convert(UInt32, interval_end))
                write_optional_fields(buf, interval, false)
                write(buf, '\0')

                item_ix += 1

                for res_try in 1:res_try_count
                    res_end = res_ends[res_try]
                    if interval_start >= res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = interval_start + res_scales[res_try]
                    end
                    while interval_end > res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = res_end + res_scales[res_try]
                    end
                end
            end

            block = takebuf_array(buf)
            max_block_size = max(max_block_size, length(block))

            if compressed
                writer = Libz.ZlibDeflateOutputStream(out)
                write(writer, block)
                flush(writer)
            else
                write(out, block)
            end

            bounds[section_ix] = BigBedBounds(block_offset, info.id,
                                              start_pos, end_pos)
        end
    end

    @assert section_ix == section_count "Incorrect number of sections when writing BigBed data"

    return max_block_size
end

function bigwig_write_blocks{T<:Number}(out::IO, intervals::IntervalCollection{T},
                                        chrom_info::Vector{BigBedChromInfo},
                                        items_per_slot, bounds::Vector{BigBedBounds},
                                        section_count, compressed::Bool,
                                        res_try_count, res_scales, res_sizes)
    section_ix = 0
    res_ends = zeros(Int, res_try_count)
    items = Array(BedGraphItem, items_per_slot)
    max_section_size = 0

    # See writeSections in BedGraphToBigWig.c
    for info in chrom_info
        tree = intervals.trees[info.name]
        it_state = start(tree)
        fill!(res_ends, 0)

        while !done(tree, it_state)
            section_ix += 1
            item_ix = 1

            while item_ix <= items_per_slot && !done(tree, it_state)
                interval, it_state = next(tree, it_state)
                interval_start = interval.first - 1
                interval_end = interval.last
                items[item_ix] = BedGraphItem(
                    interval_start, interval_end, interval.metadata)
                item_ix += 1

                for res_try in 1:res_try_count
                    res_end = res_ends[res_try]
                    if interval_start >= res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = interval_start + res_scales[res_try]
                    end
                    while interval_end > res_end
                        res_sizes[res_try] += 1
                        res_ends[res_try] = res_end = res_end + res_scales[res_try]
                    end
                end
            end

            # write a section
            chrom_id = info.id
            section_start = items[1].chrom_start
            section_end = items[end].chrom_end

            section_offset = position(out)
            header = BigWigSectionHeader(chrom_id, section_start, section_end,
                                         0, 0, BIGWIG_DATATYPE_BEDGRAPH, 0, item_ix - 1)

            if compressed
                writer = Libz.ZlibDeflateOutputStream(out)
                write(writer, header)
                for i in 1:item_ix-1
                    write(writer, items[i])
                end
                flush(writer)
            else
                write(out, header)
                for i in 1:item_ix-1
                    write(out, items[i])
                end
            end
            section_end_offset = position(out)
            max_section_size = max(max_section_size, section_end_offset - section_offset)

            bounds[section_ix] = BigBedBounds(section_offset, info.id,
                                              section_start, section_end)
        end
    end

    @assert section_ix == section_count "Incorrect number of sections when writing BigWig data"

    return max_section_size
end

type BigBedRTree
    next::Nullable{BigBedRTree}
    children::Nullable{BigBedRTree}
    parent::Nullable{BigBedRTree}
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    start_file_offset::UInt64
    end_file_offset::UInt64

    function BigBedRTree(next, children, parent, start_chrom_ix, start_base,
                         end_chrom_ix, end_base, start_file_offset, end_file_offset)
        return new(next, children, parent, start_chrom_ix, end_chrom_ix,
                   start_base, end_base, start_file_offset, end_file_offset)
    end

    function BigBedRTree()
        return new()
    end
end

function Base.copy(tree::BigBedRTree)
    return BigBedRTree(tree.next, tree.children, tree.parent,
                       tree.start_chrom_ix, tree.start_base,
                       tree.end_chrom_ix, tree.end_base,
                       tree.start_file_offset, tree.end_file_offset)
end

function Base.length(list::Nullable{BigBedRTree})
    count = 0
    while !isnull(list)
        count += 1
        list = get(list).next
    end
    return count
end

function bigbed_calc_rtree_level_sizes!(level_sizes, tree::Nullable{BigBedRTree},
                                        level, max_level)
    # See calcLevelSizes in cirTree.c

    el = tree
    while !isnull(el)
        level_sizes[level] += 1
        if level <= max_level
            bigbed_calc_rtree_level_sizes!(level_sizes, get(el).children,
                                           level + 1, max_level)
        end
        el = get(el).next
    end
end

function bigbed_rtree_write_index_level(out::IO, block_size, child_node_size,
                                        tree::BigBedRTree, offset_of_first_child,
                                        cur_level, dest_level)
    offset = offset_of_first_child
    if cur_level == dest_level
        # we've reached the right level, write out a node header
        count_one = length(tree.children)
        write(out, BigBedRTreeNode(0, 0, count_one))

        # write out elements of this node
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            write(out, BigBedRTreeInternalNode(el.start_chrom_ix, el.start_base,
                                               el.end_chrom_ix, el.end_base, offset))
            offset += child_node_size
            nel = el.next
        end

        # write out zeroes for empty slots in node
        for i in count_one:(block_size-1)
            for j in 1:sizeof(BigBedRTreeInternalNode)
                write(out, 0x00)
            end
        end
    else
        # otherwise recurse on children
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            offset = bigbed_rtree_write_index_level(out, block_size,
                                                    child_node_size, el,
                                                    cur_level + 1,
                                                    dest_level,
                                                    offset)
            nel = el.next
        end
    end
    return offset
end

function bigbed_rtree_write_leaf_level(out::IO, items_per_slot, leaf_node_size,
                                       tree::BigBedRTree, cur_level, leaf_level)
    if cur_level == leaf_level
        # we've reached the right level, write out node header
        count_one = length(tree.children)
        write(out, BigBedRTreeNode(1, 0, count_one))

        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            write(out, BigBedRTreeLeafNode(el.start_chrom_ix, el.start_base,
                                           el.end_chrom_ix, el.end_base,
                                           el.start_file_offset,
                                           el.end_file_offset - el.start_file_offset))
            nel = el.next
        end

        # write out zeroes for empty slots in node
        for i in count_one:(items_per_slot-1)
            for j in 1:sizeof(BigBedRTreeInternalNode)
                write(out, 0x00)
            end
        end
    else
        # otherwise recurse on children
        nel = tree.children
        while !isnull(nel)
            el = get(nel)
            bigbed_rtree_write_leaf_level(out, items_per_slot, leaf_node_size,
                                          el, cur_level + 1, leaf_level)
            nel = el.next
        end
    end
end

function Base.write(out::IO, tree::BigBedRTree, block_size, level_count)
    # See writeTreeToOpenFile in cirTree.c

    level_sizes = zeros(Int, level_count)
    bigbed_calc_rtree_level_sizes!(level_sizes, Nullable{BigBedRTree}(tree), 1, level_count)

    level_offsets = Array(UInt64, level_count)

    index_node_size = sizeof(BigBedRTreeNode) + block_size * sizeof(BigBedRTreeInternalNode)
    leaf_node_size = sizeof(BigBedRTreeNode) + block_size * sizeof(BigBedRTreeLeafNode)

    offset = position(out)
    for i in 1:level_count
        level_offsets[i] = offset
        offset += level_sizes[i] * index_node_size
    end

    final_level = level_count - 2
    for i in 1:final_level
        child_node_size = i == final_level ? leaf_node_size : index_node_size
        bigbed_rtree_write_index_level(out, block_size, child_node_size, tree,
                                       level_offsets[i+1], 1, i)
        @assert position(out) == level_offsets[i+1]
    end

    leaf_level = level_count - 1
    bigbed_rtree_write_leaf_level(out, block_size, leaf_node_size, tree,
                                  1, leaf_level)
end

function Base.reverse!(tree::Nullable{BigBedRTree})
    if isnull(tree)
        return tree
    end

    nodes = BigBedRTree[]
    node = tree
    while !isnull(node)
        push!(nodes, get(node))
        node = get(node).next
    end

    for i in length(nodes):-1:2
        nodes[i].next = nodes[i-1]
    end
    nodes[1].next = Nullable{BigBedRTree}()
    return Nullable{BigBedRTree}(nodes[end])
end

function bigbed_build_rtree(bounds::Vector{BigBedBounds}, block_size,
                            items_per_slot, end_file_offset)
    # See rTreeFromChromRangeArray in cirTree.c

    list  = Nullable{BigBedRTree}()
    next_offset = isempty(bounds) ? 0 : bounds[1].offset
    for i in 1:items_per_slot:length(bounds)
        final_iteration = false
        one_size = length(bounds) - i + 1
        if one_size > items_per_slot
            one_size = items_per_slot
        else
            final_iteration = true
        end

        bound = bounds[i]
        el = BigBedRTree()
        el.next = Nullable{BigBedRTree}()
        el.children = Nullable{BigBedRTree}()
        el.parent = Nullable{BigBedRTree}()
        el.start_chrom_ix = el.end_chrom_ix = bound.chrom_ix
        el.start_base = bound.start
        el.end_base = bound.stop
        el.start_file_offset = next_offset

        # Figure out end of element from offset of next element (or file size
        # for final element.)
        if final_iteration
            next_offset = end_file_offset
        else
            next_offset = bounds[i + one_size].offset
        end
        el.end_file_offset = next_offset

        # expand area spanned to include all items in block
        for j in 2:one_size
            bound = bounds[i + j - 1]
            if bound.chrom_ix < el.start_chrom_ix
                el.start_chrom_ix = bound.chrom_ix
                el.start_base = bound.start
            elseif bound.chrom_ix == el.start_chrom_ix
                el.start_base = min(el.start_base, bound.start)
            elseif bound.chrom_ix > el.end_chrom_ix
                el.end_chrom_ix = bound.chrom_ix
                el.end_base = bound.stop
            elseif bound.chrom_ix == el.end_chrom_ix
                el.end_base = max(el.end_base, bound.stop)
            end
        end

        el.next = list
        list = Nullable{BigBedRTree}(el)
    end
    list = reverse!(list)

    # Now iterate through making more and more condensed versions until have
    # just one
    tree = list
    level_count = 1
    while (!isnull(tree) && !isnull(get(tree).next)) || level_count < 2
        list = Nullable{BigBedRTree}()
        slots_used = block_size
        parent = Nullable{BigBedRTree}()
        nel = tree
        while !isnull(nel)
            el = get(nel)
            next_el = el.next
            if slots_used >= block_size
                slots_used = 1
                parent = Nullable{BigBedRTree}(copy(el))
                get(parent).children = el
                el.parent = parent
                el.next = Nullable{BigBedRTree}()

                get(parent).next = list
                list = parent
            else
                slots_used += 1
                el.next = get(parent).children
                get(parent).children = el
                el.parent = parent
                if el.start_chrom_ix < get(parent).start_chrom_ix
                    get(parent).start_chrom_ix = el.start_chrom_ix
                    get(parent).start_base = el.start_base
                elseif el.start_chrom_ix == get(parent).start_chrom_ix
                    get(parent).start_base = min(el.start_base, get(parent).start_base)
                end

                if el.end_chrom_ix > get(parent).end_chrom_ix
                    get(parent).end_chrom_ix = el.end_chrom_ix
                    get(parent).end_base = el.end_base
                elseif el.end_chrom_ix == get(parent).end_chrom_ix
                    get(parent).end_base  = max(el.end_base, get(parent).end_base)
                end
            end
            nel = next_el
        end

        list = reverse!(list)
        el = list
        while !isnull(el)
            get(el).children = reverse!(get(el).children)
            el = get(el).next
        end
        tree = list
        level_count += 1
    end

    @assert !isnull(tree)

    (get(tree), level_count)
end

function bigbed_write_index(out::IO, bounds::Vector{BigBedBounds},
                            block_size, items_per_slot, end_file_offset)
    # See cirTreeFileBulkIndexToOpenFile in cirTree.c

    rtree, level_count = bigbed_build_rtree(bounds, block_size, items_per_slot,
                                            end_file_offset)
    write(out, BigBedRTreeHeader(BIGBED_RTREE_MAGIC, block_size, length(bounds),
                                 rtree.start_chrom_ix, rtree.start_base,
                                 rtree.end_chrom_ix, rtree.end_base,
                                 end_file_offset, items_per_slot, 0))
    write(out, rtree, block_size, level_count)
end

# Write out sum to file, keeping track of minimal info on it in, and also adding
# it to second level summary
function output_one_summary_further_reduce(out::IO, sum::BigBedZoomData,
                                           twice_reduced_list::Vector{BigBedZoomData},
                                           double_reduction_size,
                                           bounds_array::Vector{BigBedBounds},
                                           bounds_array_idx,
                                           chrom_size)
    # See bbiOutputOneSummaryFurtherReduce in bbiWrite.c

    @assert bounds_array_idx <= length(bounds_array)

    bounds_array[bounds_array_idx] = BigBedBounds(position(out), sum.chrom_id,
                                                  sum.chrom_start, sum.chrom_end)
    write(out, sum)

    # Fold summary info into twice_reduced_list
    if isempty(twice_reduced_list) || twice_reduced_list[end].chrom_id != sum.chrom_id ||
       twice_reduced_list[end].chrom_start + double_reduction_size < sum.chrom_end
        push!(twice_reduced_list, sum)
    else
        twice_reduced_list[end] = BigBedZoomData(
            twice_reduced_list[end].chrom_id,
            twice_reduced_list[end].chrom_start,
            sum.chrom_end,
            twice_reduced_list[end].valid_count + sum.valid_count,
            min(sum.min_val, twice_reduced_list[end].min_val),
            max(sum.max_val, twice_reduced_list[end].max_val),
            twice_reduced_list[end].sum_data + sum.sum_data,
            twice_reduced_list[end].sum_squares + sum.sum_squares)
    end

    return bounds_array_idx + 1
end

# Write out data reduced by factor of initialReduction.  Also calculate and keep
# in memory next reduction level.  This is more work than some ways, but it
# keeps us from having to keep the first reduction entirely in memory.
function bigbed_write_reduced_once_return_reduced_twice(
                out::IO, intervals::IntervalCollection,
                interval_tree_transform::Function,
                chrom_info::Vector{BigBedChromInfo},
                initial_reduction, initial_reduced_count,
                zoom_increment, block_size, items_per_slot, compressed)

    # See writeReducedOnceReturnReducedTwice in bedToBigBed.c

    double_reduction_size = initial_reduction * zoom_increment
    bounds_array = Array(BigBedBounds, initial_reduced_count)
    bounds_array_idx = 1
    twice_reduced_list = BigBedZoomData[]

    ret_data_start = position(out)
    write(out, convert(UInt32, initial_reduced_count))

    # This gets a little complicated I'm afraid.  The strategy is to:
    #  1) Build up a range tree that represents coverage depth on that chromosome
    #     This also has the nice side effect of getting rid of overlaps.
    #  2) Stream through the range tree, outputting the initial summary level and
    #     further reducing.
    #

    valid_count = UInt64(0)
    min_val = 0.0
    max_val = 0.0
    sum_data = 0.0
    sum_squares = 0.0
    first_time = true

    sum = Nullable{BigBedZoomData}()

    for info in chrom_info
        sum = Nullable{BigBedZoomData}()
        for range in interval_tree_transform(intervals.trees[info.name])
            val = convert(Float64, range.metadata)
            chrom_start = range.first - 1
            chrom_end = range.last
            size = chrom_end - chrom_start
            @assert size >= 0 "Bad size: $chrom_start, $chrom_end"
            if first_time
                valid_count = size
                min_val = max_val = val
                sum_data = val * size
                sum_squares = val * val * size
                first_time = false
            else
                valid_count += size
                min_val = min(min_val, val)
                max_val = max(max_val, val)
                sum_data += val * size
                sum_squares += val * val * size
            end

            # If start past existing block then output it
            if !isnull(sum) && get(sum).chrom_end <= chrom_start
                bounds_array_idx = output_one_summary_further_reduce(
                        out, get(sum), twice_reduced_list,
                        double_reduction_size, bounds_array,
                        bounds_array_idx, info.size)
                sum = Nullable{BigBedZoomData}()
            end

            # If don't have a summary we're working on now, make one
            if isnull(sum)
                sum = Nullable{BigBedZoomData}(
                    BigBedZoomData(info.id, chrom_start,
                                   min(info.size, chrom_start + initial_reduction),
                                   0, val, val, 0.0, 0.0))
            end

            # Deal with case where might have to split an item between multiple
            # summaries. This loop handles all but the final affected summary in
            # that case
            s = get(sum)
            while chrom_end > s.chrom_end
                # Fold in bits that ovelap with existing summary and output
                overlap = min(chrom_end, s.chrom_end) - max(chrom_start, s.chrom_start)
                @assert overlap > 0
                @assert overlap <= size

                s = BigBedZoomData(
                        s.chrom_id, s.chrom_start, s.chrom_end,
                        s.valid_count + overlap,
                        min(s.min_val, val), max(s.max_val, val),
                        s.sum_data + val * overlap,
                        s.sum_squares + val * val * overlap)

                bounds_array_idx = output_one_summary_further_reduce(
                        out, s, twice_reduced_list, double_reduction_size,
                        bounds_array, bounds_array_idx, info.size)
                size -= overlap

                # Move summary to next part
                chrom_start = s.chrom_end
                #chrom_end = min(info.size, chrom_start + initial_reduction)
                s = BigBedZoomData(
                        s.chrom_id, chrom_start,
                        min(info.size, chrom_start + initial_reduction),
                        0, val, val, 0.0, 0.0)
            end

            # Add to summary
            s = BigBedZoomData(
                    s.chrom_id, s.chrom_start, s.chrom_end,
                    s.valid_count + size,
                    min(s.min_val, val), max(s.max_val, val),
                    s.sum_data + val * size,
                    s.sum_data + val * val * size)
            sum = Nullable(s)
        end

        if !isnull(sum)
            bounds_array_idx = output_one_summary_further_reduce(
                    out, get(sum), twice_reduced_list, double_reduction_size,
                    bounds_array, bounds_array_idx, info.size)
        end
    end

    # Write out 1st zoom index
    index_offset = position(out)
    @assert bounds_array_idx == length(bounds_array) + 1
    bigbed_write_index(out, bounds_array, block_size, items_per_slot,
                       index_offset)
    return twice_reduced_list, BigBedTotalSummary(valid_count, min_val, max_val,
                                                  sum_data, sum_squares)
end

function bigbed_write_summary_and_index(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot, compressed)
    # See: bbiWriteSummaryAndIndex in bbiWrite.c
    if compressed
        bigbed_write_summary_and_index_comp(out, summary_list, block_size, items_per_slot)
    else
        bigbed_write_summary_and_index_unc(out, summary_list, block_size, items_per_slot)
    end
end

function bigbed_write_summary_and_index_unc(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot)
    # See: bbiWriteSummaryAndIndexUnc in bbiWrite.c
    count = length(summary_list)
    write(out, convert(UInt32, count))
    bounds_array = Array(BigBedBounds, count)
    for (i, summary) in enumerate(summary_list)
        bounds_array[i] = BigBedBounds(
            position(out), summary.chrom_id, summary.chrom_start,
            summary.chrom_end)
        write(out, summary)
    end
    index_offset = position(out)
    bigbed_write_index(out, bounds_array, block_size, items_per_slot, false)
    return index_offset
end

function bigbed_write_summary_and_index_comp(
        out::IO, summary_list::Vector{BigBedZoomData}, block_size,
        items_per_slot)
    # See: bbiWriteSummaryAndIndexComp in bbiWrite.c
    count = length(summary_list)
    write(out, convert(UInt32, count))
    bounds_array = Array(BigBedBounds, count)

    items_left = count
    sum_ix = 1
    i = 1
    while items_left > 0
        items_in_slot = min(items_per_slot, items_left)
        writer = Libz.ZlibDeflateOutputStream(out)
        file_pos = position(out)
        for _ in 1:items_in_slot
            summary = summary_list[sum_ix]
            bounds_array[sum_ix] = BigBedBounds(file_pos, summary.chrom_id,
                                                summary.chrom_start, summary.chrom_end)
            write(writer, summary_list[sum_ix])
            sum_ix += 1
            if sum_ix > count
                break
            end
        end
        items_left -= items_in_slot
        flush(writer)
    end
    index_offset = position(out)
    bigbed_write_index(out, bounds_array, block_size, items_per_slot, false)
    return index_offset
end

function bigbed_summary_simple_reduce(summary_list::Vector{BigBedZoomData}, reduction)
    new_summary_list = BigBedZoomData[]
    for summary in summary_list
        if isempty(new_summary_list) ||
           new_summary_list[end].chrom_id != summary.chrom_id ||
           summary.chrom_end > new_summary_list[end].chrom_start + reduction
            push!(new_summary_list, summary)
        else
            @assert new_summary_list[end].chrom_end < summary.chrom_end
            new_summary_list[end] = BigBedZoomData(
                new_summary_list[end].chrom_id,
                new_summary_list[end].chrom_start,
                summary.chrom_end,
                new_summary_list[end].valid_count + summary.valid_count,
                min(new_summary_list[end].min_val, summary.min_val),
                max(new_summary_list[end].max_val, summary.max_val),
                new_summary_list[end].sum_data + summary.sum_data,
                new_summary_list[end].sum_squares + summary.sum_squares)
       end
    end
    return new_summary_list
end

function bed_field_count(intervals::IntervalCollection{BEDMetadata})
    num_fields = 0
    for interval in intervals
        interval_num_fields = 3 + interval.metadata.used_fields
        if num_fields == 0 || num_fields > interval_num_fields
            num_fields = interval_num_fields
        end
    end
    return num_fields
end

# Unified write function for bigwig and bigbed
function write_bigbed_bigwig(
        out::IO,
        fmt::Symbol,
        intervals::IntervalCollection,
        block_size::Int,
        items_per_slot::Int,
        compressed::Bool)

    # See the function bbFileCreate in bedToBigBed.c in kent to see the only
    # other implementation of this I'm aware of. Also bedGraphToBigWig in
    # bedGraphToBigWig.c
    @assert fmt ∈ (:bigwig, :bigbed)

    if !applicable(seek, out, 1)
        error("BigBed can only be written to seekable output streams.")
    end

    field_count = fmt == :bigbed ? bed_field_count(intervals) : 0

    # write dummy headers, come back and fill them later
    write_zeros(out, sizeof(BigBedHeader))
    write_zeros(out, BIGBED_MAX_ZOOM_LEVELS * sizeof(BigBedZoomHeader))

    # TODO: optional autoSql specification
    as_offset = 0

    total_summary_offset = position(out)
    total_summary = BigBedTotalSummary(0, 0.0, 0.0, 0.0, 0.0)
    write(out, total_summary)

    chrom_tree_offset = position(out)
    chrom_info = bigbed_write_chrom_tree(out, intervals, block_size)

    # set up to keep track of possible initial reduction levels
    total_bases = 0
    for interval in intervals
        total_bases += interval.last - interval.first + 1
    end
    ave_span = total_bases / length(intervals)

    res_try_count = 10
    res_increment = 4
    res_scales = Array(Int, res_try_count)
    res_sizes = Array(Int, res_try_count)
    min_zoom = 10
    res = max(round(Int, ave_span), min_zoom)
    for res_try in 1:res_try_count
        res_sizes[res_try] = 0
        res_scales[res_try] = res
        if res > 1000000000
            res_try_count = res_try
            break
        end
        res *= res_increment
    end

    # write out primary full resolution data in sections, collect stats to
    # use for reductions
    data_offset = position(out)
    # TODO: The BigBed/BigWig specification differs from Kent's impementation
    # here. The spec says this count is a 32bit int, but the implmentation uses
    # a 64bit int. We follow the implementation for practical reasons.
    write(out, convert(UInt64, length(intervals)))
    block_count = bigbed_count_sections_needed(chrom_info, items_per_slot)
    bounds = Array(BigBedBounds, block_count)

    if fmt == :bigbed
        max_block_size = bigbed_write_blocks(out, intervals, chrom_info, items_per_slot,
                                             bounds, block_count, compressed,
                                             res_try_count, res_scales, res_sizes)
    else
        max_block_size = bigwig_write_blocks(out, intervals, chrom_info, items_per_slot,
                                             bounds, block_count, compressed,
                                             res_try_count, res_scales, res_sizes)
    end

    # write out primary data index
    index_offset = position(out)
    bigbed_write_index(out, bounds, block_size, 1, index_offset)

    # declare arrays and vars that track the zoom levels we actually output
    zoom_levels = 0
    zoom_amounts = Array(UInt32, BIGBED_MAX_ZOOM_LEVELS)
    zoom_data_offsets = Array(UInt64, BIGBED_MAX_ZOOM_LEVELS)
    zoom_index_offsets = Array(UInt64, BIGBED_MAX_ZOOM_LEVELS)
    zoom_levels = 0

    if ave_span > 0
        data_size = index_offset - data_offset
        max_reduced_size = div(data_size, 2)
        initial_reduction = 0
        initial_reduced_count = 0

        for res_try in 1:res_try_count
            reduced_size = res_sizes[res_try] * sizeof(BigBedZoomData)
            if compressed
                reduced_size = div(reduced_size, 2) # estimate!
            end

            if reduced_size <= max_reduced_size
                initial_reduction = res_scales[res_try]
                initial_reduced_count = res_sizes[res_try]
                break
            end
        end

        if initial_reduction > 0
            zoom_increment = 4
            rezoomed_list, total_summary = bigbed_write_reduced_once_return_reduced_twice(
                    out, intervals,
                    fmt == :bigbed ? coverage : identity,
                    chrom_info, initial_reduction,
                    initial_reduced_count, zoom_increment, block_size,
                    items_per_slot, compressed)

            zoom_amounts[1] = initial_reduction
            zoom_levels = 1

            zoom_count = initial_reduced_count
            reduction = initial_reduction * zoom_increment
            while zoom_levels < BIGBED_MAX_ZOOM_LEVELS
                rezoom_count = length(rezoomed_list)
                if rezoom_count >= zoom_count
                    break
                end

                zoom_count = rezoom_count
                zoom_data_offsets[zoom_levels] = position(out)
                zoom_index_offsets[zoom_levels] =
                        bigbed_write_summary_and_index(out, rezoomed_list, block_size,
                                                       items_per_slot, compressed)
                zoom_amounts[zoom_levels] = reduction
                zoom_levels += 1
                reduction *= zoom_increment

                rezoomed_list = bigbed_summary_simple_reduce(rezoomed_list, reduction)
            end
        end
    end

    uncompress_buf_size = 0
    if compressed
        max_zoom_uncomp_size = items_per_slot * sizeof(BigBedZoomData)
        uncompress_buf_size = max(max_block_size, max_zoom_uncomp_size)
    end

    # go back and rewrite header
    seek(out, 0)
    write(out,
        BigBedHeader(fmt == :bigbed ? BIGBED_MAGIC : BIGWIG_MAGIC,
                     BIGBED_CURRENT_VERSION, zoom_levels,
                     chrom_tree_offset, data_offset, index_offset,
                     field_count, min(field_count, 12), as_offset,
                     total_summary_offset, uncompress_buf_size, 0))

    # write summary headers with data
    for i in 1:zoom_levels
        write(out, BigBedZoomHeader(zoom_amounts[i], 0, zoom_data_offsets[i],
                                    zoom_index_offsets[i]))
    end

    # write rest of summary header with no data
    for i in zoom_levels:(BIGBED_MAX_ZOOM_LEVELS-1)
        write(out, BigBedZoomHeader(0, 0, 0, 0))
    end

    # write total summary
    seek(out, total_summary_offset)
    write(out, total_summary)

    # write end signature
    seekend(out)
    write(out, convert(UInt32, fmt == :bigbed ? BIGBED_MAGIC : BIGWIG_MAGIC))
end

function write_zeros(out::IO, n::Integer)
    for i in 1:n
        write(out, 0x00)
    end
end


text:
using DelimitedFiles

puzzleinput = readdlm("inputs/day08.txt")

mutable struct CPU
    pc::Int
    a::Int
end

program = hcat(puzzleinput, fill(false, (size(puzzleinput)[1], 1)))
# The program is a 3 cols array where col 1 is the ops, col 2 the arg
# and col 3 a flag to check if the instruction has been executed

function execute(cpu, prg)
    eop = length(prg[:, 1:1])
    # Address of the last line of the program

    while true
        # Execution loop

        pc = cpu.pc
        if pc == eop + 1 # Program ended correctly
            return cpu
        elseif pc > eop + 1 # Programs jumped to a wrong address
            cpu.pc = 0 # Error code for bad address
            return cpu
        end

        ops, arg, flg = prg[pc, 1:3]

        if flg
            # If the code is going to execute an instruction already executed
            # the code is in an infinite loop then exit execution loop

            cpu.pc = -1 # Error code for infinite loop
            return cpu
        end

        if ops == "acc"
            cpu.a += arg
            cpu.pc += 1
        elseif ops == "jmp"
            cpu.pc += arg
        elseif ops == "nop"
            cpu.pc += 1
        else
            error("Unknown ops")
        end
        prg[pc, 3] = true
        # After execution the instrucion in the program is marked as executed

    end
end

cpu = execute(CPU(1, 0), program)
cpu.a

# Second Half

function fix(prg)
    eop = length(prg[:, 1:1])
    # Address of the last line of the program

    nops = [i for i = 1:eop if program[i, 1] == "nop" && program[i, 2] != 0]
    # Addresses of lines of program with ops = nop and arg != 0

    jmps = [i for i = 1:eop if program[i, 1] == "jmp"]
    # Addresses of lines of program with ops = jmp

    for i in nops
        # Try to change every nop ops into a jmp to fix the program

        newprg = copy(prg)
        newprg[i, 1] = "jmp"
        cpu = execute(CPU(1, 0), newprg)
        if cpu.pc == eop + 1 # The program ended correctly, we have fixed it!
            return cpu.a
        end
    end

    for i in jmps
        # Try to change every nop ops into a jmp to fix the program

        newprg = copy(prg)
        newprg[i, 1] = "nop"
        cpu = execute(CPU(1, 0), newprg)
        if cpu.pc == eop + 1 # The program ended correctly, we have fixed it!
            return cpu.a
        end
    end
    0
end

fix(program)


text:
using Constraints
using CompositionalNetworks
using ConstraintDomains
using Documenter

makedocs(;
    modules=[Constraints, ConstraintDomains, CompositionalNetworks],
    authors="Jean-François Baffier",
    repo="https://github.com/JuliaConstraints/Constraints.jl/blob/{commit}{path}#L{line}",
    sitename="Constraints.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://JuliaConstraints.github.io/Constraints.jl",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "Learning (ICN)" => "learning.md",
        "Dependencies" => [
            "ConstraintDomains.jl" => "domain.md",
            "CompositionalNetworks.jl" => "icn.md",
        ],
        "Library" => [
            "Public" => "public.md",
            "Internal" => "internal.md"
        ],
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/Constraints.jl.git",
    devbranch="main",
)


text:
# Configuration for nearly neutral simuation
function dfe_disadv( x::Int64 )
  dfe_deleterious( x::Int64; alpha=dfe_disadv_alpha, beta=dfe_disadv_beta )
end

const nn_simtype = 2    # nearly neutral infinite sites
const type_str = "D3"
const mu_list_flag = true
@everywhere const N_list     = [25,50,100,200,400,800,1600,3200,6400]     # popsize
const mu_list= [0.025,0.01,0.005,0.0025]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 50000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_disadv
dfe_str = "dfe_disadv alpha_disadv=0.2 beta_disadv=0.5"
dfe_disadv_alpha=0.2
dfe_disadv_beta=0.5

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using Iterators
using Combinatorics

export coalitions

type Coalitions{T}
    members::Array{T,1}
    k::Int64
    lastk::Int64
    kIter::Iterators.Binomial{T}
    invSet::Array{T,1}
end

Base.eltype(it::Coalitions) = Array{eltype(it.members),1}
function Base.length(it::Coalitions)
    length(it.members) < 40 ? 2^length(it.members) : Inf
end

"Enumerate all the coalitions with a weight for computing the Shapley value."
coalitions{T}(members::Array{T,1}) = Coalitions(members, 0, 0, subsets(members, 0), T[])

function Base.start(it::Coalitions)
    start(it.kIter)
end

function coalition_weight(q, s)
    if q > s-q-1
        factorial(s-q-1)/factorial(s,q)
    else
        factorial(q)/factorial(s,s-q-1)
    end
end

function Base.next{T}(it::Coalitions{T}, state::(Tuple{Array{Int64,1}, Bool}))
    if length(it.invSet) > 0
        nextSet = it.invSet
        it.invSet = T[]

        # check if this was a final inv set to emit
        if length(state[1]) == 0
            state = (Int64[],true)
        end
    else
        nextSet,state = next(it.kIter, state)
        it.lastk = it.k

        if 2*it.k < length(it.members)

            # we also create an inverse set if we are not
            it.invSet = setdiff(it.members, nextSet)

            # if we are done with current set of sizes move to the next one
            if state[2] && 2*it.k < length(it.members)-1
                it.k += 1
                it.kIter = subsets(it.members, it.k)
                state = start(it)
            end

            # prevent us from ending if we still have a inv set to emit
            if state[2]
                state = (Int64[],false)
            end
        end
    end

    weight = coalition_weight(it.lastk, length(it.members)+1)

    (nextSet,weight),state
end

Base.done(it::Coalitions, state::Tuple{Array{Int64,1}, Bool}) = state[2]


text:
using Test

include("trinary.jl")

@testset "trinary 1 is decimal 1" begin
    @test trinary_to_decimal("1") == 1
end

@testset "trinary 2 is decimal 2" begin
    @test trinary_to_decimal("2") == 2
end

@testset "trinary 10 is decimal 3" begin
    @test trinary_to_decimal("10") == 3
end

@testset "trinary 11 is decimal 4" begin
    @test trinary_to_decimal("11") == 4
end

@testset "trinary 100 is decimal 9" begin
    @test trinary_to_decimal("100") == 9
end

@testset "trinary 112 is decimal 14" begin
    @test trinary_to_decimal("112") == 14
end

@testset "trinary 222 is decimal 26" begin
    @test trinary_to_decimal("222") == 26
end

@testset "trinary 1122000120 is decimal 32091" begin
    @test trinary_to_decimal("1122000120") == 32091
end

@testset "invalid trinary digits returns 0" begin
    @test trinary_to_decimal("1234") == 0
end

@testset "invalid word as input returns 0" begin
    @test trinary_to_decimal("carrot") == 0
end

@testset "invalid numbers with letters as input returns 0" begin
    @test trinary_to_decimal("0a1b2c") == 0
end


text:
using ApproxFun, Test

@testset "Complex" begin
    @testset "Differentiation" begin
        f=Fun(x->exp(im*x))
        @test norm(f'-im*f) < ncoefficients(f)*100*eps()
        @test norm(integrate(f)+im*f-f.coefficients[1]*im) < 100eps()
        @test norm(real(f)-Fun(cos)) < 10eps()
        @test norm(real(f-Fun(cos))) < 10eps()
    end

    @testset "Other intervals" begin
        f=Fun(x->exp(im.*x),1..2)
        @test norm(f'-im*f) < 1000eps()
        @test norm(integrate(f)+im*f-f.coefficients[1]*im) < 100eps()
        @test norm(real(f)-Fun(cos,domain(f))) < eps()
        @test norm(real(f-Fun(cos,domain(f)))) < eps()
    end

    @testset "Complex segments" begin
        f=Fun(x->exp(im.*x),Segment(1im,2+.5im))
        @test sum(f) ≈ (0.5515167681675808 + 0.6202852564797062im)
        @test f(1im) ≈ exp(im.*im)
        @test f(1+.75im) ≈ exp(im.*(1+.75im))
        @test norm(f'-im*f) < 1000eps()
        @test norm(integrate(f)+im*f-f.coefficients[1]*im) < 100eps()
    end
end


text:
function l1_norm(a)
    return reduce(+, map(x-> abs(x),a))
end

function l2_norm(a)
    return sqrt(reduce(+, map(x -> x^2,a)))
end

function lp_norm(a,power)
    return (reduce(+, map(x -> abs(x)^power,a)))^(1/power)
end


text:
__precompile__()

module DiffEqFinancial
  using DiffEqBase, DiffEqNoiseProcess, Markdown, LinearAlgebra

  import RandomNumbers: Xorshifts

  include("problems.jl")

  export HestonProblem, BlackScholesProblem, GeneralizedBlackScholesProblem,
  ExtendedOrnsteinUhlenbeckProblem, OrnsteinUhlenbeckProblem, GeometricBrownianMotionProblem,
  MfStateProblem

end # module

max_stars_repo_name: JuliaDiffEq/FinancialModels

max_issues_repo_name: JuliaDiffEq/FinancialModels

max_forks_repo_name: JuliaDiffEq/FinancialModels


text:
"""
    chainscat(c::AbstractChains...)

Concatenate multiple chains.

By default, the chains are concatenated along the third dimension by calling
`cat(c...; dims=3)`.
"""
chainscat(c::AbstractChains...) = cat(c...; dims=3)

"""
    chainsstack(c::AbstractVector)

Stack chains in `c`.

By default, the vector of chains is returned unmodified. If `eltype(c) <: AbstractChains`,
then `reduce(chainscat, c)` is called.
"""
chainsstack(c) = c
chainsstack(c::AbstractVector{<:AbstractChains}) = reduce(chainscat, c)

"""
    bundle_samples(samples, model, sampler, state, chain_type[; kwargs...])

Bundle all `samples` that were sampled from the `model` with the given `sampler` in a chain.

The final `state` of the `sampler` can be included in the chain. The type of the chain can
be specified with the `chain_type` argument.

By default, this method returns `samples`.
"""
function bundle_samples(
    samples,
    ::AbstractModel,
    ::AbstractSampler,
    ::Any,
    ::Type;
    kwargs...
)
    return samples
end

function bundle_samples(
    samples::Vector,
    ::AbstractModel,
    ::AbstractSampler,
    ::Any,
    ::Type{Vector{T}};
    kwargs...
) where T
    return map(samples) do sample
        convert(T, sample)
    end
end

"""
    step(rng, model, sampler[, state; kwargs...])

Return a 2-tuple of the next sample and the next state of the MCMC `sampler` for `model`.

Samples describe the results of a single step of the `sampler`. As an example, a sample
might include a vector of parameters sampled from a prior distribution.

When sampling using [`sample`](@ref), every `step` call after the first has access to the
current `state` of the sampler.
"""
function step end

"""
    samples(sample, model, sampler[, N; kwargs...])

Generate a container for the samples of the MCMC `sampler` for the `model`, whose first
sample is `sample`.

The method can be called with and without a predefined number `N` of samples.
"""
function samples(
    sample,
    ::AbstractModel,
    ::AbstractSampler,
    N::Integer;
    kwargs...
)
    ts = Vector{typeof(sample)}(undef, 0)
    sizehint!(ts, N)
    return ts
end

function samples(
    sample,
    ::AbstractModel,
    ::AbstractSampler;
    kwargs...
)
    return Vector{typeof(sample)}(undef, 0)
end

"""
    save!!(samples, sample, iteration, model, sampler[, N; kwargs...])

Save the `sample` of the MCMC `sampler` at the current `iteration` in the container of
`samples`.

The function can be called with and without a predefined number `N` of samples. By default,
AbstractMCMC uses `push!!` from the Julia package
[BangBang](https://github.com/tkf/BangBang.jl) to append to the container, and widen its
type if needed.
"""
function save!!(
    samples::Vector,
    sample,
    iteration::Integer,
    ::AbstractModel,
    ::AbstractSampler,
    N::Integer;
    kwargs...
)
    s = BangBang.push!!(samples, sample)
    s !== samples && sizehint!(s, N)
    return s
end

function save!!(
    samples,
    sample,
    iteration::Integer,
    ::AbstractModel,
    ::AbstractSampler;
    kwargs...
)
    return BangBang.push!!(samples, sample)
end

# Deprecations
Base.@deprecate transitions(
    transition,
    model::AbstractModel,
    sampler::AbstractSampler,
    N::Integer;
    kwargs...
) samples(transition, model, sampler, N; kwargs...) false
Base.@deprecate transitions(
    transition,
    model::AbstractModel,
    sampler::AbstractSampler;
    kwargs...
) samples(transition, model, sampler; kwargs...) false