text:
module KernelLearning

using StatsBase
using NearestNeighbors
using MLKernels
using LinearAlgebra
using Optim
using Serialization

using JuMP, Gurobi
using Random
using Logging
using LinearAlgebra
using Distributions
using RollingFunctions

using SVDD, OneClassActiveLearning
using ValueHistories
using MLLabelUtils

using Distances

include("init_heuristics/wang_combined.jl")
include("init_heuristics/xiao_gamma.jl")
include("init_heuristics/c_on_sample.jl")

methods_dir = joinpath(@__DIR__, "methods")
for f in readdir(methods_dir)
    include(joinpath(methods_dir, f))
end

include("alkernel/nnpool.jl")
include("alkernel/alignment.jl")
include("alkernel/al_kernel.jl")

include("fixes.jl")

export
    NNPool,

    ALKernel,
    ALRandom,
    WangHeuristic,
    SupervisedAlignment,
    SupervisedGridSearch,
    DFNHeuristic,
    QMSHeuristic,

    find_parameters,

    DFNGammaStrategy,
    RandomSampleGrid
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
include("../src/jdplyr.jl")

using DataFrames, Test, .jdplyr

df_a = DataFrame(one = [1, 2, 3, 4, 5], two = [6, 7, 8, 9, 10])
df_b = DataFrame(one = [1, 1, 2, 2, 2], two = [20, 21, 22, 23, 24])
df_c = DataFrame(one = [1, 2, 3], two = [1, 2, 3], three = [1, 2, 3])

@testset "Functions" begin
    @test _select(df_a, :one, :two) == df_a

    mutate_test = copy(df_a)
    mutate_test[!, :new] = df_a.one .+ 1
    @test _mutate(df_a, new = df_a.one .+ 1) == mutate_test

    @test _filter(df_a, df_a.one .> 2, df_a.two .> 8) ==
          df_a[(df_a.one.>2).&(df_a.two.>8), :]

    @test _arrange(df_a, :one, desc(:two)) ==
          sort(df_a, (:one, order(:two, rev = true)))

    @test _summarise(df_a, sum = :one => sum).sum[1] == 15
    @test _summarize(df_a, sum = :one => sum).sum[1] == 15
    @test _summarise(_group_by(df_b, :one), sum = :two => sum).sum[1] == 41

    @test _left_join(df_a, df_b, on = :one).two_1[1] ==
          join(df_a, df_b, on = :one, kind = :left, makeunique = true).two_1[1]

    @test (
        _gather(df_c, key = :variable, factor_cols = [:one]) |>
        x -> _select(x, :variable, :value) |> x -> _arrange(x, desc(:variable))
    ) == stack(df_c, [:two, :three])[!, [:variable, :value]]

    @test _head(df_a, 2) == df_a[1:2, :]
    @test _tail(df_a, 2) == df_a[end-1:end, :]

    rename_test = copy(df_a)
    rename!(rename_test, :one => :new_one)
    @test _rename(df_a, new_one=:one) == rename_test

    @test _slice(df_a, 1:2) == df_a[1:2, :]
end

ext: jl


text:
module test_fspace
using ApproXD
using Test

@testset "testing FSpaceXD" begin
	ndims = 4

	# bounds
	lb = [-1.1,-1.5,-0.9,-1.0]
	ub = [1.2,1.6,0.9,1]

	# number of eval points and basis functions:
	# we require square basis matrices!
	npoints = [13,13,13,13]

	# number of basis funcs
	nbasis = npoints

	# splien degrees
	degs = [3,3,3,3]

	# implies a number of knots for each spline
	# remember the restriction that nknots == ncoefs
	nknots = Dict(i => nbasis[i] - degs[i] + 1 for i=1:ndims)

	# eval points
	points = Dict(i => collect(range(lb[i],stop = ub[i], length = npoints[i])) for i=1:ndims)

	# set up ApproXD
	bsp = Dict{Integer,BSpline}()
	for i in 1:ndims
		bsp[i] = BSpline(nknots[i],degs[i],lb[i],ub[i])
	end

	# set of basis functions
	d = Dict{Integer,Array{Float64,2}}()
	for i=1:ndims
		d[i] = Array(getBasis(points[i],bsp[i]))
	end

	# set of INVERSE basis functions
	id = Dict{Integer,Array{Float64,2}}()
	for k in collect(keys(d))
		id[k] = inv(d[k])
	end

	#  get a function
	function f(x,y,z,w)
		sin(sqrt(x^2+y^2)) + (z-w)^3
	end

	y = Float64[f(i,j,k,w) for i in points[1], j in points[2], k in points[3], w in points[4]]

	yvec = y[:]

	# get coefs using the function
	mycoef = getTensorCoef(id,yvec)

	# setup the FSpace
	fx = FSpaceXD(ndims,mycoef,bsp)

	rval1 = lb[1] + 0.3
	rval2 = lb[2] + 0.23
	rval3 = lb[3] + 0.111
	rval4 = lb[4] + 0.099
	# println("approx value = $(getValue(fx,[rval1,rval2,rval3,rval4]))")
	# println("true value = $(f(rval1,rval2,rval3,rval4))")
	@test isapprox(getValue(fx,[rval1,rval2,rval3,rval4]), f(rval1,rval2,rval3,rval4),atol=3e-3)

end

end


text:
module LevelSetShapes

using AbbreviatedTypes
using LinearAlgebra
using ForwardDiff: derivative, gradient, jacobian

include("util.jl")
include("shape/shape.jl")

end # module

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
# установка пакетов, выполняется один раз!
using Pkg
Pkg.add("CSV")
Pkg.add("RDatasets")
Pkg.add("Econometrics")
Pkg.add("DataFrames")
# то, что выше, делаем один раз :)

# матрички :)
a = [5 7 8 -1] # строка
b = [5 7 8 -1]' # столбец
b = [5, 7, 8, -1] # столбец
b = [5; 7; 8; -1] # столбец

y = [3 4 5 -2]' # столбец зависимой переменной
X = [1 3; 1 7; 1 5; 1 5] # матрица регрессоров (предикторов)
X'y
X'X
beta_hat = (X'X)^(-1) * X'y

y .* y # поточечное (поэлеметное умножение)

f(x) = 2x^3 + 3x - 7
f(4)

[2x + 5y for x = 2:8, y in [5, -3, 9]]

using RDatasets, Econometrics # аналог import numpy...

psid = dataset("Ecdat", "PSID") # Panel Study of Income Dynamics
psid

psid.Age
psid."Age"
psid[!, "Age"]
psid[2:5, ["Age", "Educatn"]]

model_0 = fit(EconometricModel,
    @formula(log(1 + Earnings) ~ Age + Educatn), psid)


text:
# Julia wrapper for header: /usr/include/scip/pub_bandit_ucb.h
# Automatically generated using Clang.jl wrap_c

function SCIPcreateBanditUcb(scip, ucb, priorities, alpha, nactions, initseed)
    ccall((:SCIPcreateBanditUcb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_BANDIT}}, Ptr{Cdouble}, Cdouble, Cint, UInt32), scip, ucb, priorities, alpha, nactions, initseed)
end

function SCIPgetConfidenceBoundUcb(ucb, action)
    ccall((:SCIPgetConfidenceBoundUcb, libscip), Cdouble, (Ptr{SCIP_BANDIT}, Cint), ucb, action)
end

function SCIPgetStartPermutationUcb(ucb)
    ccall((:SCIPgetStartPermutationUcb, libscip), Ptr{Cint}, (Ptr{SCIP_BANDIT},), ucb)
end


text:
using MPI
MPI.Init()
comm = MPI.COMM_WORLD
print("Hello world, I am rank $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))\n")
MPI.Barrier(comm)
MPI.Finalize()


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

@testset ExtendedTestSet "WarpAffine" begin
    tfm = WarpAffine(0.1)
    image = Image(rand(RGB, 50, 50))

    @test_nowarn apply(tfm, image)
    @test sum.(getbounds(apply(tfm, image)).rs) != (50, 500)
    buf = apply(tfm, image)
    ctfm = tfm |> CenterCrop((50, 50))
    @test_nowarn apply(ctfm, image)
    buf = apply(ctfm, image)
    @test_nowarn apply!(buf, ctfm, image)
end


text:
#=
BMM.jl

10/27/2015
Adham Beyki, odinay@gmail.com
=#

##########################################
###       Bayesian Mixture Model       ###
##########################################
immutable BMM{T} <: MixtureModel
    component::T
    KK::Int
    aa::Float64

    BMM{T}(c::T, KK::Int, aa::Float64) = new(c, KK, aa)
end
BMM{T}(c::T, KK::Int, aa::Real) = BMM{typeof(c)}(c, KK, convert(Float64, aa))

function Base.show(io::IO, bmm::BMM)
    println(io, "Finite Mixture Model with $(bmm.KK) $(typeof(bmm.component)) components")
end

function storesample(zz::Vector{Int}, sample_n::Int, filename::ASCIIString)

    if endswith(filename, "_")
        dummy_filename = string(filename, sample_n, ".h5")
    else
        dummy_filename = string(filename, "_", sample_n, ".h5")
    end

    println("\nstoring on disk...\n")
    HDF5.h5open(dummy_filename, "w") do file
        HDF5.write(file, "zz", zz)
        HDF5.write(file, "sample_n", sample_n)
    end
end

function collapsed_gibbs_sampler{T1, T2}(
    bmm::BMM{T1},
    xx::Vector{T2},
    zz::Vector{Int},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    store_every::Int=250, results_path="", filename::ASCIIString="BMM_results")

    n_iterations   = n_burnins + (n_samples)*(n_lags+1)
    NN             = length(xx)
    nn             = zeros(Int, bmm.KK)
    pp             = zeros(Float64, bmm.KK)
    log_likelihood = 0.0

    # constructing the components
    components = Array(typeof(bmm.component), bmm.KK)
    for kk = 1:bmm.KK
        components[kk] = deepcopy(bmm.component)
    end

    # initializing the model
    tic()
    for ii = 1:NN
        kk = zz[ii]
        additem!(components[kk], xx[ii])
        nn[kk] += 1
        log_likelihood += loglikelihood(components[kk], xx[ii])
    end
    elapsed_time = toq()

    # starting the MCMC chain
    for iteration = 1:n_iterations

        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, time=%.2f, likelihood=%.2f", iteration, bmm.KK, elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()
        for ii = randperm(NN)

            # 1
            # remove the datapoint
            kk = zz[ii]
            nn[kk] -= 1
            delitem!(components[kk], xx[ii])

            # 2
            # sample zz
            for kk = 1:bmm.KK
                pp[kk] = log(nn[kk] + bmm.aa) + logpredictive(components[kk], xx[ii])
            end
            lognormalize!(pp)
            kk = sample(pp)

            # 3
            # add the datapoint to the newly sampled cluster
            zz[ii] = kk
            nn[kk] += 1
            additem!(components[kk], xx[ii])
            log_likelihood += loglikelihood(components[kk], xx[ii])
        end # ii
        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(zz, sample_n, normalized_filename)
            end
        end # storesample

    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(zz, sample_n, normalized_filename)
end

function posterior{T1, T2}(bmm::BMM{T1}, xx::Vector{T2}, zz::Vector{Int})

    components = Array(typeof(bmm.component), bmm.KK)
    for kk = 1:bmm.KK
        components[kk] = deepcopy(bmm.component)
    end

    NN = length(zz)
    nn = zeros(Int, bmm.KK)

    for ii = 1:NN
        kk = zz[ii]
        additem!(components[kk], xx[ii])
        nn[kk] += 1
    end

    return([posterior(components[kk]) for kk =1:bmm.KK], nn)
end


text:
module Numkin2019

    hello( name ) = " hello $name "

end


text:
function sparse_random_factory(;
    states=10,
    state_type=Float64,
    actions=2,
    action_type=Float64,
)
    return MDPs.Examples.sprandom(state_type, action_type, states, actions)
end

"sparse random called without type arguments  should return a vector of "
"sparse Float64 transition matrices"
let
    transition, reward = MDPs.Examples.sprandom(10, 2)
    @test typeof(transition) == Vector{SparseMatrixCSC{Float64,Int64}}
end

"sparse random should return a sparse Float64 reward"
let
    transition, reward = sparse_random_factory()
    @test typeof(reward) == SparseMatrixCSC{Float64,Int64}
end

