text:
export
ruleSPBetaOutNPP,
ruleSPBetaOutNMM,
ruleSPBetaAMNM,
ruleSPBetaBMMN,
ruleVBBetaOut,
ruleVBBetaA,
ruleVBBetaB

ruleSPBetaOutNPP(msg_out::Nothing,
                 msg_a::Message{PointMass, Univariate},
                 msg_b::Message{PointMass, Univariate}) =
    Message(Univariate, Beta, a=deepcopy(msg_a.dist.params[:m]), b=deepcopy(msg_b.dist.params[:m]))

function ruleSPBetaOutNMM(msg_out::Nothing,
                          msg_a::Message{<:FactorNode, Univariate},
                          msg_b::Message{<:FactorNode, Univariate})

    n_samples = default_n_samples
    samples_a = sample(msg_a.dist,n_samples)
    samples_b = sample(msg_b.dist,n_samples)
    s_list = betainvcdf.(samples_a, samples_b, rand(n_samples))
    w_list = ones(n_samples)/n_samples

    return Message(Univariate, SampleList, s=s_list, w=w_list)
end

function ruleSPBetaAMNM(msg_out::Message{<:FactorNode, Univariate},
                        msg_a::Nothing,
                        msg_b::Message{<:FactorNode, Univariate})

    n_samples = default_n_samples
    samples_out = sample(msg_out.dist,n_samples)
    samples_b = sample(msg_b.dist,n_samples)
    logp(x,a,b) = (a-1)*log(x) + (b-1)*log(1.0-x) - loggamma(a) - loggamma(b) + loggamma(a+b)
    logp(a) = sum(logp.(samples_out,a,samples_b))/n_samples

    return Message(Univariate, Function, log_pdf=logp)
end

function ruleSPBetaBMMN(msg_out::Message{FactorNode, Univariate},
                        msg_a::Message{FactorNode, Univariate},
                        msg_b::Nothing)

    n_samples = default_n_samples
    samples_out = sample(msg_out.dist,n_samples)
    samples_a = sample(msg_a.dist,n_samples)
    logp(x,a,b) = (a-1)*log(x) + (b-1)*log(1.0-x) - loggamma(a) - loggamma(b) + loggamma(a+b)
    logp(b) = sum(logp.(samples_out,samples_a,b))/n_samples

    return Message(Univariate, Function, log_pdf = logp)
end

ruleVBBetaOut(marg_out::Any,
              dist_a::ProbabilityDistribution{Univariate},
              dist_b::ProbabilityDistribution{Univariate}) =
    Message(Univariate, Beta, a=unsafeMean(dist_a), b=unsafeMean(dist_b))

function ruleVBBetaA(dist_out::ProbabilityDistribution{Univariate},
                     dist_a::Any,
                     dist_b::ProbabilityDistribution{Univariate})

    logGamAsumB(a) = sum(loggamma.(a.+sample(dist_b,default_n_samples)))/default_n_samples
    logp(a) = (a-1)*unsafeLogMean(dist_out) + logGamAsumB(a) - loggamma(a)

    return Message(Univariate, Function, log_pdf = logp)
end

function ruleVBBetaB(dist_out::ProbabilityDistribution{Univariate},
                     dist_a::ProbabilityDistribution{Univariate}, dist_b::Any)

    logGamAsumB(b) = sum(loggamma.(b.+sample(dist_a,default_n_samples)))/default_n_samples
    logp(b) = (b-1)*unsafeMirroredLogMean(dist_out) + logGamAsumB(b) - loggamma(b)

    return Message(Univariate, Function, log_pdf = logp)
end


text:
function recurringA0(M2::Int64,m1::Int64,m2::Int64,m3::Int64,m4::Int64)
    # return (M-m1-m2-m3-m4+4)/(2M-2m1-2m2-2m3-2m4+8)/(2M-2m1-2m2-2m3-2m4+7)
    # return BigFloat(1)/BigFloat(M2-2m1-2m2-2m3-2m4+7) #BigFloat(1/(M2-2m1-2m2-2m3-2m4+7))
    return 1//BigInt(M2-2m1-2m2-2m3-2m4+7)
end

function recurringA1(n1::Int64,n2::Int64,n3::Int64,n4::Int64,m1::Int64,m2::Int64,m3::Int64,m4::Int64)

    # M = Int64((n1 + n2 + n3 + n4)/2) # integer

    if m1 == m2 == m3 == 0
       # coeff = (-1)*BigFloat(n4-2m4+2)*BigFloat(n4-2m4+1)/BigFloat(m4)*recurringA0(n1+n2+n3+n4,1,1,1,m4)
       coeff = (-1)*BigInt(n4-2m4+2)*BigInt(n4-2m4+1)//BigInt(m4)*recurringA0(n1+n2+n3+n4,1,1,1,m4)
    elseif m1 == m2 == 0
       # coeff = (-1)*BigFloat(n3-2m3+2)*BigFloat(n3-2m3+1)/BigFloat(m3)*recurringA0(n1+n2+n3+n4,1,1,m3,m4)
       coeff = (-1)*BigInt(n3-2m3+2)*BigInt(n3-2m3+1)//BigInt(m3)*recurringA0(n1+n2+n3+n4,1,1,m3,m4)
    elseif m1 == 0
       # coeff = (-1)*BigFloat(n2-2m2+2)*BigFloat(n2-2m2+1)/BigFloat(m2)*recurringA0(n1+n2+n3+n4,1,m2,m3,m4)
       coeff = (-1)*BigInt(n2-2m2+2)*BigInt(n2-2m2+1)//BigInt(m2)*recurringA0(n1+n2+n3+n4,1,m2,m3,m4)
    else
       # coeff = (-1)*BigFloat(n1-2m1+2)*BigFloat(n1-2m1+1)/BigFloat(m1)*recurringA0(n1+n2+n3+n4,m1,m2,m3,m4)
       coeff = (-1)*BigInt(n1-2m1+2)*BigInt(n1-2m1+1)//BigInt(m1)*recurringA0(n1+n2+n3+n4,m1,m2,m3,m4)
    end

    return coeff

end

function Vijkl2(n1::Int64,n2::Int64,n3::Int64,n4::Int64)

    if isodd(n1 + n2 + n3 + n4)
       return 0.0
    end

    M = Int64((n1 + n2 + n3 + n4)/2)
    maxm1 = floor(Int64,n1/2)+1
    maxm2 = floor(Int64,n2/2)+1
    maxm3 = floor(Int64,n3/2)+1
    maxm4 = floor(Int64,n4/2)+1
    matA = zeros(Rational{BigInt},2,maxm2,maxm3,maxm4) #zeros(BigFloat,2,maxm2,maxm3,maxm4) #zeros(Float64,2,maxm2,maxm3,maxm4)

    # m1=m2=m3=m4=1
    matA[1,1,1,1] = Rational(BigInt(1)) #BigFloat(1.0)
    # matA[1,1,1,1] = exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))

    # m1=m2=m3=1
    for m4 = 1:maxm4-1
        # matA[1,1,1,m4+1] = matA[1,1,1,m4]*(-1)*2*(n4-2m4+2)*(n4-2m4+1)/m4*(M-m4+1)/(2M-2m4+2)/(2M-2m4+1)
        matA[1,1,1,m4+1] = matA[1,1,1,m4]*recurringA1(n1,n2,n3,n4,0,0,0,m4)
    end

    # m1=m2=1
    for m3 = 1:maxm3-1
        # Threads.@threads for m4 = 1:maxm4 # parfor
        for m4 = 1:maxm4
            # matA[1,1,m3+1,m4] = matA[1,1,m3,m4]*(-1)*2*(n3-2m3+2)*(n3-2m3+1)/m3*(M-m3-m4+2)/(2M-2m3-2m4+4)/(2M-2m3-2m4+3)
            matA[1,1,m3+1,m4] = matA[1,1,m3,m4]*recurringA1(n1,n2,n3,n4,0,0,m3,m4)
        end
    end

    # m1=1
    for m2 = 1:maxm2-1
        # Threads.@threads for m34 = 1:maxm3*maxm4 # parfor
        for m34 = 1:maxm3*maxm4
            m3 = div(m34,maxm4)+1
            m4 = mod(m34,maxm4) # m34-div(m34,maxm4)*maxm4
            if m4 == 0
               m3 = m3-1 # div(m34,maxm4)
               m4 = maxm4 # m34-(div(m34,maxm4)-1)*maxm4 = m34-(m3-1)*maxm4 = maxm4
            end
            # matA[1,m2+1,m3,m4] = matA[1,m2,m3,m4]*(-1)*2*(n2-2m2+2)*(n2-2m2+1)/m2*(M-m2-m3-m4+3)/(2M-2m2-2m3-2m4+6)/(2M-2m2-2m3-2m4+5)
            matA[1,m2+1,m3,m4] = matA[1,m2,m3,m4]*recurringA1(n1,n2,n3,n4,0,m2,m3,m4)
        end
    end

    sumA = sum(matA[1,:,:,:])

    for m1 = 1:maxm1-1

        # Threads.@threads for m234 = 1:maxm2*maxm3*maxm4 # parfor
        for m234 = 1:maxm2*maxm3*maxm4
            m2 = div(m234,maxm3*maxm4)+1
            m3 = div(m234 - (m2-1)*maxm3*maxm4,maxm4)+1
            m4 = m234 - (m2-1)*maxm3*maxm4 - (m3-1)*maxm4
            if mod(m234,maxm3*maxm4) == 0
               m2 = m2-1
               m3 = maxm3
               m4 = maxm4
            elseif mod(m234 - (m2-1)*maxm3*maxm4,maxm4) == 0
               m3 = m3-1
               m4 = maxm4
            end
            # matA[2,m2,m3,m4] = matA[1,m2,m3,m4]*(-1)*2*(n1-2m1+2)*(n1-2m1+1)/m1*(M-m1-m2-m3-m4+4)/(2M-2m1-2m2-2m3-2m4+8)/(2M-2m1-2m2-2m3-2m4+7)
            matA[2,m2,m3,m4] = matA[1,m2,m3,m4]*recurringA1(n1,n2,n3,n4,m1,m2,m3,m4)
        end

        sumA = sumA + sum(matA[2,:,:,:])
        matA[1,:,:,:] .= matA[2,:,:,:]
        matA[2,:,:,:] .= Rational(BigInt(0)) #0

    end

    # sumA = sumA*exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))
    sumA = Float64(sumA)*exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))

    # if sumA < 0
    #    sumA = -sumA
    #    sumA = log(sumA) -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
    #    sumA = exp(sumA)
    #    sumA = -sumA
    # else
       # sumA = log(sumA) -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
       # sumA = exp(sumA)
    # end

    sumA = sumA/sqrt(2*pi)/2

    return sumA
    # return matA

end

function A0000(n1::Int64,n2::Int64,n3::Int64,n4::Int64)
   M = Int64((n1 + n2 + n3 + n4)/2)
   A0 = -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
   return exp(A0)/sqrt(2*pi)/2
end

ext: jl

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
#####
##### Duel.jl
##### Simple script to run custom benchmark duels
#####

using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "connect-four")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

SESSION_DIR = "sessions/$GAME"

AlphaZero.UserInterface.run_duel(Game, Training.Network{Game}, SESSION_DIR,
  Benchmark.Duel(
    Benchmark.NetworkOnly(),
    Benchmark.MinMaxTS(depth=5, τ=0.2),
    num_games=200,
    flip_probability=0.5,
    color_policy=CONTENDER_WHITE),
  params=Training.params)


text:
module Components
import JSON2
export Component, <|, hasproperty, is_prop_available, validate

abstract type AbstractComponent end

struct Component <: AbstractComponent
    type ::String
    namespace ::String
    props ::Dict{Symbol, Any}
    available_props ::Set{Symbol}
    wildcard_props ::Set{Symbol}
end

JSON2.@format Component begin
    available_props => (exclude = true,)
    wildcard_props => (exclude = true,)
end

function <|(comp::Component, value::Any)
    comp.props["children"] = value
    return comp
end

hasproperty(c::Component, prop::Symbol) = haskey(c.props, prop)
function is_prop_available(c::Component, prop::Symbol)
    if length(c.wildcard_props) > 0
        wild_regs = Regex("^(?<prop>$(join(c.wildcard_props, "|")))")
        if !isnothing(match(wild_regs, string(prop)))
            return true
        end
    end

    return prop in c.available_props
end

function _validate_children(children::Union{Vector, Tuple}, ids::Set{Symbol})
    for child in children
        _validate(child, ids)
    end
end
function _validate_children(children::Component, ids::Set{Symbol})
    _validate(children, ids)
end
function _validate_children(children, ids::Set{Symbol}) end

function _validate(comp::Component, ids::Set{Symbol})
    if haskey(comp.props, :id)
        id = Symbol(comp.props[:id])
        id in ids && error("Duplicate component id found in the initial layout: $(id)")
        push!(ids, id)
    end
    _validate_children(get(comp.props, :children, nothing), ids)
end
function _validate(non_comp, ids::Set{Symbol}) end

function validate(comp::Component)
    _validate(comp, Set{Symbol}())
end

end


text:
using Test
using MLJBase
include(joinpath("..", "..", "test", "_models", "models.jl"))
using .Models

@testset "_categorical" begin
    a = [1, 1, 2, 3]
    b = [3, 3, 4, 5]
    c = [missing, a...]
    d = [missing, b...]
    e = categorical(a)
    f = categorical(b)
    g = categorical(c)
    h = categorical(d)
    rhs = (Set(1:5), Set(1:5))
    @test Set.(levels.(MLJBase._categorical(a, b))) == rhs
    @test Set.(levels.(MLJBase._categorical(a, d))) == rhs
    @test Set.(levels.(MLJBase._categorical(c, b))) == rhs
    @test Set.(levels.(MLJBase._categorical(c, d))) == rhs
    @test Set.(levels.(MLJBase._categorical(a, f))) == rhs
    @test Set.(levels.(MLJBase._categorical(a, h))) == rhs
    @test Set.(levels.(MLJBase._categorical(b, a))) == rhs
    @test Set.(levels.(MLJBase._categorical(d, a))) == rhs
    @test Set.(levels.(MLJBase._categorical(b, c))) == rhs
    @test Set.(levels.(MLJBase._categorical(d, c))) == rhs
    @test Set.(levels.(MLJBase._categorical(f, a))) == rhs
    @test Set.(levels.(MLJBase._categorical(h, a))) == rhs
