text:
## AlmostBandedMatrix

abstract type AbstractAlmostBandedLayout <: MemoryLayout end
struct AlmostBandedLayout <: AbstractAlmostBandedLayout end

struct AlmostBandedMatrix{T,D,A,B,R} <: LayoutMatrix{T}
    bands::BandedMatrix{T,D,R}
    fill::LowRankMatrix{T,A,B}
    AlmostBandedMatrix{T,D,A,B,R}(bands, fill) where {T,D,A,B,R} = new{T,D,A,B,R}(bands,fill)
end

function AlmostBandedMatrix{T}(bands::BandedMatrix{T,D,R}, fill::LowRankMatrix{T,A,B}) where {T,D,A,B,R}
    if size(bands) ≠ size(fill)
        error("Data and fill must be compatible size")
    end
    AlmostBandedMatrix{T,D,A,B,R}(bands,fill)
end

AlmostBandedMatrix(bands::BandedMatrix, fill::LowRankMatrix) =
    AlmostBandedMatrix{promote_type(eltype(bands),eltype(fill))}(bands,fill)

AlmostBandedMatrix(bands::AbstractMatrix, fill::AbstractMatrix) =
    AlmostBandedMatrix(BandedMatrix(bands), LowRankMatrix(fill))

AlmostBandedMatrix{T}(::UndefInitializer, nm::NTuple{2,Integer}, lu::NTuple{2,Integer}, r::Integer) where {T} =
    AlmostBandedMatrix(BandedMatrix{T}(undef,nm,lu), LowRankMatrix{T}(undef,nm,r))

function AlmostBandedMatrix{T}(A::AlmostBandedMatrix, (l,u)::NTuple{2,Integer}) where T
    B_in,L = bandpart(A),fillpart(A)
    l_in,u_in = bandwidths(B_in)
    B = BandedMatrix(B_in, (l,u))
    for b = u_in+1:u
        B[band(b)] .= L[band(b)]
    end
    AlmostBandedMatrix{T}(B, deepcopy(L))
end

AlmostBandedMatrix(A::AlmostBandedMatrix{T}, (l,u)::NTuple{2,Integer}) where T = AlmostBandedMatrix{T}(A, (l,u))

AlmostBandedMatrix{T}(Z::Zeros, lu::NTuple{2,Integer}, r::Integer) where {T} =
    AlmostBandedMatrix(BandedMatrix{T}(Z, lu), LowRankMatrix{T}(Z, r))

AlmostBandedMatrix{T}(A::AbstractMatrix, lu::NTuple{2,Integer}=almostbandwidths(A), r::Integer=almostbandedrank(A)) where T =
    copyto!(AlmostBandedMatrix{T}(undef, size(A), lu, r), A)

AlmostBandedMatrix(A::AbstractMatrix{T}, lu::NTuple{2,Integer}=almostbandwidths(A), r::Integer=almostbandedrank(A)) where T =
    AlmostBandedMatrix{T}(A, lu, r)

MemoryLayout(::Type{<:AlmostBandedMatrix}) = AlmostBandedLayout()

for MAT in (:AlmostBandedMatrix, :AbstractMatrix, :AbstractArray)
    @eval convert(::Type{$MAT{T}}, A::AlmostBandedMatrix) where {T} =
        AlmostBandedMatrix(AbstractMatrix{T}(A.bands),AbstractMatrix{T}(A.fill))
end

bandpart(A::AlmostBandedMatrix) = A.bands
fillpart(A::AlmostBandedMatrix) = A.fill

almostbandwidths(_, A) = bandwidths(bandpart(A))
almostbandedrank(_, A) = separablerank(fillpart(A))
almostbandwidths(A) = almostbandwidths(MemoryLayout(typeof(A)), A)
almostbandedrank(A) = almostbandedrank(MemoryLayout(typeof(A)), A)

size(A::AlmostBandedMatrix) = size(A.bands)
Base.IndexStyle(::Type{ABM}) where {ABM<:AlmostBandedMatrix} =
    IndexCartesian()

function colsupport(::AbstractAlmostBandedLayout, A, j)
    l,_ = almostbandwidths(A)
    Base.OneTo(min(maximum(j)+l,size(A,1)))
end

function rowsupport(::AbstractAlmostBandedLayout, A, k)
    l,_ = almostbandwidths(A)
    max(1,minimum(k)-l):size(A,2)
end

function getindex(B::AlmostBandedMatrix, k::Integer, j::Integer)
    if j > k + bandwidth(B.bands,2)
        B.fill[k,j]
    else
        B.bands[k,j]
    end
end

# can only change the bands, not the fill
function setindex!(B::AlmostBandedMatrix, v, k::Integer, j::Integer)
    l,u = bandwidths(bandpart(B))
    if j-k ≤ u
        B.bands[k,j] = v
    else
        B.fill[k,j] = v
    end
    B
end

function triu!(A::AlmostBandedMatrix)
    triu!(bandpart(A))
    A
end

###
# SubArray
###

sublayout(::AlmostBandedLayout, ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}) =
    AlmostBandedLayout()

sub_materialize(::AbstractAlmostBandedLayout, V) = AlmostBandedMatrix(V)

bandpart(V::SubArray) = view(bandpart(parent(V)), parentindices(V)...)
fillpart(V::SubArray) = view(fillpart(parent(V)), parentindices(V)...)

###
# QR
##

_factorize(::AbstractAlmostBandedLayout, _, A) = qr(A)
_qr(::AbstractAlmostBandedLayout, _, A) = almostbanded_qr(A)
_qr!(::AlmostBandedLayout, _, A) = almostbanded_qr!(A)

almostbanded_qr(A) = _almostbanded_qr(axes(A), A)
function _almostbanded_qr(_, A)
    l,u = almostbandwidths(A)
    qr!(AlmostBandedMatrix{float(eltype(A))}(A, (l,l+u)))
end

function _almostbanded_qr!(A::AbstractMatrix{T}, τ::AbstractVector{T}) where T
    m,n = size(A)
    _almostbanded_qr!(A, τ, min(m - 1 + !(T<:Real), n))
