text:
module Mixtures
# My own implementation of mixture models (hard and soft, currently Gaussian only) with split-merge EM
using Distributions, Base.LinAlg, PDMats, StatsBase, StatsFuns, Clustering, Roots
import Base.split, Base.merge
export TMixture, em!, fsmem!, split, merge, score, assignments

"""
    verbose_iter(::Symbol)

Returns true if the verbosity level indicates we should print the log
likelihood on each iteration.
"""
verbose_iter(s::Symbol) = s == :iter

"""
    verbose_op(::Symbol)

Returns true if the verbosity level indicates we should print
information about split and merge operations.
"""
verbose_op(s::Symbol) = s == :iter || s == :op

function regularize_and_factorize!(c::Union{FullNormal,MvTDist}, λ::Real)
    Σmat = c.Σ.mat
    if λ != 0
        for i = diagind(Σmat)
            Σmat[i] += λ
        end
    end
    copy!(c.Σ.chol.factors, Σmat)
    LinAlg.chol!(c.Σ.chol.factors, Val{:U})
    c
end

function weighted_mean_cov!(mu, Σ, tmp, x, w)
    m = size(x, 1)
    n = size(x, 2)
    inv_sw = 1.0 / sum(w)
    Base.LinAlg.BLAS.gemv!('N', inv_sw, x, w, 0.0, mu)

    for j = 1:n
        cj = sqrt(w[j])
        for i = 1:m
            @inbounds tmp[i,j] = (x[i,j] - mu[i]) * cj
        end
    end
    LinAlg.copytri!(Base.LinAlg.BLAS.syrk!('U', 'N', inv_sw, tmp, 0.0, Σ), 'U')
end

type TMixture{T,M<:StridedMatrix}
    components::Vector{MvTDist}
    logπ::Vector{T}

    uniform_logpdf::T
    uniform_logπ::T

    X::M
    assignments::Vector{UInt8}
    λ::T
    fix_ν::Bool
    ll::T

    assignedu::Vector{T}
    maxtmp::Vector{T}
    mahaltmp::Vector{T}
    Xtmp::Matrix{T}
end

function compute_uniform_pdf(X)
    maxes = maximum(X, 2)
    mins = minimum(X, 2)
    -sum(log(maxes - mins))
end

function TMixture{T}(components::Vector{MvTDist}, X::StridedMatrix{T};
                     use_uniform::Bool=false, λ::T=convert(T, 0),
                     logπ::Vector{T}=fill(log(size(X, 1)/(size(X, 1)*length(components)+use_uniform)), length(components)),
                     fix_ν::Bool=false)
    length(components) < typemax(UInt8) || throw(ArgumentError("a maximum of $(typemax(UInt8)) components are supported"))
    TMixture(components, logπ,
             T(use_uniform ? compute_uniform_pdf(X) :-Inf),
             T(use_uniform ? log(1/(size(X, 1)*length(components)+use_uniform)) : -Inf),
             X, zeros(UInt8, size(X, 2)), λ, fix_ν, T(-Inf), zeros(T, size(X, 2)), zeros(T, size(X, 2)),
             zeros(T, size(X, 2)), zero(X))::TMixture{T,typeof(X)}
end

function initial_centers(k::Int, X::StridedMatrix)
    if k <= 0
        throw(ArgumentError("cannot create a mixture model with no components"))
    elseif k == 1
        mean(X, 2)
    else
        kmeans(X, k).centers
    end
end

function TMixture{T}(k::Int, X::StridedMatrix{T}; λ::T=convert(T, 0), centers::Matrix{T}=initial_centers(k, X), use_uniform::Bool=false, ν::T=T(50), fix_ν::Bool=false)
    k < typemax(UInt8) || throw(ArgumentError("a maximum of $(typemax(UInt8)) components are supported"))
    d = size(X, 1)
    sd = diagm(vec(std(X, 2)))/4
    components = MvTDist[MvTDist(ν, centers[:, i], PDMat(copy(sd), cholfact(sd))) for i = 1:k]
    TMixture(components, X; λ=λ, use_uniform=use_uniform, fix_ν=fix_ν)::TMixture{T,typeof(X)}
end

function merge_partial{T,C}(m::TMixture{T,C}, mpartial::TMixture{T,C}, original_components::Vector{Int}, partial_indices::Vector{Int})
    keepcomps = deleteat!(collect(1:length(m.components)), original_components)
    components = [deepcopy(m.components[i])::eltype(m.components) for i in keepcomps]
    append!(components, mpartial.components)
    merged = TMixture(components, m.X; λ=m.λ, use_uniform=m.uniform_logpdf != -Inf, fix_ν=m.fix_ν)
    merged.logπ = m.logπ[keepcomps]
    append!(merged.logπ, mpartial.logπ - StatsFuns.logsumexp(m.logπ[original_components]))

    old_assignments = m.assignments
    assignments = merged.assignments = similar(m.assignments)
    keepcomps_map = zeros(Int, length(m.components)+1)
    keepcomps_map[keepcomps+1] = 1:length(keepcomps)
    for i = 1:length(assignments)
        assignments[i] = keepcomps_map[old_assignments[i]+1]
    end
    partial_assignments = mpartial.assignments
    for i = 1:length(partial_indices)
        assignments[partial_indices[i]] = partial_assignments[i]
    end

    merged
end

function emstep!(m::TMixture)
    components = m.components
    logπ = m.logπ
    X = m.X
    assignments = m.assignments
    assignedu = m.assignedu
    maxtmp = m.maxtmp
    mahaltmp = m.mahaltmp
    Xtmp = m.Xtmp
    λ = m.λ
    fix_ν = m.fix_ν

    has_uniform = m.uniform_logpdf != -Inf

    # E step
    fill!(maxtmp, m.uniform_logpdf + m.uniform_logπ)
    fill!(assignments, 0)
    for k = 1:length(components)
        c = components[k]
        # Compute squared Mahalanobis distance
        broadcast!(-, Xtmp, X, c.μ)
        invquad!(mahaltmp, c.Σ, Xtmp)

        shdfhdim, v = Distributions.mvtdist_consts(c)
        lπ = logπ[k]
        for i = 1:size(X, 2)
            # Compute logpdf and see if it is better than for other components
            logpdf = v - shdfhdim * log1p(mahaltmp[i] / c.df)
            l = logpdf + lπ
            if l > maxtmp[i]
                maxtmp[i] = l
                assignments[i] = k
                assignedu[i] = (size(X, 1) + c.df)/(mahaltmp[i] + c.df)
            end
        end
    end
    m.ll = sum(maxtmp)
    isempty(components) && return

    # M step
    nuniform = 0
    if has_uniform
        for i = 1:length(assignments)
            nuniform += assignments[i] == 0
        end
        m.uniform_logπ = log((1+nuniform)/(size(X, 2)+1))
    end
    y = -digamma((size(X, 1) + components[1].df)/2)*(size(X, 2) - nuniform)
    lognsamples = log(size(X, 2)+has_uniform)
    for k = 1:length(components)
        c = components[k]
        incomponent = find(assignments .== k)

        # Update μ and Σ
        # We have actually already computed Xtmp before, so if this is
        # a computationally intensive step, we should just use more memory...
        Xcomp = X[:, incomponent]
        wt = assignedu[incomponent]
        weighted_mean_cov!(c.μ, c.Σ.mat, Xcomp, Xcomp, wt)
        regularize_and_factorize!(c, λ)

        # Update mixing proportions
        logπ[k] = log(length(incomponent)) - lognsamples

        # Update ν
        if !fix_ν
            for i = 1:length(wt)
                y -= log(2*(wt[i]/(c.dim + c.df))) - wt[i]
            end
        end
    end

    # Update degrees of freedom
    if !fix_ν
        dgconst = digamma((size(X, 1) + components[1].df)/2)
        y /= (size(X, 2) - nuniform)
        ν = fzero(ν->-y+log(ν/2)+1-digamma(ν/2), 1.0, 100000.0)
        yplogym1 = 1/(y+log(y)-1)
        # ν = min(2*yplogym1+0.0416*(1+erf(0.6594*log(2.1971*yplogym1))), 100.0)
        for k = 1:length(components)
            components[k] = MvTDist(ν, components[k].μ, components[k].Σ)
        end
    end

    m
end

function nparameters(d::Union{FullNormal,MvTDist})
    n = length(d.μ)
    div(n*(n-1), 2)+n
end

function score(m::TMixture)
    k = length(m.components)
    N = nparameters(m.components[1])
    n = size(m.X, 2)

    penalty = N/2*(k*log(n/12)+sum(m.logπ)) + k/2*log(n/12) + k*(N+1)/2
    @assert penalty > 0
    m.ll - penalty
end

function computeγ{T}(m::TMixture{T})
    components = m.components
    logπ = m.logπ
    X = m.X
    sumtmp = m.mahaltmp
    maxtmp = m.maxtmp

    γ = zeros(T, size(X, 2), length(components))

    fill!(maxtmp, -Inf)
    for k = 1:length(components)
        γ_k = sub(γ, :, k)
        # Compute log(π_k*p(x_n|θ_k))
        logpdf!(γ_k, components[k], X)
        lπ = logπ[k]
        for i = 1:size(X, 2)
            γ_k[i] += lπ
            maxtmp[i] = max(maxtmp[i], γ_k[i])
        end
    end

    # Now compute softmax over rows of γ to complete Eq. 9.23
    fill!(sumtmp, 0)
    for k = 1:length(components), i = 1:size(X, 2)
        sumtmp[i] += γ[i, k] = exp(γ[i, k] - maxtmp[i])
    end

    # Finish the softmax
    broadcast!(*, γ, γ, broadcast!(inv, sumtmp, sumtmp))

    γ
end

function logsumexp(a, b)
    u = max(a, b)
    v = min(a, b)
    u + log1p(exp(v-u))
end

"""
    merge!(out, c1, logπ1, c2, logπ2)

Merge components c1 and c2 into out and return new log mixing coefficient
"""
function merge!{C<:Union{MvTDist,FullNormal}}(out::C, c1::C, logπ1::Real, c2::C, logπ2::Real)
    π1 = exp(logπ1)
    π2 = exp(logπ2)

    # Merge μ
    for m = 1:length(c1.μ)
        out.μ[m] = (π1*c1.μ[m] + π2*c2.μ[m])/(π1+π2)
    end

    # Merge Σ
    for n = 1:size(c1.Σ.mat, 2), m = 1:size(c2.Σ.mat, 1)
        out.Σ.chol.factors[m, n] = out.Σ.mat[m, n] = (π1*c1.Σ.mat[m, n] + π2*c2.Σ.mat[m, n])/(π1+π2)
    end
    LinAlg.chol!(out.Σ.chol.factors, Val{:U})

    # New log mixing coefficient
    logsumexp(logπ1, logπ2)
end

function merge_priority(m::TMixture)
    γ = computeγ(m)
    c = Base.corzm(γ)
    ncandidates = div(size(γ, 2)*(size(γ, 2) - 1), 2)
    candidates = Array(Tuple{Int,Int}, ncandidates)
    priority = zeros(ncandidates)
    n = 1
    for i = 1:size(γ, 2), j = i+1:size(γ, 2)
        candidates[n] = (i, j)
        priority[n] = c[i, j]
        n += 1
    end
    (priority, candidates)
end

function merge{T}(m::TMixture{T}, i::Int, j::Int; tol::T=1e-3, maxiter::Int=200, verbose::Symbol=:none)
    ncomponents = length(m.components)
    mnew = deepcopy(m)
    deleteat!(mnew.logπ, j)
    deleteat!(mnew.components, j)
    assignmap = zeros(Int, ncomponents+1)
    assignmap[[2:j; j+2:ncomponents+1]] = 1:ncomponents-1
    assignmap[j+1] = i
    mnew.assignments = [assignmap[x+1] for x in m.assignments]
    mnew.logπ[i] = merge!(mnew.components[i], m.components[i], m.logπ[i], m.components[j], m.logπ[j])
    em!(mnew; tol=tol, maxiter=maxiter, verbose=verbose)
end