end

@testset "basics" begin
    yraw = ['m',     'm', 'f', 'n', missing, 'f', 'm', 'n', 'n', 'm', 'f']
    ŷraw = [missing, 'f', 'f', 'm', 'f',     'f', 'n', 'm', 'n', 'm', 'f']
    y = categorical(yraw)
    ŷ = categorical(ŷraw)
    l = levels(y) # f, m, n
    cm = MLJBase._confmat(ŷ, y; warn=false)
    ŷ_clean, y_clean = MLJBase.skipinvalid(ŷ, y)
    ee(l,i,j) = sum((ŷ_clean .== l[i]) .& (y_clean .== l[j]))
    for i in 1:3, j in 1:3
        @test cm[i,j] == ee(l,i,j)
    end

    cm2 = @test_logs (:warn, r"The classes are") MLJBase._confmat(ŷraw, yraw)
    @test cm2.mat == cm.mat

    perm = [3, 1, 2]
    l2 = l[perm]
    cm2 = @test_logs MLJBase._confmat(ŷ, y; perm=perm)
    m = ConfusionMatrix(perm=perm)
    for i in 1:3, j in 1:3
        @test cm2[i,j] == ee(l2,i,j)
    end
    @test_logs (:warn, r"The classes are un") MLJBase._confmat(ŷ, y)
    ŷc = coerce(ŷ, Union{Missing,OrderedFactor})
    yc = coerce(y, Union{Missing,OrderedFactor})
    @test MLJBase._confmat(ŷc, yc).mat == cm.mat

    y = categorical(['a','b','a','b'])
    ŷ = categorical(['b','b','a','a'])
    @test_logs (:warn, r"The classes are un") MLJBase._confmat(ŷ, y)

    # more tests for coverage
    y = categorical([1,2,3,1,2,3,1,2,3])
    ŷ = categorical([1,2,3,1,2,3,1,2,3])
    @test_throws ArgumentError MLJBase._confmat(ŷ, y, rev=true)

    # silly test for display
    ŷ = coerce(y, OrderedFactor)
    y = coerce(y, OrderedFactor)
    iob = IOBuffer()
    Base.show(iob, MIME("text/plain"), MLJBase._confmat(ŷ, y))
    siob = String(take!(iob))
    @test strip(siob) == strip("""
                         ┌─────────────────────────────────────────┐
                         │              Ground Truth               │
           ┌─────────────┼─────────────┬─────────────┬─────────────┤
           │  Predicted  │      1      │      2      │      3      │
           ├─────────────┼─────────────┼─────────────┼─────────────┤
           │      1      │      3      │      0      │      0      │
           ├─────────────┼─────────────┼─────────────┼─────────────┤
           │      2      │      0      │      3      │      0      │
           ├─────────────┼─────────────┼─────────────┼─────────────┤
           │      3      │      0      │      0      │      3      │
           └─────────────┴─────────────┴─────────────┴─────────────┘""")
end

@testset "ConfusionMatrix measure" begin

    @test info(confmat).orientation == :other
    model = DeterministicConstantClassifier()

    X = (x=rand(10),)
    long = categorical(collect("abbaacaabbbbababcbac"), ordered=true)
    y = long[1:10]
    yhat =long[11:20]

    @test confmat(yhat, y).mat == [1 2 0; 3 1 1; 1 1 0]
    @test ConfusionMatrix(perm=[2, 1, 3])(yhat, y).mat ==
        MLJBase._confmat(yhat, y, perm=[2, 1, 3]).mat

    MLJBase.value(confmat, yhat, X, y, nothing)

    e = evaluate(model, X, y,
                 measures=[misclassification_rate, confmat],
                 resampling=Holdout(fraction_train=0.5))
    cm = e.measurement[2]
    @test cm.labels == ["a", "b", "c"]
    @test cm.mat == [2 2 1; 0 0 0; 0 0 0]
end


text:
using Plots

dt = 0.01
g = 9.80665

function freefall(x_0, v_0)
    x = x_0
    v = v_0
    t = 0
    xs = []
    ts = []
    while (x > 0)
        v += g*dt
        x -= v*dt
        t += dt
        push!(xs, x)
        push!(ts, t)
    end

    (xs, ts)
end

xs_0, ts_0 = freefall(100, 0)
xs_100, ts_100 = freefall(100, -100)

plot(ts_0, xs_0, label="v_0=0",xlabel="t",ylabel="x", color=:cornflowerblue)
savefig("freefall_v0=0.png")
plot(ts_100, xs_100, label="v_0=-100",xlabel="t",ylabel="x", color=:coral)
savefig("freefall_v0=-100.png")


text:
module Compress

import Pkg.TOML
import Pkg.Types: VersionSpec, compress_versions

function load_versions(path::String)
    versions_file = joinpath(dirname(path), "Versions.toml")
    versions_dict = TOML.parsefile(versions_file)
    sort!([VersionNumber(v) for v in keys(versions_dict)])
end

function load(path::String,
    versions::Vector{VersionNumber} = load_versions(path))
    compressed = TOML.parsefile(path)
    uncompressed = Dict{VersionNumber,Dict{Any,Any}}()
    for (vers, data) in compressed
        vs = VersionSpec(vers)
        for v in versions
            v in vs || continue
            merge!(get!(uncompressed, v, Dict()), deepcopy(data))
        end
    end
    return uncompressed
end

function compress(path::String, uncompressed::Dict,
    versions::Vector{VersionNumber} = load_versions(path))
    inverted = Dict()
    for (ver, data) in uncompressed, (key, val) in data
        val isa TOML.TYPE || (val = string(val))
        push!(get!(inverted, key => val, VersionNumber[]), ver)
    end
    compressed = Dict()
    for ((k, v), vers) in inverted
        for r in compress_versions(versions, sort!(vers)).ranges
            get!(compressed, string(r), Dict{String,Any}())[k] = v
        end
    end
    return compressed
end

function save(path::String, uncompressed::Dict,
    versions::Vector{VersionNumber} = load_versions(path))
    compressed = compress(path, uncompressed)
    open(path, write=true) do io
        TOML.print(io, compressed, sorted=true)
    end
end

end # module


text:
# This file is part of InvertedFiles.jl

using InvertedFiles, SimilaritySearch, LinearAlgebra
using Test
using Random
Random.seed!(0)

@testset "WeightedInvertedFile" begin
    A = [normalize!(rand(300)) for i in 1:1000]
    B = VectorDatabase([Dict(enumerate(a)) for a in A])

    # testing with Vector container (map=nothing)
    I = append!(WeightedInvertedFile(300), B)

    k = 30
    for i in 1:10
        qid = rand(1:length(A))
        a = search(ExhaustiveSearch(NormalizedCosineDistance(), A), A[qid], KnnResult(k))
        b = search(I, B[qid], KnnResult(k))
        @test recallscore(a.res, b.res) == 1.0
    end

    # testing with Dict container (map != nothing)
    I = append!(WeightedInvertedFile(300), B)

    k = 30
    for i in 1:10
        qid = rand(1:length(A))
        a = search(ExhaustiveSearch(NormalizedCosineDistance(), A), A[qid], KnnResult(k))
        b = search(I, B[qid], KnnResult(k))
        @test recallscore(a.res, b.res) == 1.0
    end

    k = 30
    for i in 1:10
        @info i
        qid = rand(1:length(A))
        @time a = search(ExhaustiveSearch(NormalizedCosineDistance(), A), A[qid], KnnResult(k))
        @time b = search(I, B[qid], KnnResult(k))
        @test recallscore(a.res, b.res) == 1.0
    end

    ## working on sparse data
    # increasing sparsity of the arrays
    for A_ in A
        t = partialsort(A_, 7, rev=true)
        for i in eachindex(A_)
            A_[i] = A_[i] < t ? 0.0 : A_[i]
        end
        normalize!(A_)
    end

    create_sparse(A_) = Dict([i => a for (i, a) in enumerate(A_) if a > 0.0])

    #B = VectorDatabase([create_sparse(A_) for A_ in A])
    B = VectorDatabase(A)
    I = append!(WeightedInvertedFile(300), B)
    k = 1  # the aggresive cut of the attributes need a small k
    for i in 1:10
        @info i
        qid = rand(1:length(A))
        @time a = search(ExhaustiveSearch(NormalizedCosineDistance(), A), A[qid], KnnResult(k))
        @time b = search(I, B[qid], KnnResult(k))
        @test recallscore(a.res, b.res) == 1.0
        @show recallscore(a.res, b.res)
    end
end

@testset "BinaryInvertedFile" begin
    vocsize = 128
    n = 10_000
    m = 100
    len = 10
    k = 10
    db = VectorDatabase([sort!(unique(rand(1:vocsize, len))) for i in 1:n])
    queries = VectorDatabase([sort!(unique(rand(1:vocsize, len))) for i in 1:m])

    #for dist in [JaccardDistance(), DiceDistance(), CosineDistanceSet(), IntersectionDissimilarity()]
    for dist in [JaccardDistance()]
        S = ExhaustiveSearch(; dist, db)
        gI, gD = searchbatch(S, queries, k)

        IF = BinaryInvertedFile(vocsize, dist)
        append!(IF, db)
        iI, iD = searchbatch(IF, queries, k)
        @time search(IF, queries[1], SimilaritySearch.getknnresult(k))
        @time search(IF, queries[2], SimilaritySearch.getknnresult(k))

        recall = macrorecall(gI, iI)
        @show dist, recall
        @test recall > 0.95  # sets can be tricky since we can expect many similar distances
        err = 0.0
        for i in 1:m
            d = evaluate(L2Distance(), gD[:, i], iD[:, i])
            err += d
            if d > 0.1
                @info dist, i, gD[:, i], iD[:, i]
                @info dist, i, gI[:, i], iI[:, i]
                @info dist, i, queries[i]
            end
        end
        @show dist, err
        @test err < 0.01  # acc. floating point errors
    end
end


text:
import Base.@deprecate_binding

function floattype(::Type{T}) where {T <: Real}
    Base.depwarn("""
        In a future release, the fallback definition of `floattype` will throw a MethodError if it cannot return a type `<:AbstractFloat`.
        See the documentation on `floattype` for guidance on whether to define a custom `floattype(::Type{$T})` method.
        """, :floattype)
    return T
end


text:
module Classification

import Plots
import Statistics

include("confusion_matrix.jl"); export confusion_params, confusion_matrix, class_confusion
include("metrics.jl"); export visualize
include("visualization.jl"); export classification_report, condition_positive, condition_negative, predicted_positive,predicted_negative, correctly_classified, incorrectly_classified, sensitivity_score, recall_score, specificity_score, precision_score, positive_predictive_value, accuracy_score, balanced_accuracy_score, negative_predictive_value, false_negative_rate, false_positive_rate, false_discovery_rate, false_omission_rate, f1_score, prevalence_threshold, threat_score, matthews_correlation_coeff, fowlkes_mallows_index, informedness, markedness, cohen_kappa_score, hamming_loss, jaccard_score, confusion_params

end


text:
function _hash_file(fn::String)
    h0 = hash("")
    !isfile(fn) && return h0
    for line in eachline(fn)
        h0 = hash(line, h0)
    end
    return h0
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function readconfig(app, mode)

filename = app.configfile;

if mode==0 # binary
    tmp = reinterpret(Float64,read(filename));
else # text
    tmp = Main.DelimitedFiles.readdlm(filename);
end
#tmp = tmp[:];

if mode==4
    config = readLAMMPS(tmp);
    config.nconfigs = length(config.natom); # number of configurations
    config.natomall = sum(config.natom);# number of atoms for all configurations
    return config;
end

if mode<=2
    nconfigs = int64(tmp[1]);          # number of configurations
    dim = int64(tmp[2]);
    ncq = int64(tmp[3]);
    ncv = int64(tmp[4]);
    nce = int64(tmp[5]);
    ncf = int64(tmp[6]);
    ncx = dim;

    app.nconfigs = nconfigs; # number of configurations
    app.dim = dim; # physical dimension
    app.ncx = dim; # number of compoments of x
    app.ncv = ncv; # number of compoments of v
    app.ncf = ncf; # number of compoments of f
    app.ncq = ncq; # number of compoments of q
    app.nce = nce; # number of compoments of e
    config = initializeconfig(app);
end

if mode<=1
    n2 = 6;
    m  = (1+nce+dim*dim)*nconfigs;
    n1 = n2+1;
    n2 = n2+m;
    tm = reshape(tmp[n1:n2], nconfigs, (1+nce+dim*dim));
    config.natom = tm[:,1]; # number of atoms per each configuration
    config.a = tm[:,2:(1+dim)]'; # 1st principal vectors per each configuration
    config.b = tm[:,(2+dim):(1+2*dim)]'; # 2nd principal vectors per each configuration
    if dim==3
        config.c = tm[:,(2+2*dim):(1+dim*dim)]'; # 2nd principal vectors per each configuration
    else
        config.c = [];
    end
    if nce>0
        config.e = tm[:,2+dim*dim]; # potential energy per each configuration
    end

    config.natomall = sum(config.natom);# number of atoms for all configurations
    m  = (1+ncx+ncq+ncv+ncf)*config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    tm = reshape(tmp[n1:n2], config.natomall, (1+ncx+ncq+ncv+ncf));
    config.t = tm[:,1];
    config.x = tm[:,2:(1+ncx)]';
    config.q = tm[:,(2+ncx):(1+ncx+ncq)]';
    config.v = tm[:,(2+ncx+ncq):(1+ncx+ncq+ncv)]';
    config.f = tm[:,(2+ncx+ncq+ncv):(1+ncx+ncq+ncv+ncf)]';