"sparse random with Float32 types should return Float32 sparse matrices"
let
    transition, reward = sparse_random_factory(
        state_type=Float32,
        action_type=Float32,
    )
    @test typeof(transition) == Vector{SparseMatrixCSC{Float32,Int64}}
    @test typeof(reward) == SparseMatrixCSC{Float32,Int64}
end

"sparse random transition matrix should be square"
let
    transition, reward = sparse_random_factory()
    rows, columns = size(transition[1])
    @test rows == columns
end

"sparse random transition rows should sum to one"
let
    transition, reward = sparse_random_factory()
    @test_approx_eq sum(transition[1], 1) ones(size(transition[1])[1])
end

"sparse random transition probabilities should be non-negative"
let
    transition, reward = sparse_random_factory()
    @test all(nonzeros(transition[1]) .>= 0)
end

"sparse random transition matrices should be randomly generated"
let
    transition1, reward = sparse_random_factory()
    transition2, reward = sparse_random_factory()
    @test transition1[1] != transition2[1]
end


text:
__precompile__()

"""

A module for solving linear boundary value problems by the difference method using
the Thomas algorithm method.

See `Problem`, `Options`, and `solve` for more info.

"""
baremodule A1

using Base

# Load the internal logic
Base.include(A1, "Internal/Internal.jl")

# Export the second-level components
using .Internal

# Export the first-level components
export Problem, Options, solve

end


text:
"""
SimpleCachingSettings is a struct containing settings for exported macros `@scache` and
`@scachefast`.

## Fields

- `log`: a Bool indicating whether to log or not operations like computing, saving and
loading. Default is `false`.
- `log_output`: a IO containing the file descriptor the log will be written to; ignored if
`log` is `false`. Default is `stdout`.
- `date_format`: a String representing the format the DateTime will be printed when logging;
ignored if `log` is `false`. Default is `"dd/mm/yyyy HH:MM:SS"`.
- `line_started`: a String that will be placed at the beginning of the log; ignored if `log`
is `false`. Default is `"● "`.
- `create_cache_record`: a Bool indicating whether to create a tsv file containing human
readable informations about the saved file. Default is `true`.
"""
mutable struct SimpleCachingSettings
    log::Bool
    log_output::IO
    date_format::AbstractString
    line_starter::AbstractString
    create_cache_record::Bool
end

const settings = SimpleCachingSettings(false, stdout, "dd/mm/yyyy HH:MM:SS", "● ", true)

_use_serialize = false

_closelog(io::IO) = nothing
_closelog(io::IOStream) = close(io)
_closelog(io::Channel) = close(io)

_closelog() = _closelog(settings.log_output)


text:
# This file was generated, do not modify it.

using MLJ, Statistics
X, y = @load_iris;

@load DecisionTreeClassifier

tree_model = DecisionTreeClassifier()
tree = machine(tree_model, X, y)

train, test = partition(eachindex(y), 0.7, shuffle=true)
test[1:3]

fit!(tree, rows=train)

tree.fitresult

ŷ = predict(tree, rows=test)
@show ŷ[1]

ȳ = predict_mode(tree, rows=test)
@show ȳ[1]
@show mode(ŷ[1])

mce = cross_entropy(ŷ, y[test]) |> mean
round(mce, digits=4)

v = [1, 2, 3, 4]
stand_model = UnivariateStandardizer()
stand = machine(stand_model, v)

fit!(stand)
w = transform(stand, v)
@show round.(w, digits=2)
@show mean(w)
@show std(w)

vv = inverse_transform(stand, w)
sum(abs.(vv .- v))

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
export goldensection, hyperband

"""

    goldensection(f,n;kwargs) => (fmin,xmin)

Find the minimum of `f` using concurrent golden section search in `n`
dimensions. See `Knet.goldensection_demo()` for an example.

`f` is a function from a `Vector{Float64}` of length `n` to a
`Number`.  It can return `NaN` for out of range inputs.  Goldensection
will always start with a zero vector as the initial input to `f`, and
the initial step size will be 1 in each dimension.  The user should
define `f` to scale and shift this input range into a vector
meaningful for their application. For positive inputs like learning
rate or hidden size, you can use a transformation such as `x0*exp(x)`
where `x` is a value `goldensection` passes to `f` and `x0` is your
initial guess for this value. This will effectively start the search
at `x0`, then move with multiplicative steps.

I designed this algorithm combining ideas from [Golden Section
Search](http://apps.nrbook.com/empanel/index.html?pg=492) and [Hill
Climbing Search](https://en.wikipedia.org/wiki/Hill_climbing). It
essentially runs golden section search concurrently in each dimension,
picking the next step based on estimated gain.

# Keyword arguments
* `dxmin=0.1`: smallest step size.
* `accel=φ`: acceleration rate. Golden ratio `φ=1.618...` is best.
* `verbose=false`: use `true` to print individual steps.
* `history=[]`: cache of `[(x,f(x)),...]` function evaluations.

"""
function goldensection(f,n; dxmin=0.1, accel=1.618033988749895, history=[], verbose=false)

    function feval(x)           # so we don't repeat function evaluations
        for (k,v) in history
            if isapprox(x,k)
                return v
            end
        end
        fx = f(x)
        push!(history, (x,fx))
        return fx
    end

    function setindex(x,v,d)    # non-mutating setindex
        y = copy(x)
        y[d] = v
        return y
    end

    x0 = zeros(n)               # initial point
    f0 = feval(x0)              # initial value
    dx = ones(n)                # step sizes
    df = Inf * ones(n)          # expected gains
    while maximum(abs.(dx)) > dxmin
        i = argmax(df)
        x1 = setindex(x0,x0[i]+dx[i],i)
        f1 = feval(x1)
        if verbose; println((; f0,x0,f1,x1,dx,df)); end
        isnan(f1) && (f1=f0+df[i])
        if f1 < f0
            dx[i] = accel * dx[i]
            df[i] = accel * (f0-f1)
            x0,f0 = x1,f1
            for j = 1:length(df)
                if abs(dx[j]) < dxmin * accel
                    dx[j] = sign(dx[j]) * dxmin * accel
                    df[j] = max(df[j],0) # max(df[j],-1-df[j])
                end
            end
        else
            dx[i] = -dx[i] / accel
            df[i] = (f1-f0) / accel
            if abs(dx[i]) < dxmin
                df[i] = -1 # -1-df[i]
            end
        end
    end
    return (f0,x0)
end

function goldensection_demo()
    include(Knet.dir("examples","mnist.jl"))
    neval = 0

    function f(x)
        neval += 1
        winit,lr,hidden = xform(x)
        if hidden < 10000
            w = MNIST.main("--winit $winit --lr $lr --hidden $hidden --seed 1 --epochs 1 --fast")
            corr,loss = MNIST.accuracy(w,MNIST.dtst)
        else
            corr,loss = NaN,NaN # prevent huge weights
        end
        println((:neval,neval,:loss,loss,:corr,corr,:winit,winit,:lr,lr,:hidden,hidden))
        return loss
    end

    function xform(x)
        winit,lr,hidden = exp(x) .* [ 0.01, 0.01, 100.0 ]
        hidden = ceil(Int, hidden)
        (winit,lr,hidden)
    end

    goldensection(f,3)
end

"""

    hyperband(getconfig, getloss, maxresource=27, reduction=3)

Hyperparameter optimization using the hyperband algorithm from ([Lisha
et al. 2016](https://arxiv.org/abs/1603.06560)).  You can try a simple
MNIST example using `Knet.hyperband_demo()`.

# Arguments
* `getconfig()` returns random configurations with a user defined type and distribution.
* `getloss(c,n)` returns loss for configuration `c` and number of resources (e.g. epochs) `n`.
* `maxresource` is the maximum number of resources any one configuration should be given.
* `reduction` is an algorithm parameter (see paper), 3 is a good value.

"""
function hyperband(getconfig, getloss, maxresource=27, reduction=3)
    smax = floor(Int, log(maxresource)/log(reduction))
    B = (smax + 1) * maxresource
    best = (Inf,)
    for s in smax:-1:0
        n = ceil(Int, (B/maxresource)*((reduction^s)/(s+1)))
        r = maxresource / (reduction^s)
        curr = halving(getconfig, getloss, n, r, reduction, s)
        if curr[1] < best[1]; (best=curr); end
    end
    return best
end

# TODO: document Successive/Sequential Halving:
# http://www.jmlr.org/proceedings/papers/v51/jamieson16.pdf,
# http://www.jmlr.org/proceedings/papers/v28/karnin13.pdf,
# Successive  Reject:
# http://certis.enpc.fr/~audibert/Mes%20articles/COLT10.pdf
function halving(getconfig, getloss, n, r=1, reduction=3, s=round(Int, log(n)/log(reduction)))
    best = (Inf,)
    T = [ getconfig() for i=1:n ]
    for i in 0:s
        ni = floor(Int,n/(reduction^i))
        ri = r*(reduction^i)
        # println((:s,s,:n,n,:r,r,:i,i,:ni,ni,:ri,ri,:T,length(T)))
        L = [ getloss(t, ri) for t in T ]
        l = sortperm(L); l1=l[1]
        L[l1] < best[1] && (best = (L[l1],ri,T[l1])) # ;println("best1: $best"))
        T = T[l[1:floor(Int,ni/reduction)]]
    end
    # println("best2: $best")
    return best
end

function hyperband_demo()
    include(Knet.dir("examples","mnist.jl"))
    best = (Inf,)
    neval = 0

    function getloss(config,epochs)
        neval += 1
        winit,lr,hidden = config
        epochs = round(Int,epochs)
        w = MNIST.main("--winit $winit --lr $lr --hidden $hidden --epochs $epochs --seed 1 --fast")
        corr,loss = MNIST.accuracy(w,MNIST.dtst)
        println((:epochs,epochs,:corr,corr,:loss,loss,:winit,winit,:lr,lr,:hidden,hidden))
        if loss < best[1]
            best = (loss, config, epochs)
        end
        return loss
    end

    function getconfig()
        winit = 0.001^rand()
        lr = 0.001^rand()
        hidden = 16 + floor(Int, 10000^rand())
        return (winit,lr,hidden)
    end

    hyperband(getconfig, getloss)
    println((:neval,neval,:minloss,best[1],:epochs,best[3],:winit,best[2][1],:lr,best[2][2],:hidden,best[2][3]))
end


text:
"""
    render([io], tokens, view)
    render([io], tokens; kwargs...)
    (tokens::MustacheTokens)([io]; kwargs...)

Render a set of tokens with a view, using optional `io` object to print or store.

Arguments
---------

* `io::IO`: Optional `IO` object.

* `tokens`: Either Mustache tokens, or a string to parse into tokens

* `view`: A view provides a context to look up unresolved symbols
  demarcated by mustache braces. A view may be specified by a
  dictionary, a module, a composite type, a vector, a named tuple, a
  data frame, a `Tables` object, or keyword arguments.

"""
function render(io::IO, tokens::MustacheTokens, view)
    _writer = Writer()
    render(io, _writer, tokens, view)
end
function render(io::IO, tokens::MustacheTokens; kwargs...)
    render(io, tokens, kwargs)
end

render(tokens::MustacheTokens, view) = sprint(io -> render(io, tokens, view))
render(tokens::MustacheTokens; kwargs...) = sprint(io -> render(io, tokens, Dict(kwargs...)))

## make MustacheTokens callable for kwargs...
(m::MustacheTokens)(io::IO; kwargs...) = render(io, m; kwargs...)
(m::MustacheTokens)(; kwargs...) = sprint(io -> m(io; kwargs...))

## Exported call without first parsing tokens via mt"literal"
##
## @param template a string containing the template for expansion
## @param view a Dict, Module, CompositeType, DataFrame holding variables for expansion
function render(io::IO, template::AbstractString, view; tags= ("{{", "}}"))
    _writer = Writer()
    render(io, _writer, parse(template, tags), view)
end
function render(io::IO, template::AbstractString; kwargs...)
    _writer = Writer()
    render(io, _writer, parse(template), Dict(kwargs...))
end
render(template::AbstractString, view; tags=("{{", "}}")) = sprint(io -> render(io, template, view, tags=tags))
render(template::AbstractString; kwargs...) = sprint(io -> render(io, template, Dict(kwargs...)))

# Exported, but should be deprecated....
"""
    render_from_file(filepath, view)
    render_from_file(filepath; kwargs...)

Renders a template from `filepath` and `view`.
"""
render_from_file(filepath, view) = render(Mustache.load(filepath), view)
render_from_file(filepath::AbstractString; kwargs...) = render(Mustache.load(filepath); kwargs...)