function merge{T}(m::TMixture{T}; tol::T=1e-3, maxiter::Int=200, verbose::Symbol=:none)
    ncomponents = length(m.components)
    ncomponents == 1 && return m

    # Find merge candidates
    (priority, candidates) = merge_priority(m)

    # Try some of the best candidates
    maxcandidates = min(ceil(Int, ncomponents*sqrt(ncomponents)/2), length(priority))
    p = sortperm(priority)
    origscore = score(m)
    for icandidate = length(p):-1:length(p)-maxcandidates+1
        i, j = candidates[p[icandidate]]

        # Do the merge
        verbose_op(verbose) && println("Attempting to merge components $i and $j...")
        mnew = merge(m, i, j; tol=tol, maxiter=maxiter, verbose=verbose)

        newscore = score(mnew)
        if newscore > origscore
            verbose_op(verbose) && println("    score increased ($origscore -> $newscore); accepted")
            return mnew
        else
            verbose_op(verbose) && println("    score decreased ($origscore -> $newscore); rejected")
        end
    end
    return m
end

function split_priority{T}(m::TMixture{T})
    γ = computeγ(m)
    components = m.components
    priority = zeros(length(components))
    fk = zeros(T, size(m.X, 2))
    logpk = zeros(T, size(m.X, 2))
    for k = 1:length(components)
        γ_k = sub(γ, :, k)
        scale!(fk, γ_k, 1/sum(γ_k))
        logpdf!(logpk, components[k], m.X)
        Jsplit = zero(T)
        for i = 1:size(m.X, 2)
            Jsplit += fk[i]*(log(fk[i])-logpk[i])
        end
        priority[k] = Jsplit
    end
    priority
end

function split{T}(m::TMixture{T}, i::Int; tol::T=1e-3, maxiter::Int=200, verbose::Symbol=:none)
    indpartial = find(m.assignments .== i)
    if length(indpartial) < size(m.X, 1)*2
        return m
    end

    mpartial = TMixture(2, m.X[:, indpartial]; λ=m.λ, use_uniform=false, ν=m.components[i].df, fix_ν=true)
    em!(mpartial)
    if length(mpartial.components) == 1
        m
    else
        # Merge the partial EM results into a new model
        mnew = merge_partial(m, mpartial, [i], indpartial)
        em!(mnew; tol=tol, maxiter=maxiter, verbose=verbose)
    end
end

function split{T}(m::TMixture{T}; tol::T=1e-3, maxiter::Int=200, verbose::Symbol=:none)
    # Find split candidates
    priority = split_priority(m)

    # Try each candidate
    p = sortperm(priority)
    origscore = score(m)
    for icandidate = length(p):-1:1
        candidate = p[icandidate]

        # Do the split
        verbose_op(verbose) && println("Attempting to split component $candidate...")
        mnew = split(m, candidate; tol=tol, maxiter=maxiter, verbose=verbose)
        if mnew == m
            verbose_op(verbose) && println("    new component degenerate; rejected")
        else
            newscore = score(mnew)
            if newscore > origscore
                verbose_op(verbose) && println("    score increased ($origscore -> $newscore); accepted")
                return mnew
            else
                verbose_op(verbose) && println("    score decreased ($origscore -> $newscore); rejected")
            end
        end
    end
    return m
end

function drop_degenerate!(m::TMixture)
    assigned = zeros(Int, length(m.components))
    assignments = m.assignments
    for i = 1:length(assignments)
        if assignments[i] != 0
            assigned[assignments[i]] += 1
        end
    end

    assigned_ind = find(assigned .>= size(m.X, 1))
    length(assigned_ind) == length(m.components) && return
    assignmap = zeros(UInt8, length(m.components))
    assignmap[assigned_ind] = 1:length(assigned_ind)
    m.components = m.components[assigned_ind]
    m.logπ = m.logπ[assigned_ind]
    for i = 1:length(assignments)
        if assignments[i] != 0
            assignments[i] = assignmap[assignments[i]]
        end
    end
    m
end

function em!{T}(m::TMixture{T}; tol::T=1e-2, maxiter::Int=200, verbose::Symbol=:none)
    niter = 0
    m.ll = -Inf
    old_assignments = copy(m.assignments)
    while niter < maxiter
        oldll = m.ll
        old_assignments, m.assignments = m.assignments, old_assignments
        emstep!(m)
        drop_degenerate!(m)
        verbose_iter(verbose) && println("    $niter $(m.ll)")
        !isfinite(m.ll) && error("log likelihood is not finite")
        # m.assignments == old_assignments && return m
        m.ll - oldll < tol && return m
        niter += 1
    end
    warn("EM did not converge")
    m
end

function fsmem!{T}(m::TMixture{T}; tol::T=1e-2, maxiter::Int=200, verbose::Symbol=:none)
    verbose_op(verbose) && println("Performing initial EM...")
    em!(m, tol=tol, maxiter=maxiter, verbose=verbose)
    while true
        mouter = m
        while true
            minner = m
            m = merge(m, tol=tol, maxiter=maxiter, verbose=verbose)
            m === minner && break
        end
        while true
            minner = m
            m = split(m, tol=tol, maxiter=maxiter, verbose=verbose)
            m === minner && break
        end
        m === mouter && break
    end
    m
end

assignments{T}(m::TMixture{T}) = m.assignments

end # module


text:
#! /usr/bin/julia

# Rosetta Code, Number names

const stext = ["one", "two", "three", "four", "five",
               "six", "seven", "eight", "nine"]
const teentext = ["eleven", "twelve", "thirteen", "fourteen",
                  "fifteen", "sixteen", "seventeen",
                  "eighteen", "nineteen"]
const tenstext = ["ten", "twenty", "thirty", "forty", "fifty",
                  "sixty", "seventy", "eighty", "ninety"]
const ordstext = ["million", "billion", "trillion",
                  "quadrillion", "quintillion", "sextillion",
                  "septillion", "octillion", "nonillion",
                  "decillion", "undecillion", "duodecillion",
                  "tredecillion", "quattuordecillion", "quindecillion",
                  "sexdecillion", "septendecillion", "octodecillion",
                  "novemdecillion", "vigintillion"]

function normalize_digits!{T<:Integer}(a::Array{T,1})
    while  0 < length(a) && a[end] == 0
        pop!(a)
    end
    return length(a)
end

function digits2text!{T<:Integer}(d::Array{T,1}, use_short_scale=true)
    ndig = normalize_digits!(d)
    0 < ndig || return ""
    if ndig < 7
        s = ""
        if 3 < ndig
            t = digits2text!(d[1:3])
            s = digits2text!(d[4:end])*" thousand"
            0 < length(t) || return s
            if contains(t, "and")
                return s*" "*t
            else
                return s*" and "*t
            end
        end
        if ndig == 3
            s *= stext[pop!(d)]*" hundred"
            ndig = normalize_digits!(d)
            0 < ndig || return s
            s *= " and "
        end
        1 < ndig || return s*stext[pop!(d)]
        j, i = d
        j != 0 || return s*tenstext[i]
        i != 1 || return s*teentext[j]
        return s*tenstext[i]*"-"*stext[j]
    end
    s = digits2text!(d[1:6])
    d = d[7:end]
    dgrp = use_short_scale ? 3 : 6
    ord = 0
    while(dgrp < length(d))
        ord += 1
        t = digits2text!(d[1:dgrp])
        d = d[(dgrp+1):end]
        0 < length(t) || continue
        t = t*" "*ordstext[ord]
        if length(s) == 0
            s = t
        else
            s = t*" "*s
        end
    end
    ord += 1
    t = digits2text!(d)*" "*ordstext[ord]
    0 < length(s) || return t
    return t*" "*s
end

function num2text{T<:Integer}(n::T, use_short_scale=true)
    -1 < n || return "minus "*num2text(-n, use_short_scale)
    0 < n || return "zero"
    toobig = use_short_scale ? big(10)^66 : big(10)^126
    n < toobig || return "too big to say"
    return digits2text!(digits(n, 10), use_short_scale)
end

function wrap{T<:String,U<:Integer}(str::T, width::U=70, indent::U=4)
    sdent = " "^(indent-1)
    s = sdent
    llen = indent
    for t in split(str, " ")
        s *= " "*t
        llen += length(t) + 1
        width < llen || continue
        s *= "\n"*sdent
        llen = indent
    end
    return rstrip(s)
end

println("Testing num2text")
println()
println("Some easy ones to start with")
println()
tests = [-1:21, 100, 101, 10000, 10001, 1000000, 1010101]
for i in tests
    println(@sprintf("%8d is %s", i, num2text(i)))
end

println()
println("Some larger numbers")
println()
println("The largest signed literal integer (short-scale)")
i = typemax(1)
println("    ", i, " is")
println(wrap(num2text(i)))
println()
println("The largest signed literal integer (long-scale)")
println("    ", i, " is")
println(wrap(num2text(i, false)))
println()
println("The largest unsigned integer (short-scale)")
i = typemax(Uint128)
println("    ", i, " is")
println(wrap(num2text(i)))
println()
println("50! (short-scale)")
i = factorial(big(50))
println("    ", i, " is")
println(wrap(num2text(i)))
println()
println("51! (short-scale)")
i = factorial(big(51))
println("    ", i, " is")
println(wrap(num2text(i)))
println()
println("51! (long-scale)")
println("    ", i, " is")
println(wrap(num2text(i, false)))

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
type AtomicGrid{T,A<:AbstractMatrix,B<:AbstractMatrix,C<:AbstractGrid1d} <: AbstractGrid3d{T,2}
    c::A  # coords
    w::B  # weights
    r::C  # radial grid
    Ω::AngularGrid{T}  # angular grid
    R::Point3{T}  # origin of spherical coord system
    function AtomicGrid(c::AbstractMatrix{Point3{T}},w::AbstractMatrix{T},
            r::AbstractGrid1d{T,1},Ω::AngularGrid{T},R::Point3{T})
        @assert size(c) == size(w)
        new(c,w,r,Ω,R)
    end
end

AtomicGrid{T}(c::AbstractMatrix{Point3{T}},w::AbstractMatrix{T},
        r::AbstractGrid1d{T,1},Ω::AngularGrid{T},R::Point3{T}) =
    AtomicGrid{eltype(w),typeof(c),typeof(w),typeof(r)}(c,w,r,Ω,R)

function atomicgrid{T}(r::AbstractGrid1d{T,1}, Ω::AngularGrid{T}; R::Point3{Float64}=(0.,0.,0.))
    shp = (length(Ω),length(r))
    c = Array{Point3{Float64}}(shp)
    w = Array{Float64}(shp)
    for j in 1:length(r)
        for i in 1:length(Ω)
            c[i,j] = r[j]*Ω[i] + R
            w[i,j] = r[j]^2*r.w[j]*Ω.w[i]
        end
    end
    AtomicGrid(c,w,r,Ω,R)
end

# integrate_angular{T}(g::AtomicGrid{T}, f::AbstractMatrix{T}) = At_mul_B(f,g.Ω.w)

# integrate_angular{T}(g::AtomicGrid{T}, f::AbstractVector{T}) = reshape(f,size(g.w))'*g.Ω.w

# function integrate_angular{T}(g::AtomicGrid{T}, f::AbstractArray{T})
#     shp = (size(g.w,1),prod(size(f)[2:end]))
#     reshape(f,shp)'*g.Ω.w
# end

# type AtomicGrid{T,N,Tc<:AbstractArray,Tw<:AbstractArray,Tr<:AbstractGrid1d} <: AbstractGrid3d{T,N}
#     c::Tc  # coords
#     w::Tw  # weights
#     r::Tc  # radial grid
#     Ω::AngularGrid{T}  # angular grid
#     R::Point3{T}  # origin of spherical coord system
#     function AtomicGrid(c::AbstractArray{Point3{T},N},
#                         w::AbstractArray{T,N},
#                         r::AbstractGrid1d{T,1},
#                         Ω::AngularGrid{T},
#                         R::Point3{T})
#         @assert size(c) == size(w)
#         new(c,w,r,Ω,R)
#     end
# end

# function AtomicGrid{T,N}(c::AbstractArray{Point3{T},N},
#                          w::AbstractArray{T,N},
#                          r::AbstractGrid1d{T,1},
#                          Ω::AngularGrid{T},
#                          R::Point3{T})
#     AtomicGrid{eltype(w),N,typeof(c),typeof(w),typeof(r)}(c,w,r,Ω,R)
# end

max_stars_repo_name: johncsnyder/Element


text:
module M4

export SeriesValue, Series, SeriesDict, SeriesInfo, SeriesInfoDict
export read_series, read_info
export ForecastScore, evaluate, calc_stats
export NaiveForecaster
export plot_pieces

include("data_loader.jl")
include("evaluation.jl")
include("prep.jl")
include("naive.jl")
include("plotting.jl")

end


text:
@testset "Scanners" begin
  include("DummyScannerTest.jl")
  include("FlexibleScannerTest.jl")
  include("SimpleSimulatedScannerTest.jl")
end