elseif mode==2
    n2 = 6;
    m  = nconfigs;
    n1 = n2+1;
    n2 = n2+m;
    config.natom = reshape(tmp[n1:n2], 1, m); # number of atoms per each configuration
    config.natomall = sum(config.natom);# number of atoms for all configurations

    # the 1st principal vector of the simulation box
    m  = dim*nconfigs;
    n1 = n2+1;
    n2 = n2+m;
    config.a = reshape(tmp[n1:n2], dim, nconfigs); # the 1st principal vector of the simulation box

    # the 2nd principal vector of the simulation box
    n1 = n2+1;
    n2 = n2+m;
    config.b = reshape(tmp[n1:n2], dim, nconfigs); # the 2nd principal vector of the simulation box

    if dim==3
        # the 3rd principal vector of the simulation box
        n1 = n2+1;
        n2 = n2+m;
        config.c = reshape(tmp[n1:n2], dim, nconfigs); # the 3rd principal vector of the simulation box
    end

    m  = nconfigs;
    n1 = n2+1;
    n2 = n2+m;
    config.e = reshape(tmp[n1:n2], 1, m); # number of atoms per each configuration

    # read atom types for all configurations
    m  = config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    config.t = reshape(int64(tmp[n1:n2]), 1, m);

    # read atom positions for all configurations
    m  = ncx*config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    config.x = reshape(tmp[n1:n2], ncx, config.natomall);

    # read atom charges for all configurations
    m  = ncq*config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    config.q = reshape(tmp[n1:n2], ncq, config.natomall);

    # read atom velocities for all configurations
    m  = ncv*config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    config.v = reshape(tmp[n1:n2], ncv, config.natomall);

    # read atom forces for all configurations
    m  = ncf*config.natomall;
    n1 = n2+1;
    n2 = n2+m;
    config.f = reshape(tmp[n1:n2], ncf, config.natomall);
end

return config

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
EzXML.parsexml(x::ReadableFile) = parsexml(readstring(x))
EzXML.parsexml(x::Reader) = parsexml(readstring(x))

"""
    ScientificFormat

Scientific notation numbers in Excel.
If the format code is `0.00E+00` and the cell value is 12,200,000, then 1.22E+07 is
stored in the xml.
"""
struct ScientificFormat end
function Base.parse(::Type{ScientificFormat}, x::AbstractString)::Float64
    if contains(x, "E")
        x = lowercase(x)
    end
    parse(Float64, x)
end

"""
    GeneralFormat

The primary goal when a cell is using `GeneralFormat` is to render the cell content without user-specified
guidance to the best ability of the application.
It can be `String`, `Boolean`, `Error`, `Numbers`
`DateTime` do not fall into "GeneralFormat`
"""
struct GeneralFormat end
function Base.parse(::Type{GeneralFormat}, x)
    function isfloat(x)::Bool
        b = false
        if length(x) > 1
            if contains(x, ".")
                b = length(matchall(r"[.]", x)) < 2 ? true : false
            end
            if startswith(x, "-")
                b = length(matchall(r"[-]", x)) < 2 ? true : false
            end
        end
        b
    end
    function isint(x::AbstractString)::Bool
        b = false
        if !isempty(x)
            cha = matchall(r"[^0-9]", x)
            b = isempty(cha)                         ? true :
                (length(cha) == 1 && cha[1] == "-" ) ? true :
                                                     false
        end
        b
    end
    if isfloat(x)
        parse(Float64, x)
    else
        in(x, keys(EzExcel.ECMA_ERROR_VALUES))    ? missing :
        in(x, ("True", "true", "False", "false")) ? parse(Bool, lowercase(x)) :
        isint(x)                                  ? parse(Int, x) :
                                                  String(x)
    end
end

"""
    ExcelDateTime

Parse an Excel number (presumed to represent a date, a datetime or a time) to a DateTime

xldate: The Excel number
datemode: true:1904-based, false: 1900-based.

    `1904-01-01` is not regarded as a valid date in the `datemode == true`
     its "serial number" is zero.
 """
struct ExcelDateTime end
function Base.parse(::Type{ExcelDateTime}, xldate::AbstractString; args...)
    parse(ExcelDateTime, parse(ScientificFormat, xldate); args...)
end
function Base.parse(::Type{ExcelDateTime}, xldate::Real; datemode::Bool=false)::DateTime
    # Set the epoch based on the 1900/1904 datemode.
    epoch = if datemode
        DateTime(1904, 1, 1)
    else # Workaround Excel 1900 leap year bug by adjusting the epoch.
        xldate < 60 ? DateTime(1899, 12, 31) : DateTime(1899, 12, 30)
    end
    xldays, fraction = divrem(xldate, 1)
    xldays == 0 && (xldays = 1)

    sec = round(Int, fraction * 86400000.0)
    sec, milli_sec = divrem(sec, 1000)

    return epoch + Dates.Day(xldays) + Dates.Second(sec) + Dates.Millisecond(milli_sec)
end

struct SharedString end


text:
using LinearAlgebra: checksquare
using LinearAlgebra.BLAS: gemv, gemv!, gemm!, trsm!, axpy!, ger!

#####
##### `svd`
#####

function rrule(::typeof(svd), X::AbstractMatrix{<:Real})
    F = svd(X)
    function svd_pullback(Ȳ::Composite)
        # `getproperty` on `Composite`s ensures we have no thunks.
        ∂X = svd_rev(F, Ȳ.U, Ȳ.S, Ȳ.V)
        return (NO_FIELDS, ∂X)
    end
    return F, svd_pullback
end

function rrule(::typeof(getproperty), F::T, x::Symbol) where T <: SVD
    function getproperty_svd_pullback(Ȳ)
        C = Composite{T}
        ∂F = if x === :U
            C(U=Ȳ,)
        elseif x === :S
            C(S=Ȳ,)
        elseif x === :V
            C(V=Ȳ,)
        elseif x === :Vt
            # TODO: https://github.com/JuliaDiff/ChainRules.jl/issues/106
            throw(ArgumentError("Vt is unsupported; use V and transpose the result"))
        end
        return NO_FIELDS, ∂F, DoesNotExist()
    end
    return getproperty(F, x), getproperty_svd_pullback
end