text:
Wigner  = PyCall.PyNULL()
Optics  = PyCall.PyNULL()
Spin = PyCall.PyNULL()

function init_physics()
    copy!(Wigner,  PyCall.pyimport_conda("sympy.physics.wigner",       "sympy"))
    copy!(Optics,  PyCall.pyimport_conda("sympy.physics.optics",       "sympy"))
    copy!(Spin ,   PyCall.pyimport_conda("sympy.physics.quantum.spin", "sympy"))
end

module Physics

using PyCall
using SymPy

physics = [(SymPy.Wigner, [:clebsch_gordan,
                           :dot_rot_grad_Ynm,
                           :gaunt,
                           :racah,
                           :wigner_3j,
                           :wigner_6j,
                           :wigner_9j]),
           (SymPy.Optics, [:RayTransferMatrix,
                           :FreeSpace,
                           :CurvedRefraction,
                           :FlatMirror]),
           (SymPy.Spin  , [:WignerD])
]

for (m, meths) in physics
    for meth in meths
        meth_name = string(meth)
        @eval begin
            # @doc """
            # `$($meth_name)`: a SymPy function.
            #     The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
            #     """ ->
            ($meth)(args...;kwargs...) = pycall(getproperty($m, $meth_name), PyAny, map(Sym, args)...,
                                                [k=>Sym(v) for (k,v) in pairs(kwargs)]...)
        end
        eval(Expr(:export, meth))
    end
end

end


text:
using Gomen
using Base.MathConstants, Random, Statistics, Test
using StaticArrays
using LightGraphs, LightGraphs.SimpleGraphs

@testset "Game" begin
    @testset "Construction" begin
        @test_throws DomainError Game(-1, 1)
        @test_throws DomainError Game(1, 1)
        @test_throws DomainError Game(0, 0)
        @test_throws DomainError Game(0, 2)

        @test Game(0.5, 1).payoffs == @SArray [1.0 0.5; 1.0 0.0]
        @test Game(-0.3, 1.3).payoffs == @SArray [1.0 -0.3; 1.3 0.0]
    end

    @testset "Play" begin
        let g = Game(-0.3, 1.3)
            @test play(g, 1, 1) == 1.0
            @test play(g, 1, 2) == -0.3
            @test play(g, 2, 1) == 1.3
            @test play(g, 2, 2) == 0.0
        end
    end
end

@testset "Games" begin
    @testset "length" begin
        @test length(Games(0.5, 0.5)) == 9
        @test length(Games(0.1, 0.1)) == 121
        @test length(Games(0.01, 0.01)) == 10201
        @test length(Games(0.3, 0.3)) == 16
    end

    @testset "collect" begin
        @test map(g -> g.payoffs, Games(1, 1)) == [
                    SMatrix{2,2}([1.0 -0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.5 0.0])
                ]

        @test map(g -> g.payoffs, Games(0.5, 1)) == [
                    SMatrix{2,2}([1.0 -0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.0; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.0; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.5 0.0])
                ]

        @test map(g -> g.payoffs, Games(1, 0.5)) == [
                    SMatrix{2,2}([1.0 -0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.0 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.0 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.5 0.0])
                ]

        @test map(g -> g.payoffs, Games(0.5, 0.5)) == [
                    SMatrix{2,2}([1.0 -0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.0; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.0 0.0]),
                    SMatrix{2,2}([1.0  0.0; 1.0 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.0 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.0; 1.5 0.0]),
                    SMatrix{2,2}([1.0  0.5; 1.5 0.0])
                ]

        @test map(g -> g.payoffs, Games(0.3, 0.3)) == [
                    SMatrix{2,2}([1.0 -0.5; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.2; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.1; 0.5 0.0]),
                    SMatrix{2,2}([1.0  0.4; 0.5 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 0.8 0.0]),
                    SMatrix{2,2}([1.0 -0.2; 0.8 0.0]),
                    SMatrix{2,2}([1.0  0.1; 0.8 0.0]),
                    SMatrix{2,2}([1.0  0.4; 0.8 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.1 0.0]),
                    SMatrix{2,2}([1.0 -0.2; 1.1 0.0]),
                    SMatrix{2,2}([1.0  0.1; 1.1 0.0]),
                    SMatrix{2,2}([1.0  0.4; 1.1 0.0]),
                    SMatrix{2,2}([1.0 -0.5; 1.4 0.0]),
                    SMatrix{2,2}([1.0 -0.2; 1.4 0.0]),
                    SMatrix{2,2}([1.0  0.1; 1.4 0.0]),
                    SMatrix{2,2}([1.0  0.4; 1.4 0.0])
                ]
    end
end

@testset "Rule" begin
    @testset "Sigmoid" begin
        let points = [-1.0, -0.5, 0.0, 0.5, 1.0], p = sqrt(e), q = e^0.25
            @test_throws DomainError Sigmoid(-1.0)
            let rule = Sigmoid()
                expect = [1/(1 + e), 1/(1 + p), 0.5, p/(1 + p), e/(1 + e)]
                @test map(dp -> apply(rule, dp), points) ≈ expect
            end
            let rule = Sigmoid(0.5)
                expect = [1/(1 + p), 1/(1 + q), 0.5, q/(1 + q), p/(1 + p)]
                @test map(dp -> apply(rule, dp), points) ≈ expect
            end
        end
    end

    @testset "Heaviside" begin
        let points = [-1.0, -0.5, 0.0, 0.5, 1.0]
            @test_throws DomainError Heaviside(-1.0)
            let rule = Heaviside()
                @test map(dp -> apply(rule, dp), points) ≈ [0.0, 0.0, 0.5, 1.0, 1.0]
            end
            let rule = Heaviside(0.75)
                @test map(dp -> apply(rule, dp), points) ≈ [0.0, 0.5, 0.5, 0.5, 1.0]
            end
        end
    end
end

@testset "Scheme" begin
    game = Game(0.3, 1.3)
    graph = first(CycleGraphGenerator(1,3))
    @testset "Single Agent" begin
        let N = 1000
            σ(p) = sqrt(N*p*(1-p))
            cond(p) = N*p - 3σ(p), N*p + 3σ(p)

            seed = Random.rand(UInt32)
            println("Random Seed: ", seed)
            Random.seed!(seed)

            let cf = CounterFactual()
                arena = Arena(game, graph, cf)

                a, b = cond(0.5)
                @test a ≤ sum(Int(decide(cf, arena, 1, 0.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 0.0) == 1) for _ in 1:N) ≤ b

                a, b = cond(e / (1 + e))
                @test a ≤ sum(Int(decide(cf, arena, 1, 1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 1.0) == 1) for _ in 1:N) ≤ b

                a, b = cond(1 / (1 + e))
                @test a ≤ sum(Int(decide(cf, arena, 1, -1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, -1.0) == 1) for _ in 1:N) ≤ b
            end

            let cf = CounterFactual(Sigmoid(0.5))
                arena = Arena(game, graph, cf)

                a, b = cond(0.5)
                @test a ≤ sum(Int(decide(cf, arena, 1, 0.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 0.0) == 1) for _ in 1:N) ≤ b

                a, b = cond(sqrt(e) / (1 + sqrt(e)))
                @test a ≤ sum(Int(decide(cf, arena, 1, 1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 1.0) == 1) for _ in 1:N) ≤ b

                a, b = cond(1 / (1 + sqrt(e)))
                @test a ≤ sum(Int(decide(cf, arena, 1, -1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, -1.0) == 1) for _ in 1:N) ≤ b
            end

            let cf = CounterFactual(Heaviside())
                arena = Arena(game, graph, cf)

                a, b = cond(0.5)
                @test a ≤ sum(Int(decide(cf, arena, 1, 0.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 0.0) == 1) for _ in 1:N) ≤ b

                @test a ≤ sum(Int(decide(cf, arena, 1, 1e-4) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 1e-4) == 1) for _ in 1:N) ≤ b

                @test a ≤ sum(Int(decide(cf, arena, 1, -1e-4) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, -1e-4) == 1) for _ in 1:N) ≤ b

                a, b = cond(1.0)
                @test a ≤ sum(Int(decide(cf, arena, 1, 1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, 1.0) == 1) for _ in 1:N) ≤ b

                a, b = cond(0.0)
                @test a ≤ sum(Int(decide(cf, arena, 1, -1.0) == 2) for _ in 1:N) ≤ b
                @test a ≤ sum(Int(decide(cf, arena, 2, -1.0) == 1) for _ in 1:N) ≤ b
            end
        end
    end

    @testset "Multiple Agents" begin
        let cf = CounterFactual(Heaviside())
            arena = Arena(game, graph, cf)

            @test decide(cf, arena, [1, 2], [1., -1.]) == [2, 2]
            @test decide(cf, arena, [1, 1], [1., -1.]) == [2, 1]
        end
    end

    @testset "From Payoffs" begin
        let cf = CounterFactual(Heaviside())
            arena = Arena(game, graph, cf)

            @test_throws ArgumentError decide(cf, arena, [0, 2], Float64[0 0; 0 0])
            @test_throws ArgumentError decide(cf, arena, [3, 2], Float64[0 0; 0 0])
            @test_throws ArgumentError decide(cf, arena, [1, 1], Float64[0 0 0; 0 0 0])
            @test_throws ArgumentError decide(cf, arena, [1, 1, 1], Float64[0 0; 0 0])
            @test_throws ArgumentError decide(cf, arena, [1, 1], Float64[0 0])
            @test_throws ArgumentError decide(cf, arena, [1, 1], Float64[0 0; 0 0; 0 0])

            @test decide(cf, arena, [1, 1], Float64[1 1; 0 0]) == [1, 1]
            @test decide(cf, arena, [1, 1], Float64[1 0; 0 1]) == [1, 2]
            @test decide(cf, arena, [1, 2], Float64[1 1; 0 0]) == [1, 1]
            @test decide(cf, arena, [1, 2], Float64[1 0; 0 1]) == [1, 2]
            @test decide(cf, arena, [2, 1], Float64[1 1; 0 0]) == [1, 1]
            @test decide(cf, arena, [2, 1], Float64[1 0; 0 1]) == [1, 2]
            @test decide(cf, arena, [2, 2], Float64[1 1; 0 0]) == [1, 1]
            @test decide(cf, arena, [2, 2], Float64[1 0; 0 1]) == [1, 2]
        end
    end
end

@testset "Arena" begin
    @testset "Construction" begin
        let game = Game(0, 1)
            scheme = CounterFactual()
            @test_throws DomainError Arena(game, SimpleGraph(), scheme)
            @test_throws DomainError Arena(game, SimpleGraph(1), scheme)
            let graph = SimpleGraph(2)
                add_edge!(graph, 1, 1)
                add_edge!(graph, 2, 2)
                @test_throws DomainError Arena(game, graph, scheme)
            end
            let graph = SimpleGraph(2)
                add_edge!(graph, 1, 1)
                add_edge!(graph, 1, 2)
                @test_throws DomainError Arena(game, graph, scheme)
            end
            let graph = SimpleGraph(3)
                add_edge!(graph, 1, 2)
                @test_throws DomainError Arena(game, graph, scheme)
            end

            @test length(Arena(game, cycle_graph(10), scheme)) == 10
        end
    end

    @testset "payoffs" begin
        let arena = Arena(Game(0, 1), path_graph(2), CounterFactual())
            @test_throws ArgumentError payoffs(arena, Int[])
            @test_throws ArgumentError payoffs(arena, [1])
            @test_throws ArgumentError payoffs(arena, [1, 2, 1])
            @test_throws ArgumentError payoffs(arena, [0, 1])
            @test_throws ArgumentError payoffs(arena, [3, 1])

            @test payoffs(arena, [1, 1]) == [1 1; 1 1]
            @test payoffs(arena, [1, 2]) == [0 1; 0 1]
            @test payoffs(arena, [2, 1]) == [1 0; 1 0]
            @test payoffs(arena, [2, 2]) == [0 0; 0 0]
        end

        let arena = Arena(Game(0.25, 1.25), path_graph(2), CounterFactual())
            @test payoffs(arena, [1, 1]) == [1.00 1.00; 1.25 1.25]
            @test payoffs(arena, [1, 2]) == [0.25 1.00; 0.00 1.25]
        end

        let game = Game(0.25, 1.25)
            graph = first(GridGraphGenerator(; count=1, dims=[2,2]))
            scheme = CounterFactual(Heaviside())
            arena = Arena(game, graph, scheme)

            @test payoffs(arena, [1, 1, 1, 1]) == [2.00 2.00 2.00 2.00; 2.50 2.50 2.50 2.50]
            @test payoffs(arena, [1, 2, 1, 1]) == [1.25 2.00 2.00 1.25; 1.25 2.50 2.50 1.25]
            @test payoffs(arena, [2, 2, 2, 2]) == [0.50 0.50 0.50 0.50; 0.00 0.00 0.00 0.00]
        end
    end

    @testset "round" begin
        seed = Random.rand(UInt32)
        println("Random Seed: ", seed)
        Random.seed!(seed)

        let game = Game(0.25, 1.25)
            graph = first(GridGraphGenerator(; count=1, dims=[2,2]))
            scheme = CounterFactual(Heaviside())
            arena = Arena(game, graph, scheme)

            N = 1000
            σ(p) = sqrt(N*p*(1-p))
            cond(p) = N*p - 3σ(p), N*p + 3σ(p)

            @test play(arena, [1, 1, 1, 1]) == [2, 2, 2, 2]
            @test play(arena, [2, 2, 2, 2]) == [1, 1, 1, 1]

            a, b = cond(0.5)
            dist = sum(Array{Int}(play(arena, [1, 2, 1, 1]) .== [1, 2, 2, 1]) for _ in 1:N)
            @test all([a, N, N, a] .<= dist .<= [b, N, N, b])
        end
    end