text: const NullRules = OrderedRuleSet("Null", ("Unclassified", r -> true) )

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
# # Two-dimensional cellular automata
# Agents.jl provides a module (CA2D) to create and plot 2D cellular automata.

using Agents
using AgentsPlots
using Agents.CA2D

# ## 1. Define the rule
# Rules of Conway's game of life: DSR (Death, Survival, Reproduction).
# Cells die if the number of their living neighbors are <D,
# survive if the number of their living neighbors are <=S,
# come to life if their living neighbors are as many as R.
rules = (2, 3, 3)

# ## 2. Build the model
# "CA2D.build_model" creates a model where all cells are by default off ("0")
model = CA2D.build_model(rules = rules, dims = (100, 100), Moore = true)

# Let's make some random cells on
for i in 1:nv(model)
    if rand() < 0.1
        model.agents[i].status = "1"
    end
end

# ## 3. Animate the model

# Animations of 2D Cellular Automata are done through the following syntax:
# (`plot_CA2Dgif` is a function from `AgentsPlots` that creates the animation)

runs = 10
anim = CA2D.ca_run(model, runs, plot_CA2Dgif);

# We can now save the animation to a gif.

AgentsPlots.gif(anim, "game_of_life.gif")


text:
using GraphNeuralNetworks
using Statistics

"""
Super, super demo example

network will be this
nodeFeature -> number of features on node
innerSize -> how large the inner network layer will be
actionCount -> how many actions plus 1 (network importance)

GNNChain(GCNConv(nodeFeature => innerSize),
                        BatchNorm(innerSize),     # Apply batch normalization on node features (nodes dimension is batch dimension)
                        x -> relu.(x),
                        GCNConv(64 => innerSize, relu),
                        GlobalPool(mean),  # aggregate node-wise features into graph-wise features
                        Dense(innerSize, actionCount+1)) |> device;
"""
@kwdef struct SimpleGraphNetHP
  innerSize :: Int
end

"""
    SimpleGNN <: FluxGNN

    Something simple
"""

mutable struct SimpleGNN <: FluxGNN
  gspec
  hyper
  common
  vhead
  phead
end

function SimpleGNN(gspec::AbstractGameSpec, hyper::SimpleGraphNetHP)
    innerSize = hyper.innerSize
    nodeFeature = GI.state_dim(gspec)
    actionCount = GI.num_actions(gspec)
    common = GNNChain(GCNConv(nodeFeature => innerSize),
                        BatchNorm(innerSize),     # Apply batch normalization on node features (nodes dimension is batch dimension)
                        x -> relu.(x),
                        GCNConv(innerSize => innerSize, relu)
      )
      modelP = GNNChain(Dense(innerSize, 1),softmax)
      modelV = GNNChain( GlobalPool(mean),  # aggregate node-wise features into graph-wise features
                              Dense(innerSize, 1),
                              softmax);

    return SimpleGNN(gspec, hyper, common, modelV, modelP)
end

Network.HyperParams(::Type{<:SimpleGNN}) = SimpleGraphNetHP

function Base.copy(nn::SimpleGNN)
  return SimpleGNN(
    nn.gspec,
    nn.hyper,
    deepcopy(nn.common),
    deepcopy(nn.vhead),
    deepcopy(nn.phead)
  )
end

# function Network.forward(nn::SimpleGNN, state)
#   c = nn.common.(state)
#   applyV(graph) = nn.vhead(graph, graph.ndata.x)
#   resultv = applyV.(c)
#   v = first.(resultv)
#   applyP(graph) = nn.phead(graph)
#   resultp = applyP.(c)
#   @show @__LINE__
#   @show length(resultp[1].ndata.x)
#   @show state[1].num_nodes
#   p = [resultp[ind].ndata.x[indDepth] for indDepth in 1:state[1].num_nodes, ind in 1:length(state) ]
#   return (p, v)
# end

function Network.forward(nn::SimpleGNN, state)
  c = nn.common(state) # no broadcasting
  v = nn.vhead(c, c.ndata.x)   # size: num_features x num_graphs
  p = nn.phead(c, c.ndata.x)   # size: num_features x num_tot_nodes
  return p, v
end


text:
include("mesh2d.jl")
include("nektar.jl")

<<<<<<< HEAD

=======
module XLAPACK
include("xlapack.jl")
end
>>>>>>> ef030efb9656085dd4f7e05f695cbdb6bdca1eb0

module Makhno

# package code goes here
include("xlapack.jl")
include("element2d.jl")
include("h1space.jl")
include("numbering.jl")
include("matrix.jl")

end # module


text:
# set the working directory to the directory in which this file is contained
cd(@__DIR__)
println("Installing and compiling dependencies...")
println("This will take several minutes...")
# load the package manager
using Pkg
Pkg.activate("..")
# initialize project environment
Pkg.instantiate()
# precompile the packages
Pkg.precompile()


text:
module PeakshapeFunctions
	import PyPlot
	import SharedArrays
	import ..InterpolationFunctions
	import Statistics

	export findPeakIndices, calculatePeakshapes, getLocalPeakshape
	"searches an average spectrum for peaks signifficantly higher than the baseline noise and returns their interpolated indices in the average spectrum"
	  function findPeakIndices(massAxis, avgSpectrum, baseline, baselineNoise; noiseThreshold = 80, oddEven="both", signalLimit = 1)
	    totalMax = maximum(avgSpectrum)
	    peakIndices = Array{Float64,1}()
	    for j=searchsortedfirst(massAxis,10):length(massAxis)-1
	      if ( avgSpectrum[j] > baseline[j] + baselineNoise[j] * noiseThreshold * 400/(massAxis[j]^0.6)) && (avgSpectrum[j] > avgSpectrum[j-1]) && (avgSpectrum[j] > avgSpectrum[j+1]) && (avgSpectrum[j] < totalMax*signalLimit)
		  if (oddEven == "both") | ((oddEven == "even") & iseven(Int(round(massAxis[j],digits=0)))) | ((oddEven == "odd") & !iseven(Int(round(massAxis[j],digits=0))))
		      ############ put interpolated exact peak position in peak list ###
		      push!(peakIndices,InterpolationFunctions.interpolatedMax(j,avgSpectrum))
		  end
	      end
	    end
	    return peakIndices
	  end

	  function calculatePeakshapes(massAxis, baselineCorrectedAvgSpec, peakIndices; nbrMassRegions = 10, peakWindowWidth = peakWindowWidth, quantileValue = 0.05, regionStretch=1)
	    peakMasses = InterpolationFunctions.interpolate(peakIndices, massAxis)
	    peakValues = InterpolationFunctions.interpolate(peakIndices,baselineCorrectedAvgSpec)

	    peakShapesY = Array{Float64}(undef,2*peakWindowWidth+1, nbrMassRegions)
	    peakShapesCenterMass = Array{Float64}(undef,nbrMassRegions)

	    PyPlot.figure()
	    for massRegion = 1:nbrMassRegions
	      ax=PyPlot.subplot(ceil(nbrMassRegions/4),4,massRegion)

	      #peakshapeRangeStart = (massRegion-1) *peakMasses[end] / nbrMassRegions
	      #peakshapeRangeEnd = (massRegion) *peakMasses[end] / nbrMassRegions
	      peakshapeRangeStart = (massRegion-1)^regionStretch *peakMasses[end] / nbrMassRegions^regionStretch
	      peakshapeRangeEnd = (massRegion)^regionStretch *peakMasses[end] / nbrMassRegions^regionStretch

	      peakShapesCenterMass[massRegion] = (peakshapeRangeEnd + peakshapeRangeStart)/2
	      sel = (peakMasses.>peakshapeRangeStart) .& (peakMasses.<peakshapeRangeEnd) .& (peakIndices.> peakWindowWidth) .& (peakIndices .< length(baselineCorrectedAvgSpec)-peakWindowWidth)

	      peakIndicesInRegion = peakIndices[sel]
	      peakMassesInRegion = peakMasses[sel]
	      peakValuesInRegion = peakValues[sel]

	      peakWood = SharedArrays.SharedArray{Float64}(2*peakWindowWidth+1,length(peakMassesInRegion))
	      fill(peakWood, 0)
	      for i=1:length(peakMassesInRegion)
		  ############# interpolate exact peak position ###################
		lowIdx = Int64(floor(peakIndicesInRegion[i]))
		highIdx = lowIdx+1
		fract = peakIndicesInRegion[i]-lowIdx

		peakWood[:,i] = (
		(1-fract)*view(baselineCorrectedAvgSpec, lowIdx-peakWindowWidth : lowIdx+peakWindowWidth) +
		fract*view(baselineCorrectedAvgSpec, highIdx-peakWindowWidth : highIdx+peakWindowWidth))/peakValuesInRegion[i]
	      end

	      peakshapeY = Array{Float64}(undef,2*peakWindowWidth+1)
	      for i=1:2*peakWindowWidth+1
		#peakshapeY[i] = minimum(peakWood[i,:])
		values = peakWood[i,:]
		#upperThreshold = quantile(values,quantileValue*peakshapeRangeEnd/100)
		###### quantile based #############
		#upperThreshold = quantile(values,quantileValue*peakshapeRangeEnd/100)
		#peakshapeY[i] = mean(values[values.<upperThreshold])
		#lowest n perent based ############
		values = sort(values[values.>0])
		if length(values) > 0
		  peakshapeY[i] = Statistics.mean(values[1 : Int64(ceil(length(values)*quantileValue))])
		else
		  peakshapeY[i] = 0
		  if i == peakWindowWidth + 1
		    peakshapeY[i] = 1
		  end
		end
	      end
	      peakshapeY[peakshapeY.<0] .= 0
	      peakshapeY[peakWindowWidth+1] = 1
	      peakShapesY[:,massRegion] = peakshapeY / sum(peakshapeY)
	      if (size(peakWood,2) > 0)
		  PyPlot.semilogy(peakWood)
	      end
	      PyPlot.semilogy(peakshapeY, "o-", label="$peakshapeRangeStart - $peakshapeRangeEnd")
	      ax["set_title"]("$(round(peakshapeRangeStart,digits=1)) - $(round(peakshapeRangeEnd,digits=1))")
	    end
	    PyPlot.suptitle("Mass dependent Peakshapes")
	    return peakShapesCenterMass, peakShapesY
	  end

	  function getLocalPeakshape(mass, peakShapesCenterMass, peakShapesY)
	    peakShapeCenterMassIndexHigh = searchsortedfirst(peakShapesCenterMass, mass)
	    if (peakShapeCenterMassIndexHigh == 1)
	      localPeakshape = peakShapesY[:,1]
	    elseif (peakShapeCenterMassIndexHigh > length(peakShapesCenterMass))
	      localPeakshape = peakShapesY[:,end]
	    else
	      peakShapeCenterMassIndexLow = peakShapeCenterMassIndexHigh - 1
	      peakShapesCenterMassDistance = peakShapesCenterMass[peakShapeCenterMassIndexHigh] - peakShapesCenterMass[peakShapeCenterMassIndexLow]
	      contribLow =  (peakShapesCenterMass[peakShapeCenterMassIndexHigh] - mass) / peakShapesCenterMassDistance
	      contribHigh = 1 - contribLow
	      localPeakshape = peakShapesY[:,peakShapeCenterMassIndexLow]*contribLow + peakShapesY[:,peakShapeCenterMassIndexHigh]*contribHigh
	    end
	    return localPeakshape
	  end

end


text:
import Base.Iterators

abstract type Seq{T} end

struct EmptySeq{T} <: Seq{T} end

struct ConsSeq{T} <: Seq{T}
    genfun::Function
end

function Seq(::Type{T}, state, next::Function) where {T}
    ConsSeq{T}(() -> (state, next))
end

function Seq(::Type{T}, iter) where {T}
    ConsSeq{T}( () -> (start(iter),
                       state -> done(iter, state) ?
                       (Nullable{T}(), state) :
                       begin
                           v, newstate = next(iter, state)
                           Nullable{T}(v), newstate
                       end))
end

function Seq(range::Range{T}) where {T}
    Seq(T, range)
end

tryread(s::EmptySeq{T}) where {T} = Nullable{T}(), EmptySeq{T}()

function tryread(xs::ConsSeq{T}) where {T}
    state, next = xs.genfun()
    v, newstate = next(state)
    if isnull(v)
        Nullable{T}(), EmptySeq{T}()
    else
        v, ConsSeq{T}(() -> (newstate, next))
    end
end

Base.map(f::Function, xs::EmptySeq{T}, ::Type{S}) where {T, S} = EmptySeq{S}()