end
function _almostbanded_qr!(A::AbstractMatrix , τ::AbstractVector, ncols)
    B,L = bandpart(A),fillpart(A)
    l,u = bandwidths(B)
    m,n = size(A)
    k = 1
    while k ≤ ncols
        kr = k:min(k+l+u,m)
        jr1 = k:min(k+u,n)
        jr2 = k+u+1:min(kr[end]+u,n)
        jr3 = k:min(k+u,n,ncols)
        S = view(B,kr,jr1)
        τv = view(τ,jr3)

        R,_ = _banded_qr!(S, τv, length(jr3))
        Q = QRPackedQ(R,τv)

        B_right = view(B,kr,jr2)
        L_right = view(L,kr,jr2)
        # The following writes it as Q' * (B -L) + Q'*L
        # that is,  subtract out L from B, apply Q' to both
        # and add it back in again
        for j=1:size(B_right,2)
            B_right[j+1:end,j] -= L_right[j+1:end,j]
        end
        lmul!(Q', B_right)
        U,V = arguments(L_right)
        lmul!(Q', U)
        for j=1:size(B_right,2)
            B_right[j+1:end,j] += L_right[j+1:end,j]
        end
        k = jr1[end]+1
    end
    A, τ
end

function almostbanded_qr!(R::AbstractMatrix{T}, τ) where T
    _almostbanded_qr!(R, τ)
    QR(R, τ)
 end

almostbanded_qr!(R::AbstractMatrix{T}) where T = almostbanded_qr!(R, zeros(T, min(size(R)...)))

getQ(F::QR{<:Any,<:AlmostBandedMatrix}) = LinearAlgebra.QRPackedQ(bandpart(F.factors), F.τ)
function getR(F::QR{<:Any,<:AlmostBandedMatrix})
    n = min(size(F.factors,1),size(F.factors,2))
    UpperTriangular(view(F.factors,1:n,1:n))
end

function ldiv!(A::QR{<:Any,<:AlmostBandedMatrix}, B::AbstractVector)
    R = A.factors
    lmul!(adjoint(A.Q), B)
    B .= Ldiv(UpperTriangular(R), B)
    B
end

function ldiv!(A::QR{<:Any,<:AlmostBandedMatrix}, B::AbstractMatrix)
    R = A.factors
    lmul!(adjoint(A.Q), B)
    B .= Ldiv(UpperTriangular(R), B)
    B
end

# needed for adaptive QR
materialize!(M::Lmul{<:QRPackedQLayout{AlmostBandedLayout}}) =
    lmul!(QRPackedQ(bandpart(M.A.factors), M.A.τ), M.B)
function materialize!(M::Lmul{<:AdjQRPackedQLayout{AlmostBandedLayout}})
    Q = M.A'
    lmul!(QRPackedQ(bandpart(Q.factors), Q.τ)', M.B)
end

###
# UpperTriangular ldiv!
##

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:AlmostBandedLayout} = Tri{ML}()

function _almostbanded_upper_ldiv!(Tri, R::AbstractMatrix, b::AbstractVector{T}, buffer) where T
    B = bandpart(R)
    L = fillpart(R)
    U,V = arguments(L)
    fill!(buffer, zero(T))

    l,u = bandwidths(B)
    k = n = size(R,2)

    while k > 0
        kr = max(1,k-u):k
        jr1 = k+1:k+u+1
        jr2 = k+u+2:k+2u+2
        bv = view(b,kr)
        if jr2[1] < n
            muladd!(one(T),view(V,:,jr2),view(b,jr2),one(T),buffer)
            muladd!(-one(T),view(U,kr,:),buffer,one(T),bv)
        end
        if jr1[1] < n
            muladd!(-one(T),view(R,kr,jr1),view(b,jr1),one(T),bv)
        end
        materialize!(Ldiv(Tri(view(B,kr,kr)), bv))
        k = kr[1]-1
    end
    b
end

@inline function materialize!(M::MatLdivVec{TriangularLayout{'U','N',AlmostBandedLayout}})
    R,x = M.A,M.B
    A = triangulardata(R)
    r = size(arguments(fillpart(A))[1],2)
    _almostbanded_upper_ldiv!(UpperTriangular, A, x, Vector{eltype(M)}(undef, r))
    x
end

@inline function materialize!(M::MatLdivVec{TriangularLayout{'U','U',AlmostBandedLayout}})
    R,x = M.A,M.B
    A = triangulardata(R)
    r = size(arguments(fillpart(A))[1],2)
    _almostbanded_upper_ldiv!(UnitUpperTriangular, A, x, Vector{eltype(M)}(undef, r))
    x
end

@inline function materialize!(M::MatLdivVec{TriangularLayout{'L','N',AlmostBandedLayout}})
    R,x = M.A,M.B
    A = triangulardata(R)
    materialize!(Ldiv(LowerTriangular(bandpart(A)),x))
    x
end

@inline function materialize!(M::MatLdivVec{TriangularLayout{'L','U',AlmostBandedLayout}})
    R,x = M.A,M.B
    A = triangulardata(R)
    materialize!(Ldiv(UnitLowerTriangular(bandpart(A)),x))
    x
end

###
# VcatBanded
###

struct VcatAlmostBandedLayout <: AbstractAlmostBandedLayout end

# avoid ambuities
applylayout(::Type{typeof(vcat)}, ::ScalarOrBandedLayouts, ::AbstractBandedLayout) = ApplyBandedLayout{typeof(vcat)}()
applylayout(::Type{typeof(vcat)}, ::ScalarOrBandedLayouts, ::ScalarOrBandedLayouts, ::AbstractBandedLayout) = ApplyBandedLayout{typeof(vcat)}()
applylayout(::Type{typeof(vcat)}, ::ScalarOrBandedLayouts, ::ScalarOrBandedLayouts, ::ScalarOrBandedLayouts, ::AbstractBandedLayout) = ApplyBandedLayout{typeof(vcat)}()

applylayout(::Type{typeof(vcat)}, _, ::AbstractBandedLayout) = VcatAlmostBandedLayout()
applylayout(::Type{typeof(vcat)}, _, _, ::AbstractBandedLayout) = VcatAlmostBandedLayout()
applylayout(::Type{typeof(vcat)}, _, _, _, ::AbstractBandedLayout) = VcatAlmostBandedLayout()
sublayout(::VcatAlmostBandedLayout, ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}) =
    VcatAlmostBandedLayout()

__two_arguments(z, y) = (y,z)
__two_arguments(z, y, x...) = (Vcat(reverse(x)..., y), z)
_two_arguments(A) = __two_arguments(reverse(arguments(A))...)

arguments(::VcatAlmostBandedLayout, A) = arguments(ApplyLayout{typeof(vcat)}(), A)

function almostbandwidths(::VcatAlmostBandedLayout, A)
    a,b = _two_arguments(A)
    m̃ = size(a,1)
    l,u = bandwidths(b)
    (max(l+m̃,m̃-1),u-m̃)
end

almostbandedrank(::VcatAlmostBandedLayout, A) = size(first(_two_arguments(A)),1)

function _copyto!(::AlmostBandedLayout, ::AlmostBandedLayout, dest::AbstractMatrix{T}, A::AbstractMatrix) where T
    m,n = size(dest)
    (m,n) == size(A) || throw(DimensionMismatch())
    copyto!(bandpart(dest), bandpart(A))
    copyto!(fillpart(dest), fillpart(A))
    dest
end

function _copyto!(::AlmostBandedLayout, ::VcatAlmostBandedLayout, dest::AbstractMatrix{T}, A::AbstractMatrix) where T
    m,n = size(dest)
    (m,n) == size(A) || throw(DimensionMismatch())
    a,b = _two_arguments(A)
    r = size(a,1)
    bands = bandpart(dest)
    U,V = arguments(fillpart(dest))
    U[1:r,1:r] = Eye{T}(r)
    zero!(view(U,r+1:m,:))
    copyto!(V, a)

    copyto!(@views(bands[r+1:end,:]), b)

    for j = 1:min(r+bandwidth(bands,2),size(bands,2))
        kr = colsupport(bands,j) ∩ (1:r)
        if !isempty(kr)
            bands[kr,j] .= view(a,kr,j)
        end
    end

    dest
end

function resize(A::LowRankMatrix, m::Integer, n::Integer)
    U,V = arguments(A)
    LowRankMatrix([U; Zeros(m-size(U,1),size(U,2))], [V Zeros(size(V,1), n-size(V,2))])
end

function resize(A::AlmostBandedMatrix, m::Integer, n::Integer)
    B = bandpart(A)
    F = fillpart(A)
    AlmostBandedMatrix(resize(B, m,n), resize(F,m,n))
end

function _cache(::VcatAlmostBandedLayout, A::AbstractArray{T}) where T
    r = almostbandedrank(A)
    l,u = almostbandwidths(A)
    data = AlmostBandedMatrix(Zeros{T}(r,r+u), almostbandwidths(A), r)
    CachedArray(data, A, (0,0))
end

function resizedata!(::AlmostBandedLayout, ::VcatAlmostBandedLayout, C::CachedMatrix{T}, m::Integer, n::Integer) where T
    @boundscheck checkbounds(Bool, C, m, n) || throw(ArgumentError("Cannot resize beyound size of matrix"))
    A = C.array
    a,b = _two_arguments(A)
    r = almostbandedrank(C.data)

    # increase size of array if necessary
    μ,ν = C.datasize
    m,n = max(μ,m), max(ν,n)
    l,u = almostbandwidths(A)

    if (μ,ν) ≠ (m,n)
        λ,ω = almostbandwidths(C.data)
        if m ≥ size(C.data,1) || n ≥ size(C.data,2)
            N = 2*max(n,m+u,r)
            C.data = resize(C.data, N+λ, N)
        end
        if μ < r
            kr,jr = μ+1:min(r,m),1:ν
            copyto!(view(C.data,kr,jr), view(A,kr,jr))
        end
        if ν < r+ω+1
            kr,jr = 1:min(r,m),ν+1:min(r+ω+1,n)
            copyto!(view(C.data,kr,jr), view(A,kr,jr))
        end

        B,F = bandpart(C.data),fillpart(C.data)
        U,V = arguments(F)
        jr = max(ν+1,r+ω+2):n
        copyto!(view(V,:,jr), view(a,:,jr))
        zero!(view(U,max(μ+1,r+1):m,:))
        C_band = CachedArray(@view(B[r+1:end,:]), b, (max(0,μ-r),ν))
        resizedata!(C_band, max(0,m-r),n)
        C.datasize = (m,n)
    end

    C
end


text:
module MLBase

    using ArrayViews
    using StatsBase

    import Base: length, show, keys, precision, length, getindex
    import Base: start, next, done
    import Base.Order: lt, Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse
    import StatsBase: RealVector, IntegerVector, RealMatrix, IntegerMatrix, RealArray
    import StatsBase: sample

    export

    # reexport from Base.Order
    Forward, Reverse,

    # reexport some functions from StatsBase
    counts, addcounts!, countmap, proportions,
    indicatormat, sample, weights, WeightVec, wmean,

    # utils
    repeach,        # repeat each element in a vector
    repeachcol,     # repeat each column in a matrix
    repeachrow,     # repeat each row in a matrix
    counteq,        # count the number of equal pairs
    countne,        # count the number of non-equal pairs

    # datapre
    Standardize,    # the type to represent a standardizing transform

    indim,          # input dimension of a transform
    outdim,         # output dimension of a transform
    estimate,       # estimate a model or transformation
    transform,      # apply a transformation to data
    transform!,     # apply a transformation to data in place
    standardize,    # estimate and apply a standardization
    standardize!,   # estimate and apply a standardization in place

    # classification
    LabelMap,       # a type to represent a label map

    classify,       # predict class label(s) based on score values
    classify!,      # inplace version of classify
    classify_withscore,     # classify with additional output of best score
    classify_withscores,    # classify with additional output of best scores
    classify_withscores!,   # inplace version of classify_withscores
    labelmap,       # construct a label map from a list of labels
    labelencode,    # encode a sequence of discrete values using a label map
    labeldecode,    # decode the label to the associated discrete value
    groupindices,   # grouped indices based on labels

    # crossval
    CrossValGenerator,  # abstract base class for all cross-validation plans
    Kfold,              # K-fold cross validation plan
    LOOCV,              # leave-one-out cross validation plan
    RandomSub,          # repetitive random subsampling cross validation

    cross_validate,     # perform cross-validation

    # perfeval
    ROCNums,        # A class to capture ROC numbers

    correctrate,    # compute correct rate of predictions
    errorrate,      # compute error rate of predictions
    counthits,      # count the number of hits
    hitrate,        # compute hit-rate of ranked lists at a specific rank
    hitrates,       # compute hit-rate of ranked lists at multiple ranks
    roc,            # compute roc numbers from predictions
    true_positive,      # number of true positives
    true_negative,      # number of true negatives
    false_positive,     # number of false positives
    false_negative,     # number of false negatives
    true_positive_rate,     # rate of true positives
    true_negative_rate,     # rate of true negatives
    false_positive_rate,    # rate of false positives
    false_negative_rate,    # rate of false negatives
    recall,             # recall computed from ROCNums
    precision,          # precision computed from ROCNums
    f1score,            # F1-score computed from ROCNums

    # deviation
    sqL2dist,       # squared L2 distance between two arrays
    L2dist,         # L2 distance between two arrays
    L1dist,         # L1 distance between two arrays
    Linfdist,       # L-inf distance between two arrays
    gkldiv,         # (Generalized) Kullback-Leibler divergence between two vectors
    meanad,         # mean absolute deviation
    maxad,          # maximum absolute deviation
    msd,            # mean squared deviation
    rmsd,           # root mean squared deviation
    nrmsd,          # normalized rmsd
    psnr            # peak signal-to-noise ratio (in dB)

    # source files

    include("utils.jl")
    include("datapre.jl")
    include("classification.jl")
    include("crossval.jl")
    include("perfeval.jl")
    include("deviation.jl")
end

ext: jl


text:
# reroute through IndexStyle
# This allows converting indices based on IndexStyle, independent of the type tree
# May be this rerouting can be added to julia base someday
@inline to_indices(A::AbstractTiledArray, I::Tuple) = to_indices(IndexStyle(A), A, I)

# route to base definitions by default
@inline to_indices(::IndexStyle, A, I::Tuple) = to_indices(A, axes(A), I)
@inline to_indices(::IndexStyle, A, I::Tuple{Any}) =
    to_indices(A, (eachindex(IndexLinear(), A),), I)

"""
    IndexTiled <: IndexStyle

Signals that an array is best indexed into tile by tile.
"""
struct IndexTiled <: IndexStyle end

IndexStyle(::Type{<:AbstractTiledArray}) = IndexTiled()

"""
    TilingStyle

Defines how an array is tiled.
"""
abstract type TilingStyle end
TilingStyle(T::Type) = throw(MethodError("No TilingStyle defined for $T"))
TilingStyle(x) = TilingStyle(typeof(x))

@inline to_indices(::IndexTiled, A, I) = to_indices(TilingStyle(A), A, I)

"""
    Tile{T}(I::T)

an index or key used to retrieve a tile from an array.
"""
struct Tile{T}
    I::T
end

Tile(I::T) where {T} = Tile{T}(I)

getindex(A, t::Tile, i1, i...) = A[t][i1, i...]
setindex!(A, x, t::Tile, i1, i...) = (view(A, t)[i1, i...] = x)

"""
    GridTiles <: TilingStyle

A `TilingStyle`, where tiles make a grid with no offset.

If the grid size does not divide the size of the array, the last tiles along any given axis
are truncated.
"""
struct GridTiles <: TilingStyle end

@inline function to_indices(gt::GridTiles{N}, A, I::NTuple{N, Integer}) where {N}

end


text:
module NetworkAlign

using MatrixNetworks
using SparseArrays

import Printf: @printf
import DelimitedFiles: readdlm
import LinearAlgebra: triu, dot, norm, issymmetric, tril, fillstored!
import KahanSummation: sum_kbn

"""
Module ``NetworkAlign``: Documentation on the module

You can check the readme file here: \n
"https://github.com/nassarhuda/NetworkAlign.jl/blob/master/README.md"
"""
#NetworkAlign
include("../include/MatrixNetworks_bipartite_matching.jl")
include("normout.jl")
include("othermaxplus.jl")
include("othersum.jl")
include("round_messages.jl")
include("column_maxmatchsum.jl")
include("load_netalign_problem.jl")
include("netalignmr.jl")
include("isorank.jl")
include("netalignbp.jl")
include("count_overlap.jl")
include("make_squares.jl")
include("netalign_setup.jl")

export load_netalign_problem, netalignmr, isorank, netalignbp,
        count_overlap, make_squares, netalign_setup, build_Si

end # end module


text:
using GeometricIntegrators

set_config(:tab_compensated_summation, false)

const out_dir = "../../docs/src/integrators/runge_kutta/guiding_center_4d/"
const run_id  = "fast_barely_passing_firk"

include("guiding_center_4d_fast_barely_passing.jl")
include("../guiding_center_4d_firk.jl")


text:
using Dates

@testset "Crco.jl" begin
    filename = "data\\test.crco"
    df = DDR2import.Crco.read(filename)
    @test df[1, :COUNTRYNAME] == "ED"
    @test df[2, :CALLSIGN] == "N182QS"
    @test df[3, :ACTYPE] == "B737"
    @test df[4, :FLIGHTID] == 1523
    @test df[1, :LATENTRY_DEG] ≈  53.630278 atol = 0.001
    @test df[2, :LONENTRY_DEG] ≈  12.1066383 atol = 0.001
    @test df[3, :LATEXIT_DEG] ≈  50.3166667 atol = 0.001
    @test df[4, :LONEXIT_DEG] ≈  9.98833333 atol = 0.001
    @test df[1, :DISTANCE_M] ≈  396031.3143 atol = 0.1
    @test df[4, :ENTRYFL] ≈  270.0 atol = 0.1
    @test df[1, :EXITFL] ≈  270.0 atol = 0.1
    @test Dates.minute(df[2, :DATETIMEENTRY]) == 59
    @test Dates.second(df[3, :DATETIMEEXIT]) == 22
end


text:
using Plots
using Statistics
using StatsBase: sample
using Distributions
using Revise
using EvoTrees

features = rand(Int(1.25e4), 1)
# prepare a dataset
# features = rand(100, 10)
X = features
Y = randn(size(X, 1)) .* 0.1
Y[X[:,1] .< 0.2] .*= 2
Y[(X[:,1] .>= 0.4) .& (X[:,1] .< 0.6)] .*= 5
Y[(X[:,1] .>= 0.9)] .*= 5
𝑖 = collect(1:size(X,1))

Y .*= 0.01

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# train model
params1 = EvoTreeGaussian(
    loss=:gaussian, metric=:gaussian,
    nrounds=200,
    λ = 1.0, γ=1.0, η=0.5,
    max_depth = 4, min_weight = 100.0,
    rowsample=1.0, colsample=1.0, nbins=64)

@time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n = 10);
# @time model = fit_evotree(params1, X_train, Y_train, print_every_n = 10);
@time pred_train = EvoTrees.predict(model, X_train)
@time pred_train_gauss = EvoTrees.predict(params1, model, X_train)

pred_gauss = [Distributions.Normal(pred_train[i,1], pred_train[i,2]) for i in 1:size(pred_train,1)]
pred_q90 = quantile.(pred_gauss, 0.9)
pred_q10 = quantile.(pred_gauss, 0.1)

mean(Y_train .< pred_q90)
mean(Y_train .< pred_q10)

x_perm = sortperm(X_train[:,1])
plot(X_train[:, 1], Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train[x_perm, 1], color = "navy", linewidth = 1.5, label = "mu")
plot!(X_train[:,1][x_perm], pred_train[x_perm, 2], color = "blue", linewidth = 1.5, label = "sigma")
plot!(X_train[:,1][x_perm], pred_q10[x_perm, 1], color = "darkred", linewidth = 1.5, label = "q10")
plot!(X_train[:,1][x_perm], pred_q90[x_perm, 1], color = "green", linewidth = 1.5, label = "q90")
savefig("figures/regression_gaussian_v1.png")

# compare with zygote
using Zygote

pred = [0.0, log(1.0)]
target = 0.1

δ1 = (target - pred[1]) / max(1e-8, exp(2*pred[2]))
δ2 = (1 - (pred[1] - target)^2 / max(1e-8, exp(2*pred[2])))