# When not `Zero`s expect `Ū::AbstractMatrix, s̄::AbstractVector, V̄::AbstractMatrix`
function svd_rev(USV::SVD, Ū, s̄, V̄)
    # Note: assuming a thin factorization, i.e. svd(A, full=false), which is the default
    U = USV.U
    s = USV.S
    V = USV.V
    Vt = USV.Vt

    k = length(s)
    T = eltype(s)
    F = T[i == j ? 1 : inv(@inbounds s[j]^2 - s[i]^2) for i = 1:k, j = 1:k]

    # We do a lot of matrix operations here, so we'll try to be memory-friendly and do
    # as many of the computations in-place as possible. Benchmarking shows that the in-
    # place functions here are significantly faster than their out-of-place, naively
    # implemented counterparts, and allocate no additional memory.
    Ut = U'
    FUᵀŪ = _mulsubtrans!!(Ut*Ū, F)  # F .* (UᵀŪ - ŪᵀU)
    FVᵀV̄ = _mulsubtrans!!(Vt*V̄, F)  # F .* (VᵀV̄ - V̄ᵀV)
    ImUUᵀ = _eyesubx!(U*Ut)  # I - UUᵀ
    ImVVᵀ = _eyesubx!(V*Vt)  # I - VVᵀ

    S = Diagonal(s)
    S̄ = s̄ isa AbstractZero ? s̄ : Diagonal(s̄)

    # TODO: consider using MuladdMacro here
    Ā = add!!(U * FUᵀŪ * S, ImUUᵀ * (Ū / S)) * Vt
    Ā = add!!(Ā, U * S̄ * Vt)
    Ā = add!!(Ā, U * add!!(S * FVᵀV̄ * Vt, (S \ V̄') * ImVVᵀ))

    return Ā
end

#####
##### `eigen`
#####

# TODO:
# - support correct differential of phase convention when A is hermitian
# - simplify when A is diagonal
# - support degenerate matrices (see #144)

function frule((_, ΔA), ::typeof(eigen!), A::StridedMatrix{T}; kwargs...) where {T<:BlasFloat}
    F = eigen!(A; kwargs...)
    ΔA isa AbstractZero && return F, ΔA
    λ, V = F.values, F.vectors
    tmp = V \ ΔA
    ∂K = tmp * V
    ∂Kdiag = @view ∂K[diagind(∂K)]
    ∂λ = eltype(λ) <: Real ? real.(∂Kdiag) : copy(∂Kdiag)
    ∂K ./= transpose(λ) .- λ
    fill!(∂Kdiag, 0)
    ∂V = mul!(tmp, V, ∂K)
    _eigen_norm_phase_fwd!(∂V, A, V)
    ∂F = Composite{typeof(F)}(values = ∂λ, vectors = ∂V)
    return F, ∂F
end

function rrule(::typeof(eigen), A::StridedMatrix{T}; kwargs...) where {T<:Union{Real,Complex}}
    F = eigen(A; kwargs...)
    function eigen_pullback(ΔF::Composite{<:Eigen})
        λ, V = F.values, F.vectors
        Δλ, ΔV = ΔF.values, ΔF.vectors
        if ΔV isa AbstractZero
            Δλ isa AbstractZero && return (NO_FIELDS, Δλ + ΔV)
            ∂K = Diagonal(Δλ)
            ∂A = V' \ ∂K * V'
        else
            ∂V = copyto!(similar(ΔV), ΔV)
            _eigen_norm_phase_rev!(∂V, A, V)
            ∂K = V' * ∂V
            ∂K ./= λ' .- conj.(λ)
            ∂K[diagind(∂K)] .= Δλ
            ∂A = mul!(∂K, V' \ ∂K, V')
        end
        return NO_FIELDS, T <: Real ? real(∂A) : ∂A
    end
    eigen_pullback(ΔF::AbstractZero) = (NO_FIELDS, ΔF)
    return F, eigen_pullback
end

# mutate ∂V to account for the (arbitrary but consistent) normalization and phase condition
# applied to the eigenvectors.
# these implementations assume the convention used by eigen in LinearAlgebra (i.e. that of
# LAPACK.geevx!; eigenvectors have unit norm, and the element with the largest absolute
# value is real), but they can be specialized for `A`

function _eigen_norm_phase_fwd!(∂V, A, V)
    @inbounds for i in axes(V, 2)
        v, ∂v = @views V[:, i], ∂V[:, i]
        # account for unit normalization
        ∂c_norm = -real(dot(v, ∂v))
        if eltype(V) <: Real
            ∂c = ∂c_norm
        else
            # account for rotation of largest element to real
            k = _findrealmaxabs2(v)
            ∂c_phase = -imag(∂v[k]) / real(v[k])
            ∂c = complex(∂c_norm, ∂c_phase)
        end
        ∂v .+= v .* ∂c
    end
    return ∂V
end

function _eigen_norm_phase_rev!(∂V, A, V)
    @inbounds for i in axes(V, 2)
        v, ∂v = @views V[:, i], ∂V[:, i]
        ∂c = dot(v, ∂v)
        # account for unit normalization
        ∂v .-= real(∂c) .* v
        if !(eltype(V) <: Real)
            # account for rotation of largest element to real
            k = _findrealmaxabs2(v)
            @inbounds ∂v[k] -= im * (imag(∂c) / real(v[k]))
        end
    end
    return ∂V
end

# workaround for findmax not taking a mapped function
function _findrealmaxabs2(x)
    amax = abs2(first(x))
    imax = 1
    @inbounds for i in 2:length(x)
        xi = x[i]
        !isreal(xi) && continue
        a = abs2(xi)
        a < amax && continue
        amax, imax = a, i
    end
    return imax
end

#####
##### `eigvals`
#####

function frule((_, ΔA), ::typeof(eigvals!), A::StridedMatrix{T}; kwargs...) where {T<:BlasFloat}
    ΔA isa AbstractZero && return eigvals!(A; kwargs...), ΔA
    F = eigen!(A; kwargs...)
    λ, V = F.values, F.vectors
    tmp = V \ ΔA
    ∂λ = similar(λ)
    # diag(tmp * V) without computing full matrix product
    if eltype(∂λ) <: Real
        broadcast!((a, b) -> sum(real ∘ prod, zip(a, b)), ∂λ, eachrow(tmp), eachcol(V))
    else
        broadcast!((a, b) -> sum(prod, zip(a, b)), ∂λ, eachrow(tmp), eachcol(V))
    end
    return λ, ∂λ
end

function rrule(::typeof(eigvals), A::StridedMatrix{T}; kwargs...) where {T<:Union{Real,Complex}}
    F = eigen(A; kwargs...)
    λ = F.values
    function eigvals_pullback(Δλ)
        V = F.vectors
        ∂A = V' \ Diagonal(Δλ) * V'
        return NO_FIELDS, T <: Real ? real(∂A) : ∂A
    end
    eigvals_pullback(Δλ::AbstractZero) = (NO_FIELDS, Δλ)
    return λ, eigvals_pullback
end

#####
##### `cholesky`
#####

function rrule(::typeof(cholesky), A::Real, uplo::Symbol=:U)
    C = cholesky(A, uplo)
    function cholesky_pullback(ΔC::Composite)
        return NO_FIELDS, ΔC.factors[1, 1] / (2 * C.U[1, 1]), DoesNotExist()
    end
    return C, cholesky_pullback
end

function rrule(::typeof(cholesky), A::Diagonal{<:Real}, ::Val{false}; check::Bool=true)
    C = cholesky(A, Val(false); check=check)
    function cholesky_pullback(ΔC::Composite)
        Ā = Diagonal(diag(ΔC.factors) .* inv.(2 .* C.factors.diag))
        return NO_FIELDS, Ā, DoesNotExist()
    end
    return C, cholesky_pullback
end

# The appropriate cotangent is different depending upon whether A is Symmetric / Hermitian,
# or just a StridedMatrix.
# Implementation due to Seeger, Matthias, et al. "Auto-differentiating linear algebra."
function rrule(
    ::typeof(cholesky),
    A::LinearAlgebra.HermOrSym{<:LinearAlgebra.BlasReal, <:StridedMatrix},
    ::Val{false};
    check::Bool=true,
)
    C = cholesky(A, Val(false); check=check)
    function cholesky_pullback(ΔC::Composite)
        Ā, U = _cholesky_pullback_shared_code(C, ΔC)
        Ā = BLAS.trsm!('R', 'U', 'C', 'N', one(eltype(Ā)) / 2, U.data, Ā)
        return NO_FIELDS, _symhermtype(A)(Ā), DoesNotExist()
    end
    return C, cholesky_pullback
end

function rrule(
    ::typeof(cholesky),
    A::StridedMatrix{<:LinearAlgebra.BlasReal},
    ::Val{false};
    check::Bool=true,
)
    C = cholesky(A, Val(false); check=check)
    function cholesky_pullback(ΔC::Composite)
        Ā, U = _cholesky_pullback_shared_code(C, ΔC)
        Ā = BLAS.trsm!('R', 'U', 'C', 'N', one(eltype(Ā)), U.data, Ā)
        idx = diagind(Ā)
        @views Ā[idx] .= real.(Ā[idx]) ./ 2
        return (NO_FIELDS, UpperTriangular(Ā), DoesNotExist())
    end
    return C, cholesky_pullback
end

function _cholesky_pullback_shared_code(C, ΔC)
    U = C.U
    Ū = ΔC.U
    Ā = similar(U.data)
    Ā = mul!(Ā, Ū, U')
    Ā = LinearAlgebra.copytri!(Ā, 'U', true)
    Ā = ldiv!(U, Ā)
    return Ā, U
end

function rrule(::typeof(getproperty), F::T, x::Symbol) where {T <: Cholesky}
    function getproperty_cholesky_pullback(Ȳ)
        C = Composite{T}
        ∂F = if x === :U
            if F.uplo === 'U'
                C(U=UpperTriangular(Ȳ),)
            else
                C(L=LowerTriangular(Ȳ'),)
            end
        elseif x === :L
            if F.uplo === 'L'
                C(L=LowerTriangular(Ȳ),)
            else
                C(U=UpperTriangular(Ȳ'),)
            end
        end
        return NO_FIELDS, ∂F, DoesNotExist()
    end
    return getproperty(F, x), getproperty_cholesky_pullback
end


text:
abstract type JointDist{𝛭} end
struct JointPosterior{P, 𝛭} <: JointDist{𝛭}
    M::𝛭
    Θ::Vector{P}
    density::Vector{Float64}
    μ_hat::Vector{Float64}
    U::Matrix{Float64}
end
struct JointPosteriorRaw{G, 𝛭} <: JointDist{𝛭}
    M::𝛭
    grid::G
    μ_hat::Vector{Float64}
    U::Matrix{Float64}
end
function try_chol!(U::AbstractArray{<:Real,2}, Σ::AbstractArray{<:Real,2})
  @inbounds for i ∈ 1:size(U,1)
    U[i,i] = Σ[i,i]
    for j ∈ 1:i-1
      U[j,i] = Σ[j,i]
      for k ∈ 1:j-1
        U[j,i] -= U[k,i] * U[k,j]
      end
      U[j,i] /= U[j,j]
      U[i,i] -= U[j,i]^2
    end
    U[i,i] > 0 ? U[i,i] = √U[i,i] : return false
  end
  true
end
function chol!(U::AbstractArray{<:Real,2}, Σ::AbstractArray{<:Real,2})
  @inbounds for i ∈ 1:size(U,1)
    U[i,i] = Σ[i,i]
    for j ∈ 1:i-1
      U[j,i] = Σ[j,i]
      for k ∈ 1:j-1
        U[j,i] -= U[k,i] * U[k,j]
      end
      U[j,i] /= U[j,j]
      U[i,i] -= U[j,i]^2
    end
    U[i,i] = √U[i,i]
  end
end
function inv!(U_inverse::AbstractArray{<:Real,2}, U::AbstractArray{<:Real,2})
  @inbounds for i ∈ 1:size(U,1)
    U_inverse[i,i] = 1 / U[i,i]
    for j ∈ i+1:size(U,1)
      U_inverse[i,j] = U[i,j] * U_inverse[i,i]
      for k ∈ i+1:j-1
        U_inverse[i,j] += U[k,j] * U_inverse[i,k]
      end
      U_inverse[i,j] /= -U[j,j]
    end
  end
end
function inv!(U::AbstractArray{<:Real,2})
  @inbounds for i ∈ 1:size(U,1)
    U[i,i] = 1 / U[i,i]
    for j ∈ i+1:size(U,1)
      U[i,j] = U[i,j] * U[i,i]
      for k ∈ i+1:j-1
        U[i,j] += U[k,j] * U[i,k]
      end
      U[i,j] /= -U[j,j]
    end
  end
  true
end
function safe_inv_chol!(U, H)
  try_chol!(U, H) ? try_inv!(U) : false
end
function inv_chol!(U, H)
  chol!(U, H)
  inv!(U)
  U
end

function count(v::Vector, ::Type{LDR{g}}) where g
  @assert 0.0 < g < 1.0
  inv_v = 1 ./ v
  inadmissable = sum(inv_v .>= 1e11)
  ind_start = inadmissable + 1
  @inbounds for i ∈ ind_start:length(v)
    total_energy += inv_v[i]
  end
  limit = g*total_energy
  p = 0
  cumulative_energy = 0.0
  @inbounds while cumulative_energy < limit
    p += 1
    cumulative_energy += inv_v[ind_start]
    ind_start += 1
  end
  p, inadmissable
end

#reduce_dimensions assumes the eigenvalues are sorted from smallest to largest.
function reduce_dimensions!(M::Model, H::Array{Float64,2})
  ef = eigfact!(Symmetric(H))
  out = SparseQuadratureGrids.mats(M.Grid, sum(ef.values .> 1e-11))
  g0 = 0
  for (i, v) ∈ enumerate(ef.values)
    if v < 1e-11
      continue
    end
    g0 += 1
    out[:,g0] .= ef.vectors[:,i] ./ √ef.values[i]
  end
  out
end
function reduce_dimensions!(M::Model, H::Array{Float64,2}, ::Type{R}) where {g, R <: LDR{g}}
  ef = eigfact!(Symmetric(H))
  p, inadmissable = count(ef.values, R)
  out = SparseQuadratureGrids.mats(M.Grid, p)
  for (g0, i) ∈ enumerate(inadmissable+(1:p))
    out[:,g0] .= ef.vectors[:,i] ./ √ef.values[i]
  end
  out
end
function reduce_dimensions!(M::Model, H::Array{Float64,2}, ::Type{FixedRank{p}}) where p
  ef = eigfact!(Symmetric(H))
  out = SparseQuadratureGrids.mats(M.Grid, p)
  g0 = 0
  for (i, v) ∈ enumerate(ef.values)
    if v < 1e-11
      continue
    elseif g0 >= p
      break
    end
    g0 += 1
    out[:,g0] .= ef.vectors[:,i] ./ √ef.values[i]
  end
  out
end

function deduce_scale!(M::Model, H::Array{Float64,2}, ::Type{Dynamic})
  safe_inv_chol!(M.Grid.U, H) ? M.Grid.U : reduce_dimensions!(M, H)
end
function deduce_scale!(M::Model, H::Array{Float64,2}, ::Type{Full{p}}) where p
  inv_chol!(M.Grid.U, H)
end
function deduce_scale!(M::Model, H::Array{Float64,2}, ::Type{R}) where R
  reduce_dimensions!(M, H, R)
end

@inline function log_density!(Θ::ModelParam, data, neg_min::Float64 = 0.0)
  update!(Θ)
  log_density(Θ, data) + neg_min
end
@inline function log_density_cache(x::AbstractVector{<:Real}, θ::Tuple, data, neg_min::Float64 = 0.0)
    param = construct(x, θ)
    log_density(param, data) + neg_min, param
end

sample_size_order(::Any) = 1

index(U::Array{Float64,2},::Type{<:DynamicRank}, data, n::Int) = size(U,2),sample_size_order(data), n
index(U::Array{Float64,2},::Type{<:StaticRank}, data, n::Int) = sample_size_order(data), n
index(U::Array{Float64,2},::Type{<:DynamicRank}, data, seq::Vector{Int}) = size(U,2), seq
index(U::Array{Float64,2},::Type{<:StaticRank}, data, seq::Vector{Int}) = seq

function mode(M::Model{G, MP, P, R} where {G, MP, P}, data, f::Function = Main.log_density) where R
    md = ModelDiff(f, M.diff_buffer, data)
	optBFGS!(md) ? hessian!(md) : optNTR!(md)
	M.diff_buffer.state.x, deduce_scale!(M, 2M.diff_buffer.dr.derivs[2], R), M.diff_buffer.dr.value
end

function SlidingVecFun( Θ::ModelParam, data, neg_min::Float64 )
    ld = () -> log_density!( Θ, data, neg_min )
    SlidingVecFun(ld, Θ.v)
end

function fit( M::Model{G, MP, P, R} where P, data, n = default(B) ) where {q, B <: CacheBuild, G <: GridVessel{q, B}, MP, R}
    μ_hat, U, neg_min = mode( M, data )
    ldc(x::AbstractVector) = log_density_cache(x, M.ϕ, data, neg_min)
    Θ, density = eval_grid!( M.Grid, ldc, μ_hat, U, index(U, R, data, n), MP )
    JointPosterior( M, Θ, density, μ_hat, U )
end
function fit( M::Model{G, MP, P, R} where {MP, P}, data, n = default(B) ) where {q, B <: RawBuild, G <: GridVessel{q, B}, R}
    μ_hat, U, neg_min = mode( M, data )
    svf = SlidingVecFun( M.Θ, data, neg_min )
    g = eval_grid!( M.Grid, svf, μ_hat, U, index(U, R, data, n) )
    JointPosteriorRaw( M, g, μ_hat, U )
end


text:
struct UnitCell{D}
    primitive::NTuple{D,NTuple{D,Float64}}
    basis::Vector{NTuple{D,Float64}}
    interactions::Vector{Tuple{Int,Int,NTuple{D,Int},Matrix{Float64}}} #interactions specified as (basis1,basis2,offsetPrimitive,M)
    interactionsOnsite::Vector{Matrix{Float64}}
    interactionsField::Vector{Vector{Float64}}

    UnitCell(a1::NTuple{1,Float64}) = new{1}((a1,), Vector{NTuple{1,Float64}}(undef,0), Vector{Tuple{Int,Int,NTuple{1,Int},Matrix{Float64}}}(undef,0), Vector{Matrix{Float64}}(undef,0), Vector{Vector{Float64}}(undef,0))
    UnitCell(a1::NTuple{2,Float64}, a2::NTuple{2,Float64}) = new{2}((a1,a2), Vector{NTuple{2,Float64}}(undef,0), Vector{Tuple{Int,Int,NTuple{2,Int},Matrix{Float64}}}(undef,0), Vector{Matrix{Float64}}(undef,0), Vector{Vector{Float64}}(undef,0))
    UnitCell(a1::NTuple{3,Float64}, a2::NTuple{3,Float64}, a3::NTuple{3,Float64}) = new{3}((a1,a2,a3), Vector{NTuple{3,Float64}}(undef,0), Vector{Tuple{Int,Int,NTuple{3,Int},Matrix{Float64}}}(undef,0), Vector{Matrix{Float64}}(undef,0), Vector{Vector{Float64}}(undef,0))
    UnitCell(primitives...) = new{length(primitives)}(primitives, Vector{NTuple{length(primitives),Float64}}(undef,0), Vector{Tuple{Int,Int,NTuple{length(primitives),Int},Matrix{Float64}}}(undef,0), Vector{Matrix{Float64}}(undef,0), Vector{Vector{Float64}}(undef,0))
end

"""
Adds an interaction between spin1 located at basis site `b1` of the given `unitcell` and spin2 at basis site `b2` in a unit cell that is offset by `offset` lattice vectors.
The exchange energy is calculated as spin1'.M.spin2.
"""
function addInteraction!(unitcell::UnitCell{D}, b1::Int, b2::Int, M::Matrix{Float64}, offset::NTuple{D,Int}=Tuple(zeros(Int,D))) where D
    size(M) == (3,3) || error("Interaction matrix must be of size 3x3.")
    b1 == b2 && offset == Tuple(zeros(Int,D)) && error("Interaction cannot be local. Use setInteractionOnsite!() instead.")

    push!(unitcell.interactions, (b1,b2,offset,M))
end

function setInteractionOnsite!(unitcell::UnitCell{D}, b::Int, M::Matrix{Float64}) where D
    size(M) == (3,3) || error("Interaction matrix must be of size 3x3.")
    unitcell.interactionsOnsite[b] = M
end

function setField!(unitcell::UnitCell{D}, b::Int, B::Vector{Float64}) where D
    size(B) == (3,) || error("Field must be a vector of length 3.")
    unitcell.interactionsField[b] = B
end

function addBasisSite!(unitcell::UnitCell{D}, position::NTuple{D,Float64}) where D
    push!(unitcell.basis, position)
    push!(unitcell.interactionsOnsite, zeros(3,3))
    push!(unitcell.interactionsField, zeros(3))
    return length(unitcell.basis)
end


text:
instructions = split.(readlines("day8/input.txt"), " ")

visited_lines = Set{Int}()

function execute_line(line_num, instruction, num, acc)
    if instruction == "nop"
        return line_num+1, acc
    elseif instruction == "jmp"
        return line_num+num, acc
    elseif instruction == "acc"
        return line_num+1, acc+num
    end
end

num, acc = execute_line(1, instructions[1][1], parse(Int, instructions[1][2]), 0)
push!(visited_lines, num)
while true
    global acc
    global num, acc = execute_line(num, instructions[num][1], parse(Int, instructions[num][2]), acc)
    if num in visited_lines
        println(acc)
        break
    end
    push!(visited_lines, num)
end


text:
type SymbolDistr{SymbolT}
  symbols::Vector{SymbolT}
  probs::Vector{Float64}
end

probability_vector{T <: Number}(numbers::Vector{T}) = numbers ./ sum(numbers)

# Constructor from Dicts
SymbolDistr{ST, PT}(d::Dict{ST, PT}) = begin
  symbols = collect(keys(d))
  probs = convert(Vector{Float64}, [d[s] for s in symbols])
  SymbolDistr{ST}(symbols, probability_vector(probs))
end

# Return the probability of _symbol_ in distribution _d_.
function prob{ST}(d::SymbolDistr{ST}, symbol)
  index = findfirst(d.symbols, symbol)
  (index == 0) ? 0.0 : d.probs[index]
end

# Return the probabilities of _symbols_ in distribution _d_.
probs{ST}(d::SymbolDistr{ST}, symbols) = convert(Vector{Float64}, map((s) -> prob(d, s), symbols))

# Sample a probability vector and return the index to the
# probability that matched. Assumes the vector sums to 1.0.
function sample(probs::Vector{Float64})
  value = rand()
  for i in 1:length(probs)
    if value < probs[i]
      return i
    else
      value -= probs[i]
    end
  end
  throw("Probability vector was not normalized but has sum $(sum(probs))")
end

# Sample a symbol from a distribution.
sample{ST}(d::SymbolDistr{ST}) = d.symbols[sample(d.probs)]

using Distance

# Distance between two symbol distributions.
function distance{ST}(d1::SymbolDistr{ST}, d2::SymbolDistr{ST}, distanceMetric = JSDivergence())
  symbols = union(Set(d1.symbols), Set(d2.symbols))
  evaluate(distanceMetric, probs(d1, symbols), probs(d2, symbols))
end

# Count the number of samples of each type.
function count_samples(samples)
  # Count the number of samples of each kind
  counts = Dict{Any, Int64}()
  for s in samples
    counts[s] = get(counts, s, 0) + 1
  end
  counts
end

# Create a map for the empirical distribution of samples from a discrete alphabet
function empirical_distribution(samples)
  counts = count_samples(samples)
  # Normalize to frequencies by dividing by the number of samples
  distr = Dict{Any, Float64}()
  for (k, v) in counts
    distr[k] = counts[k] / length(samples)
  end
  distr
end

# A symbol stream represents an actual or a modeled stream of symbols from a
# finite alphabet.
abstract SymbolStream
alphabet(s::SymbolStream) = s.alphabet # Override if not explicitly saved in the stream

# Return the next sample from a stream if there is one or nothing otherwise.
next(s::SymbolStream) = nothing

type WhiteNoiseStream <: SymbolStream
  alphabet
end

sample(wns::WhiteNoiseStream) = wns.alphabet[rand(1:length(wns.alphabet))]

# We can always get a next sample from a WNS.
next(wns::WhiteNoiseStream) = sample(wns)

# We can also create a WhiteNoiseStream from a symbol distribution and from
# another SymbolStream.
WhiteNoiseStream{ST}(sd::SymbolDistr{ST}) = WhiteNoiseStream(sd.symbols)
WhiteNoiseStream(s::SymbolStream) = WhiteNoiseStream(alphabet(s))

# A fixed stream goes through the values in a given sequence and then stops.
type FixedStream <: SymbolStream
  alphabet
  values
  index
  max_index
  FixedStream(values) = new(unique(values), values, 1, length(values))
end

function next(fs::FixedStream)
  if fs.index <= fs.max_index
    fs.index += 1
    fs.values[fs.index-1]
  else
    nothing
  end
end

# An independent stream copy emits only the symbols that match a white noise
# version of a stream.
type IndependentStreamCopy <: SymbolStream
  base::SymbolStream      # The base stream that we are a copy of
  wns::WhiteNoiseStream   # A white noise stream for the same alphabet as the base stream
  IndependentStreamCopy(s::SymbolStream) = new(s, WhiteNoiseStream(s))
end

function next(isc::IndependentStreamCopy)
  next_base_symbol = next(base)
  while next_base_symbol != nothing && next_base_symbol != next(isc.wns)
    next_base_symbol = next(base)
  end
  return next_base_symbol # is either nothing if the base stream is empty or it is a symbol that matched the symbol from the WhiteNoiseStream
end


text:
"""
    PeriodicKernel(; r::AbstractVector=ones(Float64, 1))

Periodic kernel with parameter `r`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the periodic kernel with parameter ``r_i > 0`` is
defined[^DM] as
```math
k(x, x'; r) = \\exp\\bigg(- \\frac{1}{2} \\sum_{i=1}^d \\bigg(\\frac{\\sin\\big(\\pi(x_i - x'_i)\\big)}{r_i}\\bigg)^2\\bigg).
```

[^DM]: D. J. C. MacKay (1998). Introduction to Gaussian Processes.
"""
struct PeriodicKernel{T} <: SimpleKernel
    r::Vector{T}
    function PeriodicKernel(; r::AbstractVector{<:Real}=ones(Float64, 1))
        @check_args(PeriodicKernel, r, all(ri > zero(ri) for ri in r), "r > 0")
        return new{eltype(r)}(r)
    end
end

PeriodicKernel(dims::Int) = PeriodicKernel(Float64, dims)

"""
    PeriodicKernel([T=Float64, dims::Int=1])

Create a [`PeriodicKernel`](@ref) with parameter `r=ones(T, dims)`.
"""
PeriodicKernel(T::DataType, dims::Int=1) = PeriodicKernel(; r=ones(T, dims))

@functor PeriodicKernel

metric(κ::PeriodicKernel) = Sinus(κ.r)

kappa(::PeriodicKernel, d::Real) = exp(-0.5d)

function Base.show(io::IO, κ::PeriodicKernel)
    return print(io, "Periodic Kernel, length(r) = $(length(κ.r))")
end


text:
using GD
using HTTP

"""
    build_router(central_node::CentralNode)

Build the routes to the central node endpoints.
"""
function build_router(central_node::CentralNode)::HTTP.Router
    router = HTTP.Router()

    HTTP.@register(
        router,
        "POST",
        central_node.config.base.register_node,
        (request::HTTP.Request) -> register_client!(central_node.client_manager, request)
    )

    return router
end

"""
    start(central_node::CentralNode)

Start a `Fed.jl` HTTP server.
"""
function start(central_node::CentralNode)
    router = build_router(central_node)

    # start the HTTP server
    s = @async HTTP.serve(router, central_node.host, central_node.port)

    # wait for all the clients to join...
    wait_for(central_node.client_manager, central_node.config.base.num_total_clients)

    # federated training
    fit(central_node)

    # stop the HTTP server
    @async Base.throwto(s,  InterruptException())
end


text:
"""
    reproduce!(bact::Bacterium, model::ABM)

Simulates reproduction of bacteria. Energy is halved after subtracting a cost. Inherited
features are varied.
"""
function reproduce!(bact::Bacterium, model::ABM)
    add_agent!(
        bact.pos,
        Bacterium,
        model,
        bact.species,
        0,
        (bact.energy - model.reproduction_energy_cost) / 2.0,
        inherit(bact.sensory_radius, model.rng, model.σ_sensory_radius),
        inherit(bact.reproduction_threshold, model.rng, model.σ_reproduction_threshold),
        inherit(bact.speed, model.rng, model.σ_speed),
        bact.food_target,
    )

    bact.age = 0
    bact.energy = (bact.energy - model.reproduction_energy_cost) / 2.0

    bact.sensory_radius = inherit(bact.sensory_radius, model.rng, model.σ_sensory_radius)
    bact.reproduction_threshold =
        inherit(bact.reproduction_threshold, model.rng, model.σ_reproduction_threshold)
    bact.speed = inherit(bact.speed, model.rng, model.σ_speed)
end

function agent_step!(bact::Bacterium, model::ABM)
    # die of age or starvation
    (bact.age > model.lifetime || bact.energy <= 0) && (kill_agent!(bact, model); return)

    # grow older
    bact.age += 1

    # reproduce if possible
    bact.energy >= bact.reproduction_threshold && reproduce!(bact, model)

    # if currently on a food source
    if model.food[bact.pos...] > 0.0
        # eat
        eaten = min(model.food[bact.pos...], bact.speed * model.eat_rate_factor)
        model.food[bact.pos...] -= eaten
        bact.energy += eaten
        return
    else
        # can't eat, so look for something else
        bact.food_target = (-1, -1)
    end

    # if bact doesn't see any food
    if bact.food_target == (-1, -1)
        best_pos = (-1, -1)
        # look for the best food nearby
        for pos in nearby_positions(bact.pos, model, bact.sensory_radius)
            best_pos == (-1, -1) || model.food[best_pos...] < model.food[pos...] || continue
            best_pos = pos
        end
        best_pos != (-1, -1) && (bact.food_target = best_pos)
    end

    # if we found some food or already had eyes on it
    if bact.food_target != (-1, -1)
        # move toward target
        delta = bact.food_target .- bact.pos
        #move as many steps diagonally as possible
        diag = min(abs.(delta)..., bact.speed)

        steps = diag

        movement = sign.(delta) .* diag
        # move the rest steps axially, if possible
        delta = delta .- movement
        while steps < bact.speed && any(delta .> 0)
            movement = movement .+ sign.(delta)
            delta = delta .- sign.(delta)
            steps += 1
        end
    else
        # move randomly if there's no food nearby
        movement = rand.(model.rng, (AXIAL_DIRECTIONS, AXIAL_DIRECTIONS)) .* bact.speed
    end

    # move
    move_agent!(bact, clamp.(bact.pos .+ movement, 1, size(model.space.s)), model)

    # subtract cost of movement and looking
    bact.energy -=
        model.sensory_radius_cost * bact.sensory_radius +
        model.distance_cost * max(abs.(movement)...)
end

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text:
export completion, qAdicConj

#XXX: valuation(Q(0)) == 0 !!!!!
function newton_lift(f::fmpz_poly, r::qadic, prec::Int = parent(r).prec_max, starting_prec::Int = 2)
  Q = parent(r)
  n = prec
  i = n
  chain = [n]
  while i>starting_prec
    i = div(i+1, 2)
    push!(chain, i)
  end
  fs = derivative(f)
  qf = change_base_ring(Q, f, cached = false)
  qfs = change_base_ring(Q, fs, cached = false)
  r = setprecision(r, 1)
  o = Q(r)
  s = qf(r)
  o = inv(setprecision(qfs, 1)(o))
  @assert r.N == 1
  for p = reverse(chain)
    setprecision!(r, p)
    setprecision!(o, p)
    Q.prec_max = r.N
    if r.N > precision(Q)
      setprecision!(qf, r.N)
      setprecision!(qfs, r.N)
    end
    r = r - qf(r)*o
    if r.N >= n
      Q.prec_max = n
      return r
    end
    o = o*(2-qfs(r)*o)
  end
  return r
end

function newton_lift(f::fmpz_poly, r::LocalFieldElem, precision::Int = parent(r).prec_max, starting_prec::Int = 2)
  Q = parent(r)
  n = precision
  i = n
  chain = [n]
  while i > starting_prec
    i = div(i+1, 2)
    push!(chain, i)
  end
  fs = derivative(f)
  qf = change_base_ring(Q, f, cached = false)
  qfs = change_base_ring(Q, fs, cached = false)
  o = Q(r)
  o = setprecision(o, 1)
  s = qf(r)
  o = inv(setprecision!(qfs, 1)(o))
  for p = reverse(chain)
    r = setprecision!(r, p)
    o = setprecision!(o, p)
    setprecision!(Q, p)
    setprecision!(qf, p)
    setprecision!(qfs, p)
    r = r - qf(r)*o
    if Nemo.precision(r) >= n
      setprecision!(Q, n)
      return r
    end
    o = o*(2-qfs(r)*o)
  end
  return r
end

@doc Markdown.doc"""
    roots(f::fmpz_poly, Q::FlintQadicField; max_roots::Int = degree(f)) -> Vector{qadic}

The roots of $f$ in $Q$, $f$ has to be square-free (at least the roots have to be simple roots).
"""
function roots(f::fmpz_poly, Q::FlintQadicField; max_roots::Int = degree(f))
  k, mk = ResidueField(Q)
  rt = roots(f, k)
  RT = qadic[]
  for r = rt
    push!(RT, newton_lift(f, preimage(mk, r)))
    if length(RT) >= max_roots
      return RT
    end
  end
  return RT
end

is_splitting(C::qAdicRootCtx) = C.is_splitting

function roots(C::qAdicRootCtx, n::Int = 10)
  if isdefined(C, :R) && all(x -> x.N >= n, C.R)
    return [setprecision(x, n) for x = C.R]
  end
  lf = factor_mod_pk(Array, C.H, n)
  rt = qadic[]
  for Q = C.Q
    Q.prec_max = n
    for x = lf
      if is_splitting(C) || degree(x[1]) == degree(Q)
        append!(rt, roots(x[1], Q, max_roots = 1))
      end
    end
  end
  if isdefined(C, :R)
    st = qadic[]
    for r = C.R
      p = findfirst(x -> degree(parent(r)) == degree(parent(x)) && iszero(x-r), rt)
      push!(st, rt[p])
    end
    rt = st
  end
  C.R = rt
  return rt
end

@doc Markdown.doc"""
    qAdicConj(K::AnticNumberField, p::Int)

Creates a data structure to compute the conjugates in an unramified splitting field
over $Q_p$.
"""
mutable struct qAdicConj
  K::AnticNumberField
  C::qAdicRootCtx
  cache::Dict{nf_elem, Any}

  function qAdicConj(K::AnticNumberField, p::Int; splitting_field::Bool = false)
    if discriminant(map_coefficients(GF(p), K.pol)) == 0
      error("cannot deal with difficult primes yet")
    end
    #=
    is_index_divisor(maximal_order(K), p) && error("cannot deal with index divisors yet")
    is_ramified(maximal_order(K), p) && error("cannot deal with ramification yet")
    =#
    if splitting_field
      Zx = PolynomialRing(FlintZZ, cached = false)[1]
      C = qAdicRootCtx(Zx(K.pol), p, splitting_field = true)
      r = new()
      r.C = C
      r.K = K
      r.cache = Dict{nf_elem, Any}()
      return r
    end
    D = get_attribute!(K, :nf_conjugate_data_qAdic) do
      return Dict{Int, Tuple{qAdicRootCtx, Dict{nf_elem, Any}}}()
    end::Dict{Int, Tuple{qAdicRootCtx, Dict{nf_elem, Any}}}
    Dp = get!(D, p) do
      Zx = PolynomialRing(FlintZZ, cached = false)[1]
      d = lcm(map(denominator, coefficients(K.pol)))
      C = qAdicRootCtx(Zx(K.pol*d), p)
      return (C, Dict{nf_elem, Any}())
    end

    return new(K, Dp[1], Dp[2])
  end
end

function Base.show(io::IO, C::qAdicConj)
  println(io, "data for the $(C.C.p)-adic completions of $(C.K)")
end

#to compare to the classical conjugates
#  all = true/ false: only on of a pair of complex conjugates is returned
#  flat = true/ false: return (Re, Im) or the complex number
#TODO: not sure how this would work in the ramified, not-normal case.
@doc Markdown.doc"""
    conjugates(a::nf_elem, C::qAdicConj, n::Int = 10; flat::Bool = false, all:Bool = true) -> []

Returns an array of the $q$-adic conjugates of $a$: Let $p Z_K = \prod P_i$ for the maximal order
$Z_K$ of the parent of $a$. Then $K \otimes Q_p = \prod K_{P_i}$. For each of the $P_i$
a $q$-adic (unramifed) extension $K_{P_i}$ of $Q_p$ is computed, sth. $a$ has $\deg P_i = \deg K_{P_i}$
many cojugates in $K_{P_i}$.
If `all = true` and `flat = false`, the default, then all $n$ conjugates are returned.
If `all = false`, then for each $P_i$ only one conjugate is returned, the others could be
computed using automorphisms (the Frobenius).
If `flat = true`, then instead of the conjugates, only the $p$-adic coefficients are returned.
"""
function conjugates(a::nf_elem, C::qAdicConj, n::Int = 10; flat::Bool = false, all::Bool = true)
  if is_splitting(C.C)
    return expand(_conjugates(a, C, n, x -> x), flat = flat, all = all, degs = degrees(C.C.H))
  else
    return expand(_conjugates(a, C, n, x -> x), flat = flat, all = all)
  end
end

function expand(a::Vector{qadic}; all::Bool, flat::Bool, degs::Vector{Int}= Int[])
  re = qadic[]
  if all
    for ix = 1:length(a)
      x = a[ix]
      push!(re, x)
      y = x
      d = degree(parent(x))
      if ix <= length(degs)
        for i=2:degs[ix]
          y = frobenius(y)
          push!(re, y)
        end
      else
        for i=2:degree(parent(x))
          y = frobenius(y)
          push!(re, y)
        end
      end
    end
  else
    re = a
  end
  if flat
    r = padic[]
    for x = re
      for i=1:degree(parent(x))
        push!(r, coeff(x, i-1))
      end
    end
    return r
  else
    return re
  end
end

#TODO: implement a proper Frobenius - with caching of the frobenius_a element
function _conjugates(a::nf_elem, C::qAdicConj, n::Int, op::Function)
  R = roots(C.C, n)
  @assert parent(a) == C.K
  Zx = PolynomialRing(FlintZZ, cached = false)[1]
  d = denominator(a)
  f = Zx(d*a)
  res = qadic[]
  for x = R
    a = op(inv(parent(x)(d))*f(x))::qadic
    push!(res, a)
  end
  return res
end

function _log(a::qadic)
  q = prime(parent(a))^degree(parent(a))
  if iseven(q) # an error in flint
    return log((a^(q-1))^2)//2//(q-1)
  end
  return log(a^(q-1))//(q-1) # faster than the teichmuller stuff
  return log(a*inv(teichmuller(a)))
end

@doc Markdown.doc"""
    conjugates_log(a::nf_elem, C::qAdicConj, n::Int = 10; flat::Bool = false, all:Bool = true) -> []
    conjugates_log(a::FacElem{nf_elem, AnticNumberField}, C::qAdicConj, n::Int = 10; flat::Bool = false, all:Bool = true) -> []

Returns an array of the logarithms of the $q$-adic conjugates of $a$: Let $p Z_K = \prod P_i$ for the maximal order
$Z_K$ of the parent of $a$. Then $K \otimes Q_p = \prod K_{P_i}$. For each of the $P_i$
a $q$-adic (unramifed) extension $K_{P_i}$ of $Q_p$ is computed, sth. $a$ has $\deg P_i = \deg K_{P_i}$
many cojugates in $K_{P_i}$.
If `all = true` and `flat = false`, then all $n$ logarithms of conjugates are returned.
If `all = false`, then for each $P_i$ only one logarithm of a conjugate is returned, the others could be
computed using automorphisms (the Frobenius).
If `flat = true`, then instead of the conjugates, only the $p$-adic coefficients are returned.
"""
function conjugates_log(a::nf_elem, C::qAdicConj, n::Int = 10; all::Bool = false, flat::Bool = true)
  if haskey(C.cache, a)
    b = C.cache[a]
    if b[1,1].N == n
      return expand(b, all = all, flat = flat)
    end
  end
  C.cache[a] = b = _conjugates(a, C, n, _log)
  return expand(b, all = all, flat = flat)
end

function conjugates_log(a::FacElem{nf_elem, AnticNumberField}, C::qAdicConj, n::Int = 10; all::Bool = false, flat::Bool = true)
  first = true
  local res::Vector{qadic}
  for (k, v) = a.fac
    try
      y = conjugates_log(k, C, n, flat = false, all = false)
      if first
        res = v .* y
        first = false
      else
        res += v .* y
      end
    catch e
      if isa(e, DivideError) || isa(e, DomainError)
        lp = prime_decomposition(maximal_order(parent(k)), C.C.p)
        @assert Base.all(x -> has_2_elem_normal(x[1]), lp)
        val = map(x -> valuation(k, x[1]), lp)
        pe = prod(lp[i][1].gen_two^val[i] for i = 1:length(lp) if val[i] != 0)
        aa = k//pe
        y = conjugates_log(aa, C, n, all = false, flat = false)
        if first
          res = v .* y
          first = false
        else
          res += v .* y
        end
      else
        rethrow(e)
      end
    end
  end

  if is_splitting(C.C)
    return expand(res, flat = flat, all = all, degs = degrees(C.C.H))
  else
    return expand(res, all = all, flat = flat)
  end
end

function special_gram(m::Vector{Vector{qadic}})
  g = Vector{padic}[]
  for i = m
    r = padic[]
    for j = m
      k = 1
      S = 0
      while k <= length(i)
        s = i[k] * j[k]
        for l = 1:degree(parent(s))-1
          s += i[k+l] * j[k+l]
        end
        S += coeff(s, 0)
        @assert s == coeff(s, 0)
        k += degree(parent(s))
      end
      push!(r, S)
    end
    push!(g, r)
  end
  return g
end

function special_gram(m::Vector{Vector{padic}})
  n = transpose(matrix(m))
  n = transpose(n)*n
  return [[n[i,j] for j=1:ncols(n)] for i = 1:nrows(n)]
end

@doc Markdown.doc"""
    regulator(u::Vector{T}, C::qAdicConj, n::Int = 10; flat::Bool = true) where {T<: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}}
    regulator(K::AnticNumberField, C::qAdicConj, n::Int = 10; flat::Bool = true)
    regulator(R::NfAbsOrd, C::qAdicConj, n::Int = 10; flat::Bool = true)

Returns the determinant of $m^t m$ where the columns of $m$ are the `conjugates_log` of the units
in either the array, or the fundamental units for $K$ (the maximal order of $K$) or $R$.
If `flat = false`, then all prime ideals over $p$ need to have the same degree.
In either case, Leopold's conjecture states that the regulator is zero iff the units are dependent.
"""
function regulator(u::Vector{T}, C::qAdicConj, n::Int = 10; flat::Bool = true) where {T<: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}}
  c = map(x -> conjugates_log(x, C, n, all = !flat, flat = flat), u)
  return det(transpose(matrix(special_gram(c))))
end

function regulator(K::AnticNumberField, C::qAdicConj, n::Int = 10; flat::Bool = false)
  return regulator(maximal_order(K), C, n, flat = flat)
end

function regulator(R::NfAbsOrd{AnticNumberField, nf_elem}, C::qAdicConj, n::Int = 10; flat::Bool = false)
  u, mu = unit_group_fac_elem(R)
  return regulator([mu(u[i]) for i=2:ngens(u)], C, n, flat = flat)
end

@doc Markdown.doc"""
    regulator_iwasawa(u::Vector{T}, C::qAdicConj, n::Int = 10) where {T<: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}} -> qadic
    regulator_iwasawa(K::AnticNumberField, C::qAdicConj, n::Int = 10) -> qadic
    regulator_iwasawa(R::NfAbsOrd, C::qAdicConj, n::Int = 10) -> qadic

For a totally real field $K$, the regulator as defined by Iwasawa: the determinant of the
matrix containing the logarithms of the conjugates, supplemented by a column containing all $1$.
"""
function regulator_iwasawa(u::Vector{T}, C::qAdicConj, n::Int = 10) where {T<: Union{nf_elem, FacElem{nf_elem, AnticNumberField}}}
  k = base_ring(u[1])
  @assert is_totally_real(k)
  c = map(x -> conjugates_log(x, C, n, all = true, flat = false), u)
  m = matrix(c)
  m = hcat(m, matrix(base_ring(m), nrows(m), 1, [one(base_ring(m)) for i=1:nrows(m)]))
  return det(m)//degree(k)
end

function regulator_iwasawa(K::AnticNumberField, C::qAdicConj, n::Int = 10)
  @assert is_totally_real(K)
  return regulator_iwasawa(maximal_order(K), C, n)
end

function regulator_iwasawa(R::NfAbsOrd, C::qAdicConj, n::Int = 10)
  @assert is_totally_real(nf(R))
  u, mu = unit_group_fac_elem(R)
  return regulator_iwasawa([mu(u[i]) for i=2:ngens(u)], C, n)
end

function matrix(a::Vector{Vector{T}}) where {T}
  return matrix(permutedims(reduce(hcat, a), (2, 1)))
end

function eval_f_fs(f::PolyElem, x::RingElem)
  d = Int[]
  for i=1:degree(f)
    if !iszero(coeff(f, i))
      if i>0 && !((i-1) in d)
        push!(d, i-1)
      end
      push!(d, i)
    end
  end
  p = Dict{Int, typeof(x)}()
  p[0] = one(x)
  p[1] = x
  p[d[1]] = x^d[1]

  for i = 2:length(d)
    if haskey(p, d[i])
      continue
    end
    q, r = divrem(d[i], d[i-1])
    if haskey(p, r)
      xr = p[r]
    else
      xr = p[r] = x^r
    end
    p[d[i]] = p[d[i-1]]^q * xr
  end
  s1 = zero(x)
  s2 = zero(x)
  for i=0:degree(f)
    c = coeff(f, i)
    if !iszero(c)
      s1 += c*p[i]
      if i>0
        s2 += i*c*p[i-1]
      end
    end
  end
  return s1, s2
end

struct nf_elem_mod <: RingElem
  a::nf_elem
  p::fmpz
end
function *(a::fmpz, b::nf_elem_mod)
  c = a*b.a
  return nf_elem_mod(mod_sym(c, b.p), b.p)
end
function *(a::nf_elem_mod, b::nf_elem_mod)
  c = a.a*b.a
  return nf_elem_mod(mod_sym(c, a.p), a.p)
end
function one(a::nf_elem_mod)
  return nf_elem_mod(one(a.a), a.p)
end
function zero(a::nf_elem_mod)
  return nf_elem_mod(zero(a.a), a.p)
end
function +(a::nf_elem_mod, b::nf_elem_mod)
  return nf_elem_mod(a.a+b.a, a.p)
end
function ^(a::nf_elem_mod, i::Int)
  b = one(a)
  c = a
  while i > 0
    if isodd(i)
      b *= c
    end
    i >>= 1
    if !iszero(i)
      c *= c
    end
  end
  return b
end
function lift_root(f::fmpz_poly, a::nf_elem, o::nf_elem, p::fmpz, n::Int)
  #f(a) = 0 mod p, o*f'(a) = 1 mod p, want f(a) = 0 mod p^n
  k = 1
  while k < n
    p *= p
    k *= 2
    #TODO: here f wil be sparse (and possibly large degree), so
    #      this evaluation is bad.
    # in the calling cite: don't work in the large field, restrict
    # to working (mod p^k) in the field defined by the factor

    if false
      pa = [one(a)]
      while length(pa) <= degree(f)
        push!(pa, pa[end]*a)
        mod_sym!(pa[end], p)
      end
      fa  = sum(coeff(f, i-1) * pa[i] for i=1:length(pa))
      fsa = sum(coeff(f, i) * i * pa[i] for i=1:length(pa)-1)
    else
      _fa, _fsa = eval_f_fs(f, nf_elem_mod(a, p))
      fa = _fa.a
      fsa = _fsa.a
    end
    o = o*(2-fsa*o)
    a = a - fa*o
    mod_sym!(o, p)
    mod_sym!(a, p)
  end
  return a
end

@doc Markdown.doc"""
    completion(K::AnticNumberField, P::NfOrdIdl) -> FlintQadicField, Map{AnticNumberField -> FlintQadicField}

The completion of $K$ wrt to the topology induced by the valuation at $P$. $P$ needs
to be unramifed.
The map giving the embedding of $K$ into the completion, admits a pointwise pre-image to obtain a lift.
Note, that the map is not well defined by this data: $K$ will have $\deg P$ many embeddings.
"""
function completion(K::AnticNumberField, P::NfOrdIdl)
  #non-unique!! will have deg(P) many
  p = minimum(P)
  C = qAdicConj(K, Int(p))
  g = conjugates(P.gen_two.elem_in_nf, C)
#  @show map(x->valuation(x), g)
  i = findfirst(x->valuation(x) > 0, g)
  return completion(K, p, i[1])
end

completion(K::AnticNumberField, p::Integer, i::Int) = completion(K, fmpz(p), i)

@doc Markdown.doc"""
    completion(K::AnticNumberField, p::fmpz, i::Int) -> FlintQadicField, Map

The completion corresponding to the $i$-th conjugate in the non-canonical ordering of
`conjugates`.
"""
function completion(K::AnticNumberField, p::fmpz, i::Int)
  C = qAdicConj(K, Int(p))
  @assert 0<i<= degree(K)

  ca = conjugates(gen(K), C, all = true, flat = false)[i]
  return completion(K, ca)
end

function completion(K::AnticNumberField, ca::qadic)
  p = prime(parent(ca))
  C = qAdicConj(K, Int(p))
  r = roots(C.C, precision(ca))
  i = findfirst(x->parent(r[x]) == parent(ca) && r[x] == ca, 1:length(r))
  Zx = PolynomialRing(FlintZZ, cached = false)[1]
  function inj(a::nf_elem)
    d = denominator(a)
    pr = precision(parent(ca))
    if pr > precision(ca)
      ri = roots(C.C, precision(parent(ca)))[i]
    else
      ri = ca
    end
    return inv(parent(ca)(d))*(Zx(a*d)(ri))
  end
  # gen(K) -> conj(a, p)[i] -> a = sum a_i o^i
  # need o = sum o_i a^i
  R, mR = ResidueField(parent(ca))
  pa = [one(R), mR(ca)]
  d = degree(R)
  while length(pa) < d
    push!(pa, pa[end]*pa[2])
  end
  m = matrix(GF(p), d, d, [coeff(pa[i], j-1) for j=1:d for i=1:d])
  o = matrix(GF(p), d, 1, [coeff(gen(R), j-1) for j=1:d])
  s = solve(m, o)
  @hassert :qAdic 1 m*s == o
  a = K()
  for i=1:d
    _num_setcoeff!(a, i-1, lift(s[i,1]))
  end
  f = defining_polynomial(parent(ca), FlintZZ)
  fso = inv(derivative(f)(gen(R)))
  o = matrix(GF(p), d, 1, [coeff(fso, j-1) for j=1:d])
  s = solve(m, o)
  b = K()
  for i=1:d
    _num_setcoeff!(b, i-1, lift(s[i,1]))
  end

  #TODO: don't use f, use the factors i the HenselCtx
  #seems to be slower...
#  lf = factor_mod_pk(Array, C.C.H, Int(C.C.H.N))
#  jj = findfirst(x->iszero(x[1](ca)), lf)
#  Kjj = number_field(lf[jj][1], check = false, cached = false)[1]
#  ajj = Kjj(parent(Kjj.pol)(a))
#  bjj = Kjj(parent(Kjj.pol)(b))
#  cjj = lift_root(f, ajj, bjj, p, 10)
#  c = K(parent(K.pol)(cjj))

  c = lift_root(f, a, b, p, 10)
  pc = fmpz(10)
  function lif(x::qadic)
    if iszero(x)
      return K(0)
    end
    if precision(x) > pc
      #XXX this changes (c, pc) inplace as a cache
      #probably should be done with a new map type that can
      #store c, pc on the map.
      d = lift_root(f, a, b, p, precision(x))
#  Kjj = number_field(lf[jj][1], check = false, cached = false)[1]
#  ajj = Kjj(parent(Kjj.pol)(a))
#  bjj = Kjj(parent(Kjj.pol)(b))
#  djj = lift_root(f, ajj, bjj, p, 10)
#  d = K(parent(K.pol)(djj))
      ccall((:nf_elem_set, libantic), Nothing, (Ref{nf_elem}, Ref{nf_elem}, Ref{AnticNumberField}), c, d, K)
      ccall((:fmpz_set_si, libflint), Nothing, (Ref{fmpz}, Cint), pc, precision(x))
    elseif precision(x) < pc
      d = mod_sym(c, p^precision(x))
    else
      d = c
    end
    n = x.length
    r = K(lift(coeff(x, n-1)))
    pk = p^precision(x)
    while n > 1
      n -= 1
      r = mod_sym(r*d, pk) + lift(coeff(x, n-1))
    end
    return r#*K(p)^valuation(x)
  end
  return parent(ca), MapFromFunc(inj, lif, K, parent(ca))
end


text:
"""
Store the k-cells of a mesh in an octree.

    function octree{U,D,T}(mesh::Mesh{U,D,T}, kcells::Array{Int,2})
"""
function octree(mesh::AbstractMesh)

    nverts = dimension(mesh) + 1 # number of vertices per cell
    ncells = numcells(mesh)

    P = vertextype(mesh)
    T = coordtype(mesh)

    points = zeros(P, ncells)
    radii = zeros(T, ncells)

    for (i,cl) in enumerate(cells(mesh))
        ch = chart(mesh, cl)
        points[i] = cartesian(center(ch))
        radii[i] = maximum(norm(v-points[i]) for v in ch.vertices)
    end

    return Octree(points, radii)
end

"""
    overlap_gpredicate(γ::Mesh) -> (patch -> Bool)

Create a predicate that for a given patch determinees if it overlaps with the provided target mesh `γ`.
"""
function overlap_gpredicate(γ::AbstractMesh)

    if numcells(γ) == 0
        return simplex -> false
    end

    # build an octree with the k-cells of smallmesh
    tree = octree(γ)

    function pred(p1)

        # find all simplices of the small mesh that potentially
        # collide with this patch
        c1, s1 = boundingbox(p1)
        for box in boxes(tree, (c,s)->boxesoverlap(c,s,c1,s1))
            for i in box
                # p2 = chart(γ, γ.faces[i])
                p2 = chart(γ, cells(γ)[i])
                overlap(p1,p2) && return true
            end
        end

        return false
    end

    return pred
end

"""
Geometric predicate for determining in log(N) complexity if a the image of a chart is in
the closure of mesh `γ`.
"""
function inclosure_gpredicate(γ::AbstractMesh)

    if numcells(γ) == 0
        return simplex -> false
    end

    # build an octree with the k-cells of smallmesh
    tree = octree(γ)

    function pred(p1)

        # find all simplices of the small mesh that potentially
        # collide with this patch
        c1, s1 = boundingbox(p1)
        for box in boxes(tree, (c,s)->boxesoverlap(c,s,c1,s1))
            for i in box
                # p2 = simplex(vertices(γ, γ.faces[i]))
                p2 = chart(γ, cells(γ)[i])
                isinclosure(p2, p1) && return true
            end
        end

        return false
    end

    return pred
end


text:
"""
Computes Z4 topological invariants given an irvsp.out file, range of bands [nstart,nend] for systems with TRS and IS
"""
function compute_Z4(myFile::String, nstart::Int64, nend::Int64)
    parity = parse_parity_from_irvsp(myFile, nstart, nend)
    # Compute the Strong Index
    nₖ⁺ = sum(x->x>0, parity, dims=1) # Compute number of +1 parity eigenvalues
    nₖ⁻ = sum(x->x<0, parity, dims=1) # Compute number of -1 parity eigenvalues
    κ = mod(1/4*sum( nₖ⁺ .- nₖ⁻), 4) # Compute κ, the Z₄ invariant
    ν₀ = mod(κ, 2)

    # Compute the Weak Indices
    nₖ₁⁺ = nₖ⁺[[2,5,6,8]]
    nₖ₁⁻ = nₖ⁻[[2,5,6,8]]
    ν₁ = mod(1/4*sum(nₖ₁⁺.- nₖ₁⁻), 2)
    nₖ₂⁺ = nₖ⁺[[3,5,7,8]]
    nₖ₂⁻ = nₖ⁻[[3,5,7,8]]
    ν₂ = mod(1/4*sum(nₖ₂⁺ .- nₖ₂⁻), 2)
    nₖ₃⁺ = nₖ⁺[[4,6,7,8]]
    nₖ₃⁻ = nₖ⁻[[4,6,7,8]]
    ν₃ = mod(1/4*sum(nₖ₃⁺.- nₖ₃⁻), 2)

    Z₄ = "Z₄ = (κ;ν₁ν₂ν₃) = ($(trunc(Int,κ));$(trunc(Int,ν₁))$(trunc(Int,ν₂))$(trunc(Int,ν₃)))"
    Z₂ = "Z₂ = (ν₀;ν₁ν₂ν₃) = ($(trunc(Int,ν₀));$(trunc(Int,ν₁))$(trunc(Int,ν₂))$(trunc(Int,ν₃))), where ν₀ = κ mod 2"
    if κ ∈ (1, 3)
        result = "strong topological insulator"
    elseif κ == 0
        result = "trivial insulator"
        if (ν₁ == 1 || ν₂ == 1 || ν₃ == 1)
            result = "weak topological insulator"
        end
    elseif κ == 2
        result = "higher-order topological insulator"
    end

    println(Z₄)
    println(Z₂)
    println("System is a $result")
end

"""
Computes Z2 topological invarilnts given an irvsp.out file, range of bands [nstart,nend] for systems with TRS and IS
"""
function compute_Z2(myFile::String, nstart::Int64, nend::Int64)
    parity = parse_parity_from_irvsp(myFile, nstart, nend)
    δ = ones(8)
    for n = 1:2:(nend-nstart+1)
        δ .*= parity[n, :]
    end
    minus1_ν₀ = prod(δ)
    minus1_ν₁ = prod(δ[[2,5,6,8]])
    minus1_ν₂ = prod(δ[[3,5,7,8]])
    minus1_ν₃ = prod(δ[[4,6,7,8]])
    ν₀ = minus1_ν₀ == 1 ? 0 : 1
    ν₁ = minus1_ν₁ == 1 ? 0 : 1
    ν₂ = minus1_ν₂ == 1 ? 0 : 1
    ν₃ = minus1_ν₃ == 1 ? 0 : 1
    if ν₀ == 1
        result = "strong topological insulator"
    elseif ν₀ == 0
        result = "trivial insulator"
        if (ν₁ == 1 || ν₂ == 1 || ν₃ == 1)
            result = "weak topological insulator"
        end
    end
    Z₂ = "Z₂ = (ν₀;ν₁ν₂ν₃) = ($(trunc(Int,ν₀));$(trunc(Int,ν₁))$(trunc(Int,ν₂))$(trunc(Int,ν₃)))"
    println(Z₂)
    println("System is a $result")
end


text:
"""
get source spray weights
"""
function get_spray_weights!(weights::AbstractArray, denomI::AbstractArray,
			   ix1::AbstractArray{Int64}, ix2::AbstractArray{Int64}, iz1::AbstractArray{Int64}, iz2::AbstractArray{Int64},
			   mesh_x::AbstractArray, mesh_z::AbstractArray, xval::Float64, zval::Float64)

	ix1[1], ix2[1] = indminn(mesh_x,xval,2)
	iz1[1], iz2[1] = indminn(mesh_z,zval,2)

	denomI[1] = ((mesh_x[ix2[1]] - mesh_x[ix1[1]])*(mesh_z[iz2[1]] - mesh_z[iz1[1]]))^(-1.e0)
	"for iz1, ix1"
	weights[1] = 	(mesh_x[ix2[1]]-	xval)*	(mesh_z[iz2[1]]-	zval)*	denomI[1]
	"for iz1, ix2"
	weights[2] =   (xval-	mesh_x[ix1[1]])*	(mesh_z[iz2[1]]-	zval)*	denomI[1]
	"for iz2, ix1"
	weights[3] = 	(mesh_x[ix2[1]]-	xval)*	(zval-	mesh_z[iz1[1]])*	denomI[1]
	"for iz2, ix2"
	weights[4] = 	(xval-	mesh_x[ix1[1]])*	(zval-	mesh_z[iz1[1]])*	denomI[1]

	return weights, denomI, ix1, ix2, iz1, iz2
end

"""
get interpolation weights
"""
function get_interpolate_weights!(weights::AbstractArray, denomI::AbstractArray,
			   ix1::AbstractArray{Int64}, ix2::AbstractArray{Int64}, iz1::AbstractArray{Int64}, iz2::AbstractArray{Int64},
			   mesh_x::AbstractArray, mesh_z::AbstractArray, xval::Float64, zval::Float64)

	ix1[1], ix2[1] = indminn(mesh_x,xval,2)
	iz1[1], iz2[1] = indminn(mesh_z,zval,2)

	denomI[1] = ((mesh_x[ix2[1]] -mesh_x[ix1[1]])*(mesh_z[iz2[1]] - mesh_z[iz1[1]]))^(-1.e0)

	"iz1, ix1"
	weights[1] =	((mesh_x[ix2[1]]-	xval)*	(mesh_z[iz2[1]]-zval))*denomI[1]
	"iz1, ix2"
	weights[2] = ((xval-	mesh_x[ix1[1]])*(mesh_z[iz2[1]]-zval))*denomI[1]
	"recz2, ix1"
	weights[3] =	((mesh_x[ix2[1]]-	xval)*(zval-mesh_z[iz1[1]]))*denomI[1]
	"iz2, ix2"
	weights[4] =	((xval-	mesh_x[ix1[1]])*	(zval-	mesh_z[iz1[1]]))*	denomI[1]

end


text:
mutable struct RandomForest <: Classifier
    trees     :: Union{Nothing, Vector{DecisionTree}}
    topredict
    ntrees    :: UInt
    criterion
    attrfilter
end

RandomForest(;ntrees = 20, criterion = gini) = begin
    RandomForest(nothing, nothing, ntrees, criterion, nothing)
end

fit!(rf::RandomForest, df::DataFrame, topredict) = begin
    rf.topredict = topredict

    dfrows = [sample(1:nrow(df), nrow(df), replace = true) for i = 1:rf.ntrees]
    dfs = map(rows -> df[rows, :], dfrows)

    m = Int(ceil(1 + sqrt(size(df, 2) - 1)))
    rf.attrfilter = attrs -> sample(attrs, min(m, length(attrs)), replace = false)

    trees = Dict()
    Threads.@threads for i = 1:rf.ntrees
        id = Threads.threadid()
        t = DecisionTree(
            maxdepth = typemax(Int),
            criterion = rf.criterion,
            attrfilter = rf.attrfilter
        )
        fit!(t, dfs[i], topredict)
        trees[id] = [get(trees, id, []); [t]]
    end
    rf.trees = reduce((x,y) -> [x;y], values(trees))
    rf
end

predict(rf::RandomForest, df) = findmajority(map(t -> predict(t, df), rf.trees))


text:
module CovSelTest

using Test, Random, LinearAlgebra, Distributions
import CovSel, ProximalBase
using MathOptInterface
const MOI = MathOptInterface
using JuMP, SCS

Random.seed!(123)

@testset "random" begin
  for REP=1:10

      n = 1000
      p = 10

      ρ = 0.5
      covMat = Matrix(1.0I, p, p)
      for a = 1:p
        for b = a+1:p
          t = ρ^(b - a)
          covMat[a, b] = t
          covMat[b, a] = t
        end
      end
      precM = inv(covMat)

      sqCov = sqrt(covMat)
      X = randn(n, p) * sqCov

      S = X' * X / n

      λ = 0.5

      X = zeros(Float64, (p,p))
      Z = zeros(Float64, (p,p))
      U = zeros(Float64, (p,p))
      CovSel.covsel!(X, Z, U, S, λ; penalize_diag=true)

      # passing in verbose=0 to hide output from SCS
      problem = Model(JuMP.optimizer_with_attributes(SCS.Optimizer, "verbose" => 0))
      @variable(problem, Ω[1:p, 1:p], PSD)
      @variable(problem, lg_det)
      @variable(problem, B[1:p, 1:p])

      indOffDiag = Vector{Int64}(undef, 0)
      ind = 0
      for col = 1:p
          for row = 1:p
              ind += 1
              if row <= col
                  push!(indOffDiag, ind)
              end
          end
      end

      @constraint(problem, [lg_det; 1; Ω[indOffDiag]] in MOI.LogDetConeTriangle(p))
      @constraint(problem, Ω .<= B)
      @constraint(problem, -Ω .<= B)

      @objective(problem, Min, tr(Ω*S) - lg_det + λ * sum(B))
      optimize!(problem)

      Ωv = JuMP.value.(Ω)

      @test tr(Z*S) - logdet(Z) + λ * sum(abs, Z) ≈ tr(Ωv*S) - logdet(Ωv) + λ * sum(abs, Ωv) atol = 1e-3
      @test JuMP.value.(Ω) - Z ≈ zeros(p,p) atol = 1e-2
  end
end

@testset "non random" begin
    p = 10
    S = Matrix(1.0I, p, p)

    λ = 0.5

    X = zeros(Float64, (p,p))
    Z = zeros(Float64, (p,p))
    U = zeros(Float64, (p,p))
    CovSel.covsel!(X, Z, U, S, λ; penalize_diag=false)
    @test Z ≈ S atol=1e-3

    CovSel.covsel!(X, Z, U, S, λ; penalize_diag=false, options=CovSel.ADMMOptions(;abstol=1e-12,reltol=1e-12))
    @test Z ≈ S
end

end


text:
module Airplane

using NeuralNetworkAnalysis, LaTeXStrings
import Plots, DifferentialEquations
using Plots: plot, plot!, savefig, font, Measures.mm

# problem

Tψ = ψ -> [cos(ψ)  -sin(ψ)  0;
           sin(ψ)   cos(ψ)  0;
                0        0  1]

Tθ = θ -> [ cos(θ)  0  sin(θ);
                 0  1       0;
           -sin(θ)  0  cos(θ)]

Tϕ = ϕ -> [1       0        0;
           0  cos(ϕ)  -sin(ϕ);
           0  sin(ϕ)   cos(ϕ)]

Rϕθ = (ϕ, θ) -> [1  tan(θ) * sin(ϕ)  tan(θ) * cos(ϕ);
                 0           cos(θ)          -sin(ϕ);
                 0  sec(θ) * sin(ϕ)  sec(θ) * cos(ϕ)]

m = 1.0
g = 1.0
Ix = 1.0
Iy = 1.0
Iz = 1.0
Ixz = 0.0

@taylorize function airplane!(dx, x, p, t)
    _x, y, z, u, v, w, ϕ, θ, ψ, r, _p, q, Fx, Fy, Fz, Mx, My, Mz = x

    T_ψ = Tψ(ψ)
    T_θ = Tθ(θ)
    T_ϕ = Tϕ(ϕ)
    mat_1 = T_ψ * T_θ * T_ϕ
    xyz = mat_1 * vcat(u, v, w)

    mat_2 = Rϕθ(ϕ, θ)
    ϕθψ = mat_2 * vcat(_p, q, r)

    dx[1] = xyz[1]
    dx[2] = xyz[2]
    dx[3] = xyz[3]
    dx[4] = -g * sin(θ) + Fx / m - q * w + r * v
    dx[5] = g * cos(θ) * sin(ϕ) + Fy / m - r * u + _p * w
    dx[6] = g * cos(θ) * cos(ϕ) + Fz / m - _p * v + q * u
    dx[7] = ϕθψ[1]
    dx[8] = ϕθψ[2]
    dx[9] = ϕθψ[3]
    dx[10] = Mz
    dx[11] = My
    dx[12] = Mx
    dx[13] = zero(Fx)
    dx[14] = zero(Fy)
    dx[15] = zero(Fz)
    dx[16] = zero(Mx)
    dx[17] = zero(My)
    dx[18] = zero(Mz)
end

controller = read_nnet_sherlock(@modelpath "" "controller_Airplane_sherlock")

X₀ = Hyperrectangle(low=[ 0,       0,       0,       0.99999, 0.99999, 0.99999, 0.99999, 0.99999, 0.99999, 0,       0,       0],
                    high=[0.00001, 0.00001, 0.00001, 1,       1,       1,       1,       1,       1,       0.00001, 0.00001, 0.00001])
U₀ = ZeroSet(6)
vars_idx = Dict(:state_vars=>1:12, :control_vars=>13:18)
ivp = @ivp(x' = airplane!(x), dim: 18, x(0) ∈ X₀ × U₀)

period = 0.1
k = 20
T = k * period
T_warmup = 2 * period

prob = ControlledPlant(ivp, controller, vars_idx, period)

# simulation

println("simulation")
res = @timed simulate(prob, T=T, trajectories=3, include_vertices=false)
sim = res.value
print_timed(res)

# reachability analysis

alg = TMJets(abstol=1e-10, orderT=7, orderQ=1)
alg_nn = Ai2()

function benchmark(; T=T, silent::Bool=false)
    silent || println("flowpipe construction")
    res = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res.value
    silent || print_timed(res)

    return sol
end

benchmark(T=T_warmup, silent=true)
res = @timed benchmark()
sol = res.value

# plotting

vars = (2, 7)
fig = plot(xlab=L"y", ylab=L"\phi", leg=:topleft,
           legendfontsize=5,
           tickfont=font(25, "Times"),
           guidefontsize=25,
           xguidefont=font(25, "Times"),
           yguidefont=font(25, "Times"),
           bottom_margin=0mm, left_margin=0mm, right_margin=0mm, top_margin=0mm, size=(500, 450))
plot!(fig, sol, vars=vars, color=:orange, lab="")
plot!(project(X₀, vars), c=:blue, lab=L"X_0", alpha=0.6)
plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
savefig("airplane_yphi_juliareach.png")

end
nothing


text:
push!(LOAD_PATH,"../src/") # load module from local directory
using PolaronMobility

const T=300

#Ts,Kμs, Hμs, FHIPμs, ks, Ms, As, Bs, Cs, Fs, Taus
#effectivemass=0.12 # the bare-electron band effective-mass.
# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
println("Calculating polaron for T=$T for MAPIe params...")
MAPIe=polaronmobility(T, 4.5, 24.1, 2.25E12, 0.12)
#MAPIh=polaronmobility(T, 4.5, 24.1, 2.25E12, 0.15)

#function ImX(nurange,v,w,βred,α,ω,mb)

# Yup, this is a bit horrid.
v=MAPIe.v[1]
w=MAPIe.w[1]
βred=MAPIe.βred[1]
α=MAPIe.α[1]
ω=MAPIe.ω[1]
mb=MAPIe.mb[1]

nu=0:0.1:40
println("Integrating ImX for nu=$nu range...")
s=PolaronMobility.ImX(nu, v, w, βred,α,ω,mb)

println("Loading Plots for plotting...")
using Plots

plot( s.nu,s.ImX,label="ImX",
         markersize=3,marker=:downtriangle, xlab="nu (units Omega)",ylab="ImX")
xaxis!(:log10)
yaxis!(:log10)
savefig("MAPIe-ImX.png")
plot( s.nu,s.μ,label="mu",
         markersize=3,marker=:uptriangle, xlab="nu (units Omega)",ylab="Mob")
xaxis!(:log10)
yaxis!(:log10)

savefig("MAPIe-mu.png")

println("That's me!")


text:
"""
Extract expression sub-trees.
"""
module Match

export @match
export unkeyword!

iskw(expr::Expr) = expr.head == :kw
isref(expr::Expr) = expr.head == :ref
iscall(expr::Expr) = expr.head == :call

"""
Determines whether expression has commutative property.
"""
function iscommutative(op::Expr)
    iscall(op) && return iscommutative(op.args[1])
    isref(op) && return true
    false
end
iscommutative(op::Symbol) = _iscommutative(Val{op})

_iscommutative(::Type{Val{:(+)}}) = true
_iscommutative(::Type{T}) where T<:Val = false

"""
Output true if dictionary 'd' has a key 's' with a different value than 'v'.
Otherwise add value to dictionary and output false.
"""
function conflictadd!(d::Dict, s::Symbol, v)
    haskey(d,s) && (d[s] != v) && return true
    d[s] = v
    false
end

"""
Determine whether both expressions have the same head and arguments length.
"""
function partialmatch(expr::Expr, formula::Expr)
    samehead = (expr.head == formula.head)
    samelen = length(expr.args) == length(formula.args)
    (samehead & samelen) || return false
    true
end

"""
Match set of arguments from formula with expr.
"""
function match_args(symbols::Dict, expr::Expr, formula::Expr)
    offset = (iscall(expr) | isref(expr)) ? 1 : 0
    (!iscall(expr) | (expr.args[1] == formula.args[1])) || return false
    match_args(Val{iscommutative(expr)}, offset, symbols, expr, formula)
end
match_args(::Type{Val{false}}, kwargs...) = static_match(kwargs...)
match_args(::Type{Val{true}}, kwargs...) = commutative_match(kwargs...)

"""
Match each argument of formula with each argument of expr in an orderly fashion.
"""
function static_match(offset, symbols, expr, formula)
    eargs, margs = expr.args, formula.args
    (!isref(expr) | match(symbols, eargs[1], margs[1])) || return false
    for i in 1+offset:length(eargs)
        match(symbols, eargs[i], margs[i]) || return false
    end
    true
end

"""
Match each argument of formula with each argument of expr in any order.
"""
function commutative_match(offset, symbols, expr, formula)
    eargs, margs = expr.args, formula.args
    n = length(eargs)-offset
    (!isref(expr) | match(symbols, eargs[1], margs[1])) || return false
    indexes = collect(1+offset:length(eargs))
    perms = permutations(indexes)
    for perm in perms
        isempty(perm) && break
        success = true
        d = copy(symbols)
        for i in 1:n
            match(d, eargs[i+offset], margs[perm[i]]) || (success=false; break)
        end
        success && (merge!(symbols, d); return true)
    end
    false
end

"""
Match formula with expr. Overwrite dictionary with matched values.
"""
match(::Dict, expr, formula) = expr == formula
match(symbols::Dict, expr, s::Symbol) = !conflictadd!(symbols, s, expr)
function match(symbols::Dict, expr::Expr, formula::Expr)
    partialmatch(expr, formula) || return false
    match_args(symbols, expr, formula)
end
match(::Dict, expr, formula::Expr) = false

"""
Match formula with expr. Overwrite matched values directly to formula declared symbols.
"""
macro match(expr, formula)
    vars = getvars(formula)
    aux1 = "$formula"
    exec = quote end
    push!(exec.args, :(symbols = Dict{Symbol, Any}()))
    aux3 = esc(:($expr))
    push!(exec.args, :(success = $match(symbols, $aux3, Meta.parse($aux1))))
    for var in vars
        aux2 = "$var"
        aux3 = esc(:($var))
        push!(exec.args, :(success && ($aux3 = symbols[Meta.parse($aux2)])))
    end
    push!(exec.args, :(success))
    exec
end

"""
Get leaf symbols from expression.
"""
getvars(expr::Symbol) = _getvars(expr)
getvars(expr::Expr) = _getvars(unkeyword!(expr))

_getvars(::Any) = Set{Symbol}()
_getvars(s::Symbol) = Set{Symbol}([s])
function _getvars(expr::Expr)
    start = (iscall(expr) | iskw(expr)) ? 2 : 1
    union(map(_getvars, expr.args[start:end])...)
end

"""
Transform keyword to assignment.
"""
function unkeyword!(expr::Expr)
    iskw(expr) && (expr.head = :(=))
    expr
end

struct Permutations{T}
    a::T
    t::Int
end

Base.eltype(::Type{Permutations{T}}) where {T} = Vector{eltype(T)}

Base.length(p::Permutations) = (0 <= p.t <= length(p.a)) ? factorial(length(p.a), length(p.a)-p.t) : 0

"""
    permutations(a)
Generate all permutations of an indexable object `a` in lexicographic order. Because the number of permutations
can be very large, this function returns an iterator object.
Use `collect(permutations(a))` to get an array of all permutations.
"""
permutations(a) = Permutations(a, length(a))

"""
    permutations(a, t)
Generate all size `t` permutations of an indexable object `a`.
"""
function permutations(a, t::Integer)
    if t < 0
        t = length(a) + 1
    end
    Permutations(a, t)
end

function Base.iterate(p::Permutations, s = collect(1:length(p.a)))
    (!isempty(s) && max(s[1], p.t) > length(p.a) || (isempty(s) && p.t > 0)) && return
    nextpermutation(p.a, p.t ,s)
end

function nextpermutation(m, t, state)
    perm = [m[state[i]] for i in 1:t]
    n = length(state)
    if t <= 0
        return(perm, [n+1])
    end
    s = copy(state)
    if t < n
        j = t + 1
        while j <= n &&  s[t] >= s[j]; j+=1; end
    end
    if t < n && j <= n
        s[t], s[j] = s[j], s[t]
    else
        if t < n
            reverse!(s, t+1)
        end
        i = t - 1
        while i>=1 && s[i] >= s[i+1]; i -= 1; end
        if i > 0
            j = n
            while j>i && s[i] >= s[j]; j -= 1; end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i+1)
        else
            s[1] = n+1
        end
    end
    return (perm, s)
end

end


text:
using Base.Test
using WorldBankDataTd
using Dates

my_tests = ["ctrs_indicators.jl",
            "utils.jl",
            "wdi.jl",
            "format.jl",
            "search_wdi.jl",
            "doctests.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
function EventSlopeFilter(arg0::EventDetector, arg1::FilterType)
    return EventSlopeFilter((EventDetector, FilterType), arg0, arg1)
end

function g(obj::EventSlopeFilter, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function init(obj::EventSlopeFilter, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end


text:
#返回左右半链对应的序号，以Int的形式
function left_right_half_chain(the_state, tol_n)
    tr_bits = div(tol_n, 2)
    left_half::Int = (the_state >>> tr_bits)
    right_half::Int = (the_state - (left_half << tr_bits))
    # Julia数组从1开始
    left_half = left_half + 1
    right_half = right_half + 1
    left_half, right_half
end


text:
using DataStructure101
const DS = DataStructure101

export pre_order, in_order, pos_order

function pre_order(node::DS.TreeNode{T}, result::Vector{T}=T[]) where T
    push!(result, DS.value(node))
    pre_order(DS.left_child(node), result)
    pre_order(DS.right_child(node), result)

    return result
end

pre_order(node::DS.NullNode, result::Vector{T}=T[]) where T = result

function in_order(node::DS.TreeNode{T}, result::Vector{T}=T[]) where T
    in_order(DS.left_child(node), result)
    push!(result, DS.value(node))
    in_order(DS.right_child(node), result)

    return result
end

in_order(node::DS.NullNode, result::Vector{T}=T[]) where T = result

function pos_order(node::DS.TreeNode{T}, result::Vector{T}=T[]) where T
    pos_order(DS.left_child(node), result)
    pos_order(DS.right_child(node), result)
    push!(result, DS.value(node))

    return result
end

pos_order(node::DS.NullNode, result::Vector{T}=T[]) where T = result