Base.map(f::Function, xs::EmptySeq{T}) where {T} = EmptySeq{T}()

function Base.map(f::Function, xs::ConsSeq{T}, ::Type{S}) where {T, S}
    ConsSeq{S}(() ->
                begin
                    state, next = xs.genfun()
                    state, s -> begin
                                    v, newstate = next(s)
                                    map(f, v), newstate
                                end
                end
              )
end

function Base.map(f::Function, xs::ConsSeq{T}) where {T}
    map(f, xs, T)
end

fold(f::Function, init, xs::EmptySeq{T}) where {T} = init

function fold(f::Function, init, xs::ConsSeq{T}) where {T}
    eof = false
    acc = init
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if isnull(v)
            eof = true
        else
            acc = f(acc, get(v))
        end
    end
    acc
end

Iterators.take(xs::EmptySeq{T}, n::Integer) where {T} = EmptySeq{T}()

function Iterators.take(xs::ConsSeq{T}, n::Integer) where {T}
    ConsSeq{T}(() ->
                    begin
                        state, next = xs.genfun()
                        (state, n), s -> begin
                                             state = s[1]
                                             n = s[2]
                                             if n >= 1
                                                v, newstate = next(state)
                                                v, (newstate, n - 1)
                                             else
                                                 Nullable{T}(), state
                                             end
                                    end
                    end
             )
end

function chunkbysize(xs::EmptySeq{T}, n::Integer) where {T}
    EmptySeq{Vector{T}}()
end

function chunkbysize(xs::ConsSeq{T}, n::Integer) where {T}
    return ConsSeq{Vector{T}}(() ->
        begin
            state, next = xs.genfun()
            state, s -> begin
                            chunk = Vector{T}(n)
                            newstate = s
                            for k in 1:n
                                v, newstate = next(newstate)
                                if !isnull(v)
                                    chunk[k] = get(v)
                                else
                                    resize!(chunk, k - 1)
                                    break
                                end
                            end
                            if length(chunk) == 0
                                Nullable{Vector{T}}(), newstate
                            else
                                Nullable{Vector{T}}(chunk), newstate
                            end
                        end
        end)
end

function zip2(xs::Seq{T}, ys::Seq{S}) where {T} where {S}
    EmptySeq{Tuple{T, S}}()
end

function zip2(xs::ConsSeq{T}, ys::ConsSeq{S}) where {T} where {S}
    return ConsSeq{Tuple{T, S}}(() ->
            begin
                statex, nextx = xs.genfun()
                statey, nexty = ys.genfun()
                (statex, statey), state ->
                    begin
                        _statex, _statey = state
                        x, newstatex = nextx(_statex)
                        y, newstatey = nexty(_statey)
                        if isnull(x) || isnull(y)
                            Nullable{Tuple{T, S}}(), (newstatex, newstatey)
                        else
                            Nullable{Tuple{T, S}}((get(x), get(y))), (newstatex, newstatey)
                        end
                    end
            end
            )
end

function zip3(xs::Seq{T}, ys::Seq{S}, zs::Seq{U}) where {T} where {S} where {U}
    EmptySeq{Tuple{T, S, U}}()
end

function zip3(xs::ConsSeq{T}, ys::ConsSeq{S}, zs::ConsSeq{U}) where {T} where {S} where {U}
    return ConsSeq{Tuple{T, S, U}}(() ->
            begin
                statex, nextx = xs.genfun()
                statey, nexty = ys.genfun()
                statez, nextz = zs.genfun()
                (statex, statey, statez), state ->
                    begin
                        _statex, _statey, _statez = state
                        x, newstatex = nextx(_statex)
                        y, newstatey = nexty(_statey)
                        z, newstatez = nextz(_statez)
                        if isnull(x) || isnull(y) || isnull(z)
                            Nullable{Tuple{T, S, U}}(), (newstatex, newstatey, newstatez)
                        else
                            Nullable{Tuple{T, S, U}}((get(x), get(y), get(z))), (newstatex, newstatey, newstatez)
                        end
                    end
            end
            )
end

function zip4(xs::Seq{T}, ys::Seq{S}, zs::Seq{U}, vs::Seq{V}) where {T} where {S} where {U} where {V}
    EmptySeq{Tuple{T, S, U, V}}()
end

function zip4(xs::ConsSeq{T}, ys::ConsSeq{S}, zs::ConsSeq{U}, vs::Seq{V}) where {T} where {S} where {U} where {V}
    return ConsSeq{Tuple{T, S, U, V}}(() ->
            begin
                statex, nextx = xs.genfun()
                statey, nexty = ys.genfun()
                statez, nextz = zs.genfun()
                statev, nextv = vs.genfun()
                (statex, statey, statez, statev), state ->
                    begin
                        _statex, _statey, _statez, _statev = state
                        x, newstatex = nextx(_statex)
                        y, newstatey = nexty(_statey)
                        z, newstatez = nextz(_statez)
                        v, newstatev = nextv(_statev)
                        if isnull(x) || isnull(y) || isnull(z) || isnull(v)
                            Nullable{Tuple{T, S, U, V}}(), (newstatex, newstatey, newstatez, newstatev)
                        else
                            Nullable{Tuple{T, S, U, V}}((get(x), get(y), get(z), get(v))), (newstatex, newstatey, newstatez, newstatev)
                        end
                    end
            end
            )
end

function zip5(x1s::Seq{T1}, x2s::Seq{T2}, x3s::Seq{T3}, x4s::Seq{T4}, x5s::Seq{T5}) where {T1} where {T2} where {T3} where {T4} where {T5}
    EmptySeq{Tuple{T1, T2, T3, T4, T5}}()
end

function zip5(x1s::ConsSeq{T1}, x2s::ConsSeq{T2}, x3s::ConsSeq{T3}, x4s::ConsSeq{T4}, x5s::ConsSeq{T5}) where {T1} where {T2} where {T3} where {T4} where {T5}
    return ConsSeq{Tuple{T1, T2, T3, T4, T5}}(() ->
            begin
                statex1, nextx1 = x1s.genfun()
                statex2, nextx2 = x2s.genfun()
                statex3, nextx3 = x3s.genfun()
                statex4, nextx4 = x4s.genfun()
                statex5, nextx5 = x5s.genfun()
                (statex1, statex2, statex3, statex4, statex5), state ->
                    begin
                        _statex1, _statex2, _statex3, _statex4, _statex5 = state
                        x1, newstatex1 = nextx1(_statex1)
                        x2, newstatex2 = nextx2(_statex2)
                        x3, newstatex3 = nextx3(_statex3)
                        x4, newstatex4 = nextx4(_statex4)
                        x5, newstatex5 = nextx5(_statex5)
                        if isnull(x1) || isnull(x2) || isnull(x3) || isnull(x4) || isnull(x5)
                            Nullable{Tuple{T1, T2, T3, T4, T5}}(), (newstatex1, newstatex2, newstatex3, newstatex4, newstatex5)
                        else
                            Nullable{Tuple{T1, T2, T3, T4, T5}}((get(x1), get(x2), get(x3), get(x4), get(x5))), (newstatex1, newstatex2, newstatex3, newstatex4, newstatex5)
                        end
                    end
            end
            )
end

 function zipn(xs::Vector{<:Seq{T}}) where {T}
     EmptySeq{Vector{T}}()
 end

 function zipn(xs::Vector{ConsSeq{T}}) where {T}
     return ConsSeq{Vector{T}}(() ->
         begin
             y = [x.genfun() for x in xs]
             state = [s for (s, _) in y]
             next = [f for (_, f) in y]
             state, sarg ->
                 begin
                     res = [b(a) for (a, b) in zip(sarg, next)]
                     newstate = [x for (_, x) in res]
                     if all([!isnull(a) for (a, b) in res])
                         v = [get(x) for (x, _) in res]
                         Nullable{Vector{T}}(v), newstate
                     else
                         Nullable{Vector{T}}(), newstate
                     end
                 end
         end
     )
 end

function Base.foreach(f::Function, xs::Seq{T}) where {T}
    eof = false
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if isnull(v)
            eof = true
        else
            f(get(v))
        end
    end
end

function Base.collect(xs::Seq{T}) where {T}
    res = Vector{T}()
    eof = false
    tail = xs
    while !eof
        v, tail = tryread(tail)
        if isnull(v)
            eof = true
        else
            push!(res, get(v))
        end
    end
    res
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
md"""
# Panel Graph

AbstractPanels are related to one another through a directed graph
based on various relations.  Here we construct the inverse directed graph,
which could be one-to-many.
"""

export PanelGraph, injest, makePanelGraph
export nodes, edges, query

struct PanelGraph
    graph::DiGraph

    PanelGraph() = new(DiGraph())
end

# Delegation to the underlying DiGraph:

NahaGraphs.add_edge!(pg::PanelGraph, edge::Pair) = add_edge!(pg.graph, edge)

function Base.haskey(pg::PanelGraph, key)::Bool
    haskey(pg.graph, key)
end

Base.keys(pg::PanelGraph) = keys(pg.graph)

Base.values(pg::PanelGraph) = values(pg.graph)

NahaGraphs.nodes(pg::PanelGraph) = union(keys(pg), values(pg))

NahaGraphs.edges(pg::PanelGraph) = edges(pg.graph)

Base.getindex(pg::PanelGraph, key) =
    (p -> p.second).(filter(p -> p.first == key, edges(pg)))

NahaGraphs.query(pg::PanelGraph, from, to) = query(pg.graph, from, to)

# Building the graph from the panel structure

function injest(pg::PanelGraph, pair::Pair)
    add_edge!(pg.graph, pair)
end

function injest(pg::PanelGraph, panel::AbstractPanel)
    # No-op default method
end

function injest(pg::PanelGraph, panel::BoughtPanel)
    add_edge!(pg, (panel.was) => panel)
    injest(pg, panel.was)
end

function injest(pg::PanelGraph, panel::TerminalPanel)
    add_edge!(pg, (panel.was) => panel)
    injest(pg, panel.was)
end

function injest(pg::PanelGraph, panel::Panel)
    add_edge!(pg, (panel.cut_from) => panel)
    injest(pg, panel.cut_from)
end

function makePanelGraph(state::SearchState)::PanelGraph
    pg = PanelGraph()
    function inj(panels)
        for p in panels
            injest(pg, p)
        end
    end
    inj(state.bought)
    inj(state.finished)
    inj(state.scrapped)
    inj(state.working)
    return pg
end


text:
"""
    covellipse(μ, Σ; showaxes=false, n_std=1, n_ellipse_vertices=100)

Plot a confidence ellipse of the 2×2 covariance matrix `Σ`, centered at `μ`.
The ellipse is the contour line of a Gaussian density function with mean `μ`
and variance `Σ` at `n_std` standard deviations.
If `showaxes` is true, the two axes of the ellipse are also plotted.
"""
@userplot CovEllipse

@recipe function f(c::CovEllipse; showaxes=false, n_std=1, n_ellipse_vertices=100)
    μ, S = _covellipse_args(c.args; n_std=n_std)

    θ = range(0, 2π; length=n_ellipse_vertices)
    A = S * [cos.(θ)'; sin.(θ)']

    @series begin
        seriesalpha --> 0.3
        Shape(μ[1] .+ A[1,:], μ[2] .+ A[2,:])
    end
    showaxes && @series begin
        label := false
        linecolor --> "gray"
        ([μ[1] + S[1,1], μ[1], μ[1] + S[1,2]],
         [μ[2] + S[2,1], μ[2], μ[2] + S[2,2]])
    end
end

function _covellipse_args((μ,Σ)::Tuple{AbstractVector{<:Real}, AbstractMatrix{<:Real}}; n_std::Real)
    size(μ) == (2,) && size(Σ) == (2,2) || error(
            "covellipse requires mean of length 2 and covariance of size 2×2.")
    λ, U = eigen(Σ)
    μ, n_std * U * diagm(.√λ)
end
_covellipse_args(args; n_std) = error("Wrong inputs for covellipse: $(typeof.(args)). " *
                                      "Expected real-valued vector μ, real-valued matrix Σ.")


text:
module DiffCheckNumbers

using SpecialFunctions
import NaNMath
import Calculus

include("Constructor.jl")
include("ConversionPromotion.jl")
include("Helper.jl")
include("Arithmetic.jl")
include("Calculus.jl")

export DiffCheck, value, derivatives, StaticTuple

end # module


text:
# A DataFreeTree wraps a descendant of NNTree
# which does not contain a copy of the data
immutable DataFreeTree{N <: NNTree}
    size::Tuple{Int,Int}
    hash::UInt64
    tree::N
end