δ²1 = 1 / max(1e-8, exp(2*pred[2]))
δ²2 = 2 / max(1e-8, exp(2*pred[2])) * (pred[1] - target)^2

lpdf(x,μ,σ) = -log(σ) - log(2π)/2 - 1/2*((x-μ)/σ)^2
lpdf(0, pred[1], pred[2])

lpdf2(x,μ,lσ) = -log(exp(lσ)) - log(2π)/2 - 1/2*((x-μ)/exp(lσ))^2
lpdf2(0, pred[1], pred[2])

n1 = Normal(0, 1)
Distributions.logpdf(n1, 0)

# gradient(lpdf, target, pred[1], pred[2])[2:end]
gradient(lpdf2, target, pred[1], pred[2])[2:end]
Zygote.hessian(lpdf2, target, pred[1], pred[2])

gradient_lpdf(x,pred) = gradient(lpdf2, x, pred[1], pred[2])[3]
hessian_lpdf(x,pred) = gradient(gradient_lpdf, x, pred)[1]
gradient_lpdf(target, pred)
hessian_lpdf(target, pred)


text:
using LowRankModels
import StatsBase: sample

println("simple glrm examples")

# minimize ||A - XY||^2
function fit_pca(m,n,k)
	# matrix to encode
	A = randn(m,k)*randn(k,n)
	loss = QuadLoss()
	r = ZeroReg()
	glrm = GLRM(A,loss,r,r,k)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

# minimize_{X>=0, Y>=0} ||A - XY||^2
function fit_nnmf(m,n,k)
	# matrix to encode
	A = rand(m,k)*rand(k,n)
	loss = QuadLoss()
	r = NonNegConstraint()
	glrm = GLRM(A,loss,r,r,k)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

# minimize ||A - XY||^2 + .1||X||^2 + .1||Y||^2
function fit_pca_nucnorm(m,n,k)
	# matrix to encode
	A = randn(m,k)*randn(k,n)
	loss = QuadLoss()
	r = QuadReg(.1)
	glrm = GLRM(A,loss,r,r,k)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

# minimize_{X<=0} ||A - XY||^2
function fit_kmeans(m,n,k)
	# matrix to encode
	Y = randn(k,n)
	A = zeros(m,n)
	for i=1:m
		A[i,:] = Y[mod(i,k)+1,:]
	end
	loss = QuadLoss()
	ry = ZeroReg()
	rx = UnitOneSparseConstraint()
	glrm = GLRM(A,loss,rx,ry,k+4)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

function fit_pca_nucnorm_sparse(m,n,k,s)
	# matrix to encode
	A = randn(m,k)*randn(k,n)
	loss = QuadLoss()
	r = QuadReg(.1)
	obsx = sample(1:m,s); obsy = sample(1:n,s)
	obs = [(obsx[i],obsy[i]) for i=1:s]
	glrm = GLRM(A,loss,r,r,k, obs=obs)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

function fit_pca_nucnorm_sparse_nonuniform(m,n,k,s)
	# matrix to encode
	A = randn(m,k)*randn(k,n)
	loss = QuadLoss()
	r = QuadReg(.1)
	obsx = [sample(1:round(Int,m/4),round(Int,s/2)); sample(round(Int,m/4)+1:m,s-round(Int,s/2))]
	obsy = sample(1:n,s)
	obs = [(obsx[i],obsy[i]) for i=1:s]
	glrm = GLRM(A,loss,r,r,k, obs=obs)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

function fit_soft_kmeans(m,n,k)
	# PCA with loadings constrained to lie on unit simplex
	# constrain columns of X to lie on unit simplex
	Xreal = rand(k,m)
	Xreal ./= sum(Xreal,1)
	A = Xreal' * randn(k,n)

	loss = QuadLoss()
	rx = SimplexConstraint()
	ry = ZeroReg()
	glrm = GLRM(A,loss,rx,ry,k)
	X,Y,ch = fit!(glrm)
	println("Convergence history:",ch.objective)
	return A,X,Y,ch
end

if true
	Random.seed!(10)
	fit_pca(100,100,2)
	fit_pca_nucnorm(100,100,2)
	fit_pca_nucnorm_sparse(500,500,2,10000)
	fit_pca_nucnorm_sparse_nonuniform(1000,1000,2,20000)
	fit_kmeans(50,50,10)
	fit_nnmf(50,50,2)
end


text:
"""
    ManyElectronWavefunction

A many-electron wave function configuration can either be given as a
CSF (summed over spins) or a configuration of spin-orbitals (where all
quantum numbers are specified).
"""
const ManyElectronWavefunction{O<:AbstractOrbital} = Union{CSF{O},Configuration{<:SpinOrbital}}

"""
    Atom(radial_orbitals, orbitals, configurations, mix_coeffs, potential)

An atom constitutes a set of single-electron `orbitals` with
associated `radial_orbitals`, `configurations` which are
[`ManyElectronWavefunction`](@ref):s, comprising of anti-symmetrized
combinations of such orbitals. The expansion coefficients `mix_coeffs`
determine the linear combination of the `configurations` for
multi-configurational atoms.

The `potential` can be used to model either the nucleus by itself (a
point charge or a nucleus of finite extent) or the core orbitals
(i.e. a pseudo-potential).
"""
mutable struct Atom{T,B<:Basis,O<:AbstractOrbital,TC<:ManyElectronWavefunction,CV<:AbstractVector,P<:AbstractPotential} <: AbstractQuantumSystem
    radial_orbitals::RadialOrbitals{T,B}
    orbitals::Vector{O}
    configurations::Vector{TC}
    mix_coeffs::CV # Mixing coefficients for multi-configurational atoms
    potential::P
end

get_config(config::Configuration) = config
get_config(csf::CSF) = csf.config

function Base.copy(atom::A) where {A<:Atom}
    R,Φ = atom.radial_orbitals.args
    A(applied(*, R, copy(Φ)), copy(atom.orbitals),
      copy(atom.configurations), copy(atom.mix_coeffs),
      atom.potential)
end

"""
    outsidecoremodel(configuration::Configuration, potential::P)

Return the part of the electronic `configuration` that is not part of
the the configuration modelled by the `potential`. For a point charge,
this is the same as the `configuration` itself, but for
pseudopotential, typically only the outer shells remain.
"""
outsidecoremodel(configuration::Configuration, potential::P) where P =
    filter((orb,occ,state) -> nonrelorbital(getspatialorb(orb)) ∉ core(ground_state(potential)), configuration)
outsidecoremodel(configuration::Configuration, ::PointCharge) where P =
    configuration
outsidecoremodel(csf::CSF, potential) = outsidecoremodel(get_config(csf), potential)

"""
    Atom(undef, ::Type{T}, R::AbstractQuasiMatrix, configurations, potential, ::Type{C}[, mix_coeffs])

Create an [`Atom`](@ref) on the space spanned by `R`, from the list of
electronic `configurations`, with a nucleus modelled by `potential`,
and leave the orbitals uninitialized. `T` determines the `eltype` of
the radial orbitals and `C` the mixing coefficients, which by default,
are initialized to `[1,0,0,...]`.
"""
function Atom(::UndefInitializer, ::Type{T}, R::B, configurations::Vector{TC}, potential::P,
              ::Type{C},
              mix_coeffs::CV=vcat(one(C), zeros(C, length(configurations)-1))) where {T<:Number,B<:BasisOrRestricted,TC,C,CV<:AbstractVector{<:C},P}
    isempty(configurations) &&
        throw(ArgumentError("At least one configuration required to create an atom"))
    all(isequal(num_electrons(first(configurations))),
        map(config -> num_electrons(config), configurations)) ||
            throw(ArgumentError("All configurations need to have the same amount of electrons"))
    all(isequal(core(get_config(first(configurations)))),
        map(config -> core(get_config(config)), configurations)) ||
            throw(ArgumentError("All configurations need to share the same core orbitals"))

    pot_cfg = core(ground_state(potential))
    for cfg in configurations
        cfg_closed_orbitals = sort(unique(nonrelorbital.(getspatialorb.(core(get_config(cfg)).orbitals))))
        pot_cfg.orbitals ⊆ cfg_closed_orbitals ||
            throw(ArgumentError("Configuration modelled by nuclear potential ($(pot_cfg)) must belong to the closed set of all configurations"))
    end

    orbs = unique_orbitals(outsidecoremodel.(get_config.(configurations), Ref(potential)))

    Φ = Matrix{T}(undef, size(R,2), length(orbs))
    RΦ = applied(*, R, Φ)
    Atom(RΦ, orbs, configurations, mix_coeffs, potential)
end

"""
    Atom(init, ::Type{T}, R::AbstractQuasiMatrix, configurations, potential, ::Type{C})

Create an [`Atom`](@ref) on the space spanned by `R`, from the list of
electronic `configurations`, with a nucleus modelled by `potential`,
and initialize the orbitals according to `init`. `T` determines the
`eltype` of the radial orbitals and `C` the mixing coefficients.
"""
function Atom(init::Symbol, ::Type{T}, R::B, args...; kwargs...) where {T<:Number,B<:BasisOrRestricted,TC,C,P}
    atom = Atom(undef, T, R, args...)
    if init == :hydrogenic
        hydrogenic!(atom; kwargs...)
    elseif init == :screened_hydrogenic
        screened_hydrogenic!(atom; kwargs...)
    elseif init == :zero
        atom.radial_orbitals.args[2] .= zero(T)
    else
        throw(ArgumentError("Unknown orbital initialization mode $(init)"))
    end
    atom
end

"""
    DiracAtom

A `DiracAtom` is a specialization of [`Atom`](@ref) for the
relativistic case.
"""
const DiracAtom{T,B,TC<:RelativisticCSF,C,P} = Atom{T,B,TC,C,P}

"""
    Atom(init, R::AbstractQuasiMatrix, configurations, potential, ::Type{C})

Create an [`Atom`](@ref) on the space spanned by `R`, from the list of
electronic `configurations`, with a nucleus modelled by `potential`,
and initialize the orbitals according to `init`. `C` determines the
`eltype` of the mixing coefficients.
"""
Atom(init::Init, R::B, args...; kwargs...) where {Init,T,B<:AbstractQuasiMatrix{T}} =
    Atom(init, T, R, args...; kwargs...)

DiracAtom(init::Init, R::B, configurations::Vector{TC}, args...; kwargs...) where {Init,T,B<:AbstractQuasiMatrix{T},TC<:RelativisticCSF} =
    Atom(init, T, R, configurations, args...; kwargs...)

"""
    Atom(R::AbstractQuasiMatrix, configurations, potential[, ::Type{C}=eltype(R)])

Create an [`Atom`](@ref) on the space spanned by `R`, from the list of
electronic `configurations`, with a nucleus modelled by `potential`,
and initialize the orbitals to their hydrogenic values.
"""
Atom(R::B, configurations::Vector{<:TC}, potential::P, ::Type{C}=eltype(R), args...;
     kwargs...) where {B<:AbstractQuasiMatrix,TC<:ManyElectronWavefunction,C,P<:AbstractPotential} =
         Atom(:hydrogenic, R, configurations, potential, C, args...; kwargs...)

"""
    Atom(other_atom::Atom, configurations)

Create a new atom using the same basis and nuclear potential as
`other_atom`, but with a different set of `configurations`. The
orbitals of `other_atom` are copied over as starting guess.
"""
function Atom(other_atom::Atom{T,B,O,TC,CV,P},
              configurations::Vector{<:TC}; kwargs...) where {T,B,O,TC,CV,P}
    core(first(configurations)) == core(first(other_atom.configurations)) ||
        throw(ArgumentError("Core orbitals must be the same"))
    R = radial_basis(other_atom)
    # It is slightly wasteful to first initialize all orbitals
    # hydrogenically, and the overwrite with orbitals from
    # other_atom. Ideally, missing orbitals would be estimated by
    # diagonalizing e.g. their energy expressions without the exchange
    # terms.
    atom = Atom(R, configurations, other_atom.potential, eltype(CV); kwargs...)
    copyto!(atom, other_atom)
    atom
end

"""
    copyto!(dst::Atom, src::Atom)

Copy the radial orbitals of `src` to the corresponding ones in
`dst`. Orbitals of `src` missing from `dst` are skipped. It is assumed
that the underlying bases are compatible, i.e. that the radial
coordinate of `dst` encompasses that of `src` and grid spacing
&c. agree.
"""
function Base.copyto!(dst::Atom, src::Atom)
    Rsrc = radial_basis(src)
    Rdst = radial_basis(dst)
    axes(Rsrc,1).domain ⊆ axes(Rdst,1).domain &&
        axes(Rsrc,2) ⊆ axes(Rdst,2) ||
        throw(ArgumentError("$(Rsrc) not a subset of $(Rdst)"))

    for o in src.orbitals
        o ∈ dst.orbitals || continue
        copyto!(view(dst, o).args[2], view(src, o).args[2])
    end
end

"""
    num_electrons(atom)

Return number of electrons in `atom`.
"""
AtomicLevels.num_electrons(atom::Atom) =
    num_electrons(first(atom.configurations))

function Base.show(io::IO, atom::Atom{T,B,O,TC,C,P}) where {T,B,O,TC,C,P}
    Z = charge(atom.potential)
    N = num_electrons(atom)
    write(io, "Atom{$(T)}(R=$(radial_basis(atom)); $(atom.potential); $(N) e⁻ ⇒ Q = $(Z-N)) with ")
    write(io, "$(length(atom.configurations)) $(TC)")
    length(atom.configurations) > 1 && write(io, "s")
