text:
using Documenter, Provenance

makedocs(;
    modules = [Provenance],
    authors = "Tamas Gal",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        # assets = ["assets/logo.ico"],
    ),
    pages=[
        "Introduction" => "index.md",
        "API" => "api.md",
    ],
    repo="https://github.com/tamasgal/Provenance.jl/blob/{commit}{path}#L{line}",
    sitename="Provenance.jl",
)

deploydocs(;
    repo="github.com/tamasgal/Provenance.jl",
)


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
MLJ.color_off() # hide

KNNRegressor = @load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

ext: jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "Electrostatics Functional Tests" begin
    k = 9e9

    @testset "One particle rotates around another" begin
    # small mass with negative charge rotating around more massive object with positive charge

        r = 100.0
        q1 = 1e-3
        q2 = -1e-3
        m1 = 100.0
        m2 = 0.1
        v2 = sqrt(abs(k * q1 * q2 / m2 / r))
        t = 2 * pi * r / v2
        p1 = ChargedParticle(SVector(0.0, 0.0, 0.0), SVector(0.0, 0, 0.0), m1, q1)
        p2 = ChargedParticle(SVector(r, 0.0, 0.0), SVector(0.0, v2, 0.0), m2, q2)
        system = ChargedParticles([p1, p2], k)
        simulation = NBodySimulation(system, (0.0, t))
        sim_result = run_simulation(simulation)

        solution = sim_result.solution;
        ε = 0.1 * r
        for j = 1:2, i = 1:3
            @test solution[1][i,j] ≈ solution[end][i,j] atol = ε
        end

        (qs_act, ms_act, indxs_act, exclude_act) = NBodySimulator.obtain_data_for_electrostatic_interaction(simulation.system)
        @test qs_act[1] == q1 && qs_act[2] == q2
        @test ms_act[1] == m1 && ms_act[2] == m2
        @test length(qs_act) == length(ms_act)
    end

    @testset "Two positive charges repelling from each other" begin
    #   ("               <---⊕-----r1-----⊕--->                ")

        q1 = 1e-3 # C
        q2 = 1e-3 # C
        m1 = 1.0 # kg
        m2 = 1.0 # kg
        r1 = 1 # m

        E0 = k * q1 * q2 / r1 # initial energy
        t1 = 0.0
        t2 = 1.0
        τ = (t2 - t1) / 1000

        p1 = ChargedParticle(SVector(-r1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m1, q1)
        p2 = ChargedParticle(SVector(r1 / 2, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m2, q2)
        system = ChargedParticles([p1, p2], k)
        simulation = NBodySimulation(system, (t1, t2))
        sim_result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        r2 = get_position(sim_result, t2, 2) - get_position(sim_result, t2, 1)
        v_expected = sqrt(k * q1 * q2 / m1 * (1 / norm(r1) - 1 / norm(r2)))
        v_actual = norm(get_velocity(sim_result, t2, 2))

        ε = 0.001 * v_expected
        @test v_expected ≈ v_actual atol = ε
    end

    @testset "Constructing electorstatic potential parameters entity" begin
        default_potential_1 = ElectrostaticParameters()
        @test 9e9 == default_potential_1.k

        default_potential_2 = ElectrostaticParameters(k, 8.0)
        @test 64.0 == default_potential_2.R2

        io = IOBuffer()
        potential1 = ElectrostaticParameters()
        potential2 = ElectrostaticParameters(1.0)
        @test sprint(io -> show(io, potential1)) == "Electrostatic:\n\tk:9.0e9\n"
        @test sprint(io -> show(io, potential2)) == "Electrostatic:\n\tk:1.0\n"
    end

    @testset "Energy conservation test" begin
        n = 8
        bodies = ChargedParticle[]
        L = 1.0
        m = 1.0
        q = 1.0
        count = 1
        dL = L / (ceil(n^(1 / 3)) + 1)
        for x = dL / 2:dL:L, y = dL / 2:dL:L, z = dL / 2:dL:L
            if count > n
                break
            end
            r = SVector(x, y, z)
            v = SVector(.0, .0, .0)
            body = ChargedParticle(r, v, m, q)
            push!(bodies, body)
            count += 1
        end

        k = 9e9
        τ = 0.01 * dL / sqrt(2 * k * q * q / (dL * m))
        t1 = 0.0
        t2 = 1000 * τ

        potential = ElectrostaticParameters(k, 0.45 * L)
        system = PotentialNBodySystem(bodies, Dict(:electrostatic => potential))
        pbc = CubicPeriodicBoundaryConditions(L)
        simulation = NBodySimulation(system, (t1, t2), pbc)
        result = run_simulation(simulation, VelocityVerlet(), dt=τ)

        e_tot_1 = total_energy(result, t1)
        e_tot_2 = total_energy(result, t2)

        ε = 0.001
        @test (e_tot_2 - e_tot_1) / e_tot_1 ≈ 0.0 atol = ε
    end
end


text:
using MLIR
using Test

@testset "MLIR" begin
end

include("examples.jl")


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Convenience functions for syntactic parsing.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-07

# TODO: Each perceptron could have its own fidx dictionary, sparser?
# TODO: If so, then they would also need to be independently extended in size.

module Parse

export done, isequal, length, next, parse, start, train

require("structs.jl")
require("conllx.jl")
require("dep.jl")
require("hybrid.jl")
require("percep.jl")

using Structs
using CoNLLX
using DepGraph
using Hybrid
using Percep

function parse(sent::Sentence, predict::Function, coder)
    conf = config(sent, coder)
    while !isterminal(conf)
        trans = predict(conf)
        apply!(conf, trans)
    end

    return sentence(conf.graph, coder)
end

# TODO: Is it a good idea to override these?
# TODO: This is really a parsing specifing `isequal` implementation.
import Base: isequal
function isequal(a::Token, b::Token)
    return (a.id == b.id && a.form == b.form && a.head == b.head
        && a.deprel == b.deprel)
end

function isequal(a::Sentence, b::Sentence)
    # TODO: Prettify using map etc.
    if length(a) != length(b)
        return false
    end

    for (ae, be) in zip(a, b)
        if !isequal(ae, be)
            return false
        end
    end

    return true
end

function predict(conf, model, featids)
    legaltranss = collect(filter(t -> islegal(conf, t), model.transs))

    if length(legaltranss) == 1
        # There is only one possibly transition, use it.
        return legaltranss[1]
    end

    bestscore = -Inf
    besttrans = nothing
    for trans in legaltranss
        score = testeval(model, trans, conf, featids)
        if score > bestscore
            bestscore = score
            besttrans = trans
        end
    end

    return besttrans
end

function parse(sent, model, featids)
    conf = config(sent, model.coder)
    while !isterminal(conf)
        trans = predict(conf, model, featids)
        apply!(conf, trans)
    end

    return sentence(conf.graph, model.coder)
end

# TODO: Implement LAS as well.
export evaluate
function evaluate(model, goldsents)
    # Cache the feature id vector between parses.
    featids = Array(Uint, NUMFEATS)
    total = 0
    correct = 0
    for goldsent in goldsents
        predsent = parse(goldsent, model, featids)
        for (goldtok, predtok) in zip(goldsent, predsent)
            total += 1
            if predtok.head == goldtok.head
                correct += 1
            end
        end
    end

    if correct < 1
        return 0.0
    else
        return correct / total
    end
end

function testeval(model, trans, conf, featids)
    # Short-hands.
    weights = model.percepbytrans[trans].weights
    bias = model.percepbytrans[trans].bias
    fidbyfstruct = model.fidbyfstruct

    # Apply the transition.
    undo = apply!(conf, trans)
    # Featurise the resulting configuration.
    featids = featurise!(featids, fidbyfstruct, conf)

    score = bias
    for featid in featids
        # Only consider previously observed features.
        if featid > 0
            score += weights[featid]
        end
    end

    # Revert to the original configuration.
    undo!(conf, undo)

    return score
end

function traineval(cache, model, trans, conf)
    # Short-hands.
    percep = model.percepbytrans[trans]
    fidbyfstruct = model.fidbyfstruct
    featids = cache.fidsbytrans[trans]

    # Apply the transition.
    undo = apply!(conf, trans)
    # Featurise the resulting configuration.
    featurise!(featids, fidbyfstruct, conf)

    # Revert to the original configuration.
    undo!(conf, undo)

    # Extend the weight vector if we observed new features.
    numfeats = length(fidbyfstruct)
    if length(percep.weights) < numfeats
        grow!(model, numfeats)
    end

    # Calculate feats' * weights + bias.
    spscore(featids, percep)
end

function trainpred(train, model, conf)
    goldtrans = oracle(conf)
    legaltranss = collect(filter(t -> islegal(conf, t), model.transs))

    if length(legaltranss) == 1
        besttrans = legaltranss[1]
        # There is only one possibly transition, no need to evaluate.
        return (goldtrans, besttrans)
    end

    # TODO: Can this be made prettier and type stable?
    besttrans = nothing
    bestscore = -Inf
    for trans in legaltranss
        score = traineval(train.cache, model, trans, conf)

        if !isequal(trans, goldtrans)
            # Add the margin.
            # TODO: Extract into an argument.
            score += 1
        end

        if score > bestscore
            bestscore = score
            besttrans = trans
        end
    end

    return (goldtrans, besttrans)
end

type TrainModel
    # XXX: I dislike this... how can we fix it?
    #   type(Symbol, Label)?
    transs
    percepbytrans::Dict{Any,AvgPerceptron}
    fidbyfstruct::Identifier{FeatStruct,Uint}
    coder::Coder
end

function TrainModel()
    transs = transitions()
    dic = Dict{Any,AvgPerceptron}()
    for trans in transs
        dic[trans] = AvgPerceptron()
    end
    TrainModel(transs, dic, Identifier(FeatStruct, Uint), Coder())
end

function grow!(model::TrainModel, size)
    for (_, p) in model.percepbytrans
        extend!(p, size)
    end
    model
end

function update!(m, goldtrans, goldfids, besttrans, bestfids)
    spupdate!(m.percepbytrans[goldtrans], goldfids, true)
    spupdate!(m.percepbytrans[besttrans], bestfids, false)
end

type Model
    transs
    percepbytrans::Dict{Any,Perceptron}
    fidbyfstruct::Identities{FeatStruct,Uint}
    coder::Coder
end

# Finalise a model in training to construct the resulting model.
function Model(m)
    dic = Dict{Any,Perceptron}()
    for (trans, p) in m.percepbytrans
        dic[trans] = Perceptron(p)
    end
    # TODO: Should be a deep copy of transs, fidbyfstruct and coder?
    #   This slows down the training though...
    Model(m.transs, dic, Identities(m.fidbyfstruct), m.coder)
end

import Base: length
function length(model::Union(TrainModel,Model))
    for (_, p) in model.percepbytrans
        # We intentionally short-circuit since all perceptrons have the same
        #   number of weights.
        return length(p.weights)
    end
end

type Cache
    fidsbytrans::Dict{Any,Vector{Uint}}
end

function cache(model)
    # Pre-allocate for each transition.
    fidsbytrans = Dict{Any,Vector{Uint}}()
    for trans in model.transs
        fidsbytrans[trans] = Array(Uint, NUMFEATS)
    end

    return Cache(fidsbytrans)
end

type Train
    epoch::Uint
    epochs::Uint
    sents::Sentences
    trainmodel::TrainModel
    cache::Cache
    earlyupdates::Bool
end

function train(sents; epochs=0, earlyupdates=false)
    tmodel = TrainModel()
    Train(0, epochs, sents, tmodel, cache(tmodel), earlyupdates)
end

import Base: start
start(::Train) = nothing

import Base: next
function next(itr::Train, nada)
    itr.epoch += 1

    # Short-hands.
    tmodel = itr.trainmodel
    cache = itr.cache

    # Shuffle the sentences before for each pass.
    shuffle!(itr.sents)
    for sent in itr.sents
        # Parse the sentence.
        conf = config(sent, tmodel.coder)
        while !isterminal(conf)
            goldtrans, besttrans = trainpred(itr, tmodel, conf)

            if !isequal(besttrans, goldtrans)
                goldfids = itr.cache.fidsbytrans[goldtrans]
                bestfids = itr.cache.fidsbytrans[besttrans]
                # Update the weight vectors.
                update!(tmodel, goldtrans, goldfids, besttrans, bestfids)

                if itr.earlyupdates
                    # Skip to the next sentence.
                    break
                end
            end

            # Proceed along the gold path.
            apply!(conf, goldtrans)
        end
    end

    return (Model(tmodel), nothing)
end

# Helper function to more naturally progress the iterator by masking the
#   internals of the iterator protocol.
next(itr::Train) = next(itr, nothing)[1]

import Base: done
done(itr::Train, nada) = itr.epochs != 0 && itr.epoch >= itr.epochs

end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Directories
toastn = "2.0.2"
toastv = "v" * toastn
prefix = dirname(@__FILE__)
mkpath(prefix)

# Get source and binaries
if Sys.isapple()
  osroot = "darwin"
  oslibn = ".dylib"
end

if Sys.iswindows()
  osroot = "windows"
  oslibn = ".dll"
end

if Sys.islinux()
  osroot = "linux"
  oslibn = ".so"
end

suffix = ".zip"

srcroot = "https://github.com/toastpp/toastpp/archive/"
srcfile = toastv * suffix
srcdir = "toastpp-" * toastn

@info "Downloading Toast++ source"
download(srcroot * srcfile, joinpath(prefix, srcfile))

osroot = (Sys.isapple() ? "darwin" : Sys.islinux() ? "linux" : "windows")
dlroot = "https://github.com/toastpp/toastpp/releases/download/" * toastv * "/"
dlfile = "toast_bin_" * osroot * "64" * suffix

@info "Downloading Toast++ library binaries"
download(dlroot * dlfile, joinpath(prefix, dlfile))

# Unzip
@info "Uncompressing Toast++ source and binaries"
if Sys.islinux()
  # On linux, unzip binaries into toast subdirectory to match OS X
  run(`unzip -q -o $(joinpath(prefix, dlfile)) -d $(joinpath(prefix, "toast"))`)
else
  run(`unzip -q -o $(joinpath(prefix, dlfile)) -d $prefix`)
end
run(`unzip -q -o $(joinpath(prefix, srcfile)) -d $prefix`)

# Move the binaries into the source tree
@info "Adding binaries to source tree"
run(`cp -r $(joinpath(prefix, "toast", osroot * "64" )) $(joinpath(prefix, "toastpp-" * toastn))`)
run(`rm -rf $(joinpath(prefix, osroot * "64"))`)

# Remove temporary download
@info "Deleteing temporary files"
rm(joinpath(prefix, dlfile))
rm(joinpath(prefix, srcfile))

# Set rpath on linux
if Sys.islinux()
  libroot = joinpath(prefix, "toastpp-" * toastn, "linux64", "lib")
  libs = filter(x-> occursin(".so", x), readdir(libroot))
  for lib in libs
    run(`patchelf --set-rpath $libroot $(joinpath(libroot, lib))`)
    @info "Patched rpath for library $lib"
  end
end

# Write path
@info "Storing Toast++ installation path"
libprefix = joinpath(prefix, "toastpp-" * toastn, osroot * "64", "lib");

s = """
const _jl_toast_dir = "$(joinpath(prefix, "toastpp-" * toastn))"
const _jl_toast_libdir =  "$libprefix"
const _jl_toast_libfe = "$(joinpath(libprefix, "libfe" * oslibn))"
const _jl_toast_libmath = "$(joinpath(libprefix, "libmath" * oslibn))"
const _jl_toast_libstoast = "$(joinpath(libprefix, "libstoast" * oslibn))"
const _jl_toast_libsuperlu = "$(joinpath(libprefix, "libsuperlu" * oslibn))"
"""

f = open(joinpath(prefix, "path.jl"), "w")
write(f, s)
close(f)

# Fin
@info "Installation complete"


text:
import Tables: rows, columns, istable, rowaccess, columnaccess, schema, Schema

import Base: propertynames, getproperty

export istable

istable(::Type{JDFFile}) = true
istable(::JDFFile) = true

rowaccess(::JDFFile) = false
columnaccess(::JDFFile) = true

rowaccess(::Type{<:JDFFile}) = false
columnaccess(::Type{<:JDFFile}) = true

propertynames(jdf::JDFFile) = names(jdf)

getproperty(jdf::JDFFile, col::Symbol) = jdf[!, col]

schema(jdf::JDFFile) = begin
    meta = metadata(jdf)
    Schema(meta.names, map(x -> x.type, meta.metadatas))
end

columns(jdf::JDFFile) = jdf


text:
using MonteCarloX
using Random
using StatsBase
using Printf
using ProgressMeter
using LinearAlgebra
using DelimitedFiles
using StaticArrays

"""
# generate a stochastic trajectory
SIR dynamics is goverened by differential equation
``
    \\frac{dI}{dt} = \\sum_{i=1}^{I}\\lambda_i\\frac{S}{N} - \\mu I
``

# Arguments
* distribution for the heterogeneous rates P(lambda)
"""
function trajectory!(rng, list_time, list_I, system)
    rates = current_rates(system)
    pass_update!(rates, event) = update!(rates, event, system, rng)

    list_I[1] = system.measure_I
    time_simulation = Float64(list_time[1])
    for i in 2:length(list_time)
        if time_simulation < list_time[i]
            dtime = list_time[i] - time_simulation
            dtime_sim = advance!(KineticMonteCarlo(), rng, rates, pass_update!, dtime)
            time_simulation += dtime_sim
        end
        list_I[i] = system.measure_I
    end
end

"""
Ingreedients:
- the rate of each agent to become infected or to recover
-> the rate of infected to make a connection with any other agent (well-mixed),
   which adds to the connection matrix, and only infects with a certain
   probability
   !! Connection matrix needs to be deleted after 2 weeks -> LightGraphs plus a
   queue where we store the edges and the time when stored. (clean the graph whenever it gets updated)

- the rate to be tested for those that will develop symptoms (probability)
    -> as soon as S->I in update we draw a random number, and if the id is symptom based then add this rate to the system
- the rate to be tested randomly
-> upon positive test: add contacts to a list of tracing
- the rate to test people from tracing list
"""
mutable struct SIR
    rng::AbstractRNG
    # rates
    epsilon::Float64 # rate of susceptible individual getting infected from outside
    mu::Float64      # rate of infected individual to recover
    lambda::Float64  # rate of infected individual to spread the disease
    nu::Float64      # rate of any individual to be in contact with anyone else
    omega_r::Float64 # rate of random test on any population individual
    omega_s::Float64 # rate of symptomatic-driven test on symptomatic
    # probabilities
    p_asymptomatic::Float64
    p_infection::Float64
    p_prevalence::Float64
    # compartments
    N::Int          # population size (S+I+R)
    S::Int          # susceptible
    I::Int          # infected (=H_a+H_s+T)
    H_a::Int        # hidden asymptomatic infected
    H_s::Int        # hidden symptomatic infected
    T::Int          # traced infected
    R::Int          # recovered
    measure_S::Int
    measure_I::Int
    measure_R::Int

    function SIR(rng::AbstractRNG, mu::Number, lambda::Number, nu::Number, omega_r::Number, omega_s::Number;
                 epsilon::Number = 0,
                 p_asymptomatic::Number = 0.5, # wise knowledge of Seba, reference will follow
                 p_prevalence::Number = 3e-3, # wise knowledge of Seba, reference will follow
                 N::Int = Int(1e4),
                )
        H_a = H_s = T = 0
        if rand(rng) < p_asymptomatic
            H_a = 1
        else
            H_s = 1
        end
        R0 = rand(rng, 0:2*N*p_prevalence)
        I0 = H_a + H_s + T
        S0 = N - I0 - R0
        epsilon = Float64(epsilon)
        mu      = Float64(mu)
        lambda  = Float64(lambda)
        nu      = Float64(nu)
        omega_r = Float64(omega_r)
        omega_s = Float64(omega_s)
        p_infection = lambda/nu

        new(rng, epsilon, mu, lambda, nu, omega_r, omega_s, p_asymptomatic, p_infection, p_prevalence, N, S0, I0, H_a, H_s, T, R0, S0, I0, R0)
    end
end

"""
update rates and system according to the last event that happend:
1: recovery
2: contact internal
3: infection external
4: random test
5: symptom-driven test
"""
function update!(rates::AbstractVector, event::Int, system::SIR)
    system.measure_S = system.S
    system.measure_I = system.I
    system.measure_R = system.R
    if event == 1 # recovery (from any pool of infected individuals)
        i = sample(system.rng, [1,2,3], AbstractWeights([system.H_a, system.H_s, system.T]./system.I))
        if i==1
            H_a -= 1
        elseif i==2
            H_s -= 1
        else
            T -= 1
        end
        system.I -= 1
        system.R += 1
    elseif event == 2 # infection (internal, through contacts) these go randomly to s or a (we need ar random number)
        i = sample(system.rng, [1,2], AbstractWeights([system.p_asymptomatic, 1-system.p_asymptomatic]))
        if i==1
            H_a += 1
        else i==2
            H_s += 1
        end
        system.S -= 1
        system.I += 1
    elseif event == 3 # infection (external)
        i = sample(system.rng, [1,2], AbstractWeights([system.p_asymptomatic, 1-system.p_asymptomatic]))
        if i==1
            H_a += 1
        else i==2
            H_s += 1
        end
        system.S -= 1
        system.I += 1
    elseif event == 4 # random testing (from or(H_a,H_s) to T)
        i = sample(system.rng, [1,2], AbstractWeights([system.H_a, system.H_s]./(system.I-system.T)))
        if i==1
            H_a -= 1
        else i==2
            H_s -= 1
        end
        system.T += 1
    elseif event == 5 # symptom-driven testing (from H_s to T)
        H_s      -= 1
        system.T += 1
    elseif event == 0 # absorbing state of zero infected
        system.measure_S = system.S
        system.measure_I = system.I
        system.measure_R = system.R
    else
        throw(UndefVarError(:index))
    end

    rates .= current_rates(system)
end

"""
evaluate current rates of SIR system
"""
function current_rates(system::SIR)
    rate_recovery = system.mu * system.I
    rate_contact_internal   = system.nu * system.I
    rate_infection_external = system.epsilon * system.S
    # random tests can be restricted in simulation to only testing infected and
    # then simply to occur less often
    rate_test_random  = system.omega_r * (system.H_s + system.H_a)
    rate_test_symptom = system.omega_s * system.H_s
    return MVector{5,Float64}(rate_recovery, rate_contact_internal, rate_infection_external, rate_test_random, rate_test_symptom)
end

#function delete_from_system!(system, index)
#    update_sum!(system, -system.current_lambda[index])
#    deleteat!(system.current_lambda, index)
#end
#
#function add_to_system!(system, lambda)
#    push!(system.current_lambda, lambda)
#    update_sum!(system, lambda)
#end
#
#function update_sum!(system, change)
#    system.sum_current_lambda += change
#    system.update_current_lambda += 1
#    if system.update_current_lambda > 1e5
#        system.sum_current_lambda = sum(system.current_lambda)
#        system.update_current_lambda = 0
#    end
#end


text:
include("FourierTransforms.jl/FourierTransforms.jl")

# This is needed to flag that the fft_generic.jl file has already been loaded
const GENERIC_FFT_LOADED = true

@info("Code paths for generic floating-point types activated in DFTK. Remember to add " *
      "'using GenericLinearAlgebra' to your user script in case not yet done. " *
      "See https://docs.dftk.org/stable/examples/arbitrary_floattype/ for details.")

# Utility functions to setup FFTs for DFTK. Most functions in here
# are needed to correct for the fact that FourierTransforms is not
# yet fully compliant with the AbstractFFTs interface and has still
# various bugs we work around.

function next_working_fft_size(::Any, size)
    # TODO FourierTransforms has a bug, which is triggered
    #      only in some factorizations, see
    #      https://github.com/JuliaComputing/FourierTransforms.jl/issues/10
    # To be safe we fall back to powers of two

    adjusted = nextpow(2, size)
    if adjusted != size
        @info "Changing fft size to $adjusted (smallest working size for generic FFTs)"
    end
    adjusted
end

# Generic fallback function, Float32 and Float64 specialization in fft.jl
function build_fft_plans(T, fft_size)
    tmp = Array{Complex{T}}(undef, fft_size...)

    # Note: FourierTransforms has no support for in-place FFTs at the moment
    # ... also it's extension to multi-dimensional arrays is broken and
    #     the algo only works for some cases
    @assert all(ispow2, fft_size)

    # opFFT = FourierTransforms.plan_fft(tmp)   # TODO When multidim works
    # opBFFT = inv(opFFT).p
    opFFT  = generic_plan_fft(tmp)               # Fallback for now
    opBFFT = generic_plan_bfft(tmp)
    # TODO Can be cut once FourierTransforms supports AbstractFFTs properly
    ipFFT  = DummyInplace{typeof(opFFT)}(opFFT)
    ipBFFT = DummyInplace{typeof(opBFFT)}(opBFFT)

    ipFFT, opFFT, ipBFFT, opBFFT
end

struct GenericPlan{T}
    subplans
    factor::T
end

function generic_apply(p::GenericPlan, X::AbstractArray)
    pl1, pl2, pl3 = p.subplans
    ret = similar(X)
    for i in 1:size(X, 1), j in 1:size(X, 2)
        @views ret[i, j, :] .= pl3 * X[i, j, :]
    end
    for i in 1:size(X, 1), k in 1:size(X, 3)
        @views ret[i, :, k] .= pl2 * ret[i, :, k]
    end
    for j in 1:size(X, 2), k in 1:size(X, 3)
        @views ret[:, j, k] .= pl1 * ret[:, j, k]
    end
    p.factor .* ret
end

LinearAlgebra.mul!(Y, p::GenericPlan, X) = Y .= p * X
LinearAlgebra.ldiv!(Y, p::GenericPlan, X) = Y .= p \ X

import Base: *, \, inv, length
length(p::GenericPlan) = prod(length, p.subplans)
*(p::GenericPlan, X::AbstractArray) = generic_apply(p, X)
*(p::GenericPlan{T}, fac::Number) where T = GenericPlan{T}(p.subplans, p.factor * T(fac))
*(fac::Number, p::GenericPlan{T}) where T = p * fac
\(p::GenericPlan, X) = inv(p) * X
inv(p::GenericPlan{T}) where T = GenericPlan{T}(inv.(p.subplans), 1 / p.factor)

function generic_plan_fft(data::AbstractArray{T, 3}) where T
    GenericPlan{T}([FourierTransforms.plan_fft(data[:, 1, 1]),
                    FourierTransforms.plan_fft(data[1, :, 1]),
                    FourierTransforms.plan_fft(data[1, 1, :])], T(1))
end
function generic_plan_bfft(data::AbstractArray{T, 3}) where T
    GenericPlan{T}([FourierTransforms.plan_bfft(data[:, 1, 1]),
                    FourierTransforms.plan_bfft(data[1, :, 1]),
                    FourierTransforms.plan_bfft(data[1, 1, :])], T(1))
end

# A dummy wrapper around an out-of-place FFT plan to make it appear in-place
# This is needed for some generic FFT implementations, which do not have in-place plans
struct DummyInplace{opFFT}
    fft::opFFT
end
LinearAlgebra.mul!(Y, p::DummyInplace, X) = (Y .= mul!(similar(X), p.fft, X))
LinearAlgebra.ldiv!(Y, p::DummyInplace, X) = (Y .= ldiv!(similar(X), p.fft, X))

import Base: *, \, length
*(p::DummyInplace, X) = p.fft * X
\(p::DummyInplace, X) = p.fft \ X
length(p::DummyInplace) = length(p.fft)


text:
# Autogenerated wrapper script for OpenJpeg_jll for i686-w64-mingw32
export libopenjp2, opj_compress, opj_decompress, opj_dump

using LittleCMS_jll
using libpng_jll
using Libtiff_jll
JLLWrappers.@generate_wrapper_header("OpenJpeg")
JLLWrappers.@declare_library_product(libopenjp2, "libopenjp2.dll")
JLLWrappers.@declare_executable_product(opj_compress)
JLLWrappers.@declare_executable_product(opj_decompress)
JLLWrappers.@declare_executable_product(opj_dump)
function __init__()
    JLLWrappers.@generate_init_header(LittleCMS_jll, libpng_jll, Libtiff_jll)
    JLLWrappers.@init_library_product(
        libopenjp2,
        "bin\\libopenjp2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        opj_compress,
        "bin\\opj_compress.exe",
    )

    JLLWrappers.@init_executable_product(
        opj_decompress,
        "bin\\opj_decompress.exe",
    )

    JLLWrappers.@init_executable_product(
        opj_dump,
        "bin\\opj_dump.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# bnatt350.jl
#
# Solve a nontrivial MIP problem by injecting the optimal solution. The LP
# relaxation has the optimal objective value, so optimality is proven
# immediately. The JuMP model is populated from an .MPS file by using
# MathProgBase to construct another model and then copying the data over to
# a JuMP model. The problem comes from the MIPLIB 2010 benchmark collection.
#############################################################################

using JuMP
using MathProgBase

## Uncomment the following three lines to solve using Gurobi
# using Gurobi
# mod = Model(solver=GurobiSolver(PreCrush=1, Cuts=0, Presolve=0, DisplayInterval=1))
# m_internal = MathProgBase.LinearQuadraticModel(GurobiSolver())

using CPLEX
mod = Model(solver=CplexSolver(CPX_PARAM_PRELINEAR=0, CPX_PARAM_PREIND=0, CPX_PARAM_ADVIND=0, CPX_PARAM_MIPSEARCH=1,CPX_PARAM_MIPCBREDLP=0))
m_internal = MathProgBase.LinearQuadraticModel(CplexSolver())

## Uncomment the following three lines to solve using GLPK
# using GLPKMathProgInterface
# mod = Model(solver=GLPKSolverMIP())
# m_internal = MathProgBase.LinearQuadraticModel(GLPKSolverMIP())

# Load the model from .MPS file
MathProgBase.loadproblem!(m_internal, "data/bnatt350.mps")

# grab MathProgBase data
c = MathProgBase.getobj(m_internal)
A = MathProgBase.getconstrmatrix(m_internal)
m, n = size(A)
xlb = MathProgBase.getvarLB(m_internal)
xub = MathProgBase.getvarUB(m_internal)
l = MathProgBase.getconstrLB(m_internal)
u = MathProgBase.getconstrUB(m_internal)
vtypes = MathProgBase.getvartype(m_internal)

# populate JuMP model with data from internal model
@variable(mod, x[1:n])
for i in 1:n
    setlowerbound(x[i], xlb[i])
    setupperbound(x[i], xub[i])
    mod.colCat[x[i].col] = vtypes[i]
end
At = A' # transpose to get useful row-wise sparse representation
for i in 1:At.n
    @constraint( mod, l[i] <= sum{ At.nzval[idx]*x[At.rowval[idx]], idx = At.colptr[i]:(At.colptr[i+1]-1) } <= u[i] )
end
@objective(mod, Min, sum{ c[i]*x[i], i=1:n })

function myheuristic(cb)
    fp = open("data/bnatt350.sol", "r")
    readline(fp)
    for i in 1:n
        line = chomp(readline(fp))
        spl = filter(x->!isempty(x), split(line, " "))
        setsolutionvalue(cb, x[i], round(float(spl[2]),4))
    end
    addsolution(cb)
end  # End of callback function

addheuristiccallback(mod, myheuristic)
stat = solve(mod)
println("Solve status: ", stat)
println("Objective value: ", getobjectivevalue(mod))


text:
using Plots, Measures, HTTP, Serialization, Random

#These four convenience functions extract the state variable from the state vector
#It is assumed the layout of the vector u is u = [v_x, v_y, x, y]
state_v_x(u::Vector{Float64}) = u[1]
state_v_y(u::Vector{Float64}) = u[2]
state_x(u::Vector{Float64}) = u[3]
state_y(u::Vector{Float64}) = u[4]

"""
Modified version of the given function. Takes a Vector of tuples of (Float, Vector), first element being h and the second being the values corresponding to that h
"""
function plot_solutions( t::Int64,
    u::Vector{Tuple{Float64,Vector{Vector{Float64}}}};
    title::String = "",
    label::Union{String, Bool} = false,
    show_energy = true,
    combine_positions = false) where T

    position_plots = [] # Storing plots
    energy_values = [] # Storing values to plot later

    # Iterates through given sets of values, does mostly the same thing as the given function
    for (h, u_vec) in u
        x, y, v_x, v_y = state_x.(u_vec), state_y.(u_vec), state_v_x.(u_vec), state_v_y.(u_vec)

        r = @. sqrt(x^2 + y^2)
        E = @. 0.5*(v_x^2 + v_y^2) - 1.0/r

        # Pushes to plot list instead of plotting outright
        push!(position_plots, plot(  reduce_data_points([x], [y])..., label = label, xlabel= "X", ylabel = "Y",
        title = title*" (position): h = $h", aspectratio=1,legend=:topleft,ylim=(-7,7)))
        scatter!([0], [0], ms=15, msw=0, c=:orange, shape =:star, label="Sun")
        scatter!([x[1]], [y[1]], ms=4, msw=0, c=:blue, shape =:circle, label="Earth initial position")

        # Push values without plotting
        push!(energy_values, (h, E))
    end

    # Plot all energy values all at once
    x_values_full, y_values_full = map(x->1:x[1]:t, energy_values), map(x->x[2], energy_values)
    p4 = plot(  reduce_data_points(x_values_full, y_values_full)..., xlabel = "Time", ylabel = "Energy",
    label = transpose(map(x->x[1], energy_values)), title = title*" (energy)")

    # Exclude the energy plot if not needed
    plots_to_plot = position_plots
    if show_energy
        push!(plots_to_plot, p4)
    end

    # Collate all plots into the apex plot
    return plot(plots_to_plot..., margin = 10mm,size=(1000,1000))
end;

"""
"Simple" function for reducing the amount of data points in a graph. Both args must be vectors of vectors.
"""
function reduce_data_points(x_values_full, y_values_full)
    x_values_truncated, y_values_truncated = [], []
    for i in eachindex(x_values_full)
        len = length(x_values_full[i])
        randsub = randsubseq(1:len, min(1000, len)/len)
        push!(x_values_truncated, x_values_full[i][randsub])
        push!(y_values_truncated, y_values_full[i][randsub])
    end
    return x_values_truncated, y_values_truncated
end

# Using the same starting values as given
h_list = map(x->10^(-x), 2.0:4.0)
t_max = 200
u_0 = [0., 1, 1.5, 0]

"""
Computes the RHS for the one body problem (the t arg has been removed)
"""
function df_dt_one_body(u::Vector{Float64})::Vector{Float64}
    M, G = 1, 1
    r = sqrt(state_x(u)^2 + state_y(u)^2)
    return [-M*G*state_x(u)/r^3, -M*G*state_y(u)/r^3, state_v_x(u), state_v_y(u)]
end;

"""
Performs euler method. f is a function which takes a single arg of same type as u_0
"""
function euler_method(f, u_0, step_size, t_vec)
    values = [u_0]
    previous = u_0
    for i in 2:length(t_vec)
        next = previous .+ step_size * f(previous)
        previous = next
        push!(values, next)
    end
    return values
end

"""
Performs RK4 method. f takes a single argument of same type as u_0
"""
function rk4_method(f, u_0, step_size, t_vec)
    values = [u_0]
    prev_value = u_0
    # Could go at it recursively, but I also don't want to make my code even more unreadable
    for i in 2:length(t_vec)
        next_value = prev_value .+ rk4_calculate_section(f, prev_value, step_size)
        push!(values, next_value)
        prev_value = next_value
    end
    return values
end

"""
Helper function for the cumbersome calculation. Args must be vectors.
"""
function rk4_calculate_section(f, current_elem, step_size)
    k_1 = f(current_elem)
    k_2 = f(current_elem .+ ((step_size/2) .* k_1))
    k_3 = f(current_elem .+ ((step_size/2) .* k_2))
    k_4 = f(current_elem .+ (step_size .* k_3))
    return (step_size/6) .* (k_1 .+ (2 .* k_2) .+ (2 .* k_3) .+ k_4)
end

# The part where the graphs get drawn

display(plot_solutions(t_max, map(h -> (h, euler_method(df_dt_one_body, u_0, h, 1:h:t_max)), h_list), title="Euler's Method"))

h_best, h_worst = 0.01, 0.01 # for RK4 and Euler respectively

anim = @animate for i = 30:max(t_max÷50, 1):t_max
    plot_solutions(i, [(h_best, rk4_method(df_dt_one_body, u_0, h_best, 1:h_best:i)), (h_worst, euler_method(df_dt_one_body, u_0, h_worst, 1:h_worst:i))], title="Combined animation", show_energy = false)
end

gif(anim, "tutorial_anim_fps30.gif", fps = 30)

h_list_rk4 = [0.01, 0.1, 0.5]
display(plot_solutions(t_max, map(h -> (h, rk4_method(df_dt_one_body, u_0, h, 1:h:t_max)), h_list_rk4), title="RK4 Method"))


text:
using SCARGC, Test, PyCall

@testset "Dataset reading" begin
    dataset = "../src/datasets/synthetic/1CDT.txt"

    data = SCARGC.extractValuesFromFile(dataset, 16000, 3)

    @testset "Array size" begin
        @test size(data, 1) == 16000
        @test size(data, 2) == 3
    end

    @testset "Array content" begin
        @test data[1346, 2] == 7.287
        @test data[1905, 2] == 8.7572
        @test data[13560, 1] == 1.9553
        @test data[14193, 1] == 9.2213
        @test data[7309, 3] == 2
        @test data[10050, 3] == 2
    end
end

@testset "Dataset fit" begin
    dataset = rand(10^3, 5)

    labeledData, labeledDataLabels, streamData, streamLabels, features = SCARGC.fitData(dataset, 5.0)

    @testset "Arrays size" begin
        @testset "Row count" begin
            @test size(labeledData)[1] == 50
            @test size(labeledDataLabels)[1] == size(labeledData)[1]
            @test size(streamData)[1] == 950
            @test size(streamLabels)[1] == size(streamData)[1]
        end

        @testset "Column count" begin
            @test size(labeledData)[2] == 4
            @test_throws BoundsError size(labeledDataLabels)[2]
            @test size(streamData)[2] == 4
            @test_throws BoundsError size(streamLabels)[2]
        end
    end

    @testset "Array content" begin
        @testset "Training data" begin
            for cases = 1:100
                row = rand(1:50)
                column = rand(1:4)

                @test labeledData[row, column] == dataset[row, column]
                @test labeledDataLabels[row] == dataset[row, 5]
            end
        end

        @testset "Stream data" begin
            for cases = 1:100
                streamRow = rand(1:950)
                row = streamRow + 50
                column = rand(1:4)

                @test streamData[streamRow, column] == dataset[row, column]
                @test streamLabels[streamRow] == dataset[row, 5]
            end
        end
    end

    @testset "Feature count" begin
        @test features == 4
    end
end

@testset "KNN Classification" begin
    dataset = [
        2.8958 -7.4953 2; 8.2872 -0.68142 4; 6.3705 -1.0572 4; -10.777 1.6717 3; -3.9285 9.8147 1;
        -8.7798 -0.88944 3; -13.005 -0.092225 3; -1.8922 -7.7006 2; -1.5711 -7.7011 2; -1.3478 -7.427 2;
        -3.0136 13.741 1; -7.443 1.0671 3; -1.222 -10.88 2; 0.97767 9.6513 1; 0.42688 -12.291 2;
        0.54482 -9.2378 2; 10.375 -1.8902 4; -9.5809 0.40066 3; 10.736 1.6078 4; 14.971 2.2448 4;
        0.071359 -6.9043 2; -11.347 -1.6508 3; 3.2586 8.8188 1; 8.5209 -1.297 4; -11.167 -1.3958 3;
        -9.0726 1.4411 3; 0.25484 10.488 1; -1.0525 -9.2088 2; -0.22116 11.641 1; -0.048457 14.003 1;
        -10.575 1.8298 3; -0.0585 9.7192 1; 10.158 -0.55542 4; 0.18513 8.3579 1; -13.068 -0.46916 3;
        9.2202 2.4973 4; 9.2265 0.12338 4; 9.9718 -0.051272 4; -8.4675 -0.85861 3; -1.3546 -12.383 2;
        1.1401 -9.9543 2; 11.135 -0.622 4; 9.5399 0.0086581 4; -8.5245 -0.98091 3; 10.476 5.2178 4;
        -2.1146 -10.042 2; -5.1069 -2.0286 3; 10.077 -1.1442 4; -1.0195 -8.5242 2; -2.015 10.026 1
    ]

    nearestNeighbors = [5, 6, 9, 5, 9, 9, 2, 6, 2, 2, 10, 7, 5, 2, 7, 4, 5, 9, 5, 5, 4, 5, 2, 5, 7, 2, 2, 2, 6, 9, 1, 2, 2, 6, 2, 8, 6, 2, 9, 5]

    labeledData, labeledDataLabels, streamData, _, _ = SCARGC.fitData(dataset, 20.0)

    @testset "Output concordance" begin
        for row = 1:40
            output, _ = SCARGC.knnClassification(labeledData, labeledDataLabels, streamData[row, :])

            @test output == dataset[nearestNeighbors[row], 3]
        end
    end

    @testset "Data equivalence" begin
        for row = 1:40
            _, data = SCARGC.knnClassification(labeledData, labeledDataLabels, streamData[row, :])

            @test data == dataset[nearestNeighbors[row], 1:2]
        end
    end

end

@testset "Finding centroids" begin
    dataset = [
        2.8958 -7.4953 2; 8.2872 -0.68142 4; 6.3705 -1.0572 4; -10.777 1.6717 3; -3.9285 9.8147 1;
        -8.7798 -0.88944 3; -13.005 -0.092225 3; -1.8922 -7.7006 2; -1.5711 -7.7011 2; -1.3478 -7.427 2;
        -3.0136 13.741 1; -7.443 1.0671 3; -1.222 -10.88 2; 0.97767 9.6513 1; 0.42688 -12.291 2;
        0.54482 -9.2378 2; 10.375 -1.8902 4; -9.5809 0.40066 3; 10.736 1.6078 4; 14.971 2.2448 4;
        0.071359 -6.9043 2; -11.347 -1.6508 3; 3.2586 8.8188 1; 8.5209 -1.297 4; -11.167 -1.3958 3;
        -9.0726 1.4411 3; 0.25484 10.488 1; -1.0525 -9.2088 2; -0.22116 11.641 1; -0.048457 14.003 1;
        -10.575 1.8298 3; -0.0585 9.7192 1; 10.158 -0.55542 4; 0.18513 8.3579 1; -13.068 -0.46916 3;
        9.2202 2.4973 4; 9.2265 0.12338 4; 9.9718 -0.051272 4; -8.4675 -0.85861 3; -1.3546 -12.383 2;
        1.1401 -9.9543 2; 11.135 -0.622 4; 9.5399 0.0086581 4; -8.5245 -0.98091 3; 10.476 5.2178 4;
        -2.1146 -10.042 2; -5.1069 -2.0286 3; 10.077 -1.1442 4; -1.0195 -8.5242 2; -2.015 10.026 1
    ]

    @testset "K == C" begin
        trainingData, trainingLabels, streamData, streamLabels, features = SCARGC.fitData(dataset, 100.0)

        centroids = SCARGC.findCentroids(features, trainingData, trainingLabels, 4)

        @testset "Array size" begin
            @test size(centroids) == (4, 3)
        end

        @testset "Array content" begin
            @test centroids[1, :] == [-1.0525, -9.2088, 2.0]
            @test centroids[2, :] == [10.0244, -0.303346, 4.0]
            @test centroids[3, :] == [-9.5809, -0.46916, 3.0]
            @test centroids[4, :] == [-0.0534785, 9.92035, 1.0]
        end
    end

    @testset "K != C" begin
        py"""
            from sklearn.cluster import KMeans
            import numpy as np
            import sys, os

            def kmeans(X, k, init, centroids=[]):
                if (init == 0):
                    return KMeans(n_clusters=k).fit(X).cluster_centers_
                else:
                    centroids = np.array(centroids)
                    sys.stdout = open(os.devnull, 'w')
                    return KMeans(n_clusters=k, init=centroids).fit(X).cluster_centers_"""

        trainingData, trainingLabels, streamData, streamLabels, features = SCARGC.fitData(dataset, 100.0)

        centroids = SCARGC.findCentroids(features, trainingData, trainingLabels, 2)

        @testset "Array size" begin
            @test size(centroids) == (2, 3)
        end
    end
end

@testset "Finding labels for current centroids" begin
    pastCentroids = [-1.45945 -7.59795 1.0 1.0 2.0; 7.32885 -0.86931 1.0 1.0 4.0; -10.777 -0.092225 1.0 1.0 3.0; -3.9285 9.8147 1.0 1.0 1.0]
    currentCentroids = [
        -3.0136 13.741 1.0 1.0; -7.443 1.0671 1.0 1.0; -1.222 -10.88 1.0 1.0; 0.97767 9.6513 1.0 1.0; 0.42688 -12.291 1.0 1.0;
        0.54482 -9.2378 1.0 1.0; 10.375 -1.8902 1.0 1.0; -9.5809 0.40066 1.0 1.0; 10.736 1.6078 1.0 1.0; 14.971 2.2448 1.0 1.0;
        0.071359 -6.9043 1.0 1.0; -11.347 -1.6508 1.0 1.0; 3.2586 8.8188 1.0 1.0; 8.5209 -1.297 1.0 1.0; -11.167 -1.3958 1.0 1.0;
        -9.0726 1.4411 1.0 1.0; 0.25484 10.488 1.0 1.0; -1.0525 -9.2088 1.0 1.0; -0.22116 11.641 1.0 1.0; -0.048457 14.003 1.0 1.0
    ]

    labels = [1.0, 3.0, 2.0, 1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 4.0, 2.0, 3.0, 1.0, 4.0, 3.0, 3.0, 1.0, 2.0, 1.0, 1.0]

    _, centroidLabels = SCARGC.findLabelForCurrentCentroids(pastCentroids, currentCentroids)

    for row = 1:20
        @test centroidLabels[row] == labels[row]
    end
end

@testset "Claculating new labeled data" begin
    poolData = [
        1.1401 -9.9543 2;
        11.135 -0.622 4;
        9.5399 0.0086581 4;
        -8.5245 -0.98091 3;
        10.476 5.2178 4;
        -2.1146 -10.042 2;
        -5.1069 -2.0286 3;
        10.077 -1.1442 4;
        -1.0195 -8.5242 2;
        -2.015 10.026 1
    ]

    currentCentroids = [
        2.8958 -7.4953 2; 8.2872 -0.68142 4; 6.3705 -1.0572 4; -10.777 1.6717 3;
        -3.9285 9.8147 1; -8.7798 -0.88944 3; -13.005 -0.092225 3; -1.8922 -7.7006 2;
        -1.5711 -7.7011 2; -1.3478 -7.427 2; -3.0136 13.741 1; -7.443 1.0671 3;
        -1.222 -10.88 2; 0.97767 9.6513 1; 0.42688 -12.291 2; 0.54482 -9.2378 2;
        10.375 -1.8902 4; -9.5809 0.40066 3; 10.736 1.6078 4; 14.971 2.2448 4
    ]

    intermedCentroids = [
        0.071359 -6.9043 2; -11.347 -1.6508 3; 3.2586 8.8188 1; 8.5209 -1.297 4;
        -11.167 -1.3958 3; -9.0726 1.4411 3; 0.25484 10.488 1; -1.0525 -9.2088 2;
        -0.22116 11.641 1; -0.048457 14.003 1; -10.575 1.8298 3; -0.0585 9.7192 1;
        10.158 -0.55542 4; 0.18513 8.3579 1; -13.068 -0.46916 3; 9.2202 2.4973 4;
        9.2265 0.12338 4; 9.9718 -0.051272 4; -8.4675 -0.85861 3; -1.3546 -12.383 2
    ]

    expectedLabels = [2, 4, 4, 3, 4, 2, 3, 4, 2, 1]

    testLabeledData = SCARGC.calculateNewLabeledData(poolData, currentCentroids, intermedCentroids)

    @testset "Array size" begin
        @test size(testLabeledData) == size(expectedLabels)
    end

    @testset "Array content" begin
        @test testLabeledData == expectedLabels
    end
end

@testset "Information update" begin
    poolData = [
        -8.7798 -0.88944 1.0 1.0 3.0; -13.005 -0.092225 1.0 1.0 3.0; -1.8922 -7.7006 1.0 1.0 2.0; -1.5711 -7.7011 1.0 1.0 2.0; -1.3478 -7.427 1.0 1.0 2.0;
        -3.0136 13.741 1.0 1.0 1.0; -7.443 1.0671 1.0 1.0 3.0; -1.222 -10.88 1.0 1.0 2.0; 0.97767 9.6513 1.0 1.0 1.0; 0.42688 -12.291 1.0 1.0 2.0;
        0.54482 -9.2378 1.0 1.0 2.0; 10.375 -1.8902 1.0 1.0 4.0; -9.5809 0.40066 1.0 1.0 3.0; 10.736 1.6078 1.0 1.0 4.0; 14.971 2.2448 1.0 1.0 4.0;
        0.071359 -6.9043 1.0 1.0 2.0; -11.347 -1.6508 1.0 1.0 3.0; 3.2586 8.8188 1.0 1.0 1.0; 8.5209 -1.297 1.0 1.0 4.0; -11.167 -1.3958 1.0 1.0 3.0
    ]

    centroids = [2.8958 -7.4953 1.0 1.0 2.0; 7.32885 -0.86931 1.0 1.0 4.0; -10.777 1.6717 1.0 1.0 3.0; -3.9285 9.8147 1.0 1.0 1.0]

    labeledData = [2.8958 -7.4953 1.0 1.0; 8.2872 -0.68142 1.0 1.0; 6.3705 -1.0572 1.0 1.0; -10.777 1.6717 1.0 1.0; -3.9285 9.8147 1.0 1.0]

    labeledDataLabels = [2.0, 4.0, 4.0, 3.0, 1.0]

    newLabeledData = [3.0, 3.0, 2.0, 2.0, 2.0, 1.0, 3.0, 2.0, 1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 4.0, 2.0, 3.0, 1.0, 4.0, 3.0]

    currentCentroids = [
        -8.7798 -0.88944 1.0 1.0 3.0; -13.005 -0.092225 1.0 1.0 3.0; -1.8922 -7.7006 1.0 1.0 2.0; -1.5711 -7.7011 1.0 1.0 2.0; -1.3478 -7.427 1.0 1.0 2.0;
        -3.0136 13.741 1.0 1.0 1.0; -7.443 1.0671 1.0 1.0 3.0; -1.222 -10.88 1.0 1.0 2.0; 0.97767 9.6513 1.0 1.0 1.0; 0.42688 -12.291 1.0 1.0 2.0;
        0.54482 -9.2378 1.0 1.0 2.0; 10.375 -1.8902 1.0 1.0 4.0; -9.5809 0.40066 1.0 1.0 3.0; 10.736 1.6078 1.0 1.0 4.0; 14.971 2.2448 1.0 1.0 4.0;
        0.071359 -6.9043 1.0 1.0 2.0; -11.347 -1.6508 1.0 1.0 3.0; 3.2586 8.8188 1.0 1.0 1.0; 8.5209 -1.297 1.0 1.0 4.0; -11.167 -1.3958 1.0 1.0 3.0
    ]

    intermed = [
        -9.7784 0.39113 1.0 1.0 3.0; -11.891 0.7897375 1.0 1.0 3.0; 0.5017999999999999 -7.59795 1.0 1.0 2.0; 0.66235 -7.5982 1.0 1.0 2.0; 0.7739999999999999 -7.46115 1.0 1.0 2.0;
        -3.47105 11.77785 1.0 1.0 1.0; -9.11 1.3694 1.0 1.0 3.0; 0.8369 -9.187650000000001 1.0 1.0 2.0; -1.475415 9.733 1.0 1.0 1.0; 1.66134 -9.89315 1.0 1.0 2.0;
        1.72031 -8.36655 1.0 1.0 2.0; 8.851925 -1.379755 1.0 1.0 4.0; -10.17895 1.0361799999999999 1.0 1.0 3.0; 9.032425 0.36924499999999993 1.0 1.0 4.0; 11.149925 0.687745 1.0 1.0 4.0;
        1.4835795 -7.1998 1.0 1.0 2.0; -11.062 0.01044999999999996 1.0 1.0 3.0; -0.3349500000000001 9.316749999999999 1.0 1.0 1.0; 7.924875 -1.083155 1.0 1.0 4.0; -10.972 0.13795000000000002 1.0 1.0 3.0
    ]

    maxPoolSize = 20

    @testset "concordantLabelCount/maxPoolSize >= 1" begin
        concordantLabelCount = 20
        labeledDataLabels = [2.0, 4.0, 4.0, 3.0, 1.0, 2.0, 4.0, 4.0, 3.0, 1.0, 2.0, 4.0, 4.0, 3.0, 1.0, 2.0, 4.0, 4.0, 3.0, 1.0, 1.0]

        centroidsGot, labeledDataGot, labeledDataLabelsGot = SCARGC.updateInformation(poolData, centroids, labeledData, labeledDataLabels, newLabeledData,
                                                                                currentCentroids, intermed, concordantLabelCount, maxPoolSize)

        @testset "Array size" begin
            @test size(centroidsGot) == size(centroids)
            @test size(labeledDataGot) == size(labeledData)
            @test size(labeledDataLabelsGot) == size(labeledDataLabels)
        end

        @testset "Array content" begin
            @test centroidsGot == centroids
            @test labeledDataGot == labeledData
            @test labeledDataLabelsGot == labeledDataLabels
        end
    end

    @testset "concordantLabelCount/maxPoolSize < 1" begin
        labeledDataLabels = [2.0, 4.0, 4.0, 3.0, 1.0]
        concordantLabelCount = 6

        expectedCentroids = [
            -8.7798 -0.88944 1.0 1.0 3.0; -13.005 -0.092225 1.0 1.0 3.0; -1.8922 -7.7006 1.0 1.0 2.0; -1.5711 -7.7011 1.0 1.0 2.0; -1.3478 -7.427 1.0 1.0 2.0;
            -3.0136 13.741 1.0 1.0 1.0; -7.443 1.0671 1.0 1.0 3.0; -1.222 -10.88 1.0 1.0 2.0; 0.97767 9.6513 1.0 1.0 1.0; 0.42688 -12.291 1.0 1.0 2.0;
            0.54482 -9.2378 1.0 1.0 2.0; 10.375 -1.8902 1.0 1.0 4.0; -9.5809 0.40066 1.0 1.0 3.0; 10.736 1.6078 1.0 1.0 4.0; 14.971 2.2448 1.0 1.0 4.0;
            0.071359 -6.9043 1.0 1.0 2.0; -11.347 -1.6508 1.0 1.0 3.0; 3.2586 8.8188 1.0 1.0 1.0; 8.5209 -1.297 1.0 1.0 4.0; -11.167 -1.3958 1.0 1.0 3.0;
            -9.7784 0.39113 1.0 1.0 3.0; -11.891 0.7897375 1.0 1.0 3.0; 0.5017999999999999 -7.59795 1.0 1.0 2.0; 0.66235 -7.5982 1.0 1.0 2.0; 0.7739999999999999 -7.46115 1.0 1.0 2.0;
            -3.47105 11.77785 1.0 1.0 1.0; -9.11 1.3694 1.0 1.0 3.0; 0.8369 -9.187650000000001 1.0 1.0 2.0; -1.475415 9.733 1.0 1.0 1.0; 1.66134 -9.89315 1.0 1.0 2.0;
            1.72031 -8.36655 1.0 1.0 2.0; 8.851925 -1.379755 1.0 1.0 4.0; -10.17895 1.0361799999999999 1.0 1.0 3.0; 9.032425 0.36924499999999993 1.0 1.0 4.0; 11.149925 0.687745 1.0 1.0 4.0;
            1.4835795 -7.1998 1.0 1.0 2.0; -11.062 0.01044999999999996 1.0 1.0 3.0; -0.3349500000000001 9.316749999999999 1.0 1.0 1.0; 7.924875 -1.083155 1.0 1.0 4.0; -10.972 0.13795000000000002 1.0 1.0 3.0
        ]

        expectedLabeledData = [
            -8.7798 -0.88944 1.0 1.0; -13.005 -0.092225 1.0 1.0; -1.8922 -7.7006 1.0 1.0; -1.5711 -7.7011 1.0 1.0; -1.3478 -7.427 1.0 1.0;
            -3.0136 13.741 1.0 1.0; -7.443 1.0671 1.0 1.0; -1.222 -10.88 1.0 1.0; 0.97767 9.6513 1.0 1.0; 0.42688 -12.291 1.0 1.0;
            0.54482 -9.2378 1.0 1.0; 10.375 -1.8902 1.0 1.0; -9.5809 0.40066 1.0 1.0; 10.736 1.6078 1.0 1.0; 14.971 2.2448 1.0 1.0;
            0.071359 -6.9043 1.0 1.0; -11.347 -1.6508 1.0 1.0; 3.2586 8.8188 1.0 1.0; 8.5209 -1.297 1.0 1.0; -11.167 -1.3958 1.0 1.0
        ]

        expectedLabeledDataLabels = [3.0, 3.0, 2.0, 2.0, 2.0, 1.0, 3.0, 2.0, 1.0, 2.0, 2.0, 4.0, 3.0, 4.0, 4.0, 2.0, 3.0, 1.0, 4.0, 3.0]

        centroidsGot, labeledDataGot, labeledDataLabelsGot = SCARGC.updateInformation(poolData, centroids, labeledData, labeledDataLabels, newLabeledData,
                                                                                currentCentroids, intermed, concordantLabelCount, maxPoolSize)

        @testset "Array size" begin
            @test size(centroidsGot, 1) == size(currentCentroids, 1) + size(intermed, 1)
            @test size(labeledDataGot, 1) == size(poolData, 1)
            @test size(labeledDataLabelsGot, 1) == size(poolData, 1)

            @test size(centroidsGot, 2) == size(currentCentroids, 2) == size(intermed, 2)
            @test size(labeledDataGot, 2) == size(poolData, 2) - 1
        end

        @testset "Array content" begin
            @test centroidsGot == expectedCentroids
            @test labeledDataGot == expectedLabeledData
            @test labeledDataLabelsGot == expectedLabeledDataLabels
        end
    end
end

@testset "SCARGC with 1NN classifier" begin
    dataset = [
        2.8958 -7.4953 2; 8.2872 -0.68142 4; 6.3705 -1.0572 4;-10.777 1.6717 3;-3.9285 9.8147 1;-8.7798 -0.88944 3;-13.005 -0.092225 3;-1.8922 -7.7006 2;-1.5711 -7.7011 2;-1.3478 -7.427 2;
        -3.0136 13.741 1;-7.443 1.0671 3;-1.222 -10.88 2;0.97767 9.6513 1;0.42688 -12.291 2;0.54482 -9.2378 2;10.375 -1.8902 4;-9.5809 0.40066 3;10.736 1.6078 4;14.971 2.2448 4;
        0.071359 -6.9043 2;-11.347 -1.6508 3;3.2586 8.8188 1;8.5209 -1.297 4;-11.167 -1.3958 3;-9.0726 1.4411 3;0.25484 10.488 1;-1.0525 -9.2088 2;-0.22116 11.641 1;-0.048457 14.003 1;
        -10.575 1.8298 3;-0.0585 9.7192 1;10.158 -0.55542 4;0.18513 8.3579 1;-13.068 -0.46916 3;9.2202 2.4973 4;9.2265 0.12338 4;9.9718 -0.051272 4;-8.4675 -0.85861 3;-1.3546 -12.383 2;
        1.1401 -9.9543 2;11.135 -0.622 4;9.5399 0.0086581 4;-8.5245 -0.98091 3;10.476 5.2178 4;-2.1146 -10.042 2;-5.1069 -2.0286 3;10.077 -1.1442 4;-1.0195 -8.5242 2;-2.015 10.026 1;
        2.2047 12.298 1;-0.074896 10.403 1;0.90863 -7.3251 2;2.5506 10.718 1;-1.7999 -11.768 2;10.041 1.9354 4;0.71283 -10.392 2;-0.24533 12.384 1;-10.916 -3.0079 3;13.157 -0.087517 4;
        12.36 -2.1673 4;12.158 1.6823 4;8.3389 1.1933 4;0.51037 -11.019 2;0.39083 6.5315 1;-1.277 -12.846 2;9.0778 0.36851 4;-12.049 2.7196 3;-12.677 -2.1093 3;-10.661 0.84432 3;
        -0.36669 10.772 1;1.8032 6.9087 1;-1.3556 -8.8028 2;-0.84456 -7.1571 2;-0.55158 8.9382 1;-1.9958 -7.0433 2;-1.7224 -15.051 2;-1.406 -10.686 2;2.5305 11.46 1;1.0389 -9.7381 2;
        -11.587 1.5958 3;0.28096 8.516 1;-0.00018349 9.131 1;-10.635 2.8625 3;-2.0658 11.189 1;0.98778 10.147 1;-10.577 -0.26176 3;-1.5552 -9.5398 2;1.995 9.5663 1;-8.5554 3.3411 3;
        -1.4494 10.32 1;2.2654 12.618 1;10.908 -0.56058 4;2.5514 10.125 1;9.4298 0.82499 4;7.7632 0.24353 4;10.986 0.72106 4;-5.187 1.0313 3;-10.733 2.1273 3;8.7337 1.0812 4
    ]

    expectedLabels = [
        1, 3, 2, 1, 2, 2, 4, 3, 4, 4, 2, 3, 1, 4, 3, 3, 1, 2, 1, 1, 3, 1, 4, 1, 3, 4, 4, 4, 3, 2,
        2, 4, 4, 3, 4, 2, 3, 4, 2, 1, 1, 1, 2, 1, 2, 4, 2, 1, 3, 4, 4, 4, 4, 2, 1, 2, 4, 3, 3, 3,
        1, 1, 2, 2, 1, 2, 2, 2, 1, 2, 3, 1, 1, 3, 1, 1, 3, 2, 1, 3, 1, 1, 4, 1, 4, 4, 4, 3, 3, 4
    ]

    labelsGot, _ = SCARGC.scargc_1NN(dataset, 10.0, 20, 4)

    @test (count(label -> label == 1, labelsGot .== expectedLabels)/90) * 100 >= 90
end


text:
# XTC molecular trajectory file format

module XTC

using ..Dorothy
using ..Dorothy.Geometry
using ..Dorothy.PBC
using ..Dorothy.XDR
using Formats
using FormatStreams
using StaticArrays

export register_xtc, readxtc!

const xtc_format_code = 1995
const xtc_magicints =
		[0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
		80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
		1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003,
		16384, 20642, 26007, 32768, 41285, 52015, 65536,82570, 104031,
		131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
		832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
		4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216]
const xtc_firstidx = 9
const xtc_lastidx = length(xtc_magicints)

function register_xtc()
	Formats.addformat("trajectory/x-xtc")
	Formats.addextension("trajectory/x-xtc", ".xtc")
	Formats.addsignature("trajectory/x-xtc", [0x00,0x00,0x07,0xcb])
	Formats.addreader("trajectory/x-xtc", DorothyIO())
	Formats.addwriter("trajectory/x-xtc", DorothyIO())
	FormatStreams.addstreamer("trajectory/x-xtc", DorothyIO())
end

FormatStreams.streamf(::DorothyIO, ::MIME"trajectory/x-xtc", io::T, args...;
		kwargs...) where {T<:IO} = XTCStream{T}(io, args...; kwargs...)

Base.read(::DorothyIO, ::MIME"trajectory/x-xtc", io::IO, args...; kwargs...) =
		readxtc!(io, MolecularModel(), args...; kwargs...)

Base.read!(::DorothyIO, ::MIME"trajectory/x-xtc", io::IO, model::MolecularModel,
		args...; kwargs...) = readxtc!(io, model, args...; kwargs...)

Base.write(::DorothyIO, ::MIME"trajectory/x-xtc", io::IO, model::MolecularModel,
		args...; kwargs...) = writextc(io, model, args...; kwargs...)

struct XTCMetadata
	nparticles::Int
	precision::Float64

	function XTCMetadata(nparticles::Integer, precision::Real = 1000.0)
		@boundscheck begin
			nparticles >= 0 || error("expected positive number of particles")
			precision >= 0.0 || error("expected positive precision")
		end
		if nparticles <= 9
			precision = 0.0
		end
		new(nparticles, precision)
	end
end

function read_xtc_metadata(io::IO)
	skip(io, 4)
	nparticles = read_xdr_int(io)
	if nparticles <= 9
		precision = 0.0
	else
		skip(io, 48)
		precision = read_xdr_float(io)
	end
	XTCMetadata(nparticles, precision)
end

function peek_xtc_metadata(io::IO)
	mark(io)
	meta = read_xtc_metadata(io)
	reset(io)
	meta
end

struct XTCBuffer
	cic::Vector{UInt8}
	uic::Vector{Int}
	fic::Vector{Int}
	bytes::MVector{32,Int}
	state::MVector{3,Int}
	cell::MMatrix{3,3,Float64}

	function XTCBuffer(nparticles::Integer)
		cic = Vector{UInt8}(undef, floor(Int, 3*nparticles*1.5))
		uic = Vector{UInt8}(undef, 3*nparticles)
		fic = Vector{UInt8}(undef, 3*nparticles)
		bytes = MVector{32,Int}(undef)
		state = MVector{3,Int}(undef)
		cell = MMatrix{3,3,Float64}(undef)
		new(cic, uic, fic, bytes, state, cell)
	end
end

mutable struct XTCStream{T<:IO} <: MolecularTrajectory
	io::T
	frameindices::Vector{Int}
	pos::Int
	meta::XTCMetadata
	buffer::XTCBuffer

	function XTCStream{T}(io::T; nparticles::Integer = -1,
			precision::Real = 1000.0) where {T<:IO}
		origin = position(io)
		if eof(io)
			meta = XTCMetadata(nparticles, precision)
		else
			meta = peek_xtc_metadata(io)
		end
		new(io, [origin], 0, meta, XTCBuffer(meta.nparticles))
	end
end

Base.close(s::XTCStream) = close(s.io)

Base.eof(s::XTCStream) = eof(s.io)

Base.position(s::XTCStream) = s.pos

function Base.seek(s::XTCStream, pos::Integer)
	i = pos + 1
	@boundscheck begin
		i >= 0 || error("expected positive position")
		if i > length(s.frameindices)
			error("cannot seek to unknown position")
		end
	end
	seek(s.io, s.frameindices[i])
	s.pos = pos
	s
end

function Base.seekstart(s::XTCStream)
	seek(s.io, s.frameindices[1])
	s.pos = 0
	s
end

Base.read(s::XTCStream) = read!(s, MolecularModel())

function Base.read!(s::XTCStream, model::MolecularModel)
	readxtc!(s.io, model, s.meta, s.buffer)
	s.pos += 1
	if s.pos + 1 > length(s.frameindices)
		push!(s.frameindices, position(s.io))
	end
	model
end

#=
function Base.write(s::XTCStream, model::MolecularModel)
	@boundscheck eof(s.io) || error("writing must happen at end of trajectory")
	writextc(s.io, model, s.meta, s.buffer)
	s.pos += 1
	push!(s.frameindices, position(s.io))
	s.nframes += 1
end
=#

function Base.truncate(s::XTCStream, n::Integer)
	@boundscheck begin
		n > 0 || error("expected positive value")
		n + 1 > length(s.frameindices) &&
				error("cannot truncate to unknown position")
	end
	truncate(s.io, s.frameindices[n+1])
	deleteat!(s.frameindices, n+2:lastindex(s.frameindices))
	if s.pos > n
		s.pos = n
	end
	s
end

function readxtc!(io::IO, model::MolecularModel,
		meta::XTCMetadata = peek_xtc_metadata(io),
		buffer::XTCBuffer = XTCBuffer(meta.nparticles))
	skip(io, 8)
	resize!(model, meta.nparticles)
	model.header.step = read_xdr_int(io)
	model.header.time = read_xdr_float(io)
	for i = 1:9
		buffer.cell[i] = read_xdr_float(io) * 10.0
	end
	model.header.cell = TriclinicCell(buffer.cell)
	skip(io, 4)
	R = get!(model, :R, undef)
	if meta.nparticles > 9
		read3dfcoord!(io, R, meta, buffer)
	else
		readxdrcoord!(io, R, meta.nparticles)
	end
	model
end

function read3dfcoord!(io::IO, R::AbstractVector{Vector3D}, meta::XTCMetadata,
		buffer::XTCBuffer)

	magicints = xtc_magicints
	firstidx = xtc_firstidx
	lastidx = xtc_lastidx

	nparticles = meta.nparticles
	precision = meta.precision
	cic = buffer.cic
	uic = buffer.uic
	fic = buffer.fic
	bytes = buffer.bytes
	state = buffer.state

	skip(io, 4)
	minint = SVector{3,Int}(read_xdr_int(io), read_xdr_int(io),
			read_xdr_int(io))
	maxint = SVector{3,Int}(read_xdr_int(io), read_xdr_int(io),
			read_xdr_int(io))
	smallidx = Int(read_xdr_int(io))
	nbytes = Int(read_xdr_int(io))
	if smallidx == 0
		error("not sure what this error is...")
	end
	if nbytes == 0
		error("zero-length compressed coordinates")
	end

	sizeint = maxint - minint .+ 1
	if (sizeint[1] | sizeint[2] | sizeint[3]) > 0xffffff
		bitsizeint = size_of_xdr_int.(sizeint)
		bitsize = 0
	else
		bitsize = size_of_xdr_ints!(bytes, sizeint)
	end

	tmp = smallidx + 8
	maxidx = lastidx < tmp ? lastidx : tmp
	minidx = maxidx - 8
	tmp = smallidx - 1
	tmp = firstidx > tmp ? firstidx : tmp
	smaller = div(magicints[tmp+1], 2)
	smallnum = div(magicints[smallidx+1], 2)
	sizesmall = SVector(magicints[smallidx+1], magicints[smallidx+1],
			magicints[smallidx+1])
	larger = magicints[maxidx+1]

	resize!(cic, nbytes)
	read!(io, cic)
	remainder = nbytes % 4
	if remainder != 0
		skip(io, 4 - remainder)
	end
	for i in 1:3
		state[i] = 0
	end

	i = 0
	j = 1
	k = 0
	run = 0
	while i < nparticles
		d = 3 * i
		x = d + 1
		y = d + 2
		z = d + 3
		if bitsize == 0
			uic[x] = decode_xdr_bits!(state, cic, bitsizeint[1])
			uic[y] = decode_xdr_bits!(state, cic, bitsizeint[2])
			uic[z] = decode_xdr_bits!(state, cic, bitsizeint[3])
		else
			decode_xdr_ints!(uic, bytes, state, cic, 3, bitsize, sizeint, x)
		end
		i += 1

		uic[x] += minint[1]
		uic[y] += minint[2]
		uic[z] += minint[3]

		prevcoordx = uic[x]
		prevcoordy = uic[y]
		prevcoordz = uic[z]

		flag = decode_xdr_bits!(state, cic, 1)
		is_smaller = 0
		if flag == 1
			run = decode_xdr_bits!(state, cic, 5)
			is_smaller = run % 3
			run -= is_smaller
			is_smaller -= 1
		end
		if run > 0
			d += 3
			x = d + 1
			y = d + 2
			z = d + 3
			k = 0
			while k < run
				decode_xdr_ints!(uic, bytes, state, cic, 3, smallidx,
						sizesmall, x)
				i += 1
				uic[x] += prevcoordx - smallnum
				uic[y] += prevcoordy - smallnum
				uic[z] += prevcoordz - smallnum
				if k == 0
					tmp = uic[x]
					uic[x] = prevcoordx
					prevcoordx = tmp
					tmp = uic[y]
					uic[y] = prevcoordy
					prevcoordy = tmp
					tmp = uic[z]
					uic[z] = prevcoordz
					prevcoordz = tmp
					fic[j] = prevcoordx
					j += 1
					fic[j] = prevcoordy
					j += 1
					fic[j] = prevcoordz
					j += 1
				else
					prevcoordx = uic[x]
					prevcoordy = uic[y]
					prevcoordz = uic[z]
				end
				fic[j] = uic[x]
				j += 1
				fic[j] = uic[y]
				j += 1
				fic[j] = uic[z]
				j += 1
				k += 3
			end
		else
			fic[j] = uic[x]
			j += 1
			fic[j] = uic[y]
			j += 1
			fic[j] = uic[z]
			j += 1
		end
		smallidx += is_smaller
		if is_smaller < 0
			smallnum = smaller
			if smallidx > firstidx
				smaller = div(magicints[smallidx], 2)
			else
				smaller = 0
			end
		elseif is_smaller > 0
			smaller = smallnum
			smallnum = div(magicints[smallidx+1], 2)
		end
		sizesmall = SVector(magicints[smallidx+1], magicints[smallidx+1],
				magicints[smallidx+1])
	end

	for i in eachindex(R)
		d = (i-1) * 3
		x = (fic[d+1] * 10) / precision
		y = (fic[d+2] * 10) / precision
		z = (fic[d+3] * 10) / precision
		R[i] = Vector3D(x,y,z)
	end
	R
end

function readxdrcoord!(io::IO, R::AbstractArray{<:Real}, nparticles::Integer)
	for i = 1:nparticles
		x = read_xdr_float(io)
		y = read_xdr_float(io)
		z = read_xdr_float(io)
		R[i] = Vector3D(x,y,z)
	end
	R
end

end # module


text:
# # Two-Stage Stochastic Program

# ## Introduction
# First let's consider a standard two-stage stochastic program. Such problems
# consider 1st stage variables ``x \in X \subseteq \mathbb{R}^{n_x}`` which denote
# upfront (here-and-now) decisions made before any realization of the random
# parameters ``\xi \in \mathbb{R}^{n_\xi}`` is observed, and 2nd stage variables
# ``y(\xi) \in \mathbb{R}^{n_y}`` which denote recourse (wait-and-see) decisions
# that are made in response to realizations of ``\xi``. Moreover, the objective
# seeks to optimize 1st stage costs ``f_1(x)`` and second stage costs
# ``f_2(x, y(\xi))`` which are evaluated over the uncertain domain via a risk
# measure ``R_\xi[\cdot]`` (e.g., the expectation ``\mathbb{E}_\xi[\cdot]``).
# Putting this together, we obtain the two-stage stochastic program:
# ```math
# \begin{aligned}
#     &&\min_{x, y(\xi)} &&& f_1(x) + R_\xi[f_2(x, y(\xi))] \\
#     &&\text{s.t.} &&&  g_i(x, y(\xi)) = 0, && i \in I\\
#     &&&&& h_j(x, y(\xi)) \leq 0, && j \in J\\
#     &&&&&  x \in X\\
# \end{aligned}
# ```
# where ``g_i(x, y(\xi)), \ i \in I,`` denote 2nd stage equality constraints,
# ``h_j(x, y(\xi)), \ j \in J,`` are 2nd stage inequality constraints, and ``X``
# denotes the set of feasible 1st stage decisions.

# ## Formulation

# For an example, we consider the classic farmer problem. Here the farmer must
# allocate farmland ``x_c`` for each crop ``c \in C`` with random yields per acre
# ``\xi_c`` such that he minimizes expenses (i.e., maximizes profit) while fulfilling
# contractual demand ``d_c``. If needed he can purchase crops from other farmers
# to satisfy his contracts. He can also sell extra crop yield that exceeds his
# contractual obligations. Thus, here we have 1st stage variables ``x_c`` and
# 2nd stage variables of crops sold ``w_c(\xi)`` and crops purchased ``y_c(\xi)``.
# Putting this together using the expectation ``\mathbb{E}_\xi[\cdot]`` as our
# risk measure we obtain:
# ```math
# \begin{aligned}
#     &&\underset{x, y(\xi), w(\xi)}{\text{min}} &&& \sum_{c \in C} \alpha_c x_c + \mathbb{E}_{\xi}\left[\sum_{c \in C}\beta_c y_c(\xi) - \lambda_c w_c(\xi)\right] \\
#     &&\text{s.t.} &&&  \sum_{c \in C} x_c \leq \bar{x}\\
#     &&&&& \xi_c x_c + y_c(\xi) - w_c(\xi) \geq d_c, && c \in C \\
#     &&&&& 0 \leq x_c \leq \bar{x}, && c \in C \\
#     &&&&& 0 \leq y_c(\xi) \leq \bar{y}_c, && c \in C \\
#     &&&&& 0 \leq w_c(\xi) \leq \bar{w}_c, && c \in C \\
#     &&&&& \xi_c \in \Xi_c, && c \in C
# \end{aligned}
# ```
# where ``\alpha_c`` are production costs, ``\beta_c`` are the purchase prices,
# ``\lambda_c`` are the selling prices, ``\bar{x}`` is the total acreage,
# ``\bar{y}_c`` are purchases limits, ``\bar{w}_c`` are selling limits, and
# ``\Xi_c`` are the underlying distributions.

# ## Problem Setup

# First let's import the necessary packages:
using InfiniteOpt, Distributions, Ipopt

# Next let's specify the problem data:
num_scenarios = 10 # small amount for example
C = 1:3
α = [150, 230, 260] # land cost
β = [238, 210, 0]   # purchasing cost
λ = [170, 150, 36]  # selling price
d = [200, 240, 0]   # contract demand
xbar = 500          # total land
wbar3 = 6000        # no upper bound on the other crops
ybar3 = 0           # no upper bound on the other crops
Ξ = [Uniform(0, 5), Uniform(0, 5), Uniform(10, 30)]; # the distributions

# ## Problem Definition

# Let's start by setting up the infinite model that uses Ipopt as the optimizer
# that will ultimately be used to solve the transcribed variant:
model = InfiniteModel(Ipopt.Optimizer)
set_optimizer_attribute(model, "print_level", 0);

# Now let's define the infinite parameters using [`@infinite_parameter`](@ref):
@infinite_parameter(model, ξ[c in C] ~ Ξ[c], num_supports = num_scenarios)

# Now let's define all of the decision variables using `@variables`:
@variables(model, begin
    ## 1st stage variables
    0 <= x[C] <= xbar
    ## 2nd stage variables
    0 <= y[C], Infinite(ξ)
    0 <= w[C], Infinite(ξ)
end)

# Next, the objective is defined using `@objective` and [`𝔼`](@ref):
@objective(model, Min, sum(α[c] * x[c] for c in C) +
                       𝔼(sum(β[c] * y[c] - λ[c] * w[c] for c in C), ξ))

# Finally, all we need to do is define the constraints using `@constraints`:
@constraints(model, begin
    ## capacity constraint
    sum(x[c] for c in C) <= xbar
    ## balances
    [c in C], ξ[c] * x[c] + y[c] - w[c] >= d[c]
    ## crop limits
    w[3] <= wbar3
    y[3] <= ybar3
end)

# ## Problem Solution

# With the model defined, let's optimize and get the results
optimize!(model)
x_opt = value.(x)
profit = -objective_value(model)

println("Land Allocations: ", [round(x_opt[k], digits = 2) for k in keys(x_opt)])
println("Expected Profit: \$", round(profit, digits = 2))

# We did it!

# ## CVaR Objective

# An interesting modification to the above problem would be to use a ``CVaR``
# risk measure instead of an expectation. This also can be readily achieved via
# `InfiniteOpt`. The ``CVaR`` measure is defined:
# ```math
# CVaR_\epsilon(X) = \underset{t \in \mathbb{R}}{\text{inf}}\left\{t + \frac{1}{1-\epsilon} \mathbb{E}[\text{max}(0, X - t)] \right\}
# ```
# where ``\epsilon`` is the confidence level. Inserting this into the formulation,
# we now obtain:
# ```math
# \begin{aligned}
#     &&\underset{x, y(\xi), w(\xi), t, q(\xi)}{\text{min}} &&& \sum_{c \in C} \alpha_c x_c + t + \frac{1}{1-\epsilon} \mathbb{E}_{\xi}[q(\xi)] \\
#     &&\text{s.t.} &&& \sum_{c \in C} x_c \leq \bar{x}\\
#     &&&&& \xi_c x_c + y_c(\xi) - w_c(\xi) \geq d_c, && c \in C \\
#     &&&&& 0 \leq x_c \leq \bar{x}, && c \in C \\
#     &&&&& 0 \leq y_c(\xi) \leq \bar{y}_c, && c \in C \\
#     &&&&& 0 \leq w_c(\xi) \leq \bar{w}_c, && c \in C \\
#     &&&&& \xi_c \in \Xi_c, && c \in C \\
#     &&&&& q(\xi) \geq \sum_{c \in C}\beta_c y_c(\xi) - \lambda_c w_c(\xi) - t \\
#     &&&&& q(\xi) \geq 0
# \end{aligned}
# ```
# where ``q(\xi)`` is introduced to handle the max operator. Let's update and
# resolve our `InfiniteOpt` model using ``\epsilon = 0.95``:

# Define the additional variables:
@variables(model, begin
    t
    q >= 0, Infinite(ξ)
end)

# Redefine the objective:
@objective(model, Min, sum(α[c] * x[c] for c in C) + t + 1 / (1 - 0.95) * 𝔼(q, ξ))

# Add the max constraint:
@constraint(model, q >= sum(β[c] * y[c] - λ[c] * w[c] for c in C) - t)

# Optimize and get the results
optimize!(model)
x_opt = value.(x)
y_opt = value.(y)
w_opt = value.(w)
profit = -sum(α[c] * x_opt[c] for c in C) - 1 / num_scenarios *
            sum(β[c] * y_opt[c][k] - λ[c] * w_opt[c][k] for c in C, k in 1:num_scenarios)

println("Land Allocations: ", [round(x_opt[k], digits = 2) for k in keys(x_opt)])
println("Expected Profit: \$", round(profit, digits = 2))

# That's it!

# ### Maintenance Tests
# These are here to ensure this example stays up to date.
using Test
@test termination_status(model) == MOI.LOCALLY_SOLVED
@test x_opt isa JuMPC.DenseAxisArray{<:Real}
@test profit isa Real


text:
function convert_arg(arr::AbstractArray) pointer(arr) end
function prepare_arg(arr::AbstractArray) HSA.memory_register(arr, sizeof(arr)) end
function cleanup_arg(arr::AbstractArray) HSA.memory_deregister(arr, sizeof(arr)) end

function convert_arg(x) x end
function prepare_arg(x) end
function cleanup_arg(x) end

export prepare_args
function prepare_args(args)
    kernel_args = [ convert_arg(a) for a in args ]

    for a in args
        prepare_arg(a)
    end

    return kernel_args
end

export cleanup_args
function cleanup_args(args)
    for a in args
        cleanup_arg(a)
    end
end

export allocate_args
function allocate_args(agent, kernel_info, kernel_args)
    karg_size = kernel_info.kernarg_size

    debug_print("allocate_args: Allocate $karg_size bytes of kernarg memory for $(string(kernel_info.func))")

    karg_region = find_kernarg_region(agent)

    karg_mem = HSA.memory_allocate(karg_region, karg_size)

    arg_offset = 0
    for i = 1:size(kernel_args, 1)
        a = kernel_args[i]
        arg_size = sizeof(a)

        if (arg_offset + arg_size) > karg_size
            error("Kernel Argument overflow at argument index $i, offset $arg_offset with kernarg_size $karg_size")
        end

        typed_ptr = convert(Ptr{typeof(a)}, karg_mem.ptr + arg_offset);

        unsafe_store!(typed_ptr, a)

        arg_offset += arg_size
    end

    return karg_mem
end

function find_kernarg_region(agent)
    regions = HSA.regions(agent)

    karg_region_idx = findfirst(r -> begin
        segment = HSA.region_info_segment(r)
        if (HSA.HSA_REGION_SEGMENT_GLOBAL != segment)
            return false
        end

        flags = HSA.region_info_global_flags(r)
        if (flags & HSA.HSA_REGION_GLOBAL_FLAG_KERNARG) != 0
            return true
        end

        return false
    end, regions)
    if karg_region_idx == 0
        error("No Kernarg Region found")
    end

    return regions[karg_region_idx]
end


text:
# %%
push!(LOAD_PATH, ".");
# %%
@time using CellularBase
#@time using Plots
#@time using LsqFit
@time using Random
# %%
mutable struct SandPile <: AbstractGrid{Int8, 2}
    state::Array{Int8,2}
    h_c::Int8
    neighborhood::AbstractNeighborhood
    bc::BoundaryCondition
    function SandPile(state::Array{Int8,2}, h_c::Int8)
        new(state, h_c, VonNeumannNeighborhood(1, 2), FixedMin)
    end
end

CellularBase.possible_states(sp::SandPile) = 0:sp.h_c

@inline iscritical(val::Int8, g::SandPile) = val >= g.h_c
@inline iscritical(g::SandPile) = p -> iscritical(p, g)

function (_g::SandPile)(current::Int8, neighbors::Array{Int8}; kwargs...)
    h = iscritical(_g)
    if h(current) kwargs[:topple][end] += 1 end
    return current + (h(current) ? -4 : 0) + count(h, neighbors)
end

function ready_next(grid::SandPile, step::Int; topple::Array{Int64}, max_steps)
    if count(iscritical(grid), state(grid)) == 0
        # No critical left
        if length(topple) == max_steps
            # Enough, my god
            return :interrupt
        else
            # Moar, moar, MOAR!
            xrand, yrand = rand.(Base.OneTo.(size(grid)))
            # Potentially destabilize - HAHAHAA
            grid.state[xrand, yrand] += 1
            # Make new counter
            push!(topple, 0)
        end
    end
    return :continue
end
# %%
Random.seed!(0)
h_c = Int8(4)
topple_count = Int[0]
w = 20
rand_arr = rand(Int8.(1:h_c-1), w, w)
zero_arr = zeros(Int8, w, w)

grid = SandPile(rand_arr, h_c)
@time simulate!(
    grid, typemax(Int),
    postrunhook=ready_next, topple=topple_count,
    max_steps=100_000, store_results=false);

# %%

function freqtable(k::Array{T}) where T
    d = Dict{T,Int}()
    for i ∈ k
        d[i] = get(d, i, 0) + 1
    end
    return d
end

topple_freq = freqtable(topple_count)

topple_freq = filter!((p) -> p[1] * p[2] != 0, topple_freq)

x = log.(keys(topple_freq))
y = log.(values(topple_freq))

plot(x,y, st=:scatter)
# %%
@. model(x, p) = p[1] * x + p[2]
fit = curve_fit(model, x, y, [0.,0.])
println(fit.param)
xp = 0:0.01:maximum(x)
yp = model(xp, fit.param)
plot!(xp, yp)

# %%
#=
anim = @animate for frame ∈ results
    wireframe(frame, zlims=(0,4))
end
gif(anim, fps=30)
=#
#%%


text:
using Hypatia
using Random
using Empirikos
using LinearAlgebra
using Test
using JuMP
using StatsBase

G = Normal(1,0.3)
F_G = marginalize(StandardNormalSample(), G)
Random.seed!(1)

Zs = StandardNormalSample.(rand(F_G, 1000))

gcal = MixturePriorClass(Normal.(-3:0.5:3, 0.2))

ghat = StatsBase.fit( NPMLE(;convexclass=gcal, solver=Hypatia.Optimizer),  Zs)

floc = StatsBase.fit(DvoretzkyKieferWolfowitz(;α=0.01), Zs)

discr = interval_discretizer(-3:0.2:3)

amari_withF  = AMARI(;convexclass=gcal, flocalization=floc,solver=Hypatia.Optimizer, discretizer=discr,
              plugin_G = ghat.prior, modulus_model=Empirikos.ModulusModelWithF)

amari_withoutF  = AMARI(;convexclass=gcal, flocalization=floc,solver=Hypatia.Optimizer, discretizer=discr,
              plugin_G = ghat.prior, modulus_model=Empirikos.ModulusModelWithoutF)

target = MarginalDensity(NormalSample(1.0, 0.5))
target(G)

fit_withoutF = StatsBase.fit(amari_withoutF, target, Zs)
fit_withF = StatsBase.fit(amari_withF, target, Zs)

@test fit_withoutF.Q(StandardNormalSample(2.0)) ≈ fit_withF.Q(StandardNormalSample(2.0)) atol = 0.001

ci_withF = confint(amari_withF, target, Zs)
ci_withoutF = confint(amari_withoutF, target, Zs)

@test ci_withF.lower ≈ ci_withoutF.lower rtol =0.001
@test ci_withF.upper ≈ ci_withoutF.upper rtol =0.001
@test ci_withF.halflength ≈ ci_withoutF.halflength rtol =0.001
@test ci_withF.estimate ≈ ci_withoutF.estimate rtol =0.001
@test ci_withF.maxbias ≈ ci_withoutF.maxbias rtol =0.001

_init_with_F = Empirikos.initialize_method(amari_withF, target, Zs)
Empirikos.set_δ!(_init_with_F.modulus_model,  0.01);

_init_without_F = Empirikos.initialize_method(amari_withoutF, target, Zs)
Empirikos.set_δ!(_init_without_F.modulus_model,  0.01);

@test JuMP.dual(_init_with_F.modulus_model.bound_delta) ≈ JuMP.dual(_init_without_F.modulus_model.bound_delta) rtol=0.0001

Random.seed!(1)
G = Normal(0, 1.0)
σs = rand(1:10,1000)./5
zs = rand(G, 1000) .+ randn(1000) .* σs
Zs = NormalSample.(zs, σs)
@test length(Empirikos.heteroskedastic(Zs).vec) == 10

gcal_mix = MixturePriorClass(Normal.(-4:0.2:4, 0.5))

discr = interval_discretizer(-4:0.1:4)

amari  = AMARI(;convexclass=gcal, flocalization=DvoretzkyKieferWolfowitz(;α=0.01),
              solver=Hypatia.Optimizer, discretizer=discr,
              plugin_G = NPMLE(;convexclass=gcal, solver=Hypatia.Optimizer))

postmean_target = PosteriorMean( NormalSample(1.0, 0.5))
targets = [MarginalDensity( NormalSample(1.0, 0.3)), Empirikos.PriorDensity(  1.5)]
for target in targets
    amari_fit = StatsBase.fit(amari, target, Zs)
    amari_ci = confint(amari_fit, target, Zs)
    lp_biases = Empirikos.worst_case_bias_lp(amari_fit.method, amari_fit.Q, target)
    @show target
    @show lp_biases
    @test lp_biases.maxbias ≈ amari_ci.maxbias atol = 0.00001
    @test lp_biases.maxbias ≈ -lp_biases.minbias atol = 0.00001
    # high tolerance for this test
    @test sqrt(lp_biases.expected_var / 1000) ≈ amari_ci.se rtol=0.1
    @test_throws ErrorException confint(amari_fit, postmean_target, Zs)
end

postmean_ci = confint(amari, postmean_target, Zs)

postmean_ci_01 = confint(amari, postmean_target, Zs; level=0.9)
@test postmean_ci_01.upper - postmean_ci_01.lower < postmean_ci.upper - postmean_ci.lower

#Qs = amari_fit.Q.(Zs)

#X = [fill(1, length(Qs)) Qs]
#X_hats = diag(X*inv(X'X)*X')
#scatter( std.(Zs), X_hats)


text:
# Different cooldown schedules. Internally, a cooldown function
# accepts the current temperature as its first argument and the
# current energy as its second argument and return the new temperature
# of the system.

@doc raw"""
    exponential_cooldown(λ :: Float64)

Create an exponential cooldown schedule. Temperature on step $n + 1$
is calculated as $T_{n + 1} = \lambda T_{n}$.
"""
exponential_cooldown(λ :: Float64 = 0.999999) = (T :: Float64, :: Float64) -> λ * T

"""
    aarts_korst_cooldown(;n = 15, λ = 0.01)

Make the Aarts-Korst cooldown schedule. The temperature decreases
each `n` steps of annealing algorithm. The higer `λ` is the faster
decreases the temperature.

For more information, see Aarts, E.H.L. and Korst, J.H.M. (1989)
Simulated Annealing and Boltzmann Machines: A Stochastic Approach to
Combinatorial Optimization and Neural Computing. John Wiley & Sons,
Chichester.
"""
function aarts_korst_cooldown(;n :: Integer = 15, λ :: Float64 = 0.01)
    costs = Vector{Float64}(undef, n)
    counter = 1

    function f(T :: Float64, cost :: Float64)
        costs[counter] = cost
        counter = mod1(counter + 1, n)

        newT = T
        if counter == 1
            σ = std(costs)
            newT = T * σ / (σ + λ*T)
        end

        return newT
    end

    return f
end

"""
    frost_heineman_cooldown(;n = 250, λ = 0.01)

Make the Frost-Heineman cooldown schedule. The temperature decreases
each `n` steps of annealing algorithm or until some "target" energy is
reached. The target energy is based on standard deviation of the cost
function and a parameter `λ`.

For more information, see R. Frost, P. Heineman "Simulated Annealing:
A Heuristic for Parallel Stochastic Optimization" (1997)
"""
function frost_heineman_cooldown(;n :: Integer = 250, λ :: Float64 = 0.01)
    costs = Vector{Float64}(undef, 0)
    sizehint!(costs, n)
    target = Inf

    function f(T :: Float64, cost :: Float64)
        push!(costs, cost)
        newT = T

        μ = mean(costs)
        if length(costs) == n && μ < target
            σ = std(costs)
            empty!(costs)

            target_prev = isinf(target) ? μ : target
            target = μ - λ*σ
            newT = T + (target - target_prev) * (T/σ)^2
        end

        return newT
    end

    return f
end


text:
"""
    cluster_loads!(grid,nload)

Execute aglomerative clustering of all load buses into 'nload' substations
using geographic distances as metric.
"""
function cluster_loads!(grid::Grid, nload)
    loads = [bus for bus in buses(grid) if isa(bus, LoadBus)]
    for b in loads # creates 1 substation for each load bus
        push!(
            substations(grid),
            Substation(
                length(substations(grid)) + 1,
                b.coords,
                [b.voltage],
                b.load,
                0,
                b.population
            )
        )
        substations(grid)[end].grouping = [b]
    end
    # Now let's avoid computing some things several times
    load_dist = distance(loads)
    all_ids = [l.id for l in loads]
    id2ind = Dict{Int, Int}()
    for id in all_ids
        id2ind[id] = findfirst(all_ids, id)
    end
    m = length(loads)
    sum_pops = fill(0, (m, m))
    for s1 in 1:(m - 1), s2 in (s1 + 1):m
        sum_pops[s2, s1] = loads[s1].population + loads[s2].population
    end
    sum_pops = Symmetric(sum_pops, :L)
    sub_dist = Matrix(load_sub_dist(substations(grid), load_dist, id2ind, sum_pops))
    const POPLIMIT = 200000 # population limit for a substation (=> 400MW load)

    while length(substations(grid)) > nload
        n = length(substations(grid))
        ok = false
        i = j = -1
        while true
            m = indmin(sub_dist)
            if sub_dist[m] == Inf # Can't cluster anymore without going over POPLIMIT
                warn("Can't cluster anymore without going over population limit!")
                warn("Stopping clustering process before target number is reached")
                nload = Inf
                break
            end
            i, j = ind2sub(size(sub_dist), m)
            if (substations(grid)[i].population + substations(grid)[j].population) < POPLIMIT
                break
            else # this keeps loads from clustering beyond the population limit
                sub_dist[i, j] = Inf
                sub_dist[j, i] = Inf
            end
        end
        if nload == Inf
            break
        end
        merge!(grid, i, j)
        # Update substation distances
        mask = fill(true, (n, n))
        mask[j, :] = false # Eliminate rows and columns corresponding to the
        mask[:, j] = false # deleted substation
        sub_dist = reshape(sub_dist[mask], (n - 1, n - 1)) # Remove old dists
        i = i > j ? i - 1 : i # Update i to the new length
        update_sub_dist!(sub_dist, i, substations(grid), load_dist, id2ind, sum_pops)
    end
end

"""
    fill_gen!(sub, gens, MW)

Group generators ('gens') into substations ('subs') by geographic proximity
until generation threshold 'MW' is reached.
"""
function fill_gen!(sub, gens, MW)
    dists = [haversine(sub, g) for g in gens]
    while sub.generation < MW
        m = indmin(dists)
        sub.generation += gens[m].generation
        sub.voltages = [v for v in Set(vcat(sub.voltages, gens[m].voltage))]
        push!(sub.grouping, gens[m])
        deleteat!(dists, m)
        deleteat!(gens, m)
    end
end

"""
    cluster_load_gen!(grid, nboth)

Create 'nboth' substations with both load and generation, chosen randomly.
"""
function cluster_load_gen!(grid::Grid, nboth)
    gens = [bus for bus in buses(grid) if isa(bus, GenBus)]
    nb = 0
    srand(grid.seed + 66) # Just so we don't keep returning to the same point.
    s = rand(1:length(substations(grid)))
    rrange = e^(-1):PREC:e
    while nb < nboth
        while substations(grid)[s].generation > 0 # ensures we are drawing a
            s = rand(1:length(substations(grid))) # load-only substation
        end
        MW = rand(rrange) * substations(grid)[s].load
        fill_gen!(substations(grid)[s], gens, MW)
        nb += 1
    end
end

"""
    cluster_gens!(grid, ngen)

Execute aglomerative clustering of generators into 'ngen' substations.
"""
function cluster_gens!(grid::Grid, ngen)
    gen_both = GenBus[]
    for sub in substations(grid), bus in sub.grouping
        if isa(bus, GenBus)
            push!(gen_both, bus)
        end
    end
    gens = filter(buses(grid)) do bus
        isa(bus, GenBus) && !(bus in gen_both)
    end
    g_subs = Array{Substation}(length(gens))
    st = length(substations(grid))
    for (i, g) in enumerate(gens) # creates substations for the remaining generators
        substation = Substation(
            length(substations(grid)) + 1,
            g.coords, g.voltage,
            0,
            g.generation,
            0
        )
        substation.grouping = [g]
        push!(substations(grid), substation)
        g_subs[i] = substation
    end
    dists = distance(g_subs)
    while length(g_subs) > ngen
        m = indmin(dists)
        n = length(g_subs)
        ii, jj = ind2sub(size(dists), m)
        merge!(grid, ii + st, jj + st)
        deleteat!(g_subs, jj)
        # Update substation distances
        mask = fill(true, (n, n))
        mask[jj, :] = false # Eliminate rows and columns corresponding to the
        mask[:, jj] = false # deleted substation
        dists = reshape(dists[mask], (n - 1, n - 1)) # Remove old dists
        ii = ii > jj ? ii - 1 : ii # Update i to the new length
        update_distance!(dists, g_subs, ii)
    end
end

"""
    cluster!(grid::Grid, nloads, nboth, ngens)

Cluster all nodes of a grid into 'nloads' load substations, 'ngens' generation
substations and 'nboth' substations with both load and generation.

REFERENCE: Birchfield, Adam B., et al. "Grid Structural Characteristics as
Validation Criteria for Synthetic Networks."
IEEE Transactions on Power Systems (2016).
"""
function cluster!(grid::Grid, nloads, nboth, ngens)
    nloads > 0 ? cluster_loads!(grid, nloads) : nothing
    nboth > 0 ? cluster_load_gen!(grid, nboth) : nothing
    ngens > 0 ? cluster_gens!(grid, ngens) : nothing
    (nloads + nboth + ngens) > 0 ? connect_subs!(grid) : nothing
    # The processes above should mess substation ids, so we need to correct them
    for i in length(substations(grid))
        substations(grid)[i].id = i
    end
end


text:
abstract type HamiltonianModel end

const σx = Float64[0 1; 1 0]
const σy = ComplexF64[0 -1im; 1im 0]
const σz = Float64[1 0; 0 -1]
const id2 = Float64[1 0; 0 1]

@doc raw"
    hamiltonian(model<:HamiltonianModel)

return the hamiltonian of the `model` as a two-site tensor operator.
"
function hamiltonian end

struct Ising <: HamiltonianModel end

struct diaglocal{T<:Vector} <: HamiltonianModel
    diag::T
end

"""
    diaglocal(diag::Vector)

return the 2-site Hamiltonian with single-body terms given
by the diagonal `diag`.
"""
function hamiltonian(model::diaglocal)
    diag = model.diag
    n = length(diag)
    h = ein"i -> ii"(diag)
    id = Matrix(I,n,n)
    reshape(h,n,n,1,1) .* reshape(id,1,1,n,n) .+ reshape(h,1,1,n,n) .* reshape(id,n,n,1,1)
end

@doc raw"
    TFIsing(hx::Real)

return a struct representing the transverse field ising model with magnetisation `hx`.
"
struct TFIsing{T<:Real} <: HamiltonianModel
    hx::T
end

"""
    hamiltonian(model::TFIsing)

return the transverse field ising hamiltonian for the provided `model` as a
two-site operator.
"""
function hamiltonian(model::TFIsing)
    hx = model.hx
    -2 * ein"ij,kl -> ijkl"(σz,σz) -
        hx/2 * ein"ij,kl -> ijkl"(σx, id2) -
        hx/2 * ein"ij,kl -> ijkl"(id2, σx)
end

@doc raw"
    Heisenberg(Jz::T,Jx::T,Jy::T) where {T<:Real}

return a struct representing the heisenberg model with magnetisation fields
`Jz`, `Jx` and `Jy`..
"
struct Heisenberg{T<:Real} <: HamiltonianModel
    Jz::T
    Jx::T
    Jy::T
end
Heisenberg() = Heisenberg(1.0,1.0,1.0)

"""
    hamiltonian(model::Heisenberg)

return the heisenberg hamiltonian for the `model` as a two-site operator.
"""
function hamiltonian(model::Heisenberg)
    h = model.Jz * ein"ij,kl -> ijkl"(σz, σz) -
        model.Jx * ein"ij,kl -> ijkl"(σx, σx) -
        model.Jy * ein"ij,kl -> ijkl"(σy, σy)
    h = ein"ijcd,kc,ld -> ijkl"(h,σx,σx')
    real(h ./ 2)
end


text:
using PointProcessDistance
using Test
using Distances

@testset "PointProcessDistance.jl" begin
    n = 100
    x = rand(2, n)
    dist = PeriodicEuclidean(ones(2))

    rs = 0.:0.01:sqrt(2) / 2
    k = K(rs, x; dist)
    d = D(rs, x, n - 1; dist)
    t = T(rs, x; dist)
    h = H(rs, x; dist)

end


text:
__precompile__()

module DevIL

import Libdl

const libIL = Libdl.find_library(["libIL", "DevIL"],
                                 ["/usr/lib/x86_64-linux-gnu", joinpath(dirname(pathof(DevIL)), "..", "deps", string(Sys.ARCH))])

# Stolen from getCFun macro
macro ilFunc(cFun)
    arguments = map(function (arg)
                        if isa(arg, Symbol)
                            arg = Expr(:(::), arg)
                        end
                        return arg
                    end, cFun.args[1].args[2:end])

    # Get info out of arguments of `cFun`
    argumentNames = map(arg->arg.args[1], arguments)
    returnType    = cFun.args[2]
    inputTypes    = map(arg->arg.args[2], arguments)

    # Construct the result.
	cName     = cFun.args[1].args[1]
	cSym      = Expr(:quote, cName)
	symAndLib = :($cSym, $libIL)
#=
    body       = Expr(:ccall, symAndLib, returnType, Expr(:tuple, inputTypes...), argumentNames...)
    func       = Expr(:function, Expr(:call, cName, argumentNames...), body)
	exportExpr = Expr(:export, cName)
	ret        = Expr(:block, func, exportExpr)
=#
    ret = quote
        $cName($(argumentNames...)) = ccall($symAndLib, $returnType, ($(inputTypes...),), $(argumentNames...))
        export $cName
    end

    return esc(ret)
end

macro ilConst(assignment)
	constExpr  = Expr(:const, assignment)
	exportExpr = Expr(:export, assignment.args[1])
	ret        = Expr(:block, constExpr, exportExpr)

	return esc(ret)
end

include("IL.jl")

end # module