function get_points_dim(data)
    if eltype(data) <: AbstractVector
        ndim = eltype(eltype(data))
        npoints = length(data)
    elseif typeof(data) <: Matrix
        ndim = size(data, 1)
        npoints = size(data, 2)
    else
        error("Unknown input data format")
    end
    return ndim, npoints
end

"""
    DataFreeTree(treetype, data[, reorderbufffer = similar(data), indicesfor = :data, kargs...]) -> datafreetree

Creates a `DataFreeTree` which wraps a `KDTree` or `BallTree`. Keywords arguments are passed
to their respective constructors.

The `KDTree` or `BallTree` will be stored without a reference to the underlaying data. `injectdata`
has to be used to re-link them to a data array before use.

By default the `reorder` feature of `KDTree`/`BallTree` is turned off. In case a `reorderbuffer`
is provided, reordering is performed and the contents of `reorderbuffer` have to be later provided to
`injectdata`.

`indicesfor` controlls whether the indices returned by the query functions should refer to `data` or the `reorderbuffer`. Valid values are `:data` and `:reordered`.
"""
function DataFreeTree{T<:NNTree}(::Type{T}, data, args...; reorderbuffer = data[:, 1:0], kargs...)
    tree = T(data, args...; storedata = false, reorderbuffer = reorderbuffer, kargs...)
    ndim, npoints = get_points_dim(data)
    DataFreeTree((ndim, npoints), hash(tree.reordered ? reorderbuffer : data), tree)
end

"""
    injectdata(datafreetree, data) -> tree

Returns the `KDTree`/`BallTree` wrapped by `datafreetree`, set up to use `data` for the points data.
"""
function injectdata{T}(datafreetree::DataFreeTree, data::Matrix{T})
    dim = size(data, 1)
    npoints = size(data, 2)
     if isbits(T)
        new_data = reinterpret(SVector{dim, T}, data, (npoints, ))
    else
        new_data = SVector{dim, T}[SVector{dim, T}(data[:, i]) for i in 1:npoints]
    end
    new_hash = hash(data)
    injectdata(datafreetree, new_data, new_hash)
end

function injectdata{V <: AbstractVector}(datafreetree::DataFreeTree, data::Vector{V}, new_hash::UInt64=0)
    if new_hash == 0
        new_hash = hash(data)
    end
    if length(V) != datafreetree.size[1] || length(data) != datafreetree.size[2]
         throw(DimensionMismatch("NearestNeighbors:injectdata: The size of 'data' $(length(data)) × $(length(V)) does not match the data array used to construct the tree $(datafreetree.size)."))
     end
    if new_hash != datafreetree.hash
        throw(ArgumentError("NearestNeighbors:injectdata: The hash of 'data' does not match the hash of the data array used to construct the tree."))
    end

    typ = typeof(datafreetree.tree)
    fields = map(x-> getfield(datafreetree.tree, x), fieldnames(datafreetree.tree))[2:end]
    typ(data, fields...)
end


text:
using Sunny

## Make lattice
crystal = Sunny.diamond_crystal(a=1.0)

## Interactions -- units of K
J = 28.28
interactions = [
    heisenberg(J, Bond{3}(3, 6, [0,0,0])),
]

## Spin system
extent = (4,4,4)
sys = SpinSystem(crystal, interactions, extent)
sys_size = length(sys)
rand!(sys)

## Wang-Landau sampler
N_bins = 1000
bs = 2*J / N_bins
E_min = -53.0
E_max = 0.0
ps = true

WL = WangLandau(
    system=sys,
    bin_size = bs,
    hcheck_interval=10_000,
    bounds=[E_min, E_max],
    mc_step_size=0.2,
    flatness=0.6,
    ln_f_final=1e-6,
    per_spin=ps,
    max_mcsweeps=10_000_000_000
)

## Run WL simulation
@time run!(WL)


text:
# auxiliary.jl

# getindex, fistindex, lastindex
for TM in (:TaylorModel1, :RTaylorModel1, :TaylorModelN)
    @eval begin
        copy(f::$TM) = $TM(copy(f.pol), f.rem, f.x0, f.dom)
        @inline firstindex(a::$TM) = 0
        @inline lastindex(a::$TM) = get_order(a)

        getindex(a::$TM, n::Integer) = a.pol[n]
        getindex(a::$TM, u::UnitRange) = a.pol[u]
        getindex(a::$TM, c::Colon) = a.pol[c]
        # getindex(a::$TM, u::StepRange{Int,Int}) = a.pol[u[:]]

        constant_term(a::$TM) = constant_term(polynomial(a))
        linear_polynomial(a::$TM) = linear_polynomial(polynomial(a))

        # norm
        norm(x::$TM, p::Real=2) = norm(polynomial(x), p)
    end
end

# setindex, iscontained
for TM in tupleTMs
    @eval begin
        setindex!(a::$TM{T,S}, x::T, n::Integer) where {T<:Number, S} = a.pol[n] = x
        setindex!(a::$TM{T,S}, x::T, u::UnitRange) where {T<:Number, S} = a.pol[u] .= x
        function setindex!(a::$TM{T,S}, x::Array{T,1}, u::UnitRange) where {T<:Number, S}
            @assert length(u) == length(x)
            for ind in eachindex(x)
                a.pol[u[ind]] = x[ind]
            end
        end
        setindex!(a::$TM{T,S}, x::T, c::Colon) where {T<:Number, S} = a[c] .= x
        setindex!(a::$TM{T,S}, x::Array{T,1}, c::Colon) where {T<:Number, S} = a[c] .= x

        iscontained(a, tm::$TM) = a ∈ domain(tm)-tm.x0
        iscontained(a::Interval, tm::$TM) = a ⊆ domain(tm)-tm.x0
    end
end
iscontained(a, tm::TaylorModelN) = a ∈ domain(tm)-tm.x0
iscontained(a::IntervalBox, tm::TaylorModelN) = a ⊆ domain(tm)-tm.x0

# fixorder and bound_truncation
for TM in tupleTMs
    @eval begin
        function fixorder(a::$TM, b::$TM)
            @assert tmdata(a) == tmdata(b)
            a.pol.order == b.pol.order && return a, b

            order = min(a.pol.order, b.pol.order)
            apol0, bpol0 = polynomial.((a, b))
            apol, bpol = TaylorSeries.fixorder(apol0, bpol0)

            # Bound for the neglected part of the polynomial
            dom = a.dom - a.x0
            Δa = bound_truncation($TM, apol0, dom, order) + remainder(a)
            Δb = bound_truncation($TM, bpol0, dom, order) + remainder(b)

            return $TM(apol, Δa, a.x0, a.dom), $TM(bpol, Δb, b.x0, b.dom)
        end

        function bound_truncation(::Type{$TM}, a::Taylor1, aux::Interval,
                order::Int)
            order ≥ get_order(a) && return zero(aux)
            if $TM == TaylorModel1
                res = Taylor1(copy(a.coeffs))
                res[0:order] .= zero(res[0])
            else
                res = Taylor1(copy(a.coeffs[order+2:end]), get_order(a)-order )
            end
            return res(aux)
        end

    end
end

function fixorder(a::TaylorModelN, b::TaylorModelN)
    @assert tmdata(a) == tmdata(b)
    a.pol.order == b.pol.order && return a, b

    order = min(a.pol.order, b.pol.order)
    apol0, bpol0 = polynomial.((a, b))
    apol, bpol = TaylorSeries.fixorder(apol0, bpol0)

    # Bound for the neglected part of the polynomial
    dom = a.dom - a.x0
    Δa = bound_truncation(TaylorModelN, apol0, dom, order) + remainder(a)
    Δb = bound_truncation(TaylorModelN, bpol0, dom, order) + remainder(b)

    return TaylorModelN(apol, Δa, a.x0, a.dom), TaylorModelN(bpol, Δb, b.x0, b.dom)
end

function bound_truncation(::Type{TaylorModelN}, a::TaylorN, aux::IntervalBox,
        order::Int)
    order ≥ get_order(a) && return zero(aux[1])
    res = deepcopy(a)
    res[0:order] .= zero(res[0])
    return res(aux)
end

# auxiliary constructors for special intervals and boxes
@inline zero_interval(T) = zero(Interval{T})
@inline zero_box(N, T) = IntervalBox(zero_interval(T), Val(N))
@inline symmetric_box(N, T) = IntervalBox(Interval{T}(-1, 1), Val(N))


text:
using CUDA
using BenchmarkTools

N = 2^20
x_d = CUDA.fill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

function add_broadcast!(y, x)
    CUDA.@sync y .+= x
        return
end

@btime add_broadcast!($y_d, $x_d)

max_stars_repo_name: neversakura/OSQAT-GPU

max_issues_repo_name: neversakura/OSQAT-GPU

max_forks_repo_name: neversakura/OSQAT-GPU