end

function Base.show(io::IO, ::MIME"text/plain", atom::Atom{T,B,O,TC,C,P}) where {T,B,O,TC,C,P}
    show(io, atom)
    if length(atom.configurations) > 1
        write(io, ":\n")
        show(io, "text/plain", atom.configurations)
    else
        write(io, ": ")
        show(io, atom.configurations[1])
    end
end

radial_basis(atom::Atom) = first(atom.radial_orbitals.args)

function orbital_index(atom::Atom{T,B,O}, orb::O) where {T,B,O}
    i = findfirst(isequal(orb),atom.orbitals)
    i === nothing && throw(BoundsError("$(orb) not present among $(atom.orbitals)"))
    i
end

"""
    getindex(atom, j)

Returns a copy of the `j`:th radial orbital.
"""
Base.getindex(atom::Atom, j::I) where {I<:Integer} =
    applied(*, radial_basis(atom), atom.radial_orbitals.args[2][:,j])

"""
    getindex(atom, orb)

Returns a copy of the radial orbital corresponding to `orb`.
"""
Base.getindex(atom::Atom{T,B,O}, orb::O) where {T,B,O} =
    atom[orbital_index(atom, orb)]

"""
    view(atom, j)

Returns a `view` of the `j`:th radial orbital.
"""
Base.view(atom::Atom, j::I) where {I<:Integer} =
    applied(*, radial_basis(atom), view(atom.radial_orbitals.args[2], :, j))

"""
    view(atom, orb)

Returns a `view` of the radial orbital corresponding to `orb`.
"""
Base.view(atom::Atom{T,B,O}, orb::O) where {T,B,O} =
    view(atom, orbital_index(atom, orb))

"""
    SCF.coefficients(atom)

Returns a `view` of the mixing coefficients.
"""
SCF.coefficients(atom::A) where {A<:Atom} =
    view(atom.mix_coeffs, :)

"""
    SCF.orbitals(atom)

Returns a `view` of the radial orbital coefficients (NB, it does _not_
return the `MulQuasiMatrix`, but the actual underlying expansion
coefficients, since `SCF` operates on them in the self-consistent
iteration).
"""
SCF.orbitals(atom::A) where {A<:Atom} =
    view(atom.radial_orbitals.args[2], :, :)

"""
    norm(atom[, p=2; configuration=1])

This calculates the _amplitude_ norm of the `atom`, i.e. ᵖ√N where N
is the number electrons. By default, it uses the first `configuration`
of the `atom` to weight the individual orbital norms.
"""
function LinearAlgebra.norm(atom::Atom{T}, p::Real=2; configuration::Int=1) where T
    RT = real(T)
    n = zero(RT)
    for (orb,occ,state) in outsidecoremodel(get_config(atom.configurations[configuration]),
                                            atom.potential)
        n += occ*norm(view(atom, orb), p)^p
    end
    n^(one(RT)/p) # Unsure why you'd ever want anything but the 2-norm, but hey
end

LinearAlgebra.normalize!(atom::A, v::V) where {A<:Atom,V<:AbstractVector} =
    normalize!(applied(*, radial_basis(atom), v))

function SCF.norm_rot!(ro::RO) where {RO<:RadialOrbital}
    normalize!(ro)
    SCF.rotate_first_lobe!(ro.args[2])
    ro
end

export Atom, DiracAtom, radial_basis


text:
using GLMakie
using BenchmarkTools
GLMakie.activate!()

##

NUM_BOIDS = 100
VIEW_DISTANCE = 70
VIEW_ANGLE = 4π/3
SIDE_MARGIN = 50
MARGIN_TURN = 1

# Boid Modifiers

ALIGNMENT_MODIFIER = 0.1
COHESION_MODIFIER = 0.1
SEPERATION_MODIFIER = 0.02
MIN_DISTANCE = VIEW_DISTANCE / 10

MAX_VELOCITY = 120
MIN_VELOCITY = 10
MAX_ACCELERATION = 10

# Animation settings
FRAMERATE = 30
SECONDS = 30

# Plot Settings
X_WIDTH = 1000
Y_WIDTH = 1000

mutable struct Boids
    x::Vector{Float64} # m
    y::Vector{Float64} # m
    dx::Vector{Float64} # m/s
    dy::Vector{Float64} # m/s
end

Boids(x) = Boids(rand(x) * X_WIDTH, rand(x) * Y_WIDTH, rand(x)*(MAX_VELOCITY) .- MIN_VELOCITY, rand(x)*(MAX_VELOCITY) .- MIN_VELOCITY)

##

function iterateBoids!(boidsArr)
    # Iterate positions

    for i in 1:NUM_BOIDS
        dist = .√(((boidsArr.x .- boidsArr.x[i]) .^ 2) .+ ((boidsArr.y .- boidsArr.y[i]) .^ 2))

        xs = []
        ys = []
        dxttl = 0
        dyttl = 0
        moveX = 0
        moveY = 0
        numInRange = 0

        for inner in 1:NUM_BOIDS
            if dist[inner] > 0.1 && dist[inner] < VIEW_DISTANCE
                #Cohesion
                push!(ys, boidsArr.y[inner])
                push!(xs, boidsArr.x[inner])

                dxttl += boidsArr.dx[inner]
                dyttl += boidsArr.dy[inner]
                #Seperation
                moveX += boidsArr.x[i] - boidsArr.x[inner]
                moveY += boidsArr.y[i] - boidsArr.y[inner]

                numInRange += 1
            end

            if numInRange > 0
                # COHESION
                avgX = sum(xs) / numInRange
                avgY = sum(ys) / numInRange

                boidsArr.dx[i] += (avgX - boidsArr.x[i]) * COHESION_MODIFIER
                boidsArr.dy[i] += (avgY - boidsArr.y[i]) * COHESION_MODIFIER

                #SEPERATION
                boidsArr.dx[i] += moveX * SEPERATION_MODIFIER
                boidsArr.dy[i] += moveY * SEPERATION_MODIFIER

                #ALIGNMENT
                avgDX = dxttl / numInRange
                avgDY = dyttl / numInRange

                boidsArr.dx[i] += (avgDX - boidsArr.dx[i]) * ALIGNMENT_MODIFIER
                boidsArr.dy[i] += (avgDY - boidsArr.dy[i]) * ALIGNMENT_MODIFIER
            end

            #Limit Velocity
            mag = getMagnitude(boidsArr.dx[i], boidsArr.dy[i])
            if mag > MAX_VELOCITY
                boidsArr.dx[i] = (boidsArr.dx[i] / mag) * MAX_VELOCITY
                boidsArr.dy[i] = (boidsArr.dy[i] / mag) * MAX_VELOCITY
            end

            #Stay in bounds
            if boidsArr.x[i] < SIDE_MARGIN
                boidsArr.dx[i] += MARGIN_TURN
            elseif boidsArr.x[i] > X_WIDTH - SIDE_MARGIN
                boidsArr.dx[i] -= MARGIN_TURN
            end
            if boidsArr.y[i] < SIDE_MARGIN
                boidsArr.dy[i] += MARGIN_TURN
            elseif boidsArr.y[i] > Y_WIDTH - SIDE_MARGIN
                boidsArr.dy[i] -= MARGIN_TURN
            end
        end
    end

    boidsArr.x = boidsArr.x .+ (boidsArr.dx ./ FRAMERATE)
    boidsArr.y = boidsArr.y .+ (boidsArr.dy ./ FRAMERATE)
end

function getMagnitude(dx, dy)
    return √(dx^2 + dy^2)
end

##

# Define Figure

fig = Figure(resolution=(X_WIDTH, Y_WIDTH))

ax = fig[1, 1] = Axis(fig)

hidespines!(ax)
hidedecorations!(ax)

limits!(ax, 0, X_WIDTH, 0, Y_WIDTH)

##
boids = Boids(NUM_BOIDS)

timestamps = range(0, SECONDS, step=1 / FRAMERATE)

record(fig, "animation.mp4", timestamps; framerate = FRAMERATE) do i
    empty!(ax)
    iterateBoids!(boids)
    # tans = atan.(boids.yv ./ boids.xv)
    scatter!(ax, boids.x, boids.y, color = :black)
end

##

@benchmark iterateBoids!(b) setup=(b = Boids(NUM_BOIDS)) seconds=10

max_stars_repo_name: MarshallVielmetti/JuliaBoids

max_issues_repo_name: MarshallVielmetti/JuliaBoids

max_forks_repo_name: MarshallVielmetti/JuliaBoids


text:
function predictionDPR(x,degree,maxiter,burnin,seed=1,T=0,filename = "/results")
  Random.seed!(seed)
  savelocation = string(pwd(), filename, "/dpr_seed$seed")
  mkpath(savelocation)

  n = length(x)
  if T > 0
    append!(x, zeros(T))
    xpred = fill(Float64[],T)
    for t in 1:1:T
      xpred[t] = zeros(maxiter-burnin)
    end
    n_ = n + T
  else
    n_ = n
  end

  degree += 1
  theta = zeros(degree) # polynomial coeeficients vector
  sampledTheta = zeros(maxiter-burnin,degree) # matrix to store sample thetas
  thetaprev = copy(theta)
  epsj = 20.
  x0 = 0.0
  sampledX0 = zeros(maxiter-burnin)
  zp = zeros(maxiter-burnin)
  clusters = zeros(Int64,maxiter)
  C = zeros(maxiter)
  # MCMC OPTIONS

  ### initialize concentration parameter of the DP and specify prior hyperparameters
  c = 1.
  ac,bc= 1.,1.

  ### shape and rate for the gamma prior on precisions
  at, bt = 0.001,0.001

  ### DP meausre latent variables
  d = collect(1:1:n_) # clustering variables
  u = zeros(n_)       # slice variables

  ### begin MCMC
  for its in 1:1:maxiter

    M = maximum(d)
    v = zeros(M)
    w = zeros(M)
    tau = zeros(M)
    for j in 1:1:M
      flag1 = 0.0
      flag2 = 0.0
      for i in 1:1:n_
        if d[i].==j
          flag1 += 1.
        end
        if d[i].>j
          flag2 += 1.
        end
      end
      v[j] = rand(Beta(1 + flag1, c + flag2))
    end
    w[1] = v[1]
    for j in 2:1:M
      w[j] = w[j-1]*(1-v[j-1])*v[j]/(v[j-1])
    end

    # sample slice variables
    for i in 1:1:n_
      u[i] = rand(Uniform(0.0, w[d[i]]))
    end

    # sample precisions
    for j in 1:1:M
      term,counts = 0.,0.
      if d[1].==j
        counts += 1
        term += (x[1] - polyMap(theta,x0))^2
      end
      for i in 2:1:n_
        if d[i].==j
          counts += 1
          term += (x[i] - polyMap(theta,x[i-1]))^2
        end
      end
      tau[j] = rand(Gamma(at + 0.5counts,(bt + 0.5term)^(-1)))
    end

    # sample aditional number of weights
    sm, Mflag, iflag = 0, 0, 0
    sm = sum(w)
    ustar = minimum(u)
    while (sm < 1. - ustar)
      iflag = iflag + 1
      Mflag = M + iflag
      resize!(v, Mflag)
      resize!(w, Mflag)
      v[Mflag] = rand(Beta(1.0, c))
      w[Mflag] = w[Mflag - 1] * v[Mflag] * (1 - v[Mflag - 1]) / v[Mflag - 1]
      sm = sm + w[Mflag]
    end
    if Mflag > M
      append!(tau, rand(Gamma(at,1/bt), Mflag-M))
    end

    # sample clustering variables
    nc0 = 0.
    A1 = findall(u[1] .< w)
    for j in 1:1:length(A1)
      nc0 += tau[A1[j]]^(0.5) * exp(-0.5*tau[A1[j]]*(x[1]-polyMap(theta,x0))^2)
    end
    rd0 = rand()
    prob0 = 0.
    for j in 1:1:length(A1)
      prob0 += tau[A1[j]]^(0.5) * exp(-0.5*tau[A1[j]]*(x[1]-polyMap(theta,x0))^2) / nc0
      if rd0.<prob0
        d[1] = A1[j]
        break
      end
    end

    for i in 2:1:n_
      Ai = findall(u[i] .< w)
      nc = 0.
      for j in 1:1:length(Ai)
        nc += tau[Ai[j]]^(0.5) * exp(-0.5*tau[Ai[j]]*(x[i]-polyMap(theta,x[i-1]))^2)
      end
      prob = 0.
      rd = rand()
      for j in 1:1:length(Ai)
        prob += tau[Ai[j]]^(0.5) * exp(-0.5*tau[Ai[j]]*(x[i]-polyMap(theta,x[i-1]))^2) / nc
        if rd.<prob
          d[i] = Ai[j]
          break
        end
      end
    end

    clusters[its] = length(unique(d))

    # update concentration parameter
    c = learnc(ac,bc,c,n_,length(unique(d)))
    C[its] = c

    if its > burnin
      # sample from noise predictive
      W = cumsum(w)
      rp = rand()
      if rp.>W[end]
        zp[its-burnin] = rand(Normal(0., sqrt(1/rand(Gamma(at,1/bt)))))
      else
        for i in 1:1:length(W)
          if rp.<W[i]
            zp[its-burnin] = rand(Normal(0., sqrt(1/tau[i])))
            break
          end
        end
      end

      # sample the vector with the coefficients
      for j in 1:1:length(theta)
        thetaj = copy(theta)
        thetaj[j] = 0.

        tauj = tau[d[1]]*x0^(2(j-1))
        meanj = tau[d[1]]*x0^(j-1) * (x[1]-polyMap(thetaj,x0))
        for i in 2:1:n_
          tauj += tau[d[i]]*x[i-1]^(2(j-1))
          meanj += tau[d[i]]*x[i-1]^(j-1) * (x[i] - polyMap(thetaj,x[i-1]))
        end
        meanj /= tauj
        vj = -(2/tauj)*log(rand()) + (thetaprev[j]-meanj)^2
        theta[j] = rand(Uniform(max(-sqrt(vj)+meanj,-epsj), min(sqrt(vj)+meanj,epsj)))
      end
      thetaprev = copy(theta)
      sampledTheta[its-burnin,:] = thetaprev
      # sample x0
      aux = -(2.0/tau[d[1]])*log(rand()) + (x[1]-polyMap(theta,x0))^2
      poly_right = copy(theta)
      poly_right[1] = poly_right[1] - (x[1] - sqrt(aux))
      poly_right = poly_right./poly_right[end]
      roots_right = Polynomials.roots(Poly(poly_right))
      poly_left = copy(theta)
      poly_left[1] = poly_left[1] - (x[1] + sqrt(aux))
      poly_left = poly_left./poly_left[end]
      roots_left = Polynomials.roots(Poly(poly_left))
      allroots = [roots_left roots_right]
      idx = abs.(imag(allroots)).<1e-06
      realroots = sort(real(allroots[idx]))
      intervals = range_intersection(realroots,[-20. 20.])
      x0 = unifmixrnd(intervals)
      sampledX0[its-burnin] = x0

      # predict future unobserved observations directly
      # good results only for the noise process f₂
      if T > 0
        for j in 1:1:T-1
          ax1 = -(2/tau[d[n+j]])*log(rand()) + (x[n+j]-polyMap(theta,x[n+j-1]))^2
          ax2 = -(2/tau[d[n+j+1]])*log(rand()) + (x[n+j+1]-polyMap(theta,x[n+j]))^2

          xnm = polyMap(theta,x[n+j-1]) - sqrt(ax1)
          xnp = polyMap(theta,x[n+j-1]) + sqrt(ax1)

          poly_right_x = copy(theta)
          poly_right_x[1] = poly_right_x[1] - (x[n+j+1]-sqrt(ax2))
          poly_right_x = poly_right_x./poly_right_x[end]
          roots_right_x = Polynomials.roots(Poly(poly_right_x))

          poly_left_x = copy(theta)
          poly_left_x[1] = poly_left_x[1] - (x[n+j+1]+sqrt(ax2))
          poly_left_x = poly_left_x./poly_left_x[end]
          roots_left_x = Polynomials.roots(Poly(poly_left_x))

          allrootsx = [roots_left_x roots_right_x]
          idxx = abs.(imag(allrootsx)).<1e-06
          realrootsx = sort(real(allrootsx[idxx]))
          intervalsx = range_intersection(realrootsx,[xnm xnp])
          x[n+j] = unifmixrnd(intervalsx)
          xpred[j][its-burnin] = x[n+j]
        end
        new_mu = polyMap(theta,x[n+T-1])
        new_var = 1/tau[d[n+T]]
        x[n+T] = rand(Normal(new_mu, sqrt(new_var)))
        xpred[T][its-burnin] = x[n+T]
      end

    end

    if mod(its,10000)==0
      println("MCMC iterations: $its")
    end
  # end gibbs
  end

  ### write outputs to data files
  writedlm(string(savelocation, "/thetas.txt"), sampledTheta)
  writedlm(string(savelocation, "/X0.txt"), sampledX0)
  writedlm(string(savelocation, "/clusters.txt"), clusters)
  writedlm(string(savelocation, "/C.txt"), C)
  writedlm(string(savelocation, "/noise.txt"), zp)

  if T > 0
    for i in 1:1:size(xpred)[1]
      writedlm(string(savelocation, "/xpred$i.txt"), xpred[i])
    end
    return sampledTheta, zp, sampledX0, xpred
  else
    return sampledTheta, zp, sampledX0
  end