end

@testset "ROC" begin
    @testset "Construction" begin
        @test fpr(ROC()) == [0.0, 1.0]
        @test tpr(ROC()) == [0.0, 1.0]

        @test_throws DimensionMismatch ROC([0.5], Float64[])
        @test_throws DimensionMismatch ROC(Float64[], [0.5])

        @test_throws ArgumentError ROC([0.1, 0.2], [0.5, 0.4])
        @test_throws ArgumentError ROC([0.2, 0.1], [0.4, 0.5])

        @test_throws ArgumentError ROC([-0.5, 0.5], [0.5, 0.8])
        @test_throws ArgumentError ROC([0.5, -0.5], [0.8, 0.5])
        @test_throws ArgumentError ROC([0.5, 0.8], [-0.5, 0.5])
        @test_throws ArgumentError ROC([0.8, 0.5], [0.5, -0.5])

        @test_throws ArgumentError ROC([0.5, 2.0], [0.5, 0.8])
        @test_throws ArgumentError ROC([2.0, 0.5], [0.8, 0.5])
        @test_throws ArgumentError ROC([0.5, 0.8], [0.5, 2.0])
        @test_throws ArgumentError ROC([0.8, 0.5], [2.0, 0.5])

        let r = ROC([0.3, 0.8], [0.5, 0.9])
            @test length(fpr(r)) == length(tpr(r)) == 4
            @test fpr(r)[1] ≈ tpr(r)[1] ≈ 0.0
            @test fpr(r)[end] ≈ tpr(r)[end] ≈ 1.0
            @test issorted(fpr(r))
            @test issorted(tpr(r))
        end

        let r = ROC([0.8, 0.3], [0.9, 0.5])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 0.8, 0.3], [0.5, 0.9, 0.5])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 0.8, 0.3, 0.0], [0.5, 0.9, 0.5, 0.2])
            @test fpr(r) ≈ [0.0, 0.0, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.2, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 0.8, 0.3, 0.2], [0.5, 0.9, 0.5, 0.0])
            @test fpr(r) ≈ [0.0, 0.2, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.0, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 0.8, 0.3, 0.0], [0.5, 0.9, 0.5, 0.0])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 1.0, 0.8, 0.3], [0.5, 0.95, 0.9, 0.5])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 1.0, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 0.95, 1.0]
        end

        let r = ROC([0.3, 0.95, 0.8, 0.3], [0.5, 1.0, 0.9, 0.5])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 0.95, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 1.0, 1.0]
        end

        let r = ROC([0.3, 1.0, 0.8, 0.3], [0.5, 1.0, 0.9, 0.5])
            @test fpr(r) ≈ [0.0, 0.3, 0.8, 1.0]
            @test tpr(r) ≈ [0.0, 0.5, 0.9, 1.0]
        end

        let r = ROC([0.3, 0.3, 0.4, 0.4], [0.4, 0.2, 0.8, 0.6])
            @test fpr(r) ≈ [0.0, 0.3, 0.3, 0.4, 0.4, 1.0]
            @test tpr(r) ≈ [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
        end
    end

    @testset "mean" begin
        let r = ROC(), s = ROC()
            μ = mean([r, s])
            @test fpr(μ) ≈ [0.0, 1.0]
            @test tpr(μ) ≈ [0.0, 1.0]
        end

        let r = ROC(), s = ROC([0.5], [0.5])
            μ1, μ2 = mean([r, s]), mean([s, r])
            @test fpr(μ1) ≈ fpr(μ2) ≈ [0.0, 0.5, 1.0]
            @test tpr(μ1) ≈ tpr(μ2) ≈ [0.0, 0.5, 1.0]
        end

        let r = ROC(), s = ROC([0.5], [0.0])
            μ1, μ2 = mean([r, s]), mean([s, r])
            @test fpr(μ1) ≈ fpr(μ2) ≈ [0.0, 0.5, 1.0]
            @test tpr(μ1) ≈ tpr(μ2) ≈ [0.0, 0.25, 1.0]
        end

        let r = ROC([0.5], [0.0]), s = ROC([0.5], [0.5])
            μ1, μ2 = mean([r, s]), mean([s, r])
            @test fpr(μ1) ≈ fpr(μ2) ≈ [0.0, 0.5, 1.0]
            @test tpr(μ1) ≈ tpr(μ2) ≈ [0.0, 0.25, 1.0]
        end

        let r = ROC([0.25, 0.75], [0.3, 0.8]), s = ROC([0.5], [0.25])
            μ1, μ2 = mean([r, s]), mean([s, r])
            @test fpr(μ1) ≈ fpr(μ2) ≈ [0.0, 0.25, 0.5, 0.75, 1.0]
            @test tpr(μ1) ≈ tpr(μ2) ≈ [0.0, 0.2125, 0.4, 0.7125, 1.0]
        end
    end

    @testset "AUC" begin
        @test auc(ROC()) ≈ 0.5
        @test auc(ROC([0.5], [0.5])) ≈ 0.5
        @test auc(ROC([0.25, 0.6], [0.25, 0.6])) ≈ 0.5

        @test auc(ROC([0.5], [0.0])) ≈ 0.25
        @test auc(ROC([0.5], [1.0])) ≈ 0.75

        @test auc(ROC([0.2, 0.3, 0.8, 0.9], [0.1, 0.4, 0.85, 0.95])) ≈ 0.5350
    end
end


text:
struct BoundUpdate{S<:Union{LT,GT}}
    cvi::CVI
    s::S
end

mutable struct Node
    lt_bounds::Vector{BoundUpdate{LT}}
    gt_bounds::Vector{BoundUpdate{GT}}
    lt_general_constrs::Vector{AffineConstraint{LT}}
    gt_general_constrs::Vector{AffineConstraint{GT}}
    depth::Int
    dual_bound::Float64
end

function Node()
    return Node(
        BoundUpdate{LT}[],
        BoundUpdate{GT}[],
        AffineConstraint{LT}[],
        AffineConstraint{GT}[],
        0,
        -Inf,
    )
end

function Node(
    lt_bounds::Vector{BoundUpdate{LT}},
    gt_bounds::Vector{BoundUpdate{GT}},
    depth::Int,
    dual_bound::Float64 = -Inf,
)
    return Node(lt_bounds, gt_bounds, [], [], depth, dual_bound)
end

function Node(
    lt_bounds::Vector{BoundUpdate{LT}},
    gt_bounds::Vector{BoundUpdate{GT}},
    lt_general_constrs::Vector{AffineConstraint{LT}},
    gt_general_constrs::Vector{AffineConstraint{GT}},
    depth::Int,
)
    return Node(
        lt_bounds,
        gt_bounds,
        lt_general_constrs,
        gt_general_constrs,
        depth,
        -Inf,
    )
end

function Base.copy(node::Node)
    return Node(
        copy(node.lt_bounds),
        copy(node.gt_bounds),
        copy(node.lt_general_constrs),
        copy(node.gt_general_constrs),
        node.depth,
        node.dual_bound,
    )
end

mutable struct Tree
    open_nodes::DataStructures.Stack{Node}

    Tree() = new(DataStructures.Stack{Node}())
end

Base.isempty(tree::Tree) = Base.isempty(tree.open_nodes)
Base.length(tree::Tree) = Base.length(tree.open_nodes)
push_node!(tree::Tree, node::Node) = push!(tree.open_nodes, node)
pop_node!(tree::Tree) = pop!(tree.open_nodes)
num_open_nodes(tree::Tree) = length(tree.open_nodes)


text:
export VrPCA,
        reduceToRank,
        getindex,
        SvdApprox,
        isConverged

"""
    # enum SvdApprox

. Lrapprox encodes for the initialization of vrpca with the LowRankApprox package

. Ifsvd  encodes for the initialization of vrpca with incremental svd of Vempala-Kunnan

"""
@enum SvdApprox begin
    lrapprox
    ifsvd
end

"""
# VrPCA

returns an orthonormal matrix W with columns vectors the first
k left singular vectors  spanning the range of X.
**k is supposed to be small compared to d**.

If X = U S Vt, the methode computes in W the matrix U[:, 1:k]
This structure gathers data to do stochastic svd as described in :
***Fast Stochastic Algorithms for SVD and PCA : convergences Properties and Convexity
Ohad Shamir 2015***

FIELDS
------

 - X  the (d,n) matrix we want to reduce.
 - k     the dimension into which we want to reducte the data Matrix
         ***Note :  If the obtained rank is less than asked for value of k is reset!***

 - W     orthogonal (d,k) matrix
 - eta   step size
 - m     epoch length

 - U : Union{ Some{Matrix{Float64}}, Nothing} the matrix of left singular vectors if any
 - S : Union{ Some{Vector{Float64}}, Nothing} the vecor of singular values if any
 - Vt: Union{ Some{Matrix{Float64}}, Nothing} the transposed matrix of right singular vectors if any

"""
mutable struct VrPCA
    X::Matrix{Float64}
    k::Int64
    # orthogonal matrix (d,k)
    W::Matrix{Float64}
    converged::Bool
    relerr::Float64
    #
    U::Union{ Some{Matrix{Float64}}, Nothing}
    S::Union{ Some{Vector{Float64}}, Nothing}
    Vt::Union{ Some{Matrix{Float64}}, Nothing}
end

"""
# function VrPCA(data::Matrix{Float64}, k::Int64)

This function initialize a VrPCA problem with the LowRankApprox package

## Args
. data : the data matrix with data vectors in column
. k : the rank we want as output
"""
function VrPCA(data::Matrix{Float64}, k::Int64)
        @debug  "VrPCA initialized by LowRankApprox and rank target"
        # get an initial W0, we do not need very good precision
        opts = LowRankApprox.LRAOptions(rank=k, rtol = 0.001)
        # get W such that Mdata ~ W*W'*Mdata  in Froebenius norm
        W = prange(data, opts)
        if size(W)[2] < k
            @warn "VrPCA initialization got less than asked rank, got : " size(W)[2]
        end
        @debug "got rank " size(W)[2]
        VrPCA(data, size(W)[2] , W, false, 1., nothing,nothing,nothing)
end

"""
# function VrPCA(data::Matrix{Float64}, epsil::Float64)

This function initialize a VrPCA problem with the LowRankApprox package

## Args
 - data : the data matrix with data vectors in column
 - epsil: the precision we ask at initialisation, the VrPCA iteration will further reduce error
         when iterating with the rank obtained at initialization
"""
function VrPCA(data::Matrix{Float64}, epsil::Float64)
    @debug "VrPCA initialized by Vempala's LowRankApproxMc and epsil target"
    lrmc = LowRankApproxMc(data)
    samplingSize = 100
    reduceToEpsil(lrmc, samplingSize, epsil)
    U = PartialSvdStoch.getindex(lrmc, :U)
    if U === nothing
        throw("initialization of VrPCA with epsil failed")
    end
    W = something(U)
    @debug "initialization got rank " size(W)[2]
    VrPCA(data, size(W)[2], W, false, 1. , nothing, nothing, nothing)
end

"""
# function iterateVrPCA(pcaPb::VrPCA, eta::Float64, m::Int64)

## Args
 - eta : the gradient step. Should be sufficiently small to ensure inversibility of the (k,k) matrix
         W*transpose(W) obtained during iterations on W(d,k)

 - m   : the batch size. Roughly inversely proportional to eta.

    eta = 0.01 and m = 20 is a good compromises. m does not need to be large
        as sampling is weighted by L2-norm of data column.

 Output : The function does not return any output.
    It fills the fields U, S, Vt of the structure VrPCA.
    The fields can be retrieved by the function getIndex.
"""
function iterateVrPCA(pcaPb::VrPCA, eta::Float64, m::Int64)
    @debug "entering VrPCA reduce eta m block version" eta m
    small = 1.E-5
    d,n = size(pcaPb.X)
    k = pcaPb.k
    # allocates once and for all 3 temproray matrices
    Us = zeros(Float64, d, k)  # is (d,k)
    B_t = zeros(Float64, k, k)   # Bt_1 is Wt*W so it is a (k,k) matrix
    W_tmp = zeros(Float64, d, k)
    pcaPb.converged = false
    s = 1
    W_s = pcaPb.W
    W_t = pcaPb.W
    deltaWIter = Vector{Float64}()
    # get L2 proba per column
    probaC = Vector{Float64}(undef,n)
    normC = Vector{Float64}(undef,n)
    normX = 0
    Threads.@threads for i = 1:n
        normC[i] = norm(pcaPb.X[:,i]) * norm(pcaPb.X[:,i])
    end
    for i = 1:n
        probaC[i] = i > 1 ? probaC[i-1] + normC[i] : normC[i]
        normX += normC[i]
    end
    if probaC[end] <= 0
        @debug "\n non positive proba for last slot=  :  " probaC[end]
        error("cannot sample column")
    end
    map!(x -> x/probaC[end], probaC, probaC)
    #
    @debug " W_s dim" size(W_s)
    more = true
    while more
        fill!(Us,0.)
        # split in block the update of Us to get speed from BLAS.gemm
        blockSize = 1000
        y = zeros(k, blockSize)
        nbblocs = floor(Int64, n / blockSize)
        for numbloc = 1:nbblocs
            blocfirst = 1 + (numbloc -1) * blockSize
            bloclast = min(blockSize * numbloc, n)
            # get  W_s' * pcaPb.X[:,i] in y
            BLAS.gemm!('T', 'N', 1., W_s , pcaPb.X[:,blocfirst:bloclast], 0., y)
            # Us = Us + pcaPb.X[:,i] * y'  (= transpose(pcaPb.X[:,i]) * W_s)
            for i in 1:blockSize
                j = (numbloc - 1) * blockSize + i
                BLAS.ger!(1.,pcaPb.X[:,j], y[:,i], Us)
            end
            # Us = Us + pcaPb.X[:,i] * (transpose(pcaPb.X[:,i]) * W_s)
        end
        # do not forget the residual part of block splitting!!
        if n % blockSize > 0
            y = zeros(k)
            for i in 1 + nbblocs * blockSize:n
                BLAS.gemm!('T', 'N', 1., W_s , pcaPb.X[:,i], 0., y)
                BLAS.ger!(1.,pcaPb.X[:,i], y, Us)
            end
        end
        Us = Us / n
        # stochastic update pass
        # we need to store Wt W_{s-1} (Ws1 in the code) W_{t-1} (Wt in the code)
        for t = 1:m
            # compute in place to avoid allocator transpose(W_{t_1}) * W_{s-1}, B_t is (k,k)
            BLAS.gemm!('T', 'N', 1., W_t, W_s, 0., B_t)
            # compute svd of transpose(W_{t_1}) * W_{s-1}    i.e do a svd of a (k,k) matrix
            U,S,Vt = LAPACK.gesvd!('A', 'A', B_t)
            # compute B_t as V*Ut
            BLAS.gemm!('T', 'T', 1., Vt, U, 0., B_t)
            #
            # biased sampling instead of it = rand(1:n)
            xsi = rand(Float64)
            it = searchsortedfirst(probaC, xsi)
            s_weight = 1. / (n * probaC[it])
            # update Wt, vaux1 is  (1,k)
            vaux1 = transpose(pcaPb.X[:,it]) * W_t - (transpose(pcaPb.X[:,it]) * W_s) * B_t
            # vaux2 is (d,k)
            vaux2 = s_weight * pcaPb.X[:,it] * vaux1 + BLAS.gemm('N', 'N', Us, B_t)
            # @debug "norm vaux2" norm(vaux2)
            # W_s1 * B_t1 is (d,k), so W_t1 and W_s1. W_tmp _s (d,k)
            W_tmp = W_t + eta * vaux2
            #
            # final update of W_t, compute sqrt(transpose(W_t1)*W_{t}),  (k,k)
            # in fact we want to orthonormalize U and Vt are (k,k)
            U,S,Vt = LAPACK.gesvd!('A', 'A', transpose(W_tmp) * W_tmp)
            # @debug "\n size U , Vt" size(U), size(Vt)
            # check for null values of S
            nbsmall = count( x -> x < small, S)
            if nbsmall > 0
                @debug "too small values " S
                @warn "could not do W update, ||W|| " norm(W_tmp, 2)
            end
            map!(x -> 1/sqrt(x) , S, S)
            # compute the inverse, BLAS.gemm does not seem to spped up things here
            Wnorminv = transpose(Vt) * diagm(0 => S) * transpose(U)
            # orthonormalize W. compute:  W_t / sqrt(transpose(W_t)*W_t)
            # @debug "\n size W_tmp Wnorminv" size(W_tmp) size(Wnorminv)
            W_t = BLAS.gemm('N', 'N', W_tmp, Wnorminv)
            # Cshould check W_t orthogonality : OK!
            # @debug "\n W orthogonality check " norm(transpose(W_t) * W_t - Matrix{Float64}(I,k,k))
        end
        #
        # convergence detection
        normWt = norm(W_t) * norm(W_t)
        deltaW = norm(W_t - W_s) * norm(W_t - W_s)
        lastDeltaW = 0
        if length(deltaWIter) > 0
            lastDeltaW = deltaWIter[end]
        end
        push!(deltaWIter, deltaW)
        @debug " iter ,  ||W_t - W_s||_2 , ||W_t||_2" s deltaW normWt
        W_s = W_t
        #
        # do a convergence test
        if normWt > 0 && deltaW/normWt < 1.E-3
            pcaPb.W = W_t
            pcaPb.relerr = deltaW/normWt
            more = false
            pcaPb.converged = true
        else
            if s >= 3
                if deltaW >= 0.8 * lastDeltaW || deltaW <= 0.05 * deltaWIter[1]
                    @warn "\n reached stationarity  s deltaW  ||W|| ,  exiting :  " s  deltaW  normWt
                    pcaPb.relerr =  deltaW/normWt
                    pcaPb.W = W_t
                    pcaPb.converged = true
                    more  = false
                elseif s > 5
                    # we stop without setting pcaPb.converged...
                    @warn "\n stopping iteration" s  deltaW  normWt
                    pcaPb.relerr =  deltaW/normWt
                    pcaPb.W = W_t
                    more = false
                end
            end
            s = s+1
        end
    end # end of while
    # now we have a converged matrix W (d,k) consisting in k left singular vectors
    # compute singular values and right singular vectors
    #                 = (k,d) * (d,n)
    #      Vt = transpose(W_t) * pcaPb.X
    Vt = Matrix{Float64}(undef, k, n)
    BLAS.gemm!('T', 'N', 1., W_t, pcaPb.X, 0., Vt)
    S = map(i -> norm(Vt[i,:]), (1:size(Vt)[1]))
    for i in 1:length(S)
        Vt[i,:] /= S[i]
    end
    pcaPb.U = Some(W_t)
    pcaPb.S = Some(S)
    pcaPb.Vt = Some(Vt)
end

"""
# function getindex(VrPCA, Symbol)

returns U,S,Vt,V or rank according to Symbol :S , :U, :Vt , V or :k

"""
function getindex(pcaPb::VrPCA, d::Symbol)
    @match d begin
    :S    =>  pcaPb.S
    :U    =>  pcaPb.U
    :V    =>  pcaPb.Vt'
    :Vt   =>  pcaPb.Vt
    :k    =>  length(pcaPb.S)
    _     =>  throw(KeyError(d))
    end
end

"""
# function getResidualError(vrpca::VrPCA)

if X is the data matrix and U computation succeeded
    return the 2-uple (error, relerr) with:
```math
        error = norm(X - U * transpose(U) * X)
```
and
```math
        relerror = norm(X - U * transpose(U) * X)/norm(X)
```

else throws an error
"""
function getResidualError(vrpca::VrPCA)
    U = something(vrpca.U)
    resErr = norm(vrpca.X - U* BLAS.gemm('T', 'N', U, vrpca.X))
    return resErr,resErr/norm(vrpca.X)
end

"""
# function isConverged(pcaPb::VrPCA)

return true if stationarity criteria were satisfied during iterations.
"""
function isConverged(pcaPb::VrPCA)
    return pcaPb.converged
end


text:
"""
    Orihime.@plot(vega_lite_spec...)

A `@vlplot` wrapper for extra functionalities:

* `data` is trimmed down using `Orihime.cleanup`
* `extra_facets` (a list of column names) can be passed as the top-level named
   argument to group by data by these columns. The plots for these groups are
   concatenated vertically.
"""
macro plot(vega_lite_spec...)
    vlexprs = []
    options = []
    for ex in vega_lite_spec
        if Meta.isexpr(ex, :(=), 2) && ex.args[1] in (:extra_facets, :data, :title)
            push!(options, Expr(:kw, ex.args...))
        else
            push!(vlexprs, insert_cleanup_expr(ex))
        end
    end
    @gensym orihime_plot_body
    quote
        function $orihime_plot_body(data, title)
            return $VegaLite.@vlplot(data = data, title = title, $(vlexprs...))
        end
        $plot_impl($orihime_plot_body; $(options...))
    end |> esc
end

function plot_impl(
    @nospecialize(orihime_plot_body);
    data = nothing,
    extra_facets = nothing,
    title = "",
)
    if extra_facets === nothing || isempty(extra_facets)
        return orihime_plot_body(data, title)
    end
    if data === nothing
        error("`extra_facets` requires `data` as the top-level argument")
    end
    df = Orihime.cleanup(data)
    plots = mapfoldl(
        push!,
        pairs(DataFrames.groupby(df, extra_facets; sort = true));
        init = VLSpec[],
    ) do (groupkey, subdata)
        title_values = join(
            ("$k=$(Orihime.readablestring(v))" for (k, v) in pairs(groupkey)),
            ", ",
        )
        if title == ""
            subtitle = title_values
        else
            subtitle = "$title $title_values"
        end
        return orihime_plot_body(subdata, subtitle)
    end
    # `reduce(vcat, plots)` does not work with vlplot?
    return vcat(plots...)
end

insert_cleanup_expr(@nospecialize(x)) = x
function insert_cleanup_expr(ex::Expr)
    if Meta.isexpr(ex, :meta) && Meta.isexpr(ex, :call) && Meta.isexpr(ex, :ref)
        return ex
    elseif Meta.isexpr(ex, :(=), 2) && ex.args[1] == :data
        # Convert `data = $rhs` to `data = cleanup($rhs)` when it's not at the
        # top level (i.e., not handled by `plot_impl`)
        return Expr(ex.head, ex.args[1], :($Orihime.cleanup($(ex.args[2]))))
    else
        return Expr(ex.head, Iterators.map(insert_cleanup_expr, ex.args)...)
    end
end


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module FirstOrderLp

# Built-in Julia packages.
import LinearAlgebra
# We use the Logging package to control the output level of QPSReader, not for
# internal logging.
# TODO: Consider switching internal logging to Logging.
import Logging
import Printf
import SparseArrays
import Random

# Third-party Julia packages.
import GZip
import QPSReader
import Statistics
import StatsBase
import StructTypes

const Diagonal = LinearAlgebra.Diagonal
const diag = LinearAlgebra.diag
const dot = LinearAlgebra.dot
const cholesky = LinearAlgebra.cholesky
const ldlt = LinearAlgebra.ldlt
const lu = LinearAlgebra.lu
const norm = LinearAlgebra.norm
const mean = Statistics.mean
const median = Statistics.median
const nzrange = SparseArrays.nzrange
const nnz = SparseArrays.nnz
const nonzeros = SparseArrays.nonzeros
const rowvals = SparseArrays.rowvals
const sparse = SparseArrays.sparse
const SparseMatrixCSC = SparseArrays.SparseMatrixCSC
const spdiagm = SparseArrays.spdiagm
const spzeros = SparseArrays.spzeros
const quantile = Statistics.quantile
const sample = StatsBase.sample

include("quadratic_programming.jl")
include("trust_region_utils.jl")
include("solve_log.jl")
include("quadratic_programming_io.jl")
include("preprocess.jl")
include("saddle_point.jl")
include("termination.jl")
include("iteration_stats_utils.jl")
include("mirror_prox.jl")
include("primal_dual_hybrid_gradient.jl")

end  # module FirstOrderLp


text: ## correlation to calculate flux


text:
module SDFReader

export filereader, Field, particle_variable, build_vector, electric_field

using PyCall
using Unitful: Unitful, Units, @u_str
using GeometryTypes
# using RecursiveArrayTools
# needs broadcast for VectorOfArray

const sdf = PyNULL()
const sh = PyNULL()

function __init__()
    copy!(sdf, pyimport("sdf"))
    copy!(sh, pyimport("sdf_helper"))
end

include("utils.jl")
include("fields.jl")
include("particles.jl")

function build_vector(v::Vararg{T,N}) where {N,T}
    Array(Point{N}.(v...))
end

"""
    filereader(path; convert=false)

Return a file reader for the EPOCH dump file (.sdf) at the given path.
The values can be optionally converted to 32 bits with the `convert` argument.
"""
function filereader(path; convert=false)
    sdf.read(path, convert=convert)
end

function electric_field(fr)
    py_ex = fr.Electric_Field_Ex
    py_ey = fr.Electric_Field_Ey
    py_ez = fr.Electric_Field_Ez

    ex = convert_it(py_ex.data)
    ey = convert_it(py_ey.data)
    ez = convert_it(py_ez.data)

    units = get_units(py_ex.units)
    grid_mid = convert_it.(py_ex.grid_mid.data)
    Field(ex, ey, ez, grid_mid, units)
end

end # module


text:
using TinySegmenter

f = open("timemachineu8j.txt")
TinySegmenter.tokenize("プレコンパイル")
s = join(readlines(f))

@time for i in 1:100
  TinySegmenter.tokenize(s)
end


text:
"""
    fm = parse_source(filename::AbstractString, mod::Module)

Parse the source `filename`, returning a [`FileModules`](@ref) `fm`.
`mod` is the "parent" module for the file (i.e., the one that `include`d the file);
if `filename` defines more module(s) then these will all have separate entries in `fm`.

If parsing `filename` fails, `nothing` is returned.
"""
parse_source(filename::AbstractString, mod::Module) =
    parse_source!(FileModules(mod), filename, mod)

"""
    parse_source!(fm::FileModules, filename, mod::Module)

Top-level parsing of `filename` as included into module
`mod`. Successfully-parsed expressions will be added to `fm`. Returns
`fm` if parsing finished successfully, otherwise `nothing` is returned.

See also [`parse_source`](@ref).
"""
function parse_source!(fm::FileModules, filename::AbstractString, mod::Module)
    if !isfile(filename)
        @warn "$filename is not a file, omitting from revision tracking"
        return nothing
    end
    parse_source!(fm, read(filename, String), Symbol(filename), 1, mod)
end

"""
    success = parse_source!(fm::FileModules, src::AbstractString, file::Symbol, pos::Integer, mod::Module)

Parse a string `src` obtained by reading `file` as a single
string. `pos` is the 1-based byte offset from which to begin parsing `src`.

See also [`parse_source`](@ref).
"""
function parse_source!(fm::FileModules, src::AbstractString, file::Symbol, pos::Integer, mod::Module)
    local ex, oldpos
    # Since `parse` doesn't keep track of line numbers (it works
    # expression-by-expression), to ensure good backtraces we have to
    # keep track of them here. For each expression we parse, we count
    # the number of linefeed characters that occurred between the
    # beginning and end of the portion of the string consumed to parse
    # the expression.
    line_offset = 0
    while pos < lastindex(src)
        try
            oldpos = pos
            ex, pos = Meta.parse(src, pos; greedy=true)
        catch err
            ex, posfail = Meta.parse(src, pos; greedy=true, raise=false)
            warnline = line_offset + count(c->c=='\n', SubString(src, oldpos, posfail)) + 1
            with_logger(SimpleLogger(stderr)) do
                @error "omitting file $file due to parsing error near line $warnline"
            end
            showerror(stderr, err)
            println(stderr)
            return nothing
        end
        if isa(ex, Expr)
            ex = ex::Expr
            fix_line_statements!(ex, file, line_offset)  # fixes the backtraces
            parse_expr!(fm, ex, file, mod)
        end
        # Update the number of lines
        line_offset += count(c->c=='\n', SubString(src, oldpos, pos-1))
    end
    fm
end

"""
    success = parse_source!(fm::FileModules, ex::Expr, file, mod::Module)

For a `file` that defines a sub-module, parse the body `ex` of the
sub-module.  `mod` will be the module into which this sub-module is
evaluated (i.e., included). Successfully-parsed expressions will be
added to `fm`. Returns `true` if parsing finished successfully.

See also [`parse_source`](@ref).
"""
function parse_source!(fm::FileModules, ex::Expr, file::Symbol, mod::Module)
    @assert ex.head == :block
    for a in ex.args
        if isa(a, Expr)
            parse_expr!(fm, a::Expr, file, mod)
        end
    end
    fm
end

"""
    parse_expr!(fm::FileModules, ex::Expr, file::Symbol, mod::Module)

Recursively parse the expressions in `ex`, iterating over blocks and
sub-module definitions. Successfully parsed
expressions are added to `fm` with key `mod`, and any sub-modules will
be stored in `fm` using appropriate new keys. This accomplishes two main
tasks:

* add parsed expressions to the source-code cache (so that later we can detect changes)
* determine the module into which each parsed expression is `eval`uated into
"""
function parse_expr!(fm::FileModules, ex::Expr, file::Symbol, mod::Module)
    if ex.head == :block
        for a in ex.args
            a isa Expr || continue
            parse_expr!(fm, a, file, mod)
        end
        return fm
    end
    macroreplace!(ex, String(file))
    if ex.head == :line
        return fm
    elseif ex.head == :module
        parse_module!(fm, ex, file, mod)
    elseif isdocexpr(ex) && isa(ex.args[nargs_docexpr], Expr)
        exex = ex.args[nargs_docexpr]
        while is_trivial_block_wrapper(exex)
            exex = exex.args[end]
        end
        if isa(exex, Expr) && exex.head == :module
            # Module with a docstring (issue #8)
            # Split into two expressions, a module definition followed by
            # `"docstring" newmodule`
            newmod = parse_module!(fm, ex.args[nargs_docexpr], file, mod)
            ex.args[nargs_docexpr] = Symbol(newmod)
            fm[mod].defmap[convert(RelocatableExpr, ex)] = nothing
        elseif isa(exex, Expr) && (exex.head == :function || exex.head == :(=) || exex.head == :macrocall)
            # Generate the doc expressions
            local mex
            try
                mex = macroexpand(mod, ex)
            catch
                return nothing
            end
            while is_trivial_block_wrapper(mex)
                mex = mex.args[end]
            end
            if mex.head != :block
                @warn "expected block expression got $mex from $ex"
                return nothing
            end
            if mex.args[1] isa Expr && (mex.args[1].head == :(=) || mex.args[1].head == :function)
                # Separate the function definition from the expression
                # defining the docstring
                parse_expr!(fm, exex, file, mod)
                mex.args[1] = nothing
            end
            for arg in mex.args
                arg isa ExLike || continue
                is_linenumber(arg) && continue
                fm[mod].defmap[convert(RelocatableExpr, arg)] = nothing
            end
        else
            fm[mod].defmap[convert(RelocatableExpr, ex)] = nothing
        end
    elseif ex.head == :call && ex.args[1] == :include
        # skip include statements
    else
        # Any expression that *doesn't* define line numbers, new
        # modules, or include new files must be "real code."
        # Handle macros
        exorig = ex0 = ex
        if ex isa Expr && ex.head == :macrocall
            if ex.args[1] ∉ (Symbol("@warn"), Symbol("@info"), Symbol("@debug"), Symbol("@error"), Symbol("@logmsg"))  # issue #208
                # To get the signature, we have to expand any unrecognized macro because
                # the macro may change the signature
                try
                    ex0, ex = macexpand(mod, ex)
                catch
                end
            end
        end
        ex isa Expr || return fm
        ex.head == :tuple && isempty(ex.args) && return fm
        if ex.head == :block
            return parse_expr!(fm, ex, file, mod)
        end
        # Add any method definitions to the cache
        sig = ex.head == :macrocall ? nothing : get_signature(convert(RelocatableExpr, ex))
        # However, we have to store the original unexpanded expression if
        # `revise(mod)` can be expected to work (issue #174).
        rex = convert(RelocatableExpr, exorig)
        if isa(sig, ExLike)
            fm[mod].defmap[rex] = sig  # we can't safely `eval` the types because they may not yet exist
        else
            fm[mod].defmap[rex] = nothing
        end
    end
    fm
end

function macexpand(mod::Module, ex::Expr)
    ex0 = ex
    if ex.args[1] ∈ poppable_macro
        ex = ex.args[end]
        if ex isa Expr && ex.head == :macrocall
            ex0.args[end], ex = macexpand(mod, ex)
        end
    else
        ex0 = ex = macroexpand(mod, ex)
    end
    return ex0, ex
end

const nargs_docexpr = 4
isdocexpr(ex) = ex.head == :macrocall && ex.args[1] == GlobalRef(Core, Symbol("@doc")) &&
           length(ex.args) >= nargs_docexpr

"""
    newmod = parse_module!(fm::FileModules, ex::Expr, file, mod::Module)

Parse an expression `ex` that defines a new module `newmod`. This
module is "parented" by `mod`. Source-code expressions are added to
`fm` under the appropriate module name.
"""
function parse_module!(fm::FileModules, ex::Expr, file::Symbol, mod::Module)
    newname = _module_name(ex)
    if isdefined(mod, newname)
        newmod = getfield(mod, newname)
    else
        id = Base.identify_package(mod, String(newname))
        if id === nothing
            newmod = eval_module_expr(mod, ex, newname)
        else
            newmod = Base.root_module(id)
            if !isa(newmod, Module)
                newmod = eval_module_expr(mod, ex, newname)
            end
        end
    end
    fm[newmod] = FMMaps()
    parse_source!(fm, ex.args[3], file, newmod)  # recurse into the body of the module
    newmod
end

function eval_module_expr(mod, ex, newname)
    with_logger(_debug_logger) do
        @debug "parse_module" _group="Parsing" activemodule=fullname(mod) newname
    end
    try
        Core.eval(mod, ex) # support creating new submodules
    catch
        @warn "Error evaluating expression in $mod:\n$ex"
        rethrow()
    end
    return getfield(mod, newname)
end

_module_name(ex::Expr) = ex.args[2]


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

module Gaston

export closefigure, closeall, figure,
       plot, plot!, histogram, imagesc, surf,
       printfigure, set

import Base.show

# before doing anything else, verify gnuplot is present on this system
try
  success(`gnuplot --version`)
catch
  error("Gaston cannot be loaded: gnuplot is not available on this system.")
end

# load files
include("gaston_types.jl")
include("gaston_aux.jl")
include("gaston_llplot.jl")
include("gaston_hilvl.jl")
include("gaston_config.jl")

# set up global variables
# global variable that stores gnuplot's state
gnuplot_state = GnuplotState(false,nothing,[],"","","",false,
                             Figure[],displayable("image/png"))

# when gnuplot_state goes out of scope, exit gnuplot
if VERSION >= v"0.7-"
  finalizer(gnuplot_exit, gnuplot_state)
else
  finalizer(gnuplot_state, gnuplot_exit)
end

# global variable that stores Gaston's configuration
gaston_config = GastonConfig()

# get started
gnuplot_init()

end


text:
using Plots

function runsimulation_1D(input::JuSwalbe.Inputconstants)
    lx = input.lx
    ly = input.ly
    maxruntime = input.maxruntime

    h0 = 2 .+ sin.(4*π.*[1:lx;]./lx)
    if ly == 0
        println("Starting a one dimensional simulation")
        # Define an initial fluid state
        mom = JuSwalbe.Macroquant(height = h0,
                                  velocity = zeros(lx),
                                  pressure = zeros(lx),
                                  energy = zeros(lx))

        # Define an struct containing all forces
        forces = JuSwalbe.Forces(slip = zeros(lx),
                                 h∇p = zeros(lx),
                                 bathymetry = zeros(lx),
                                 thermal = zeros(lx))

        # Set up an equilibrium distribtution from the chosen fluid state
        equilibrium = calc_equilibrium_distribution(mom)
        temp_dist = deepcopy(equilibrium)

        # Next we iterate until maxruntime is reached
        for time in 1:maxruntime
            # Compute the updated moments
            mom = calculatemoments(temp_dist)

            # Compute all forces
            computeslip(mom, forces, input)
            pressure(mom, γ=input.γ, θ=ones(Float64,1)*1/9)
            ∇p(mom, forces)
            computethermalcapillarywaves(mom, forces, input)

            # Add boundary conditions here

            # Performe the collision
            out_dist = collisionBGK(mom, forces, temp_dist, input)

            # Add obstacles and bounce back here

            # Stream the distribtutions
            temp_dist = streamdistperiodic(out_dist)

            if time%100 == 0
                println("Heigth: ", mom.height[5], " Velocity: ", mom.velocity[5], " Pressure: ", mom.pressure[5], " Time step: ", time)
                println("Forces -> Pressure-grad: ", forces.h∇p[5], " slippage: ", forces.slip[5])
            end
        end
        plot(h0)
        plot!(mom.height)
        gui()
        return mom
    else
        println("Please use `runsimulation_2D`")
    end
end

function runsimulation_2D(input::JuSwalbe.Inputconstants)
    lx = input.lx
    ly = input.ly
    maxruntime = input.maxruntime

    h0 = [2 + (sin(4*π*i/lx) * sin(4*π*j/ly)) for i in 1:lx, j in 1:ly]
    if ly == 0
        println("Please use `runsimulation_1D`")
    else
        println("Starting Simulation in two dimensions")
        # Define an initial fluid state
        mom = JuSwalbe.Macroquant(height = h0,
                                  velocity = zeroTwovector(lx, ly),
                                  pressure = zeros(lx, ly),
                                  energy = zeros(lx, ly))

        # Define an struct containing all forces
        forces = JuSwalbe.Forces(slip = zeroTwovector(lx, ly),
                                 h∇p = zeroTwovector(lx, ly),
                                 bathymetry = zeroTwovector(lx, ly),
                                 thermal = zeroTwovector(lx, ly))

        # Set up an equilibrium distribtution from the chosen fluid state
        equilibrium = calc_equilibrium_distribution(mom)
        temp_dist = deepcopy(equilibrium)

        # Next we iterate until maxruntime is reached
        for time in 1:maxruntime
            # Compute the updated moments
            mom = @spawn calculatemoments(temp_dist)

            # Compute all forces
            computeslip(fetch(mom), forces, input)
            pressure(fetch(mom), γ=input.γ, θ=ones(Float64,(1,1))*1/9)
            ∇p(fetch(mom), forces)
            computethermalcapillarywaves(fetch(mom), forces, input)

            # Add boundary conditions here

            # Performe the collision
            out_dist = @spawn collisionBGK(fetch(mom), forces, temp_dist, input)

            # Add obstacles and bounce back here

            # Stream the distribtutions
            temp_dist = streamdistperiodic(fetch(out_dist))
            if time % 1000 == 0
                println("Time step: $time")
            end

        end
        s1 = heatmap(h0)
        s2 = heatmap(fetch(mom).height)
        plot(s1, s2)
        gui()
        return mom
    end

end

function runsimulation_2DGPU(input::JuSwalbe.Inputconstants)
    lx = input.lx
    ly = input.ly
    maxruntime = input.maxruntime

    h0 = [2 + (sin(4*π*i/lx) * sin(4*π*j/ly)) for i in 1:lx, j in 1:ly]
    h0gpu = Float32.(h0)
    hnew = CUDA.CuArray(h0gpu)
    if ly == 0
        println("Please use `runsimulation_1D`")
    else
        println("Starting Simulation in two dimensions")
        # Define an initial fluid state
        mom = JuSwalbe.Macroquant(height = hnew,
                                  velocity = zeroTwovectorCU(lx, ly),
                                  pressure = CUDA.zeros(lx, ly),
                                  energy = CUDA.zeros(lx, ly))

        # Define an struct containing all forces
        forces = JuSwalbe.Forces(slip = zeroTwovectorCU(lx, ly),
                                 h∇p = zeroTwovectorCU(lx, ly),
                                 bathymetry = zeroTwovectorCU(lx, ly),
                                 thermal = zeroTwovectorCU(lx, ly))

        # Set up an equilibrium distribtution from the chosen fluid state
        equilibrium = calc_equilibrium_distribution(mom)
        temp_dist = deepcopy(equilibrium)

        # Next we iterate until maxruntime is reached
        for time in 1:maxruntime
            # Compute the updated moments
            mom = calculatemoments(temp_dist)

            # Compute all forces
            computeslip(mom, forces, input)
            pressure(mom, γ=input.γ, θ=ones(Float32,(1,1))*1/9)
            ∇p(mom, forces)
            computethermalcapillarywaves(mom, forces, input)

            # Add boundary conditions here

            # Performe the collision
            out_dist = collisionBGK(mom, forces, temp_dist, input)

            # Add obstacles and bounce back here

            # Stream the distribtutions
            temp_dist = streamdistperiodic(out_dist)
            if time % 1000 == 0
                println("Time step: $time")
            end

        end
        s1 = heatmap(h0)
        s2 = heatmap(mom.height)
        plot(s1, s2)
        gui()
        return mom
    end

end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
##
#
# Draw various histograms
#

using ChristensenSet

DATA_FILE = "ex13_Float64.dat"

# set up images
images = (
  (RootsImage(-2.0-2.0im, 2.0+2.0im, 2000, 2000), "img13_1.png", "img13_1.tex"),
  (RootsImage(-1.0-0.0im, 0.0+1.0im, 2000, 2000), "img13_2.png", "img13_2.tex"),
  (RootsImage(-0.78-1.8im, 0.22-0.8im, 2000, 2000), "img13_3.png", "img13_3.tex"),
  (RootsImage(-0.78-1.6im, -0.28-1.1im, 2000, 2000), "img13_4.png", "img13_4.tex"),
)

# fill images with roots
fill_images!(collect(map(first, images)), DATA_FILE)

# save images to disk
for img in images
  save_image(img[1], img[2], latex_filename=img[3])
  run(`pdflatex $(img[3])`)
end

@info "Done."


text:
__precompile__(true)

"""
Module `Approximations.jl` -- polygonal approximation of convex sets through
support vectors.
"""
module Approximations

using LazySets, LazySets.Arrays, Requires, LinearAlgebra, SparseArrays

using LazySets: _isapprox, _leq, _geq, _rtol, _normal_Vector, isapproxzero,
                default_lp_solver, _isbounded_stiemke, require, dim, linprog,
                is_lp_optimal

import LazySets: project

using ..Assertions: @assert, activate_assertions
# activate assertions by default
activate_assertions(Approximations)

export approximate,
       ballinf_approximation,
       box_approximation, interval_hull,
       decompose,
       overapproximate,
       underapproximate,
       box_approximation_symmetric, symmetric_interval_hull,
       BoxDirections,
       DiagDirections,
       BoxDiagDirections,
       OctDirections,
       PolarDirections,
       SphericalDirections,
       CustomDirections,
       isbounding

include("box_approximation.jl")
include("iterative_refinement.jl")
include("symmetric_interval_hull.jl")
include("ballinf_approximation.jl")
include("template_directions.jl")
include("overapproximate.jl")
include("underapproximate.jl")
include("approximate.jl")
include("decompositions.jl")
include("hausdorff_distance.jl")
include("init.jl")

end # module


text:
# weights(lattice::D2Q9) = [4/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9]'
# abscissaes(lattice::D2Q9) = [0  0; -1  1; -1  0; -1 -1; 0 -1; 1 -1; 1  0; 1  1; 0  1]'
# opposite(lattice::D2Q9) = [1 2 3 4 5 6 7 8 9]
# opposite(lattice::D2Q9, idx) = idx + (9 - 1) / 2
# dimension(lattice::D2Q9) = 2
# order(lattice::D2Q9) = [1, 7, 9, 3, 5, 8, 2, 4, 6]

# D2Q9
# const original_order = [1, 7, 9, 3, 5, 8, 2, 4, 6]

struct D2Q9{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end
D2Q9() = D2Q9(
    [
        0 -1 -1 -1 0 1 1 1 0
        0 1 0 -1 -1 -1 0 1 1
    ],
    [4 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9],
    3.0,
)

order(q::D2Q9) = 5
function opposite(q::D2Q9, idx::Int64)
    if idx == 1
        return 1
    end
    if idx <= 5
        return idx + 4
    end
    return idx - 4
end

Base.show(io::IO, q::D2Q9) = show(io, "D2Q9")
Base.string(q::D2Q9) = "D2Q9"


text:
export WishartGaussian, GaussianDiagonal, DirichletMultinomial, GammaNormal, NormalNormal, BetaBernoulli, ConjugatePostDistribution

abstract type ConjugatePostDistribution end

abstract type UnivariateConjugatePostDistribution <: ConjugatePostDistribution end
abstract type DiscreteUnivariateConjugatePostDistribution <: UnivariateConjugatePostDistribution end
abstract type ContinuousUnivariateConjugatePostDistribution <: UnivariateConjugatePostDistribution end

abstract type MultivariateConjugatePostDistribution <: ConjugatePostDistribution end
abstract type DiscreteMultivariateConjugatePostDistribution <: MultivariateConjugatePostDistribution end
abstract type ContinuousMultivariateConjugatePostDistribution <: MultivariateConjugatePostDistribution end

# Gaussian with Normal Inverse Wishart Prior
mutable struct WishartGaussian <: ContinuousMultivariateConjugatePostDistribution

  D::Int

  # sufficient statistics
  n::Int
  sums::Vector{Float64}
  ssums::Array{Float64}

  # base model parameters
  μ0::Vector{Float64}
  κ0::Float64
  ν0::Float64
  Σ0::Array{Float64}

end

"""
  WishartGaussian(μ0, κ0, ν0, Σ0)

## Gaussian-inverse-Wishart distribution
A Gaussian-inverse-Wishart distribution is the conjugate prior of a multivariate normal distribution with unknown mean and covariance matrix.

## Parameters
* `μ0, Dx1`: location
* `κ0 > 0`: number of pseudo-observations
* `ν0 > D-1`: degrees of freedom
* `Σ0 > 0, DxD`: scale matrix

## Example
```julia-repl
julia> (N, D) = size(X)
julia> μ0 = mean(X, dims = 1)
julia> d = WishartGaussian(μ0, 1.0, 2*D, cov(x))
```

"""
function WishartGaussian(μ0::Vector{Float64}, κ0::Float64,
                         ν0::Float64, Σ0::Array{Float64})

  d = length(μ0)
  (D1, D2) = size(Σ0)
  @assert D1 == D2
  @assert D1 == d

  WishartGaussian(d, 0, zeros(d), zeros(d, d), μ0, κ0, ν0, Σ0)
end

# Normal with Gamma prior
mutable struct GammaNormal <: ContinuousUnivariateConjugatePostDistribution

  # sufficient statistics
  n::Int
  sums::Float64
  ssums::Float64

  # model parameters
  μ0::Float64
  λ0::Float64
  α0::Float64
  β0::Float64

end

"""
GammaNormal(; μ0 = 0.0, λ0 = 1.0, α0 = 1.0, β0 = 1.0)

## Normal-Gamma distribution
A Normal-Gamma distribution is the conjugate prior of a Normal distribution
with unknown mean and precision.

## Paramters
* `μ0`: location
* `λ0 > 0`: number of pseudo-observations
* `α0 > 0`
* `β0 > 0`

Example:
```julia
d = GammaNormal()
```
"""
function GammaNormal(;μ0 = 0.0, λ0 = 1.0, α0 = 1.0, β0 = 1.0)
  GammaNormal(0, 0.0, 0.0, μ0, λ0, α0, β0)
end

# Normal with Normal prior
mutable struct NormalNormal <: ContinuousUnivariateConjugatePostDistribution

  # sufficient statistics
  n::Int
  sums::Float64
  ssums::Float64

  # model parameters
  μ0::Float64
  σ0::Float64

  function NormalNormal(;μ0 = 0.0, σ0 = 1.0)
    new(0, 0.0, 0.0, μ0, σ0)
  end

end

# Gaussian with Diagonal Covariance
mutable struct GaussianDiagonal{T <: ContinuousUnivariateConjugatePostDistribution} <: ContinuousMultivariateConjugatePostDistribution

  # sufficient statistics
  dists::Vector{T}

  # isn't the default constructor sufficient here?
  #function GaussianDiagonal(dists::Vector{T})
  #    new(dists)
  #end

end

# Multinomial with Dirichlet Prior
mutable struct DirichletMultinomial <: DiscreteMultivariateConjugatePostDistribution

  D::Int

  # sufficient statistics
  n::Int
  counts::SparseMatrixCSC{Int,Int}

  # base model parameters
  α0::Float64

  # cache
  dirty::Bool
  Z2::Float64
  Z3::Array{Float64}

  function DirichletMultinomial(D::Int, α0::Float64)
    new(D, 0, sparsevec(zeros(D)), α0, true, 0.0, Array{Float64}(0))
  end

end

# Categorical with Dirichlet Prior
mutable struct DirichletCategorical <: DiscreteUnivariateConjugatePostDistribution

  # sufficient statistics
  n::Int
  counts::SparseMatrixCSC{Int,Int}

  # base model parameters
  α0::Float64

  # cache
  dirty::Bool
  Z2::Float64
  Z3::Array{Float64}

  function DirichletMultinomial(D::Int, α0::Float64)
    new(D, 0, sparsevec(zeros(D)), α0, true, 0.0, Array{Float64}(0))
  end

end

# Bernoulli with Beta Prior
mutable struct BetaBernoulli <: DiscreteUnivariateConjugatePostDistribution

  # sufficient statistics
  successes::Int
  n::Int

  # beta distribution parameters
  α0::Float64
  β0::Float64

  function BetaBernoulli(;α0 = 1.0, β0 = 1.0)
    new(0, 0, α0, β0)
  end

end


text:
#Implementation of UPGMA in julia to output a Newick format string. Either with or without distances. upgma() is the main function

#modules required
using DataFrames
using CSV

#convert dataframe to "lower left" must be run before using with upgma
#replace uppper values with missing
function add_missing(df::DataFrame)

    #length of dataframe
    n = nrow(df)

    #allows adding of NA (missing) values
    allowmissing!(df)

    #sets upper left and diaganol to missing
    for i in 1:n
        for j in 1:n
            if j >= i
                df[i,j] = missing
            end
        end
    end
    return df
end

#function to locate minimum cell in the table
function minimum_cell(df::DataFrame)
    #set starting value to infinity
    low_value = float(Inf)

    #set starting x,y values, non valid ints
    x = -1
    y = -1

    #get length of dataframe
    n = nrow(df)

    #evaluate each cell to see if it is lower
    for i in 1:n
        for j in 1:n
            if isless(df[i,j], low_value)
                low_value = df[i,j]
                x,y = i,j
            end
        end
    end

    #return the coordinates of the lowest value
    return x, y
end

#function to assemble the string to make output
function combine_strings(strings::Array{String,1}, a::Int64, b::Int64)

    #ensure that the order is consistant, always join on the lower one
    if b < a
        a, b = b, a
    end

    #join the values at the position of the first value
    strings[a] = string("(", strings[a], ",", strings[b], ")")

    #remove the second (redundant) value from strings
    setdiff!(strings, [strings[b]])

    return strings
end

#function to assemble the string to make output WITH DISTANCES
function combine_strings_dist(strings::Array{String,1}, a::Int64, b::Int64, dist_df::DataFrame)

    #ensure that the order is consistant, always join on the lower one
    if b < a
        a, b = b, a
    end

    #retrieve the distance value for the appropriate string from the distance dataframe
    tdist_a = dist_df[dist_df.string .== strings[a],:][2][1]
    tdist_b = dist_df[dist_df.string .== strings[b],:][2][1]

    #get the length of the data frame
    dist_df_len = nrow(dist_df)

    #Check if the string a is already part of a node. If so, save the longest distance that it is associated with
    #subtract two, because we don't want to include itself, or its pair
    for i in 1:(dist_df_len-2)
        if occursin(string("(", dist_df.string[i], ":"), strings[a]) || occursin(string(",", dist_df.string[i], ":"), strings[a])
            global diff_a = dist_df.distance[i]
        end
    end

    #same as above, but for b
    for i in 1:(dist_df_len-2)
        if occursin(string("(", dist_df.string[i], ":"), strings[b]) || occursin(string(",", dist_df.string[i], ":"), strings[b])
            global diff_b = dist_df.distance[i]
        end
    end

    #check to see if it is part of a node, if so, subtract previous distance (a)
    if occursin(",", strings[a])
        dist_a = tdist_a - diff_a
    else
        dist_a = tdist_a
    end

    #same as above, but for b
    if occursin(",", strings[b])
        dist_b = tdist_b - diff_b
    else
        dist_b = tdist_b
    end

    #join the values at the position of the first value
    strings[a] = string("(", strings[a], ":", dist_a, ",", strings[b], ":", dist_b, ")")

    #remove the second (redundant) value from strings
    setdiff!(strings, [strings[b]])

    return strings
end

#function to combine df on the column
#comma count is a way to track the weight of each entry
#If it is already an average, it is weighted as two values (or more)
function combine_df(df::DataFrame, a::Int64, b::Int64, comma_count::Array{Int64,1})

    #The number of values that are factored into an entry are one more than the number of commas
    b_den = comma_count[b]+1
    a_den = comma_count[a]+1

    #ensure that order is consistant, always join on the lower one.
    if b < a
        a, b = b, a
        a_den, b_den = b_den, a_den
    end

    #rows in dataframe
    n = nrow(df)

    #for the lower index, reconstruct the entire row
    for i in 1:a
        df[a,i] = ((df[a,i]*a_den + df[b,i]*b_den)/(a_den+b_den))
    end

    #for the lower index, reconstruct the section of the column between the indeces
    #don't include first value, because it was taken care of in the row.
    for i in a+1:b
        df[i,a] = (df[i,a]*a_den + df[b,i]*b_den)/(a_den+b_den)
    end

    #Construct the rest of the row for the lower index, getting values from the ith row.
    for i in b+1:n
        df[i,a] = (df[i,a]*a_den + df[i,b]*b_den)/(a_den+b_den)
    end

    #eliminate redundant rows and columns
    deleterows!(df, b)
    deletecols!(df, b)
    return df
end

#The function for outputting newick file (UPGMA)
#takes a data frame (distance matrix) and an array of strings which are the names of the samples in order
#verbose (if true) prints all dataframes along the way
function upgma(df::DataFrame, strings::Array{String,1} = fill("",1); output::String ="file.newick",
    header::Bool=true, verbose::Bool=false, distances::Bool=true)
    #loop through, eliminated one row/column/string each time
    #print lines are optional, for more information

    #print if verbose is true
    if verbose == true
        println("Initial dataframe:")
        println(df)
        println(strings)
        println(string("number of strings: ", length(strings)))
        println(string("number of rows in DataFrame: ", nrow(df)))
    end

    #get the column names as the strings to use in the combine function if header is true
    if header == true
        strings = fill("", nrow(df))
        colnames = names(df)
        n_cols = length(colnames)
        for i in 1:n_cols
            strings[i] = string(colnames[i])
        end
    end

    #print strings if verbose is true
    if verbose == true
        println(strings)
    end

    #error if the data frame length is not the same as the string length
    if length(strings) != nrow(df)
        error = "Error: matrix size does not match sample size, exiting"
        return error
    end

    #make distance dataframe to store values
    dist_df = DataFrame(string = String[], distance = Float64[])

    #loop through until all the columns/strings have been combined
    while length(strings) > 1

        #find the coordinates of the minimum entry
        x, y = minimum_cell(df)

        #count the number of commas in each string. This (+1) corresponds to how many values it represents
        #make the comma counter array
        comma_count = fill(0, length(strings))
        #println(comma_count)

        #actually count the commas
        for j in 1:length(strings)
            global comma_count[j] = count(i->(i==','), strings[j])
        end
        #println(comma_count)

        #calculate distance
        tdist = df[x,y]/2
        push!(dist_df, (strings[x], tdist))
        push!(dist_df, (strings[y], tdist))

        #combine the strings on the lowest coordinate with or without distances
        if distances == true
            combine_strings_dist(strings, x, y, dist_df)
        else
            combine_strings(strings, x, y)
        end

        #print strings if verbose is true
        if verbose == true
            println(strings)
        end

        #combine the matrix on the lowest coordinate
        combine_df(df, x, y, comma_count)

        #rename df column names based on join
        rename!(df, strings)

        #print the df if true
        if verbose == true
            println(df)
        end

    end

    #write the newick string to a text file specified by argument
    open(output, "w") do f
        write(f, string(strings[1], ";"))
    end

    #return newick as string
    return(string(strings[1], ";"))
end

#function to create distance matrix if in csv form
function create_dm_df(dm_csv::String)
    distance_matrix_dataframe = add_missing(CSV.read(dm_csv))
    return distance_matrix_dataframe
end


text:
using .DataFrames

"""

# convert_a3d

# Convert the output file(s) created by cmdstan to a DataFrame.

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:dataframe})
  # Inital DataFrame
  df = DataFrame(a3d_array[:, :, 1], Symbol.(cnames))

  # Append the other chains
  for j in 2:size(a3d_array, 3)
    df = vcat(df, DataFrame(a3d_array[:, :, j], Symbol.(cnames)))
  end
  df
end

"""

# convert_a3d

# Convert the output file(s) created by cmdstan to a Vector{DataFrame).

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:dataframes})

  dfa = Vector{DataFrame}(undef, size(a3d_array, 3))
  for j in 1:size(a3d_array, 3)
    dfa[j] = DataFrame(a3d_array[:, :, j], Symbol.(cnames))
  end

  dfa