text:
function batch_problem()
    m = Model()
    @variable(m, 0 <= x1 <= 1.38629436111989)
    @variable(m, 0 <= x2 <= 1.38629436111989)
    @variable(m, 0 <= x3 <= 1.38629436111989)
    @variable(m, 0 <= x4 <= 1.38629436111989)
    @variable(m, 0 <= x5 <= 1.38629436111989)
    @variable(m, 0 <= x6 <= 1.38629436111989)
    @variable(m, 5.7037824746562 <= x7 <= 8.00636756765025)
    @variable(m, 5.7037824746562 <= x8 <= 8.00636756765025)
    @variable(m, 5.7037824746562 <= x9 <= 8.00636756765025)
    @variable(m, 5.7037824746562 <= x10 <= 8.00636756765025)
    @variable(m, 5.7037824746562 <= x11 <= 8.00636756765025)
    @variable(m, 5.7037824746562 <= x12 <= 8.00636756765025)
    @variable(m, 4.45966 <= x13 <= 397.747)
    @variable(m, 3.7495 <= x14 <= 882.353)
    @variable(m, 4.49144 <= x15 <= 833.333)
    @variable(m, 3.14988 <= x16 <= 638.298)
    @variable(m, 3.04452 <= x17 <= 666.667)
    @variable(m, 0.729961 <= x18 <= 2.11626)
    @variable(m, 0.530628 <= x19 <= 1.91626)
    @variable(m, 1.09024 <= x20 <= 2.47654)
    @variable(m, -0.133531 <= x21 <= 1.25276)
    @variable(m, 0.0487901 <= x22 <= 1.43508)
    @variable(m, b23, Bin)
    @variable(m, b24, Bin)
    @variable(m, b25, Bin)
    @variable(m, b26, Bin)
    @variable(m, b27, Bin)
    @variable(m, b28, Bin)
    @variable(m, b29, Bin)
    @variable(m, b30, Bin)
    @variable(m, b31, Bin)
    @variable(m, b32, Bin)
    @variable(m, b33, Bin)
    @variable(m, b34, Bin)
    @variable(m, b35, Bin)
    @variable(m, b36, Bin)
    @variable(m, b37, Bin)
    @variable(m, b38, Bin)
    @variable(m, b39, Bin)
    @variable(m, b40, Bin)
    @variable(m, b41, Bin)
    @variable(m, b42, Bin)
    @variable(m, b43, Bin)
    @variable(m, b44, Bin)
    @variable(m, b45, Bin)
    @variable(m, b46, Bin)

    @NLobjective(m, Min, 250*exp(x1 + 0.6*x7) + 250*exp(x2 + 0.6*x8) + 250*exp(x3 + 0.6*x9)
     + 250*exp(x4 + 0.6*x10) + 250*exp(x5 + 0.6*x11) + 250*exp(x6 + 0.6*x12))

    @NLconstraint(m, 250*exp(x1 + 0.6*x7) + 250*exp(x2 + 0.6*x8) + 250*exp(x3 + 0.6*x9)
    + 250*exp(x4 + 0.6*x10) + 250*exp(x5 + 0.6*x11) + 250*exp(x6 + 0.6*x12) <= 290000)

    @constraint(m, x7 - x13 >= 2.06686275947298)
    @constraint(m, x8 - x13 >= 0.693147180559945)
    @constraint(m, x9 - x13 >= 1.64865862558738)
    @constraint(m, x10 - x13 >= 1.58923520511658)
    @constraint(m, x11 - x13 >= 1.80828877117927)
    @constraint(m, x12 - x13 >= 1.43508452528932)
    @constraint(m, x7 - x14 >= -0.356674943938732)
    @constraint(m, x8 - x14 >= -0.22314355131421)
    @constraint(m, x9 - x14 >= -0.105360515657826)
    @constraint(m, x10 - x14 >= 1.22377543162212)
    @constraint(m, x11 - x14 >= 0.741937344729377)
    @constraint(m, x12 - x14 >= 0.916290731874155)
    @constraint(m, x7 - x15 >= -0.356674943938732)
    @constraint(m, x8 - x15 >= 0.955511445027436)
    @constraint(m, x9 - x15 >= 0.470003629245736)
    @constraint(m, x10 - x15 >= 1.28093384546206)
    @constraint(m, x11 - x15 >= 1.16315080980568)
    @constraint(m, x12 - x15 >= 1.06471073699243)
    @constraint(m, x7 - x16 >= 1.54756250871601)
    @constraint(m, x8 - x16 >= 0.832909122935104)
    @constraint(m, x9 - x16 >= 0.470003629245736)
    @constraint(m, x10 - x16 >= 0.993251773010283)
    @constraint(m, x11 - x16 >= 0.182321556793955)
    @constraint(m, x12 - x16 >= 0.916290731874155)
    @constraint(m, x7 - x17 >= 0.182321556793955)
    @constraint(m, x8 - x17 >= 1.28093384546206)
    @constraint(m, x9 - x17 >= 0.8754687373539)
    @constraint(m, x10 - x17 >= 1.50407739677627)
    @constraint(m, x11 - x17 >= 0.470003629245736)
    @constraint(m, x12 - x17 >= 0.741937344729377)
    @constraint(m, x1 + x18 >= 1.85629799036563)
    @constraint(m, x2 + x18 >= 1.54756250871601)
    @constraint(m, x3 + x18 >= 2.11625551480255)
    @constraint(m, x4 + x18 >= 1.3609765531356)
    @constraint(m, x5 + x18 >= 0.741937344729377)
    @constraint(m, x6 + x18 >= 0.182321556793955)
    @constraint(m, x1 + x19 >= 1.91692261218206)
    @constraint(m, x2 + x19 >= 1.85629799036563)
    @constraint(m, x3 + x19 >= 1.87180217690159)
    @constraint(m, x4 + x19 >= 1.48160454092422)
    @constraint(m, x5 + x19 >= 0.832909122935104)
    @constraint(m, x6 + x19 >= 1.16315080980568)
    @constraint(m, x1 + x20 >= 0)
    @constraint(m, x2 + x20 >= 1.84054963339749)
    @constraint(m, x3 + x20 >= 1.68639895357023)
    @constraint(m, x4 + x20 >= 2.47653840011748)
    @constraint(m, x5 + x20 >= 1.7404661748405)
    @constraint(m, x6 + x20 >= 1.82454929205105)
    @constraint(m, x1 + x21 >= 1.16315080980568)
    @constraint(m, x2 + x21 >= 1.09861228866811)
    @constraint(m, x3 + x21 >= 1.25276296849537)
    @constraint(m, x4 + x21 >= 1.19392246847243)
    @constraint(m, x5 + x21 >= 1.02961941718116)
    @constraint(m, x6 + x21 >= 1.22377543162212)
    @constraint(m, x1 + x22 >= 0.741937344729377)
    @constraint(m, x2 + x22 >= 0.916290731874155)
    @constraint(m, x3 + x22 >= 1.43508452528932)
    @constraint(m, x4 + x22 >= 1.28093384546206)
    @constraint(m, x5 + x22 >= 1.30833281965018)
    @constraint(m, x6 + x22 >= 0.78845736036427)

    @NLconstraint(m, 250000*exp(x18 - x13) + 150000*exp(x19 - x14) + 180000*exp(x20 - x15) +
    160000*exp(x21 - x16) + 120000*exp(x22 - x17) <= 6000)

    @NLconstraint(m,x1 - 0.693147180559945*b29 - 1.09861228866811*b35
    - 1.38629436111989*b41 == 0)

    @NLconstraint(m,x2 - 0.693147180559945*b30 - 1.09861228866811*b36
    - 1.38629436111989*b42 == 0)

    @NLconstraint(m, x3 - 0.693147180559945*b31 - 1.09861228866811*b37
    - 1.38629436111989*b43 == 0)

    @NLconstraint(m, x4 - 0.693147180559945*b32 - 1.09861228866811*b38
    - 1.38629436111989*b44 == 0)

    @NLconstraint(m, x5 - 0.693147180559945*b33 - 1.09861228866811*b39
    - 1.38629436111989*b45 == 0)

    @NLconstraint(m,  x6 - 0.693147180559945*b34 - 1.09861228866811*b40
    - 1.38629436111989*b46 == 0)

    @constraint(m,b23 + b29 + b35 + b41 == 1)
    @constraint(m,b24 + b30 + b36 + b42 == 1)
    @constraint(m,b25 + b31 + b37 + b43 == 1)
    @constraint(m,b26 + b32 + b38 + b44 == 1)
    @constraint(m,b27 + b33 + b39 + b45 == 1)
    @constraint(m,b28 + b34 + b40 + b46 == 1)
    return m
end

function cvxnonsep_nsig20r_problem()
    m = Model()
    @variable(m, 1 <= i1 <= 10, Int)
    @variable(m, 1 <= i2 <= 10, Int)
    @variable(m, 1 <= i3 <= 10, Int)
    @variable(m, 1 <= i4 <= 10, Int)
    @variable(m, 1 <= i5 <= 10, Int)
    @variable(m, 1 <= i6 <= 10, Int)
    @variable(m, 1 <= i7 <= 10, Int)
    @variable(m, 1 <= i8 <= 10, Int)
    @variable(m, 1 <= i9 <= 10, Int)
    @variable(m, 1 <= i10 <= 10, Int)
    @variable(m, 1E-5 <= x11 <= 10)
    @variable(m, 1E-5 <= x12 <= 10)
    @variable(m, 1E-5 <= x13 <= 10)
    @variable(m, 1E-5 <= x14 <= 10)
    @variable(m, 1E-5 <= x15 <= 10)
    @variable(m, 1E-5 <= x16 <= 10)
    @variable(m, 1E-5 <= x17 <= 10)
    @variable(m, 1E-5 <= x18 <= 10)
    @variable(m, 1E-5 <= x19 <= 10)
    @variable(m, 1E-5 <= x20 <= 10)
    @variable(m, x21)
    @variable(m, x22)
    @variable(m, x23)
    @variable(m, x24)
    @variable(m, x25)
    @variable(m, x26)
    @variable(m, x27)
    @variable(m, x28)
    @variable(m, x29)
    @variable(m, x30)
    @variable(m, x31)
    @variable(m, x32)
    @variable(m, x33)
    @variable(m, x34)
    @variable(m, x35)
    @variable(m, x36)
    @variable(m, x37)
    @variable(m, x38)
    @variable(m, x39)
    @variable(m, x40)
    @constraint(m,  x21 + x22 + x23 + x24 + x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 + x40 <= -1.6094379124341)
    @NLconstraint(m, -0.065*log(i1) - x21 <= 0)
    @NLconstraint(m, -0.004*log(i2) - x22 <= 0)
    @NLconstraint(m, -0.084*log(i3) - x23 <= 0)
    @NLconstraint(m, -0.093*log(i4) - x24 <= 0)
    @NLconstraint(m, -0.06*log(i5) - x25 <= 0)
    @NLconstraint(m, -0.075*log(i6) - x26 <= 0)
    @NLconstraint(m, -0.074*log(i7) - x27 <= 0)
    @NLconstraint(m, -0.039*log(i8) - x28 <= 0)
    @NLconstraint(m, -0.065*log(i9) - x29 <= 0)
    @NLconstraint(m, -0.017*log(i10) - x30 <= 0)
    @NLconstraint(m, -0.07*log(x11) - x31 <= 0)
    @NLconstraint(m, -0.03*log(x12) - x32 <= 0)
    @NLconstraint(m, -0.028*log(x13) - x33 <= 0)
    @NLconstraint(m, -0.005*log(x14) - x34 <= 0)
    @NLconstraint(m, -0.01*log(x15) - x35 <= 0)
    @NLconstraint(m, -0.082*log(x16) - x36 <= 0)
    @NLconstraint(m, -0.069*log(x17) - x37 <= 0)
    @NLconstraint(m, -0.032*log(x18) - x38 <= 0)
    @NLconstraint(m, -0.095*log(x19) - x39 <= 0)
    @NLconstraint(m, -0.003*log(x20) - x40 <= 0)
    @NLobjective(m, Min,  1.5*i1 + 0.51*i2 + 1.01*i3 + 1.4*i4 + 1.78*i5 + 1.92*i6
    + 1.09*i7 + 0.48*i8 + 0.67*i9 + 0.52*i10 + 1.68*x11 + 0.51*x12 + 1.63*x13
    + 0.49*x14 + 1.86*x15 + 0.7*x16 + 0.39*x17 + 0.5*x18 + 1.23*x19 + 0.95*x20)

    return m
end


text:
"""
Merkle's Puzzles
================

A toy implementation of Merkle's puzzles, an early, impractical key exchange algorithm devised by
Ralph Merkle.

# References

- https://en.wikipedia.org/wiki/Merkle%27s_Puzzles
- https://medium.com/100-days-of-algorithms/day-75-merkles-puzzles-d9f0e8f9c9d0
"""
module MerklePuzzles

using Random: shuffle
using Nettle: digest

#====
The message format mostly follows Tomáš Bouda's Medium article.

    secret  = [ KEY_SIZE | INDEX_SIZE] | DIGEST_SIZE
    otp_key =       DIGEST_SIZE        | DIGEST_SIZE
====#

sha1(s::Vector{UInt8}) = digest("sha1", s)

# Set it to 1 if Eve takes too long to solve all the puzzles.
const INDEX_SIZE = 2
const N_PUZZLES = 1 << (8 * INDEX_SIZE)
const DIGEST_SIZE = length(sha1(UInt8[]))
# may be larger than the actual symmetric key size
const KEY_SIZE = DIGEST_SIZE - INDEX_SIZE
const OTP_SEED_SIZE = 10

@assert KEY_SIZE > 0
@assert OTP_SEED_SIZE > INDEX_SIZE

"Represents each of Merkle's puzzles."
struct Puzzle
    "The encrypted message, length = DIGEST_SIZE * 2"
    encrypted::Vector{UInt8}

    "Partial OTP key, length = OTP_SEED_SIZE - INDEX_SIZE"
    partial_otp_seed::Vector{UInt8}
end

# Helpers
@inline function pack_bigendian(a::AbstractVector{UInt8}, x::Int)
    x = x - 1
    for i in length(a):-1:1
        a[i] = x & 0xFF
        x >>= 8
    end
end

@inline function unpack_bigendian(a::AbstractVector{UInt8})::Int
    x = 0
    for i in 1:length(a)
        x <<= 8
        x += a[i]
    end
    x + 1
end

# TODO: refactor the API.
"Returns `otp_key`. Randomly sets `otp_seed` if it is empty."
function make_otp_key!(otp_seed::Vector{UInt8})::Vector{UInt8}
    if isempty(otp_seed)
        resize!(otp_seed, OTP_SEED_SIZE)
        otp_seed .= rand(UInt8, OTP_SEED_SIZE)
    end

    otp_key = Vector{UInt8}(undef, DIGEST_SIZE * 2)

    otp_digest = sha1(otp_seed)
    otp_key[1:DIGEST_SIZE] = otp_digest

    otp_digest = sha1(otp_digest)
    otp_key[(DIGEST_SIZE + 1):end] = otp_digest

    otp_key
end

function make_puzzle(key::Vector{UInt8}, index::Int)::Puzzle
    @assert length(key) == KEY_SIZE
    @assert 1 <= index <= N_PUZZLES

    secret = Vector{UInt8}(undef, DIGEST_SIZE * 2)
    secret[1:KEY_SIZE] = key
    index_view = view(secret, (KEY_SIZE + 1):DIGEST_SIZE)
    pack_bigendian(index_view, index)
    # NOTE:
    # 1. Actually index can be excluded from the digestion.
    # 2. The digest must be encrypted too, or it will be insecure.
    secret[(DIGEST_SIZE + 1):end] = sha1(secret[1:DIGEST_SIZE])

    otp_seed = UInt8[]
    otp_key = make_otp_key!(otp_seed)
    @assert length(otp_seed) == OTP_SEED_SIZE
    encrypted = map(xor, secret, otp_key)
    partial_otp_seed = otp_seed[(INDEX_SIZE + 1):end]

    Puzzle(encrypted, partial_otp_seed)
end

function generate_puzzles(n::Int, keys::Matrix{UInt8})
    [make_puzzle(keys[i, :], i) for i in 1:n]
end