end

max_stars_repo_name: cmerkatas/GSBR

max_issues_repo_name: cmerkatas/GSBR

max_forks_repo_name: cmerkatas/GSBR


text:
module Ecosystem

using StatsBase
using EcosystemCore

include("./ecosystem_macros.jl")
include("./ecosystem_agents.jl")

export World
export agent_step!, agent_count, world_step!, simulate!, every_nth

function simulate!(world::World, iters::Int; cb=()->())
    for i in 1:iters
        world_step!(world)
        cb()
    end
end

agent_count(p::Plant) = size(p)/EcosystemCore.max_size(p)
agent_count(::Animal) = 1
agent_count(as::Vector{<:Agent}) = sum(agent_count,as)

function agent_count(w::World)
    function op(d::Dict,a::Agent{S}) where S<:Species
        n = nameof(S)
        if n in keys(d)
            d[n] += agent_count(a)
        else
            d[n] = agent_count(a)
        end
        return d
    end
    foldl(op, w.agents |> values |> collect, init=Dict{Symbol,Real}())
end

function every_nth(f::Function, n::Int)
    i = 1
    function callback(args...)
        # display(i) # comment this out to see out the counter increases
        if i == n
            f(args...)
            i = 1
        else
            i += 1
        end
    end
end

end


text:
@static if VERSION < v"1.1"
    const PackageEntry = Vector{Dict{String,Any}}
else
    using Pkg.Types: PackageEntry
end

"""
    manifest(c::Pkg.Types.Context)
Retrieves the UUID -> PackageEntry map from the manifest of a Context.
"""
function manifest(c::Pkg.Types.Context)
    m = c.env.manifest
    if m isa Dict
        return m
    else
        return m.deps
    end
end

"""
    read_manifest(manifest_filename)

Read the manifest from the path and return the UUID -> PackageEntry map.
If the file can't be read, return `nothing`.
"""
function read_manifest(manifest_filename)
    try
        m = Pkg.API.read_manifest(manifest_filename)
        if m isa Dict
            return m
        else
            return m.deps
        end
    catch err
        @warn "Could not load manifest." exception=(err, catch_backtrace())
        return nothing
    end
end

"""
    project(c::Pkg.Types.Context)
Retrieves the project of a Context.
"""
project(c::Pkg.Types.Context) = c.env.project

"""
    isinproject(context, package::Union{String,UUID})
Checks whether a package is in the dependencies of a given context, e.g. is directly loadable.
"""
function isinproject end

"""
    isinmanifest(context, package::Union{String,UUID})
Checks whether a package is in the manifest of a given context, e.g. is either directly loadable or is a dependency of an loadable package.
"""
function isinmanifest end

@static if VERSION < v"1.1"
    isinmanifest(context::Pkg.Types.Context, module_name::String) = module_name in keys(manifest(context))
    isinmanifest(context::Pkg.Types.Context, uuid::UUID) = any(get(p[1], "uuid", "") == string(uuid) for (u, p) in manifest(context))
    isinmanifest(manifest::Dict{String,Any}, uuid::AbstractString) = any(get(p[1], "uuid", "") == uuid for (u, p) in manifest)
    isinmanifest(manifest::Dict{String,Any}, uuid::UUID) = isinmanifest(manifest, string(uuid))

    isinproject(context::Pkg.Types.Context, package_name::String) = haskey(deps(project(context)), package_name)
    isinproject(context::Pkg.Types.Context, package_uuid::UUID) = any(u == package_uuid for (n, u) in deps(project(context)))

    function packageuuid(c::Pkg.Types.Context, name::String)
        for pkg in manifest(c)
            if first(pkg) == name
                return UUID(last(pkg)[1]["uuid"])
            end
        end
    end
    packageuuid(pkg::Pair{Any,Any}) = last(pkg) isa String ? UUID(last(pkg)) : UUID(first(last(pkg))["uuid"])
    packageuuid(pkg::Pair{String,Any}) = last(pkg) isa String ? UUID(last(pkg)) : UUID(first(last(pkg))["uuid"])

    packagename(pkg::Pair{String,Any})::String = first(pkg)
    function packagename(c::Pkg.Types.Context, uuid)
        for (n, p) in manifest(c)
            if get(first(p), "uuid", "") == string(uuid)
                return n
            end
        end
        return nothing
    end
    function packagename(manifest::Dict{String,Any}, uuid::String)
        for (n, p) in manifest
            if get(first(p), "uuid", "") == string(uuid)
                return n
            end
        end
        return nothing
    end
    packagename(manifest::Dict{String,Any}, uuid::UUID) = packagename(manifest, string(uuid))

    function deps(uuid::UUID, c::Pkg.Types.Context)
        if any(p[1]["uuid"] == string(uuid) for (n, p) in manifest(c))
            return manifest(c)[string(uuid)][1].deps
        else
            return Dict{Any,Any}()
        end
    end
    deps(d::Dict{String,Any}) = get(d, "deps", Dict{String,Any}())
    deps(pe::PackageEntry) = get(pe[1], "deps", Dict{String,Any}())
    path(pe::PackageEntry) = get(pe[1], "path", nothing)
    version(pe::PackageEntry) = get(pe[1], "version", nothing)
    tree_hash(pe) = get(pe[1], "git-tree-sha1", nothing)

    frommanifest(c::Pkg.Types.Context, uuid) = frommanifest(manifest(c), uuid)

    function frommanifest(manifest::Dict{String,Any}, uuid)
        for p in values(manifest)
            if get(first(p), "uuid", "") == string(uuid)
                return (p)
            end
        end
        return nothing
    end
    is_package_deved(manifest, uuid) = get(first([p[2][1] for p in manifest if get(p[2][1], "uuid", "") == string(uuid)]), "path", "") != ""
else
    isinmanifest(context::Pkg.Types.Context, module_name::String) = any(p.name == module_name for (u, p) in manifest(context))
    isinmanifest(context::Pkg.Types.Context, uuid::UUID) = haskey(manifest(context), uuid)
    isinmanifest(manifest::Dict{UUID,PackageEntry}, uuid::UUID) = haskey(manifest, uuid)

    isinproject(context::Pkg.Types.Context, package_name::String) = haskey(deps(project(context)), package_name)
    isinproject(context::Pkg.Types.Context, package_uuid::UUID) = any(u == package_uuid for (n, u) in deps(project(context)))

    function packageuuid(c::Pkg.Types.Context, name::String)
        for pkg in manifest(c)
            if last(pkg).name == name
                return first(pkg)
            end
        end
    end
    packageuuid(pkg::Pair{String,UUID}) = last(pkg)
    packageuuid(pkg::Pair{UUID,PackageEntry}) = first(pkg)

    packagename(pkg::Pair{UUID,PackageEntry})::String = last(pkg).name
    packagename(c::Pkg.Types.Context, uuid::UUID) = manifest(c)[uuid].name
    packagename(manifest::Dict{UUID,PackageEntry}, uuid::UUID) = manifest[uuid].name

    function deps(uuid::UUID, c::Pkg.Types.Context)
        if haskey(manifest(c), uuid)
            return deps(manifest(c)[uuid])
        else
            return Dict{String,Base.UUID}()
        end
    end
    deps(pe::PackageEntry) = pe.deps
    deps(proj::Pkg.Types.Project) = proj.deps
    deps(pkg::Pair{String,UUID}, c::Pkg.Types.Context) = deps(packageuuid(pkg), c)
    path(pe::PackageEntry) = pe.path
    version(pe::PackageEntry) = pe.version
    version(pe::Pair{UUID,PackageEntry}) = last(pe).version
    frommanifest(c::Pkg.Types.Context, uuid) = manifest(c)[uuid]
    frommanifest(manifest::Dict{UUID,PackageEntry}, uuid) = manifest[uuid]
    tree_hash(pe::PackageEntry) = VERSION >= v"1.3" ? pe.tree_hash : get(pe.other, "git-tree-sha1", nothing)

    is_package_deved(manifest, uuid) = manifest[uuid].path !== nothing
end

function sha2_256_dir(path, sha=zeros(UInt8, 32))
    (uperm(path) & 0x04) != 0x04 && return
    startswith(path, ".") && return
    if isfile(path) && endswith(path, ".jl")
        s1 = open(path) do f
            sha2_256(f)
        end
        sha .+= s1
    elseif isdir(path)
        for f in readdir(path)
            sha = sha2_256_dir(joinpath(path, f), sha)
        end
    end
    return sha
end

function sha_pkg(pe::PackageEntry)
    path(pe) isa String && isdir(path(pe)) && isdir(joinpath(path(pe), "src")) ? sha2_256_dir(joinpath(path(pe), "src")) : nothing
end