end


text:
# Autogenerated wrapper script for Graphite2_jll for x86_64-w64-mingw32
export libgraphite2

JLLWrappers.@generate_wrapper_header("Graphite2")
JLLWrappers.@declare_library_product(libgraphite2, "libgraphite2.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgraphite2,
        "bin\\libgraphite2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export test_problem, tridiagonal

#
# test the linear program solver on a simple problem
#
using Ipopt, KNITRO
#simple_tests();

# 54
# 59
#

function solve_with_JuMP_sqrt(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, solver=IpoptSolver(max_iter=300))
	model = Model(solver=solver)

	n, k = size(A)

	@defVar(model, x[1:k] >= 0)
	@setNLObjective( model, Min, sum{log(x[i]+1), i=1:k})

	#@addConstraint(model, A*x .== b)
	@defConstrRef constr[1:n]
	for i = 1:n
		constr[i] = @addConstraint(model, sum{A[i,j]*x[j], j=1:k} == b[i])
	end

	status = solve(model)

	if status == :Optimal
		return 1, getValue(x)[:], getDual(constr)[:]
	elseif status == :Infeasible
		return 2
	elseif status == :Unbounded
		return 3
	else
		return 0
	end
end

function test_problem_nlo(name::String)
	file_name = name * ".mat";

  dir = dirname(@__FILE__) * "/Problems/"

  A, b, c = get_netlib_problem(dir, file_name);
  Q = tridiagonal(length(c),0.0, 0.0);
  #Q = tridiagonal(length(c),-1.0,0.0);

	println("Solving ", file_name, " with the homogeneous algorithm")
	println(size(A,2), " variables and ", size(A,1), " constraints")
	println("Non-zeros: ", length(nonzeros(A)))
	settings = class_settings();

	println("=================== Linear system solver is MUMPS ==================")
	settings.linear_system_solver = linear_solver_MUMPS();

	settings.newton_solver = class_homogeneous_newton();

  qp = class_quadratic_program(A, b, c, Q);
  nlo = class_nlo(A, b, nl_log(length(c)))
  homogeneous_algorithm(nlo,settings)
  #trivial_test(A, b, c, Q, 1, file_name, settings, true);

	#settings.newton_solver = class_newton_ip();
	#trivial_test(A, b, c, Q, 1, file_name, settings, true);

	println("=========================================================================")
	println("Calls IPOPT")
	solve_with_JuMP_sqrt(A, b, c, IpoptSolver(max_iter=1000,print_level=3));
end

max_forks_repo_name: ohinder/Everest


text:
@testset "Utils" begin
    mags = [0.0,0.049, 0.051, 0.10]
    ΔM = 0.1
    x, y = histmag(mags, ΔM)
    @test round.(x,digits=3) == [0.05,0.15]
    @test y == [3,1]
end