"Returns `(index, key)` pair, or `nothing` if the puzzle cannot be solved."
function solve_puzzle(encrypted::Vector{UInt8},
                      otp_seed_template::Vector{UInt8},
                      key_index::Int)::Union{Tuple{Int, Vector{UInt8}}, Nothing}
    key_index_view = view(otp_seed_template, 1:INDEX_SIZE)
    pack_bigendian(key_index_view, key_index)

    otp_key = make_otp_key!(otp_seed_template)
    decrypted = map(xor, encrypted, otp_key)
    if sha1(decrypted[1:DIGEST_SIZE]) == decrypted[(DIGEST_SIZE + 1):end]
        decrypted = decrypted[1:DIGEST_SIZE]
        index_view = view(decrypted, (KEY_SIZE + 1):DIGEST_SIZE)
        index = unpack_bigendian(index_view)
        (index, decrypted[1:KEY_SIZE])
    else
        nothing
    end
end

# The puzzles created by Bob should be serialized before being sent to Alice.
# We assume the (de)serialization code already exists and just pass around Julia objects.

function alice(puzzles::Vector{Puzzle})::Union{Tuple{Int, Vector{UInt8}}, Nothing}
    println("""
        Alice: Puzzles received.
        Alice: Pick a random one and solve it by brute force.
        """)
    i = rand(1:N_PUZZLES)
    puzzle = puzzles[i]
    otp_seed_template = [zeros(UInt8, INDEX_SIZE); puzzle.partial_otp_seed]

    for j in 1:N_PUZZLES
        maybe_index_key = solve_puzzle(puzzle.encrypted, otp_seed_template, j)
        if maybe_index_key != nothing
            return maybe_index_key
        end
    end

    println("Alice: Fsck, I can't solve it. Bob must be a liar.")
    nothing
end

function bob()::Tuple{Matrix{UInt8}, Vector{Puzzle}}
    keys = rand(UInt8, (N_PUZZLES, KEY_SIZE))
    puzzles = generate_puzzles(N_PUZZLES, keys)
    shuffled_puzzles = shuffle(puzzles)
    println("""
        Bob: Puzzles generated.
        Bob: Send puzzles to Alice via insecure channel.
        """)
    (keys, shuffled_puzzles)
end

function eve(puzzles::Vector{Puzzle})::Union{Vector{Vector{UInt8}}, Nothing}
    println("""
        Eve: Puzzles eavesdropped.
        Eve: Start to crack every puzzle.
        """)
    solved_puzzles = Vector{Union{Vector{UInt8}, Nothing}}(nothing, N_PUZZLES)

    for i in 1:N_PUZZLES
        puzzle = puzzles[i]
        otp_seed_template = [zeros(UInt8, INDEX_SIZE); puzzle.partial_otp_seed]

        for j in 1:N_PUZZLES
            maybe_index_key = solve_puzzle(puzzle.encrypted, otp_seed_template, j)
            if maybe_index_key != nothing
                (index, key) = maybe_index_key

                if solved_puzzles[index] != nothing
                    error("Duplicate indices")
                end

                solved_puzzles[index] = key
                @goto succeed
            end
        end

        return nothing
        @label succeed
    end

    solved_puzzles
end

function run()
    println("""
        ======================
        =  Merkle's Puzzles  =
        ======================
        """)
    (bob_keys, shuffled_puzzles) = bob()

    alice_resp = alice(shuffled_puzzles)
    if alice_resp == nothing
        error("Alice cannot solve the puzzle.")
    end
    (index, dec_key) = alice_resp
    println("*: Alice sends index $index back to Bob.")
    if bob_keys[index, :] == dec_key
        println("""
            *: Alice and Bob have agreed on a shared key to use
            *: ...and they communicate happily ever after.
            """)
    else
        error("Alice and Bob have not reached consensus.")
    end

    # Eve the eavesdropper knows the index too.
    eve_resp = eve(shuffled_puzzles)
    if eve_resp == nothing
        println("Eve cannot crack every puzzle.")
    end
    dec_keys = eve_resp
    if bob_keys[index, :] == dec_keys[index]
        println("*: Eve has successfully cracked the key.")
    else
        println("*: Eve has NOT cracked the key.")
    end

    println("\nEnd.")
end

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
import Base.*
import Base.convert
import Base.promote_rule

# Hierarchia typow

function type_hierarchy(t)
    if(t != Any)
        type_hierarchy(supertype(t))
        print(" -> ")
    end

    print(t)
end

type_hierarchy(Float32)

# Struktura

struct G{N} <: Integer
    x::Int

    function G(m)
        if(m >= N)
            m = m % N
        end

        if(m < 0 || gcd(m, N) != 1)
            throw(DomainError())
        end

        new(m)
    end
end

# Mnozenie

function *{N}(a::G{N}, b::G{N})
    return G{N}((a.x * b.x) % N)
end

function *{N}(a::G{N}, b::Integer)
    return a * G{N}(b)
end

function *{N}(a::Integer, b::G{N})
    return G{N}(a) * b
end

# Konwersja

function convert{N}(::Type{G{N}}, v::Int64)
    return G{N}(v)
end

function convert{N}(::Type{Int64}, v::G{N})
    return v.x
end

# Promocja

promote_rule(::Type{G{N}}, ::Type{T}) where {T <: Integer, N} = T

# Potegowanie

function ^{N}(a::G{N}, b::Integer)
    ret::G{N} = 1
    for i = 1:b
       ret = ret * a
    end
    return ret
end

# Okres

function period{N}(a::G{N})
    for i = 1:N
        if a ^ i == 1
            return i
        end
    end

    return -1
end

# Element odwrotny

function inverse{N}(a::G{N})
    for i = 1:N
        try
            if a * i == 1
                return G{N}(i)
            end
        end
    end

    throw(DomainError())
end

# Moc zbioru elementów

function card{N}(::Type{G{N}})
    count = 0

    for i = 0:N
        try
            G{N}(i)
            count += 1
        end
    end

    return count
end

# Testowanie

println(G{13}(7) * G{13}(3))

println(promote_type(G{5}, Int8))

println(G{7}(5)^5)

println(period(G{1023}(1022)))

println(inverse(G{1023}(1019)))

println(card(G{1029}))

# RSA

N = 55
c = G{N}(17)
b = G{N}(4)

r = period(b)
println(r)

d = inverse(c)
println(d)

a = b ^ d
println(a)

println(b, " ", (a ^ c))


text:
using Documenter, Weave, Econometrics
using CSV, RDatasets, StatsBase

for file ∈ readdir(joinpath(dirname(pathof(Econometrics)), "..", "docs", "jmd"))
      weave(joinpath(dirname(pathof(Econometrics)), "..", "docs", "jmd", file),
            out_path = joinpath(dirname(pathof(Econometrics)), "..", "docs", "src"),
            doctype = "github")
end

makedocs(format = Documenter.HTML(assets = ["assets/custom.css"]),
         modules = [Econometrics, StatsBase, CSV, RDatasets],
         sitename = "Econometrics.jl",
         pages = ["Introduction" => "index.md",
                  "Getting Started" => "getting_started.md",
                  "Estimators" => "estimators.md",
                  "API" => "api.md"]
    )

deploydocs(repo = "github.com/Nosferican/Econometrics.jl.git")


text:
#using FITSIO
#using Pkg
#using Dates

function make_ccf_fits_header(metadata::AbstractDict; line_list_filename::Union{Nothing,String} = "unknown",
            v::Union{Nothing,Float64} = nothing, e_v::Union{Nothing,Float64} = nothing,
            blue_ord::Union{Nothing,Int64} = 118, red_ord::Union{Nothing,Int64} = 89,  # TODO: WARNING: Hardwired
            mask_shape::Union{Nothing,String} = "tophat", barycor::Union{Nothing,String} = "chromatic", # TODO: WARNING: Hardwired
            ccord::Union{Nothing,Bool} = false, blazenorm::Union{Nothing,Bool} = false, contnorm::Union{Nothing,Bool} = true,
            div_tell::Union{Nothing,Bool} = false, rawnorm::Union{Nothing,Bool} = nothing, plate_scale::Union{Nothing,Float64} = nothing,
            )
    # TODO: WARNING: several parameters, e.g., BLUE_ORD, RED_ORD, WINDOW, BARYCOR, PLATESCL, given default value above instead of defaulting to outputing nothing
    header_dict = OrderedDict{String,Tuple{Any,String}}("VERSION"=>("penn.state." * string(Pkg.project().version), string(Pkg.project().name) * " code version"))
    # Add version info for packages used
    deps = Pkg.dependencies()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        header_dict["VERSION_" * dep.name] = (string(dep.version), "version of " * dep.name)
    end
    # header fields similar to EXPRES CCFs
    header_dict["DATE-CCF"] = (string(Dates.now()), "Date and time of CCF calculation")
    if haskey(metadata,:Filename)
        header_dict["FILENAME"] = (last(splitpath(metadata[:Filename])), "Original filename")
    end
    if haskey(metadata,:bjd)
        header_dict["MJD"] = (metadata[:bjd], "MJD of geometric midpoint")
    end
    if !isnothing(line_list_filename)
        header_dict["MASK"] = (line_list_filename, "CCF mask linelist used")
    end
    if haskey(metadata,:expcount)
        header_dict["EXPCOUNT"] = (metadata[:expcount], "Exposure meter counts")
    end
    # TODO: Skipped: SNR, CHI2 from CCF FITS headers
    if !isnothing(v)
        header_dict["V"] = (v, "Fitted velocity in cm/s")
    end
    if !isnothing(e_v)
        header_dict["E_V"] = (e_v, "Error of velocity in cm/s")
    end
    if !isnothing(blue_ord)
        header_dict["BLUE_ORD"] = (blue_ord, "Bluest echelle order in CCF")
    end
    if !isnothing(red_ord)
        header_dict["RED_ORD"] = (red_ord, "Reddest echelle order in CCF")
    end
    if !isnothing(mask_shape)
        header_dict["WINDOW"] = (mask_shape, "Shape of CCF mask used in CCF")
    end
    if !isnothing(barycor)
        header_dict["BARYCOR"] = (barycor, "Type of barycentric correction")
    end
    if haskey(metadata,:airmass)
        header_dict["AIRMASS"] = (metadata[:airmass], "Airmass of exposure")
    end
    if haskey(metadata,:exposure_time)
        header_dict["EXPTIME"] = (metadata[:exposure_time], "Exposure time (s)")
    end
    if haskey(metadata,:expres_epoch)
        header_dict["RVEPOCH"] = (metadata[:expres_epoch], "RV calibration epoch")
    end
    if haskey(metadata,:wavecal)
        header_dict["WAVE_CAL"] = (metadata[:wavecal], "Wavelength calibration used when computing CCF")
    end
    if !isnothing(ccord)
        header_dict["CCOR"] = (ccord, "HARPS_style colour correction")
    end
    if !isnothing(blazenorm)
        header_dict["BLAZE"] = (blazenorm, "Pixels weighted by blaze in CCF")
    end
    if !isnothing(contnorm)
        header_dict["CONTNORM"] = (contnorm, "Continuum normalized spectrum in CCF")
    end
    # TODO: Warning: Hardcoded SUB_CONT
    header_dict["SUB_CONT"] = (false, "Continuum subtracted before running CCF")
    if !isnothing(div_tell)
        header_dict["DIV_TELL"] = (div_tell, "Tellurics divided before CCF")
    end
    # TODO: WARNING: PLATESCL hard-wired for now
    if !isnothing(plate_scale)
        header_dict["PLATESCL"] = (plate_scale, "CCF platescale")
    end
    #=
    # Added some that I thought might be useful to have
    if haskey(metadata,:normalization)
        header_dict["NORMALIZATION"] = (string(metadata[:normalization]), "How spectrum was normalized before computing CCF")
    end
    if haskey(metadata,:snr_prelim)
        header_dict["SNR_PRELIM"] = (string(metadata[:snr_prelim]), "Something related to SNR, but different from what EXPRES provided")
    end
    if haskey(metadata,:pwv)
        header_dict["PWV"] = (metadata[:pwv], "precip. water vabor")
    end
    if haskey(metadata,:sundist)
        header_dict["SUN_DIST"] = (metadata[:sundist], "Distance to sun")
    end
    if haskey(metadata,:moondist)
        header_dict["MOON_DIST"] = (metadata[:moondist], "Distance to moon")
    end
    =#
    fits_header = FITSHeader(collect(keys(header_dict)),first.(values(header_dict)),last.(values(header_dict)))
end