function _doc(binding::Base.Docs.Binding)
    try
        sig = Union{}
        if Base.Docs.defined(binding)
            result = Base.Docs.getdoc(Base.Docs.resolve(binding), sig)
            result === nothing || return string(result)
        end
        results, groups = Base.Docs.DocStr[], Base.Docs.MultiDoc[]
    # Lookup `binding` and `sig` for matches in all modules of the docsystem.
        for mod in Base.Docs.modules
            dict = Base.Docs.meta(mod)::IdDict{Any,Any}
            if haskey(dict, binding)
                multidoc = dict[binding]
                push!(groups, multidoc)
                for msig in multidoc.order
                    sig <: msig && push!(results, multidoc.docs[msig])
                end
            end
        end
        if isempty(results)
            for group in groups, each in group.order
                push!(results, group.docs[each])
            end
        end
        md = try
            Base.Docs.catdoc(map(Base.Docs.parsedoc, results)...)
        catch
            nothing
        end
        return md === nothing ? "" : string(md)
    catch
        return ""
    end
end

_lookup(vr::FakeUnion, depot::EnvStore, cont=false) = nothing
_lookup(vr::FakeTypeName, depot::EnvStore, cont=false) = _lookup(vr.name, depot, cont)
_lookup(vr::FakeUnionAll, depot::EnvStore, cont=false) = _lookup(vr.body, depot, cont)
function _lookup(vr::VarRef, depot::EnvStore, cont=false)
    if vr.parent === nothing
        if haskey(depot, vr.name)
            val = depot[vr.name]
            if cont && val isa VarRef
                return _lookup(val, depot, cont)
            else
                return val
            end
        else
            return nothing
        end
    else
        par = _lookup(vr.parent, depot, cont)
        if par !== nothing && par isa ModuleStore && haskey(par, vr.name)
            val = par[vr.name]
            if cont && val isa VarRef
                return _lookup(val, depot, cont)
            else
                return val
            end
        else
            return nothing
        end
    end
end

maybe_lookup(x, env) = x isa VarRef ? _lookup(x, env, true) : x

"""
    maybe_getfield(k::Symbol , m::ModuleStore, server)

Try to get `k` from `m`. This includes: unexported variables, and variables
exported by modules used within `m`.
"""
function maybe_getfield(k::Symbol, m::ModuleStore, envstore)
    if haskey(m.vals, k)
        return m.vals[k]
    else
        for v in m.used_modules
            !haskey(m.vals, v) && continue
            submod = m.vals[v]
            if submod isa ModuleStore && k in submod.exportednames && haskey(submod.vals, k)
                return submod.vals[k]
            elseif submod isa VarRef
                submod = _lookup(submod, envstore, true)
                if submod isa ModuleStore && k in submod.exportednames && haskey(submod.vals, k)
                    return submod.vals[k]
                end
            end
        end
    end
end

function issubmodof(m::Module, M::Module)
    if m == M
        return true
    elseif parentmodule(m) === m
        return false
    elseif parentmodule(m) == M
        return true
    else
        return issubmodof(parentmodule(m), M)
    end
end

function Base.print(io::IO, f::FunctionStore)
    println(io, f.name, " is a Function.")
    nm = length(f.methods)
    println(io, "# $nm method", nm == 1 ? "" : "s", " for function ", f.name)
    for i = 1:nm
        print(io, "[$i] ")
        println(io, f.methods[i])
    end
end

const JULIA_DIR = normpath(joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia"))

function Base.print(io::IO, m::MethodStore)
    print(io, m.name, "(")
    for i = 1:length(m.sig)
        if m.sig[i][1] != Symbol("#unused#")
            print(io, m.sig[i][1])
        end
        print(io, "::", m.sig[i][2])
        i != length(m.sig) && print(io, ", ")
    end
    print(io, ")")
    path = replace(m.file, JULIA_DIR => "")
    print(io, " in ", m.mod, " at ", path, ':', m.line)
end

function Base.print(io::IO, t::DataTypeStore)
    print(io, t.name, " <: ", t.super)
    for i = 1:length(t.fieldnames)
        print(io, "\n  ", t.fieldnames[i], "::", t.types[i])
    end
end

Base.print(io::IO, m::ModuleStore) = print(io, m.name)
Base.print(io::IO, x::GenericStore) = print(io, x.name, "::", x.typ)

extends_methods(f) = false
extends_methods(f::FunctionStore) = f.name != f.extends
get_top_module(vr::VarRef) = vr.parent === nothing ? vr.name : get_top_module(vr.parent)

# Sorting is the main performance of calling `names`
unsorted_names(m::Module; all::Bool=false, imported::Bool=false) =
    ccall(:jl_module_names, Array{Symbol,1}, (Any, Cint, Cint), m, all, imported)

## recursive_copy
#
# `deepcopy` is reliable but incredibly slow. Its slowness comes from two factors:
# - generically iterating over, e.g., `fieldnames(typeof(x))` rather than having a method
#   optimized for each struct type
# - its care to protect against circular depenency graphs
# When you don't need to worry about cycles, you can do much better by defining your own function.

recursive_copy(x) = deepcopy(x)

recursive_copy(::Nothing) = nothing

recursive_copy(s::Symbol) = s

recursive_copy(c::Char) = c

recursive_copy(str::String) = str

recursive_copy(x::Number) = x

recursive_copy(p::Pair) = typeof(p)(recursive_copy(p.first), recursive_copy(p.second))

recursive_copy(A::Array) = eltype(A)[recursive_copy(a) for a in A]

recursive_copy(d::Dict) = typeof(d)(recursive_copy(p) for p in d)

recursive_copy(ref::VarRef) = VarRef(recursive_copy(ref.parent), ref.name)

recursive_copy(tn::FakeTypeName) = FakeTypeName(recursive_copy(tn.name), recursive_copy(tn.parameters))

recursive_copy(tb::FakeTypeofBottom) = tb

recursive_copy(u::FakeUnion) = FakeUnion(recursive_copy(u.a), recursive_copy(u.b))

recursive_copy(tv::FakeTypeVar) = FakeTypeVar(tv.name, recursive_copy(tv.lb), recursive_copy(tv.ub))

recursive_copy(ua::FakeUnionAll) = FakeUnionAll(recursive_copy(ua.var), recursive_copy(ua.body))

@static if !(Vararg isa Type)
    function recursive_copy(va::FakeTypeofVararg)
        if isdefined(va, :N)
            FakeTypeofVararg(recursive_copy(va.T), va.N)
        elseif isdefined(va, :T)
            FakeTypeofVararg(recursive_copy(va.T))
        else
            FakeTypeofVararg()
        end
    end
end

recursive_copy(m::ModuleStore) = ModuleStore(recursive_copy(m.name), recursive_copy(m.vals), m.doc,
                                             m.exported, copy(m.exportednames), copy(m.used_modules))

recursive_copy(p::Package) = Package(p.name,
                                     recursive_copy(p.val),
                                     p.uuid,
                                     recursive_copy(p.sha))

recursive_copy(ms::MethodStore) = MethodStore(ms.name,
                                              ms.mod,
                                              ms.file,
                                              ms.line,
                                              recursive_copy(ms.sig),
                                              copy(ms.kws),
                                              recursive_copy(ms.rt))

recursive_copy(dts::DataTypeStore) = DataTypeStore(recursive_copy(dts.name),
                                                   recursive_copy(dts.super),
                                                   recursive_copy(dts.parameters),
                                                   recursive_copy(dts.types),
                                                   recursive_copy(dts.fieldnames),
                                                   recursive_copy(dts.methods),
                                                   dts.doc,
                                                   dts.exported)

recursive_copy(fs::FunctionStore) = FunctionStore(recursive_copy(fs.name),
                                                  recursive_copy(fs.methods),
                                                  fs.doc,
                                                  recursive_copy(fs.extends),
                                                  fs.exported)

recursive_copy(gs::GenericStore) = GenericStore(recursive_copy(gs.name),
                                                recursive_copy(gs.typ),
                                                gs.doc,
                                                gs.exported)

# Tools for modifying source location
# env = getenvtree([:somepackage])
# symbols(env, somepackage)
# m = env[:somepackage]
# To strip actual src path:
# modify_dirs(m, f -> modify_dir(f, pkg_src_dir(somepackage), "PLACEHOLDER"))
# To replace the placeholder:
# modify_dirs(m, f -> modify_dir(f, "PLACEHOLDER", new_src_dir))
function modify_dirs(m::ModuleStore, f)
    for (k, v) in m.vals
        if v isa FunctionStore
            m.vals[k] = FunctionStore(v.name, MethodStore[MethodStore(m.name, m.mod, f(m.file), m.line, m.sig, m.kws, m.rt) for m in v.methods], v.doc, v.extends, v.exported)
        elseif v isa DataTypeStore
            m.vals[k] = DataTypeStore(v.name, v.super, v.parameters, v.types, v.fieldnames, MethodStore[MethodStore(m.name, m.mod, f(m.file), m.line, m.sig, m.kws, m.rt) for m in v.methods], v.doc, v.exported)
        elseif v isa ModuleStore
            modify_dirs(v, f)
        end
    end
end

pkg_src_dir(m::Module) = dirname(pathof(m))

# replace s1 with s2 at the start of a string
function modify_dir(f, s1, s2)
    # @assert startswith(f, s1)
    # Removed assertion because of Enums issue
    string(s2, f[length(s1)+1:end])
end

# tools to retrieve cache from the cloud

function get_file_from_cloud(manifest, uuid, environment_path, depot_dir, cache_dir = "../cache", download_dir = "../downloads/", symbolcache_upstream = "https://www.julia-vscode.org/symbolcache")
    paths = get_cache_path(manifest, uuid)
    name = packagename(manifest, uuid)
    link = string(first(splitext(join([symbolcache_upstream, "store/v1/packages", paths...], '/'))), ".tar.gz")

    dest_filepath = joinpath(cache_dir, paths...)
    dest_filepath_unavailable = string(first(splitext(dest_filepath)), ".unavailable")
    download_dir = joinpath(download_dir, first(splitext(last(paths))))
    download_filepath = joinpath(download_dir, last(paths))
    download_filepath_unavailable = string(first(splitext(download_filepath)), ".unavailable")

    @debug "Downloading cache file for $name."
    if isfile(dest_filepath_unavailable)
        @info "Cloud was unable to cache $name in the past, we won't try to retrieve it again."
        return false
    end
    file = try
        if Pkg.PlatformEngines.download_verify_unpack(link, nothing, download_dir)
            mkpath(dirname(dest_filepath))
            if !isfile(download_filepath) && isfile(download_filepath_unavailable)
                mv(download_filepath_unavailable, dest_filepath_unavailable)
                @info "Cloud is unable to cache $name, we won't try to retrieve it again."
                return false
            end
            mv(download_filepath, dest_filepath)
            dest_filepath
        else
            nothing
        end
    catch err
        @info "Couldn't retrieve cache file for $name." exception = err
        return false
    end

    if file === nothing
        @info "Couldn't retrieve cache file for $name."
        return false
    end

    cache = try
        open(file, "r") do io
            CacheStore.read(io)
        end
    catch
        @info "Couldn't read cache file for $name, deleting."
        rm(file)
        return false
    end

    pkg_path = Base.locate_package(Base.PkgId(uuid, name))
    if pkg_path === nothing || !isfile(pkg_path)
        pkg_path = get_pkg_path(Base.PkgId(uuid, name), environment_path, depot_dir)
    end
    if pkg_path === nothing
        @info "Successfully downloaded and saved $(name), but with placeholder paths"
        return false
    end

    @debug "dirname" dirname(pkg_path)
    modify_dirs(cache.val, f -> modify_dir(f, "PLACEHOLDER", joinpath(pkg_path, "src")))
    open(file, "w") do io
        CacheStore.write(io, cache)
    end

    @info "Successfully downloaded, scrubbed and saved $(name)"
    return true
end

"""
    validate_disc_store(store_path, manifest)

This returns a list of non-jll packages in the manifest that don't have caches on disc.
"""
function validate_disc_store(store_path, manifest)
    filter(manifest) do pkg
        uuid = packageuuid(pkg)
        endswith(packagename(manifest, uuid), "_jll") && return false

        file_name = joinpath(get_cache_path(manifest, uuid)...)
        yield()
        return !isfile(joinpath(store_path, file_name))
    end
end

function find_project_file(env)
    isdir(env) || return false
    for filename in ("Project.toml", "JuliaProject.toml")
        maybe_project_file = joinpath(env, filename)
        if isfile(maybe_project_file)
            return maybe_project_file
        end
    end
    return false
end

"""
    get_pkg_path(pkg::Base.PkgId, env, depot_path)

Find out where a package is installed without having to load it.
"""
function get_pkg_path(pkg::Base.PkgId, env, depot_path)
    project_file = find_project_file(env)
    project_file isa Bool && return nothing
    manifest_file = Base.project_file_manifest_path(project_file)

    d = Base.parsed_toml(manifest_file)
    entries = get(d, pkg.name, nothing)::Union{Nothing, Vector{Any}}
    entries === nothing && return nothing # TODO: allow name to mismatch?
    for entry in entries
        entry = entry::Dict{String, Any}
        uuid = get(entry, "uuid", nothing)::Union{Nothing, String}
        uuid === nothing && continue
        if UUID(uuid) === pkg.uuid
            path = get(entry, "path", nothing)::Union{Nothing, String}
            # this can only be true for explicitly dev'ed packages
            if path !== nothing
                path = normpath(abspath(dirname(manifest_file), path))
                return path
            end
            hash = get(entry, "git-tree-sha1", nothing)::Union{Nothing, String}
            hash === nothing && return nothing
            hash = Base.SHA1(hash)
            # empty default path probably means that we should use the default Julia depots
            if depot_path == ""
                depot_paths = []
                Base.append_default_depot_path!(depot_paths)
            else
                depot_paths = [depot_path]
            end
            for depot in depot_paths
                # Keep the 4 since it used to be the default
                for slug in (Base.version_slug(pkg.uuid, hash, 4), Base.version_slug(pkg.uuid, hash))
                    path = abspath(depot, "packages", pkg.name, slug)
                    ispath(path) && return path
                end
            end
            return nothing
        end
    end
    return nothing
end

function load_package(c::Pkg.Types.Context, uuid, conn, loadingbay, percentage = missing)
    isinmanifest(c, uuid isa String ? Base.UUID(uuid) : uuid) || return
    pe_name = packagename(c, uuid)

    pid = Base.PkgId(uuid isa String ? Base.UUID(uuid) : uuid, pe_name)
    if pid in keys(Base.loaded_modules)
        conn !== nothing && println(conn, "PROCESSPKG;$pe_name;$uuid;noversion;$percentage")
        loadingbay.eval(:($(Symbol(pe_name)) = $(Base.loaded_modules[pid])))
        m = getfield(loadingbay, Symbol(pe_name))
    else
        m = try
            conn !== nothing && println(conn, "STARTLOAD;$pe_name;$uuid;noversion;$percentage")
            loadingbay.eval(:(import $(Symbol(pe_name))))
            conn !== nothing && println(conn, "STOPLOAD;$pe_name")
            m = getfield(loadingbay, Symbol(pe_name))
        catch
            return
        end
    end
end

function write_cache(uuid, pkg::Package, outpath)
    mkpath(dirname(outpath))
    @info "Now writing to disc $uuid"
    open(outpath, "w") do io
        CacheStore.write(io, pkg)
    end
    outpath
end

"""
    get_cache_path(manifest, uuid)

Returns a vector containing the cache storage path for a package structured: [folder, folder, file].
"""
function get_cache_path(manifest, uuid)
    name = packagename(manifest, uuid)
    pkg_info = frommanifest(manifest, uuid)
    ver = version(pkg_info)
    ver = ver === nothing ? "nothing" : ver
    ver = replace(string(ver), '+'=>'_')
    th = tree_hash(pkg_info)
    th = th === nothing ? "nothing" : th

    [
        string(uppercase(string(name)[1]))
        string(name, "_", uuid)
        string("v", ver, "_", th, ".jstore")
    ]
end

function write_depot(server::Server, ctx, written_caches)
    for (uuid, pkg) in server.depot
        cache_paths = get_cache_path(manifest(ctx), uuid)
        outpath = joinpath(server.storedir, cache_paths...)
        outpath in written_caches && continue

        written_path = write_cache(uuid, pkg, outpath)
        !isempty(written_path) && push!(written_caches, written_path)
    end
end


text:
import DataFrames.describe

export acceptance, describe

# Compute acceptance or rejection rate of MCMCChain
function acceptance(c::MCMCChain; lags::Ranges=1:nrow(c.samples), reject::Bool=false)
  rlen = lags.len
  @assert lags[end] <= nrow(c.samples) "Range of acceptance rate not within post-burnin range of MCMC chain"

  if reject
    return (rlen-sum(c.diagnostics["accept"][lags]))*100/rlen
  else
    return sum(c.diagnostics["accept"][lags])*100/rlen
  end
end

# TODO 1: Compute MCMC quantiles based on
# Flegal J.M, Galin L.J, Neath R.C. Markov Chain Monte Carlo Estimation of Quantiles. arXiv, 2013
# TODO 2: Include these MCMC estimates of quantiles in describe()

# describe() provides summary statistics for MCMCChain objects, similarly to dataframes' describe()
describe(c::MCMCChain) = describe(STDOUT, c)

function describe(io, c::MCMCChain)
  for i in 1:ncol(c.samples)

    col = c.samples[i]
    println(io, colnames(c.samples)[i])

    if all( isna(col) )
      println(col, " * All NA * ")
        return
    end

    filtered = float(removeNA(col))
    qs = quantile(filtered, [0, 1])
    nrows = nrow(c.samples)
    varimse = mcvar_imse(filtered)
    variid = var(filtered)/nrows
    mcerror = sqrt(varimse)
    ss = nrows*variid./varimse
    act = varimse./variid

    statNames = ["Min", "Mean", "Max", "MC Error", "ESS", "AC Time"]
    statVals = [qs[1], mean(filtered), qs[2], mcerror, ss, act]
    for i = 1:6
        println(io, string(rpad(statNames[i], 10, " "), " ", string(statVals[i])))
    end
    nas = sum(isna(col))
    println(io, "NAs        $nas")
    println(io, "NA%        $(round(nas*100/length(col), 2))%")

    println(io, )
  end
end


text:
#=

  Euler #12 in Julia.

  Problem 12
  """
  The sequence of triangle numbers is generated by adding the natural numbers.
  So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.
  The first ten terms would be:

  1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

  Let us list the factors of the first seven triangle numbers:

       1: 1
       3: 1,3
       6: 1,2,3,6
      10: 1,2,5,10
      15: 1,3,5,15
      21: 1,3,7,21
      28: 1,2,4,7,14,28

  We can see that the 7th triangle number, 28, is the first triangle number
  to have over five divisors.

  Which is the first triangle number to have over five-hundred divisors?")
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

#
# This is basically factors + collect together so it's not general
# enough for placing in Euler.jl
#
function factors_map(n::Int64)
    if n == one(n)
        return Dict{Int64,Int64}(1=>1)
    end
    m = Dict{Int64,Int64}();
    while n % 2 === zero(n)
        if !haskey(m,2)
            m[2] = zero(0)
        end
        m[2] += one(n);
        # n = round(Int, n/2);
        n = n//2
    end
    t = 3
    while n > 1 && t < ceil(Int,sqrt(n))
        while n % t == 0
            if !haskey(m,t)
                m[t] = 0
            end
            m[t] += 1
            n //= t
        end
        t += 2
    end
    if n > 1
        if !haskey(m,n)
            m[n] = 0
        end
        m[n] += 1
    end
    return m;
end

# Brute force
# 459.90159305s !
 function euler12a_dont_run()
    len = 0
    i = 0
    tnum = 0
    while len <= 500
        i+=1
        tnum += i
        divs = all_divisors(tnum)
        len = length(divs)
        if i % 1000 === 0
            println("i:$i tnum:$tnum len: $len")
        end
    end
    return tnum;
end

# 0.13744243s
function euler12b()
    len = 0
    i = 0
    tnum = 0
    while len <= 500
        i+=1
        tnum += i
        f = factors_map(tnum)
        len = prod(values(f).|>(i->i+1))
    end
    return tnum
end

# 0.16336305s
function euler12c()
    len = 0
    i = 0
    tnum = 0
    while len <= 500
        i+=1
        tnum += i
        len = prod(values(collect_list(factors(tnum))).|>j->j+1)
    end
    return tnum
end

# run_euler(euler12a); # DON'T RUN! Too slow.
run_euler(euler12b);
# run_euler(euler12c);

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using SymPy
using Test

@testset "ODes" begin
    ## ODEs
    x, a = Sym("x, a")
    F = SymFunction("F")
    ex = diff(F(x), x) - a*F(x)
    ex1 = dsolve(ex)
    ex2 = rhs(ex1) |> subs(Sym(:C1), 1) |> subs(a, 2)
    @test ex2 == exp(2x)

    t, = @syms t
    X, Y = SymFunction("X, Y")
    eq = [Eq(diff(X(t),t), 12*t*X(t) + 8*Y(t)), Eq(diff(Y(t),t), 21*X(t) + 7*t*Y(t))]
    dsolve(eq)

    ## version 0.4+ allow use of u'(x) in lieu of diff(u(x), x) and `ivpsolve`
    u = SymFunction("u")
    a, x, y, y0, y1 = symbols("a, x, y, y0, y1")

    @test dsolve(u'(x) - a*u(x), x, (u, 0, 1)) == Eq(u(x), exp(a*x))
    @test dsolve(u'(x) - a*u(x), x, (u, 0, y1)) == Eq(u(x), y1*exp(a*x))
    dsolve(u'(x) - a*u(x), x, (u, y0, y1)) # == Eq(u(x), y1 * exp(a*(x - y0)))
    dsolve(x*u'(x) + x*u(x) + 1, x, (u, 1, 1))
    dsolve((u'(x))^2 - a*u(x), x, (u, 0, 1))
    dsolve(u''(x) - a * u(x), x, (u, 0, 1), (u', 0, 0))

    F, G, K = SymFunction("F, G, K")
    eqn = F(x)*u'(y)*y + G(x)*u(y) + K(x)
    dsolve(eqn, y, (u, 1, 0))

    ## dsolve eqn has two answers, but we want to eliminate one
    # based on initial condition
    dsolve(u'(x) - (u(x)-1)*u(x)*(u(x)+1), x, (u, 0, 1//2))
end


text:
import Base.copy
export MPProblem, MPSolution
export clearsamples!

type MPSolution{T}
    status::Symbol
    cost::T
    elapsed::Float64
    metadata::Dict
end

type MPProblem{T<:AbstractFloat}
    SS::StateSpace
    init::State
    goal::Goal
    CC::CollisionChecker
    V::SampleSet
    config_name::String
    status::Symbol
    solution::MPSolution{T}

    function MPProblem(SS::StateSpace,
                       init::State,
                       goal::Goal,
                       CC::CollisionChecker,
                       V::SampleSet,
                       config_name::String="$(dim(SS))D $(typeof(SS))")
        new(SS, init, goal, CC, V, config_name, "not yet solved")
    end
end
changeprecision{T<:AbstractFloat}(::Type{T}, P::MPProblem) =
    MPProblem(map(x -> changeprecision(T,x), (P.SS, P.init, P.goal, P.CC))...)  # TODO: changeprecision for `SampleSet`s

MPProblem{S}(SS::StateSpace{S}, init::State, goal::Goal, CC::CollisionChecker) =
    MPProblem{eltype(S)}(SS, S(init), goal, CC, defaultNN(SS, S(init)))
function copy(P::MPProblem)
    Pcopy = MPProblem(P.SS, P.init, P.goal, P.CC, P.V, P.config_name)
    Pcopy.status = P.status
    Pcopy.solution = P.solution
    Pcopy
end
function clearsamples!(P::MPProblem)
    P.V = defaultNN(P.SS, P.init)
end

plot_path(SS::StateSpace, V::SampleSet, sol; kwargs...) = plot_path(V[sol], SS; kwargs...)
plot_tree(SS::StateSpace, V::SampleSet, A; kwargs...) = plot_tree(V.V, A, SS; kwargs...)

function plot(P::MPProblem; SS=true, CC=true, goal=true, meta=false, sol=true, smoothed=false)
    SS && plot(P.SS)
    CC && plot(P.CC, P.SS.lo, P.SS.hi)
    goal && plot(P.goal, P.SS)
    if isdefined(P, :solution)
        S = P.solution
        if meta
            haskey(S.metadata, "tree") && plot_tree(P.SS, P.V, S.metadata["tree"], color="gray", alpha=0.5)
            # TODO: graph (PRM)
        end
        sol && plot_path(P.SS, P.V, S.metadata["path"], color="blue")
        smoothed && haskey(S.metadata, "smoothed_path") && plot_path(S.metadata["smoothed_path"], color="orange")
        P.status, P.status == :solved ? P.solution.cost : Inf
    else
        "unsolved"
    end
end


text:
#=
    The file contains the implementation of the FGLM algorithm,
    an algorithm for ordering conversion
=#

mutable struct NextMonomials
    # monomials to check
    monoms::Vector{Int}
    load::Int
    done::Dict{Int, Int}
    ord::Symbol
end

Base.isempty(m::NextMonomials) = m.load == 0

function initialize_nextmonomials(ht, ord)
    vidx = insert_in_hash_table!(ht, zero(ht.exponents[1]))
    monoms = Vector{Int}(undef, 2^3)
    monoms[1] = vidx
    load = 1
    NextMonomials(monoms, load, Dict{Int, Int}(vidx => 1), ord)
end

function insertnexts!(
            m::NextMonomials,
            ht::MonomialHashtable, monom::Int)
    while m.load + ht.nvars >= length(m.monoms)
        resize!(m.monoms, length(m.monoms) * 2)
    end

    emonom = ht.exponents[monom]
    for i in 1:ht.nvars
        eprod = copy(emonom)
        eprod[i] += 1
        eprod[end] += 1
        vidx = insert_in_hash_table!(ht, eprod)

        if !haskey(m.done, vidx)
            m.load += 1
            m.monoms[m.load] = vidx
            m.done[vidx] = 1
        end
    end

    sort_monoms_decreasing!(m.monoms, m.load, ht, m.ord)
end

function nextmonomial!(m::NextMonomials)
    # assuming m.monoms is sorted (reversed)
    monom = m.monoms[m.load]
    m.load -= 1
    monom
end

function add_generator!(basis::Basis{C}, matrix, relation, ht, ord) where {C<:Coeff}
    rexps, rcoeffs, _ = extract_sparse_row(relation)

    if debug()
        @warn "add generator"
        println(relation)
        println(rexps, " ", rcoeffs)
    end

    for i in 1:length(rexps)
        rexps[i] = matrix.rightcol2hash[rexps[i]]
    end

    sort_terms_decreasing!(rexps, rcoeffs, ht, ord)

    if debug()
        @warn "extracted"
        println(rexps)
        println(rcoeffs)
    end

    check_enlarge_basis!(basis, 1)
    basis.ndone += 1
    basis.nlead += 1
    basis.nonred[basis.nlead] = basis.nlead
    basis.gens[basis.ndone] = rexps
    basis.coeffs[basis.ndone] = rcoeffs
end

function divides_staircase(monom, staircase, ht)
    for m in staircase
        if debug()
            @warn "uwu"
            println("is $(ht.exponents[m]) divisible by $(ht.exponents[monom])")
        end
        if is_monom_divisible(monom, m, ht)
            return true
        end
    end
    false
end

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

function fglm_f4!(
            ring::PolyRing,
            basis::Basis{C},
            ht::MonomialHashtable,
            ord::Symbol) where {C<:Coeff}

    newbasis = initialize_basis(ring, basis.ntotal, C)
    nextmonoms = initialize_nextmonomials(ht, ord)
    matrix = initialize_double_matrix(basis)
    staircase = Int[]

    while !isempty(nextmonoms)
        monom = nextmonomial!(nextmonoms)

        if debug()
            @warn "new iteration" ht.exponents[monom]
            @warn "" ht
            println("#############")
            println(nextmonoms)
            println(newbasis.gens, " ", newbasis.coeffs)
        end

        if divides_staircase(monom, staircase, ht)
            continue
        end

        tobereduced    = initialize_basis(ring, [[monom]], [C[1]])
        tobereduced.ntotal = 1

        # compute normal form
        normal_form_f4!(ring, basis, ht, tobereduced)

        if debug()
            println("Normal form ", tobereduced.gens[1], " " ,tobereduced.coeffs[1])
        end

        # matrix left rows can express tobereduced?
        # reduces monom and tobereduced
        exists, relation = linear_relation!(matrix, monom, tobereduced, ht)

        if debug()
            @warn "dump"
            dump(matrix, maxdepth=2)
            println(exists)
        end

        # if linear relation between basis elements exists
        if exists
            lead = ht.exponents[monom]

            if debug()
                @error "produced element" lead
            end

            add_generator!(newbasis, matrix, relation, ht, ord)
            push!(staircase, monom)
        else
            insertnexts!(nextmonoms, ht, monom)
        end
    end

    standardize_basis!(newbasis, ht, ord)

    linbasis = extract_linear_basis(ring, matrix)

    # println(linbasis)
    # println(newbasis)

    newbasis, linbasis, ht
end

function fglm_f4(
            ring::PolyRing,
            basisexps::Vector{Vector{ExponentVector}},
            basiscoeffs::Vector{Vector{C}},
            metainfo) where {Rng, C<:Coeff}

    tablesize = select_tablesize(ring, basisexps)

    basis, ht = initialize_structures(ring, basisexps,
                                        basiscoeffs, metainfo.rng, tablesize)

    basis, linbasis, ht = fglm_f4!(ring, basis, ht, metainfo.computeord)

    export_basis_data(basis, ht)
end

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

function extract_linear_basis(ring, matrix::DoubleMacaulayMatrix{C}) where {C}
    exps = Vector{Vector{Int}}(undef, matrix.nrrows)
    coeffs = Vector{Vector{C}}(undef, matrix.nrrows)

    # println("Extracting")
    # dump(matrix, maxdepth=2)

    for i in 1:matrix.nrrows
        exps[i] = matrix.rightrows[i]
        coeffs[i] = matrix.rightcoeffs[i]
        for j in 1:length(exps[i])
            exps[i][j] = matrix.rightcol2hash[exps[i][j]]
        end
    end

    linbasis = initialize_basis(ring, exps, coeffs)

    linbasis.ndone = length(exps)
    linbasis.nlead = length(exps)
    linbasis.nonred = collect(1:linbasis.ndone)

    linbasis
end

function kbase_f4(
            ring::PolyRing,
            basisexps::Vector{Vector{ExponentVector}},
            basiscoeffs::Vector{Vector{C}},
            metainfo) where {Rng, C<:Coeff}

    tablesize = select_tablesize(ring, basisexps)

    basis, ht = initialize_structures(ring, basisexps,
                                        basiscoeffs, metainfo.rng, tablesize)

    basis, linbasis, ht = fglm_f4!(ring, basis, ht, metainfo.computeord)

    export_basis_data(linbasis, ht)
end


text:
"""
    MOSUM(data, G[, var_est_method="mosum", alpha=0.1, criterion="eta", eta=0.4, epsilon=0.2])

Runs the MOSUM procedure for the univariate array `data` with bandwidth `G`, and returns the number and position of found changepoints.

Optionally, `var_est_method` specifies the variance estimator to normalise by; this can be the average `mosum` (default) or minimum `mosum.min` across windows.
`alpha` determines the signicance level (default 0.1).
`criterion` determines whether to use the `eta` (default) or `epsilon` location procedure (see references).
`eta` and `epsilon` are tuning parameters for the mentioned procedures (default 0.4 and 0.2).

See also: [`@BS`](@ref), [`MOSUM_multi_scale`](@ref), [`mosum_plot`](@ref)

# Returns
* `out::Dict{String,Union{Bool, Float64, Int64, Array}}`:
  * `out["number"]` : Integer number of detected changes
  * `out["Reject"]` : Boolean hypothesis test outcome
  * `out["threshold"]` : Float of testing threshold
  * `out["detector"]` : vector of mosum detector values
  * `out["var.estimation"]` : vector of estimated variance values

# Example
```julia-repl
# Sample Normal time series with changing mean
n = 1000
λ = 300
μ, σ = Normal(0.0, 10.0), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)
G = 150
# Run MOSUM procedure
MOSUM_out = MOSUM(data, G)
# Plot MOSUM detector
mosum_plot(MOSUM_out)
```

# References
Eichinger, Birte, and Claudia Kirch. "A MOSUM procedure for the estimation of multiple random change points." Bernoulli 24.1 (2018): 526-564.
Meier, Aledataander, Claudia Kirch, and Haeran Cho. "mosum: A package for moving sums in change point analysis." (2018).
"""
function MOSUM( data::Array{Float64}, G::Int64; var_est_method="mosum", alpha=0.1, criterion="eta", eta=0.4, epsilon=0.2)
    var_est_method in ("mosum", "mosum.min") || error("Keyword argument var_est_method must be either \"mosum\" or \"mosum.min\"")

    CP = Array{Int64}(undef,0)
    q = 0
    Reject = false
    n = length(data)

    csum = cumsum(data)
    sqsum = cumsum(data.^2)
    variance = zeros(Float64, n)
    lvar = zeros(Float64, n)
    rvar = zeros(Float64, n)
    lsqsums = zeros(Float64, n)
    rsqsums =  zeros(Float64, n)
    detector = zeros(Float64, n)
    #jump_size = zeros(Float64, n)
    lsums = zeros(Float64, n)
    rsums = zeros(Float64, n)

    #populate
    for t in 1:(n-G)
        rsums[t] = csum[t+G] - csum[t]
        lsums[t+G] = rsums[t]

        rsqsums[t]  = sqsum[t+G] - sqsum[t]
        lsqsums[t+G] = rsqsums[t]

        rlmean = rsums[t]/G
        rvar[t] = rsqsums[t]/G -  (1/G - 2)*(rlmean^2)
        lvar[t+G] = rvar[t]
    end

    # normalise
    for t in (G+1):(n-G)
        if var_est_method == "mosum"
            variance[t] = (rvar[t] + lvar[t])/ 2
        elseif var_est_method == "mosum.min"
            variance[t] = min( rvar[t], lvar[t])
        end
        #jump_size[t] =
        detector[t] = abs(rsums[t] - lsums[t])/ sqrt(2*G*variance[t])
    end

    #threshold
    a = sqrt(2*n/G); b = 2log(n/G) + log(log(n/G))/2 + log(3/2) - log(pi)/2
    c = - log(log( (1- alpha)^(-.5) ))
    D = (b+c)/a

    mosum_stat, k = findmax(detector)

    if mosum_stat > D
        got_cps = get_cps_mosum(detector, D, G, criterion, eta, epsilon)
        CP = got_cps[1]
        q = got_cps[2]
        Reject = got_cps[3]
    end

#    cost = 0.0
#    CP = sort(CP)
#    for j in 1:(length(CP)-1)
#        cost = cost + segment_cost(CP[j]+1,CP[j+1]) + pen
#    end
#    cost = cost + segment_cost(CP[end]+1,n)
    # organise output into a dictionary
    out = Dict{String, Union{Array, Bool, Int64, Float64} }()
    out["number"] = q
    out["detector"] = detector
    out["Reject"] = Reject
    out["changepoints"] = CP
    out["threshold"] = D
    out["var.estimation"] = variance
    #multi-scale
    #if(return_jump_size) out["jump_size"] = jump_size
    return out
end

function get_cps_mosum(stat::Array{Float64} , D_n::Float64, G::Int64,  criterion::String="epsilon",
                       eta::Float64 = 0.4, epsilon::Float64 = 0.2 )
    cps = Array{Int64}(undef,0) ##assign empty cps
    Reject = false
    if criterion == "epsilon"
        n = length(stat)
        rshift = stat[2:end]
        push!(rshift, 0.0)
        lshift = stat[1:(end-1)]
        pushfirst!(lshift, 0.0)
        over =  stat .> D_n #findall(data -> data > D_n , stat ) #indices greater than D_n
        vv = over .& (lshift .< D_n )#& lshift >0 ) #lowers
        v = findall(vv .== 1)
        ww = over .& (rshift .< D_n)
        w = findall(ww .== 1)
        nu_remove = findall(w-v .>= epsilon * G) #nu(epsilon) test for distance between
        v_nu = v[nu_remove]; w_nu = w[nu_remove] #c(w[nu_remove],n)
        sub_pairs = [v_nu w_nu]

        q = size(sub_pairs)[1]
        if q>0
            for ii in 1:q
                interval = sub_pairs[ii,1]:sub_pairs[ii,2]
                kk = findmadata(stat[interval])[2] #internal cp location
                push!(cps, kk + 1 + sub_pairs[ii,1]) #- G-p
            end
        end
    end

    if criterion == "eta"
        n = length(stat)
        window = Int(floor(eta*G))
        for t in (G+1):(n-G)
            if (stat[t] > D_n) & (stat[t] == maximum(stat[(t .- window):(t .+ window)]) )
                push!(cps, t) ##add to cps
            end
        end
    end

    q = length(cps)
    if q>0
        Reject = true
    end

    return cps, q, Reject
end

"""
    MOSUM_multi_scale(data, Gset[, var_est_method="mosum", alpha=0.1, criterion="eta", eta=0.4, epsilon=0.2])

Runs the Multiple Filtre MOSUM procedure for the univariate array `data` with multiple bandwidths `Gset`, and returns the position of found changepoints.
For optional arguments, see `?MOSUM`.

See also: [`MOSUM`](@ref)

# Returns
* `cps` : Integer array of estimated change points

# Edataample
```julia-repl
# Sample Normal time series with changing mean
n = 1000
λ = 300
μ, σ = Normal(0.0, 10.0), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)
Gset = [25, 50, 100, 150]
# Run MOSUM procedure
multi_scale_out = MOSUM_multi_scale(data, Gset; alpha = 0.05)
# Plot change points
changepoint_plot(data, multi_scale_out)
```

# References
Messer M, Kirchner M, Schiemann J, Roeper J, Neininger R, Schneider G (2014). “A Multiple Filter Test for the Detection of Rate Changes in Renewal Processes with Varying Variance.” The Annals of Applied Statistics, 8(4), 2027–2067.
"""
function MOSUM_multi_scale(data::Array{Float64}, Gset::Array{Int64}; kwargs...)

    sort!(Gset)
    cps = Array{Int64}[]
    for G in Gset
        mosum_G = MOSUM(data, G; kwargs...)
        cps_G = mosum_G["changepoints"]
        if G == Gset[1]
            cps = cps_G
        else
            for k in cps_G
                min_dist = minimum(abs.(cps .- k))
                if min_dist >= G
                    push!(cps, k)
                end
            end # for
        end
    end

#    #step 1
#    PairSet = Array{Union{Float64, Int64}}[]
#    for G in Gset
#        mosum_G = MOSUM(data, G; default_args)
#        cps = mosum_G["changepoints"]
#        for k in cps
#            h = mosum_G["detector"][k]
#            push!(PairSet, [h, k, G])
#        end # for
#    end
#    sort!(Pairset)
#    #step 2
#    D = Array{Union{Float64, Int64}}[] #conflicting
#    Cset = PairSet
#    while length(PairSet) > 0
#        candidate = pop!(Pairset) #current max jump
#        for i in 1:length(PairSet)
#            if  abs(PairSet[i,2] -. candidate[2]) < min(candidate[3], PairSet[i,3])
#                take = popat(PairSet, i)
#                push!(D, take)
#            end
#        end
#    end
    return(cps)
    #2
end