function write_ccf_fits(filename::String, v_grid::AbstractArray{T1,1}, ccf_total::AbstractArray{T2,1}, ccf_total_var::AbstractArray{T3,1};
            orders::AbstractArray{T4,1}=zeros(Int64,0), ccf_orders::AbstractArray{T5,2}=zeros(T2,0,0), ccf_orders_var::AbstractArray{T6,2}=zeros(T3,0,0),
            fits_hdr::FITSHeader,
            #total_vel::Union{Nothing,Real} = nothing, sigma_total_vel::Union{Nothing,Real} = nothing,
            order_vels::Union{Nothing,AbstractArray{T7,1}} = nothing, sigma_order_vels::Union{Nothing,AbstractArray{T8,1}} = nothing ) where {
                T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real, T7<:Real, T8<:Real   }
    f = FITS(filename, "w")
    #if isnothing(total_vel) || isnothing(sigma_total_vel)
    write(f,UInt8[],header=fits_hdr,name="metadata")
    #else
    #    write(f,UInt8[],header=fits_hdr,name="metadata",v=total_vel,e_V=sigma_total_vel)
    #end
    t2 = OrderedDict("V_grid"=>collect(v_grid),"ccf"=>ccf_total,"e_ccf"=>sqrt.(ccf_total_var) )
    write(f,collect(keys(t2)),collect(values(t2)),name="ccf_total")
    if length(orders) > 0 && length(ccf_orders) > 0 && length(ccf_orders_var) > 0
        if isnothing(order_vels) || isnothing(sigma_order_vels)
            t3 = OrderedDict("orders"=>convert.(Int64,orders),"ccfs"=>ccf_orders,"errs"=>sqrt.(ccf_orders_var) )
        else
            t3 = OrderedDict("orders"=>convert.(Int64,orders),"ccfs"=>ccf_orders,"errs"=>sqrt.(ccf_orders_var),
                                "v" => order_vels,"e_v" => sigma_order_vels)
        end
        write(f,collect(keys(t3)),collect(values(t3)),name="ccf_order")
    end
    close(f)
end

function write_each_ccf_fits(metadata::AbstractArray{Dict{Symbol,Any},1}, v_grid::AbstractArray{T1,1}, ccf_total::AbstractArray{T2,2}, ccf_total_var::AbstractArray{T3,2};
            orders::AbstractArray{T4,1}=zeros(Int64,0), ccf_orders::AbstractArray{T5,3}=zeros(T2,0,0,size(ccf_total,2)), ccf_orders_var::AbstractArray{T6,3}=zeros(T3,0,0,size(ccf_total,2)),
            fits_hdr::FITSHeader, line_list_filename::Union{Nothing,String} = "unknown", #  )
            total_vel::Union{Nothing,AbstractArray{T7,1}} = nothing, sigma_total_vel::Union{Nothing,AbstractArray{T7,1}} = nothing,
            order_vels::Union{Nothing,AbstractArray{T8,2}} = nothing, sigma_order_vels::Union{Nothing,AbstractArray{T8,2}} = nothing ) where {
                T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real, T7<:Real, T8<:Real }
    @assert length(metadata) == size(ccf_total,2) == size(ccf_total_var,2)
    num_files = length(metadata)
    for i in 1:num_files
        fits_fn = last(splitpath(metadata[i][:Filename]))
        ccf_fn = replace(fits_fn, ".fits" => "_ccf.fits")
        if isnothing(total_vel) || isnothing(sigma_total_vel)
            fits_hdr = make_ccf_fits_header(metadata[i], line_list_filename=line_list_filename)
        else
            fits_hdr = make_ccf_fits_header(metadata[i],v=total_vel[i],e_v=sigma_total_vel[i], line_list_filename=line_list_filename)
        end
        if isnothing(order_vels) || isnothing(sigma_order_vels)
            write_ccf_fits(ccf_fn, v_grid, ccf_total[:,i], ccf_total_var[:,i],
                orders=orders, ccf_orders=ccf_orders[:,:,i], ccf_orders_var=ccf_orders_var[:,:,i],
                fits_hdr=fits_hdr)
        else
            @assert size(order_vels) == size(sigma_order_vels)
            @assert size(order_vels,1) == num_files
            @assert size(order_vels,2) == length(orders)
            write_ccf_fits(ccf_fn, v_grid, ccf_total[:,i], ccf_total_var[:,i],
                orders=orders, ccf_orders=ccf_orders[:,:,i], ccf_orders_var=ccf_orders_var[:,:,i],
                fits_hdr=fits_hdr,
                order_vels=order_vels[i,:], sigma_order_vels=sigma_order_vels[i,:] )
        end
        # break  # WARNING: Aborts after one file for testing purposes
    end
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
@dict_readable struct InitializedEventArguments <: Outbound
end

@dict_readable struct StoppedEventArguments <: Outbound
    reason::String
    description::Union{Missing,String}
    threadId::Union{Missing,Int}
    preserveFocusHint::Union{Missing,Bool}
    text::Union{Missing,String}
    allThreadsStopped::Union{Missing,Bool}
end

@dict_readable struct ContinuedEventArguments <: Outbound
    threadId::Int
    allThreadsContinued::Union{Missing,Bool}
end

@dict_readable struct ExitedEventArguments <: Outbound
    exitCode::Int
end

@dict_readable struct TerminatedEventArguments <: Outbound
    restart::Union{Missing,Any}
end

@dict_readable struct ThreadEventArguments <: Outbound
    reason::String
    threadId::Int
end

@dict_readable struct OutputEventArguments <: Outbound
    category::Union{Missing,String}
    output::String
    variablesReference::Union{Missing,Int}
    source::Union{Missing,Source}
    line::Union{Missing,Int}
    column::Union{Missing,Int}
    data::Union{Missing,Any}
end

@dict_readable struct BreakpointEventArguments <: Outbound
    reason::String
    breakpont::Breakpoint
end

@dict_readable struct ModuleEventArguments <: Outbound
    reason::String
    mod::DAModule
end

@dict_readable struct LoadedSourceEventArguments <: Outbound
    reason::String
    source::Source
end

@dict_readable struct ProcessEventArguments <: Outbound
    name::String
    systemProcessId::Union{Missing,Int}
    isLocalProcess::Union{Missing,Bool}
    startMethod::Union{Missing,String}
    pointerSize::Union{Missing,Int}
end

@dict_readable struct CapabilitiesEventArguments <: Outbound
    capabilities::Capabilities
end


text:
"""
    EVA

Abstract type containing the extreme value model types.

- BlockMaxima
- ThresholdExceedance

"""
abstract type EVA end

struct paramfun
    covariate::Vector{<:DataItem}
    fun::Function
end

Base.Broadcast.broadcastable(obj::EVA) = Ref(obj)

"""
    computeparamfunction(covariates::Vector{Variable})

Establish the parameter as function of the corresponding covariates.

"""
function computeparamfunction(covariates::Vector{<:DataItem})::Function

    fun =
    if isempty(covariates)
        function(β::Vector{<:Real})
            return identity(β)
        end
    else
        X = ones(length(covariates[1].value))

        for cov in covariates
            X = hcat(X, cov.value)
        end
        function(β::Vector{<:Real})
            return X*β
        end
    end
    return fun

end

"""
    getcovariatenumber(model::EVA)::Int

Return the number of covariates.

"""
function getcovariatenumber end

"""
    getdistribution(model::EVA, θ::Vector{<:Real})
    getdistribution(fm::fittedEVA)

Return the distributions corresponding to the model or the fitted model.

If an extreme value model is provided, the distributions corresponding to the
parameter vector `θ` are returned. If a fitted extreme value model is provident,
the distributions corresponding to the parameter estimates are returned.

# Implementation

In the stationary case, a single extreme value distribution is returned.

In the non-stationary case, a vector of extreme value distributions is returned,
one for each data value.

In the Bayesian fitted model case, a array of distributions is returned where
each column corresponds to a MCMC iteration.

"""
function getdistribution end

"""
    getinitialvalue(model::EVA)

Get an initial estimates of the model parameters.

"""
function getinitialvalue end

"""
    paramindex(model::EVA)

Return the postitions corresponding to the location, scale and shape parameter.
"""
function paramindex end

"""
    loglike(model::EVA, θ::Vector{<:Real})

Compute the model loglikelihood evaluated at θ.

"""
function loglike(model::EVA, θ::Vector{<:Real})::Real

    y = model.data.value

    pd = getdistribution(model, θ)

    ll = sum(logpdf.(pd, y))

    return ll

end

"""
    quantile(model::EVA, θ::Vector{<:Real}, p::Real)

Compute the quantile of level `p` from the model evaluated at `θ"".

If the model is non-stationary, then the *effective quantiles* are returned,
*i.e.* one for each covariate value.

"""
function quantile(model::EVA, θ::Vector{<:Real}, p::Real)::Vector{<:Real}

    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    pd = getdistribution(model, θ)

    q = quantile.(pd, p)

    return q

end

"""
    validatestationarity(model::T)::T where T<:EVA

Throw warning if the model is nonstationary.

"""
function validatestationarity(model::T)::T where T<:EVA

    if getcovariatenumber(model) > 0
        @warn "covariates cannot be included in the model when estimating the
            paramters by the probability weighted moment parameter estimation.
            The estimates for the stationary model is returned."

        return T(model.data)
    end

    return model

end

"""
    Base.show(io::IO, obj::EVA)

Override of the show function for the objects of type EVA.

"""
function Base.show(io::IO, obj::EVA)

    showEVA(io, obj)

end

"""
    showparamfun(name::String, param::paramfun)::String

Constructs a string describing a parameter `param` with name `name`.

"""
function showparamfun(name::String, param::paramfun)::String

    covariate = [" + $(x.name)" for x in param.covariate]
    res = string("$name ~ 1", covariate...)

    return res

end

"""
    validatelength(length::Real, explvariables::Vector{Variable})

Validate that the explanatory variables are of length `length`.
"""
function validatelength(n::Real, explvariables::Vector{<:DataItem})

    for explvariable in explvariables
        @assert length(explvariable.value) == n "The explanatory variable length should match data length."
    end

end

include(joinpath("eva", "blockmaxima.jl"))
include(joinpath("eva", "thresholdexceedance.jl"))


text:
# test DrakeVisualizer

# for more low level interaction with the DraekVisualizer/Director, please see
# www.github.com/rdeits/DrakeVisualizer.jl

using Caesar, RoME, CoordinateTransformations, Rotations
using GeometryTypes
using ColorTypes: RGB
# for illustration only
using DrakeVisualizer
import DrakeVisualizer: Triad

vis = startdefaultvisualization(draworigin=true)
# Launch the viewer application if it isn't running already:
# DrakeVisualizer.any_open_windows() || DrakeVisualizer.new_window();

setgeometry!(vis[:basics][:triad], Triad())
settransform!(vis[:basics][:triad], Translation(0.0,0,0) ∘ LinearMap(Rotations.AngleAxis(pi/4,0,0,1.)  ))

setgeometry!(vis[:basics][:box], HyperRectangle(Vec(0.,0,0), Vec(1.,1,1)) )
settransform!(vis[:basics][:box], Translation(2.0,0,0) ∘ LinearMap(Rotations.AngleAxis(0.0,0,0,1.)  ))

# and a point cloud

pointcloud = PointCloud([[x, 0, 0.5] for x in linspace(-1, 1)])
pointcloud.channels[:rgb] = [RGB(g, (1-g), 0) for g in linspace(0, 1)]
setgeometry!(vis[:pointcloud], pointcloud)
delete!(vis[:pointcloud])

drawLine!(vis[:testing][:lines2], 3*randn(3), 1+randn(3) , color=RGBA(0.0,1,0.0,0.5) )

#


text:
@testset "Comparison" begin
    for mod in (StructProg, SDDP)
        @test mod._isapproxzero(0, 1e-10)
        @test !mod._isapproxzero(1, 1e-10)
        @test mod._isapproxzero(1e-5, 1e-4)
        @test !mod._isapproxzero(1e-3, 1e-4)
        @test mod._isapproxzero([1e-5, 1e-5], 1e-4)
        @test !mod._isapproxzero([1e-5, 1e-3], 1e-4)
        @test mod._isapprox(0, 0, 1e-3)
        @test mod._isapprox(1, 1, 1e-3)
        @test !mod._isapprox(0, 1, 1e-3)
        @test mod._lt(0, 1, 1e-3)
        @test !mod._lt(1, 1, 1e-3)
        @test !mod._lt(1, 0, 1e-3)
        @test !mod._lt(1, BigInt(1), 1e-3)
    end
end