text:
using JLD2
using VortexDistributions
using GLMakie

## Load data
@load "src/VortexDetection3D/sim_examples/box_rings.jld2" psi_ring1 psi_ring2 psi_ring3 X

## Plot iso surface for select ring
function plot_iso(psi)
    density = abs2.(psi)
    pmax = maximum(density)
    density = density/pmax
    volume(density, algorithm = :iso, show_axis=true)
end
plot_iso(psi_ring3)

psi = psi_ring3 # Use "psi" as current sim you're working on
grad = gradient_3D_cent(psi, X)
grad_i = grad_itp(grad, X)
x = X[1]; y = X[2]; z = X[3]

dz = z[2] - z[1]

vorts_3d = find_vortices3D_v2(psi, X)
label = 1;
for zidx in 2:length(z)-2
    vorts_slice = vorts_3d[zidx];
    num_vorts = length(vorts_slice) # Number of vorts on z-slice
    if num_vorts != 0
        num_vorts = length(vorts_slice[:, 1])
        for v in 1:num_vorts
            if (vorts_slice[v, 5] == 0) # Check hasn't been labeled
                vorts_3d[zidx][v, 5] = label;
                label += 1;
            end

            vorts_sliceup = vorts_3d[zidx+1];

            if length(vorts_sliceup) != 0
                wps = wps_int(grad_i, [vorts_slice[v, 1], vorts_slice[v, 2], z[zidx]]);
                wps = wps .* (dz/abs(wps[3]));

                if vorts_slice[v, 4] > 0
                    approx_vort = vorts_slice[v, 1:2] + wps[1:2]
                else
                    approx_vort = vorts_slice[v, 1:2] - wps[1:2]
                end
                vort_index = find_closest_tuple(vorts_sliceup[:, 1:2], [approx_vort[1], approx_vort[2]]);
                if vorts_sliceup[vort_index, 4] == vorts_slice[v, 4]
                    vorts_sliceup[vort_index, 5] = vorts_slice[v, 5] # Assign same label
                end

            end

        end
    end
end
vorts_3d

# To see which layers have vortices found
# for i in 1:64
#     if length(vorts_3d[i]) > 0
#         println(i)
#     end
# end

vorts_label = Array{Float64, 2}[];
for i in 1:label-1
    current_vort = [1000 0 0 0 0];

    for zidx in 1:length(z)

        if length(vorts_3d[zidx]) != 0
            for v in 1:length(vorts_3d[zidx][:, 1])

                if vorts_3d[zidx][v, 5] == i
                    if current_vort[1] == 1000
                        current_vort = vorts_3d[zidx][v, :]';
                    else
                        current_vort = vcat(current_vort, vorts_3d[zidx][v, :]');
                    end
                end
            end
        end

    end
    push!(vorts_label, current_vort)
end

NUM_VORTS = length(vorts_label)
vorts_label[1]

dx = x[2] - x[1];
dy = y[2] - y[1];

vort_linked = Array{Float64, 2}[];
while (length(vorts_label) > 0)

    v = vorts_label[1]
    if !vortex_boundary_bottom(v, X, dx, dy, dz)
        vortex_link_bottom(1, vorts_label, X)
    end
    if !vortex_boundary_top(v, X, dx, dy, dz)
        vortex_link_top(1, vorts_label, X)
    end

    vorts_label = sort_vorts_label(vorts_label);
    reverse!(vorts_label);
    label = vorts_label[end][1, 5];
    while ((length(vorts_label) > 0) && (vorts_label[end][1, 5] == label))
        pop = pop!(vorts_label);
        push!(vort_linked, pop);
    end
    reverse!(vorts_label)
    println(length(vorts_label))
    # Pop all v in vorts_label with same label as v
end

vort_linked[1]
vort_linked[2]
using VortexDistributions


text:
using AGRN
using Base.Test

info("Chromosome tests")
include("test/chromosome.jl")
info("Mutation tests")
include("test/mutation.jl")
info("Crossover tests")
include("test/crossover.jl")
info("Evolution tests")
include("test/evolution.jl")

ext: jl


text:
# Work with the equation of state data files

using JLD  # for loading files

# Functions to read in the JLD files for quick EOS access
# Generates JLD files if they're not yet present

"Load a file from the data directory"
function loaddata(filename)
    load(joinpath(config.datadir, filename))
end

"Get the tabular equations of state from storage"
load_tabular_eoses() = loaddata("eos-tabular.jld")
"Get phase boundary information from storage"
load_phase_boundaries() = loaddata("phase-boundaries.jld")
"Get functional equations of state from storage"
load_functional_eoses() = loaddata("eos-functional.jld")
"Get equations of state defined piecewise from storage"
load_piecewise_eoses() = loaddata("eos-piecewise.jld")
"Get the complete stitched equation of state and thermal expansivity from storage"
load_full_eos() = loaddata("eos-full.jld")
"Get the heat capacity from storage"
load_heat_capacity() = loaddata("heatcapacity.jld")

function generate_jlds!()
    info("Generating data/phase-boundaries.jld")
    save_phase_boundaries!()
    info("Generating data/eos-functional.jld")
    save_functional_eoses!()
    info("Generating data/eos-tabular.jld")
    save_tabular_eoses!()
    info("Generating data/eos-piecewise.jld")
    save_piecewise_eoses!()
    info("Generating data/heatcapacity.jld")
    save_heat_capacity!()
    info("Generating data/eos-full.jld")
    save_full_eos!()

    nothing
end


text:
function speedstats(xy, Δt)
    if length(xy) < 2
        return missing
    end
    l = norm.(diff(xy))
    # L = cumsum(l)
    s = l/Δt
    μ = mean(s)
    σ = std(s, mean = μ)
    μ ± σ
end
speedstats(track) = (; (Symbol(k, :speed) => speedstats(getproperty(track, k), step(track.t)) for k in (:homing, :searching))...)

const intervals = sort([5, 10, 30, 60, 90, 120])
const nintervals = length(intervals) + 1

function coordinate2group(xy)
    l = norm(xy)
    for (i, L) in pairs(intervals)
        if l ≤ L
            return i
        end
    end
    return nintervals
end

function directionstats(xy, dropoff)
    if length(xy) < 2
        return missing
    end
    ss = [[Mean(), Mean()] for _ in 1:nintervals]
    for i in 2:length(xy)
        g = coordinate2group(xy[i] - dropoff)
        fit!.(ss[g], LinearAlgebra.normalize(xy[i] - xy[i-1]))
    end
    # [OnlineStats.nobs(s[1]) > 0 ? atand(reverse(value.(s))...) : missing for s in ss]
    [OnlineStats.nobs(s[1]) > 0 ? rad2deg(angular_diff_from_pos_y_axis(value.(s))) : missing for s in ss]
end
directionstats(track::Track, dropoff) = (; (Symbol(k, :direction) => directionstats(getproperty(track, k), dropoff) for k in (:homing, :searching))...)

function angular_diff_from_pos_y_axis(u)
    α = π/2 - atan(reverse(u)...)
    return α > π ? α - 2π : α
end

dropoff2tp(track::Track, dropoff) = (; vector_length = norm(track.turning_point - dropoff))

function finddiscrete(xy, dropoff)
    g = [coordinate2group(p - dropoff) for p in xy]
    Δ = diff(g)
    findall(!iszero, Δ)
end

function discretedirection(xy, dropoff)
    if length(xy) < 2
        return missing
    end
    map(finddiscrete(xy, dropoff)) do i
        isnothing(i) ? missing : rad2deg(angular_diff_from_pos_y_axis(LinearAlgebra.normalize(xy[i + 1] - xy[i])))
    end
end

discretedirection(track::Track, dropoff) = (; (Symbol(k, :discretedirection) => discretedirection(getproperty(track, k), dropoff) for k in (:homing, :searching))...)

tp_discretedirection(track::Track) = (; tpdiscretedirection = track.tp == 1 ? missing : rad2deg(angular_diff_from_pos_y_axis(LinearAlgebra.normalize(track.coords[track.tp] - track.coords[track.tp - 1]))))

path_length(s) = (; path_length = length(s.track.homing) > 1 ? sum(norm, diff(s.track.homing)) : missing, nest_corrected_vector_2 = -s.nest2feeder - last(s.track.turning_point), nest_corrected_vector_1 = last(s.track.turning_point))


text:
function delaunaytest(n::Int, dim::Int, reps::Int)
    for i = 1:reps
        pts = rand(n, dim)
        Simplices.Delaunay.delaunayn(pts)
    end

    return true
end

@testset "Delaunay" begin
    @test delaunaytest(10, 3, 10) == true
end


text:
# This file is a part of SimilaritySearch.jl

export rebuild

"""
    rebuild(g::SearchGraph; neighborhood=Neighborhood(), callbacks=SearchGraphCallbacks(), pools=getpools(index))

Rebuilds the `SearchGraph` index but seeing the whole dataset for the incremental construction, i.e.,
it can connect the i-th vertex to its knn in the 1..n possible vertices instead of its knn among 1..(i-1) as in the original algorithm.

# Arguments

- `g`: The search index to be rebuild.
- `neighborhood`: The neighborhood strategy to follow in the rebuild, it can differ from the original one.
- `callbacks`: The set of callbacks
- `pools`: The set of caches for the indexes

"""
function rebuild(g::SearchGraph; neighborhood=Neighborhood(), callbacks=SearchGraphCallbacks(), pools=getpools(g))
    n = length(g)
    @assert n > 0
    direct = Vector{Vector{Int32}}(undef, n)  # this separated links version needs has easier multithreading/locking needs
    reverse = Vector{Vector{Int32}}(undef, n)

    Threads.@threads for i in 1:n
        @inbounds direct[i] = find_neighborhood(g, g[i], neighborhood, pools, hints=g.links[i][1])
        reverse[i] = Vector{Int32}(undef, 0)
    end

    G = copy(g; links=direct, locks=copy(g.locks), hints=copy(g.hints), search_algo=copy(g.search_algo))
    _connect_reverse_links_neg(G.links, reverse, G.locks, 1, length(G))
    execute_callbacks(callbacks, G, force=true)
    G
end

function _connect_reverse_links_neg(direct, reverse, locks, sp, ep)
    Threads.@threads for i in sp:ep
        j = 0
        D = direct[i]
        @inbounds while j < length(D)
            j += 1
            id = D[j]
            if i == id
                D[j] = D[end]
                pop!(D)
                j -= 1
                continue
            end

            lock(locks[id])
            try
                push!(reverse[id], i)
            finally
                unlock(locks[id])
            end
        end
    end

    Threads.@threads for i in sp:ep
        @inbounds begin
            if length(direct[i]) >= length(reverse[i])
                append!(direct[i], reverse[i])
            else
                append!(reverse[i], direct[i])
                direct[i] = reverse[i]
            end
        end
    end
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
# *******************************************
#  Elementary Functions
# *******************************************
const ELEMFUNC = (
    :sin, :cos, :tan, :cot, :sec, :csc,
    :sinh, :cosh, :tanh, :coth, :sech, :csch,
    :asin, :acos, :atan, :acot, :asec, :acsc,
    :asinh, :acosh, :atanh, :acoth, :asech, :acsch,
    :sinpi, :cospi, :sinc, :cosc,
    :exp, :log, :sqrt, :cbrt,
    :abs
)

for fn = ELEMFUNC
    @eval begin
        Base.$fn(f::XFunc) = XFunc($fn.(f.vals), f.config)
        Base.$fn(f::KFunc) = KFunc($fn.(f.vals), f.config)
    end
end

# *******************************************
#  Analysis Tools
# *******************************************
function dVgen(config::ConfigFFT)
    xlens = (x -> -(.-(x...))).(config.xranges)
    return prod(xlens ./ config.ngrids)
end

function integ_X(f::XFunc)
    sum(f) * dVgen(f.config)
end

∫(f) = integ_X(f)

function norm_X(f::XFunc, p::Real=2)

    if p == Inf
        return max(abs(f)...)
    else
        return ( ∫( abs(f)^p ) )^(1/p)
    end

end

function l2inpr_X_X(f::XFunc{T,N}, g::XFunc{T,N}) where {T,N}
    check_config_consistency(f.config, g.config)
    return ∫(f * g)
end

# *******************************************
#  Supplimental Tools
# *******************************************
function xlens_xranges(xranges)
    (x -> -(-(x...))).(xranges)
end

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
using VegaLite, IndexedTables
function plot(br::BenchmarkResult)
    res = vec(br.sizedresults.results)
    brsizes = br.sizedresults.sizes
    sizes = Vector{eltype(brsizes)}(undef, length(res))
    ntests = length(br.tests) - 1
    for i ∈ 0:length(brsizes)-1
        si = brsizes[i+1]
        for j ∈ 1:ntests
            sizes[j + i*ntests] = si
        end
    end
    tests = vcat((@view(br.tests[2:end]) for _ ∈ eachindex(brsizes))...)
    t = table((GFLOPS = res, Size = sizes, Method = tests))
    t |> @vlplot(
        :line,
        x = :Size,
        y = :GFLOPS,
        color = :Method
    )
end


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "gradients" begin
    M = Circle()
    N = PowerManifold(M, 4)
    x = [0.1, 0.2, 0.3, 0.5]
    tvTestξ = [-1.0, 0.0, 0.0, 1.0]
    @test ∇TV(N, x) == tvTestξ
    @test ∇TV(M, (x[1], x[1])) ==
          (zero_tangent_vector(M, x[1]), zero_tangent_vector(M, x[1]))
    @test norm(N, x, ∇TV(N, x, 2) - tvTestξ) ≈ 0
    tv2Testξ = [0.0, 0.5, -1.0, 0.5]
    @test ∇TV2(N, x) == tv2Testξ
    @test norm(N, x, forward_logs(N, x) - [0.1, 0.1, 0.2, 0.0]) ≈ 0 atol = 10^(-16)
    @test norm(
        N,
        x,
        ∇intrinsic_infimal_convolution_TV12(N, x, x, x, 1.0, 1.0)[1] -
        [-1.0, 0.0, 0.0, 1.0],
    ) ≈ 0
    @test norm(N, x, ∇intrinsic_infimal_convolution_TV12(N, x, x, x, 1.0, 1.0)[2]) ≈ 0
    x2 = [0.1, 0.2, 0.3]
    N2 = PowerManifold(M, size(x2)...)
    @test ∇TV2(N2, x2) == zeros(3)
    @test ∇TV2(N2, x2, 2) == zeros(3)
    @test ∇TV(M, (0.0, 0.0), 2) == (0.0, 0.0)
    # 2d forward logs
    N3 = PowerManifold(M, 2, 2)
    N3C = PowerManifold(M, 2, 2, 2)
    x3 = [0.1 0.2; 0.3 0.5]
    x3C = cat(x3, x3; dims=3)
    tC = cat([0.2 0.3; 0.0 0.0], [0.1 0.0; 0.2 0.0]; dims=3)
    @test norm(N3C, x3C, forward_logs(N3, x3) - tC) ≈ 0 atol = 10^(-16)

    M = Circle()
    p = 0
    q = π / 4
    @test ∇distance(M, p, q) == q - p
end


text:
using ReachabilityAnalysis, Plots

# canonical direction along x₄
const e4 = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
const T = 20.0

@taylorize function laubloomis!(dx, x, params, t)
    dx[1] = 1.4*x[3] - 0.9*x[1]
    dx[2] = 2.5*x[5] - 1.5*x[2]
    dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    dx[4] = 2 - 1.3*(x[3]*x[4])
    dx[5] = 0.7*x[1] - (x[4]*x[5])
    dx[6] = 0.3*x[1] - 3.1*x[6]
    dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    return dx
end

function laubloomis(; W=0.01)

    # initial states
    X0c = [1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45]
    X0 = Hyperrectangle(X0c, fill(W, 7))

    # initial-value problem
    prob = @ivp(x' = laubloomis!(x), dim: 7, x(0) ∈ X0)

    return prob
end


text:
#
# `replace`
#
# Replace the values for a variable in the DataFrame if the `filter` condition is met.
# Operates row-by-row, therefore operators on scalars are expected. `replace` does not change the
# data type of the assigned variable.
#
# See also the notes for `generate`.
#
# this is the general form of the command
macro replace(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            replace <expression> [if] <expression>
        or
            bysort <varlist> (<varlist>): replace <expression> [if] <expression>
    """)
end

# short form. This is a copy of the above for `gen`.
macro rep(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            @replace(t, by, sort, arguments, filter, use, options)
        end
    )
end

# version with by
macro replace(t::Symbol,
    by::Vector{Symbol},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@replace_byrow!($t, $by, $sort, $arguments, $filter, $options )
        end
    )
end

# version without by
macro replace(t::Symbol,
    by::Nothing,
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@replace_byrow!($t, $sort, $arguments, $filter )
        end
    )
end

# version with `by`
macro replace_byrow!(t::Symbol, varlist_by::Vector{Symbol}, varlist_sort::Union{Vector{Symbol}, Nothing}, arguments::Expr, filter::Union{Expr, Nothing}, options::Union{Dict{String,Any}, Nothing})
    # assert that `arguments` is an assignment
    (arguments.head == :(=)) || error("`generate` expects an assignment operation, e.g. :x = :y + :z")
    # replace :varname by :varname[i] if not referenced
    Douglass.ref_quotenodes!(arguments)
    !isnothing(filter) && Douglass.ref_quotenodes!(filter)
    # get the assigned var symbol (note that it's in a QuoteNode)
    assigned_var = arguments.args[1].args[1].value
    # and the QuoteNode
    assigned_var_qn = arguments.args[1].args[1]

    # if the RHS of the assignment expression is currently a symbol, make it an Expr
    transformation = arguments.args[2]
    # return `nothing`s in case the index is < 1
    isexpr(transformation) && replace_invalid_indices!(transformation)

    return esc(
        quote
            # check variable is not present
            ($(assigned_var_qn) ∈ propertynames($t)) || error("Variable $($(assigned_var_qn)) not present in DataFrame.")

            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end
            #determine type of resulting column from the type of the first element
            #assigned_var_type = eltype($t[!,$assigned_var_qn])
            #$t[!,$(assigned_var)] = missings(assigned_var_type,size($t,1))

            # this is the function that maps every sub-df into its transformed df
            my_f = _df -> Douglass.@with _df begin
                # define _N
                _N = size(_df, 1)
                # fill the new variable, row by row
                for _n in 1:size(_df,1)
                    if (isnothing($filter) ? true : $filter)  # if condition is not satisfied, leave with missing
                        $(assigned_var_qn)[_n] = $(transformation)
                    end
                end
                _df
            end
            # execute, and copy it to another dataframe,
            # otherwise we get copies of the group variables in there as well
            t2 = combine(my_f, groupby($t,$varlist_by, sort=false, skipmissing = false ))
            Douglass.@with $t begin
                for _n = 1:size($t,1)
                    if (isnothing($filter) ? true : $filter)
                        $(assigned_var_qn)[_n] = t2[_n,^($(assigned_var_qn))]
                    end
                end
            end
            $t
        end
    )
end

# version without `by`
macro replace_byrow!(t::Symbol,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing})

    # assert that `arguments` is an assignment
    (arguments.head == :(=)) || error("`replace` expects an assignment operation, e.g. :x = :y + :z")
    # replace :varname by :varname[i] if not referenced
    Douglass.ref_quotenodes!(arguments)
    !isnothing(filter) && Douglass.ref_quotenodes!(filter)
    # get the assigned var symbol and QuoteNode
    if isa(arguments.args[1], Symbol)
        assigned_var_qn = QuoteNode(arguments.args[1])
    elseif isexpr(arguments.args[1]) && (arguments.args[1].head == :$)
        # interpolate
        # Now this becomes a Symbol
        assigned_var_qn = arguments.args[1].args[1] # ")  #QuoteNode(arguments.args[1].args[1]) #
    elseif isa(arguments.args[1].args[1], QuoteNode)
        # and the QuoteNode
        assigned_var_qn = arguments.args[1].args[1]
    else
        error("Assigned variable must be referenced using the colon syntax, e.g. :x.")
    end

    transformation = arguments.args[2]
    # return `nothing`s in case the index is < 1
    isexpr(transformation) && replace_invalid_indices!(transformation)

    return esc(
        quote
            # check variable is present
            ($(assigned_var_qn) ∉ propertynames($t)) && error("Variable $($(assigned_var_qn)) not present in DataFrame.")

            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, $varlist_sort)
            end

            Douglass.@with $t begin
                # define _N
                _N = size($t, 1)
                # fill the new variable, row by row
                for _n in 1:size($t,1)
                    if (isnothing($filter) ? true : $filter)  # if condition is not satisfied, leave with missing
                        $(assigned_var_qn)[_n] = $(transformation)
                        # This is what needs to be used if the assigned_var_qn is a Symbol
                        #$t[_n, $(assigned_var_qn)] = $(transformation)
                    end
                end
                $t
            end
            $t
        end
    )
end


text:
using Elemental, MPI
using MPIClusterManagers
const MCM = MPIClusterManagers

export ElRegressor, @mb_mpi_do

macro mb_mpi_do(mngr, ex)
    esc(:(if isa($mngr, MCM.MPIManager)
    MCM.@mpi_do($mngr, $ex)
    else
        @eval Main $ex
    end))
end

@with_kw mutable struct ElRegressor <: BaseEstimator
    W::Vector{Float32} = zeros(Float32, 0)
    name::String = "lstsq"
    alpha::Float32 = 0f0
end

is_classifier(::ElRegressor) = false

function gridparams(m::ElRegressor)
    grid = [
        "name" => ["lstsq", "ridge", "lasso"],
        "alpha" => [0.001, 0.01, 0.1, 1, 10],
    ]
    filter(gridparams(grid)) do d
        d["name"] != "lstsq" || d["alpha"] == 0.001
    end
end

function fit!(m::ElRegressor, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    @unpack name, alpha = m
    mngr = try Main.mngr catch end
    yr = reshape(y, 1, :)
    xr = reshape(x, size(x, 1), :)
    @eval Main xr, yr = $xr, $yr
    @mb_mpi_do mngr using TimerOutputs
    @mb_mpi_do mngr begin
        using Elemental, Distributed, SparseArrays, LinearAlgebra
        reset_timer!()
        if myid() != 1
            xr = spzeros(Float32, $(size(xr))...)
            yr = spzeros(Float32, $(size(yr))...)
        end
        B = Elemental.DistMatrix(Float32)
        X = Elemental.DistMatrix(Float32)
        @timeit "copy B" copy!(B, permutedims(yr))
        if !@isdefined(A) || size(A, 1) != size(B, 1)
            A = Elemental.DistMatrix(Float32)
            myid() == 1 && @info("copying to A...")
            @timeit "copy A" copy!(A, transpose(xr))
            Δ = Elemental.DistMatrix(Float32)
            Elemental.gaussian!(Δ, size(A)...)
            @timeit "A += 0.01Δ" BLAS.axpy!(5f-3, Δ, A)
        end
        if $name == "lstsq"
            Elemental.leastSquares!(A, B, X)
        elseif $name == "ridge"
            Elemental.ridge!(A, B, $alpha, X)
        elseif $name == "lasso"
            Elemental.bpdn!(A, B, $alpha, X)
        end
        X_array = Array(X)
        print_timer()
    end
    m.W = vec(Main.X_array)
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

predict(m::ElRegressor, x) = vec(transpose(m.W) * reshape(x, size(x, 1), :))

visualize(m::ElRegressor, columns) = write_feaimpt(m.W, columns)

modelhash(m::ElRegressor) = hash(m.W)


text:
include("benchmarks/vs_sort_code.jl")

using SortingLab;
using BenchmarkTools, DataFrames, CSV, Plots, StatPlots
using RCall;
N = 1_000_000;
K = 100;
randomseed = 1;

tic()
string_sort_perf(N, 100)
read_string_sort_perf_plot(N)
toc()

for N in Int.(10.^(7:8))
    string_sort_perf.(N, 100)
    read_string_sort_perf_plot(N)
end


text:
type Game
    input_tape
    output_tape
    gold_tape
    task
    head
    timestep
    prev_actions
    is_done

    function Game(x,y,task="copy")
        # make input tpae
        input_tape = nothing
        if in(task, ("copy","reverse"))
            x = convert(Array{Int64}, x)
            input_tape = reshape(x, 1, length(x))
        elseif in(task, ("add","mul","radd"))
            input_tape = make_grid(x)
        elseif task == "walk"
            input_tape = x
        end

        # make output tape
        output_tape = Int64[]

        # make gold tape
        gold_tape = nothing
        if in(task, ("copy","reverse","walk"))
            gold_tape = y
        else
            gold_tape = digits(y)
        end

        # previous action
        prev_actions = [STOP_ACTION]

        # head
        head = init_head(input_tape,task)
        timestep = 1

        new(input_tape, output_tape, gold_tape,
            task, head, timestep, prev_actions, false)
    end
end

function init_head(grid,task)
    return get_origin(grid,task)
end

function init_head(g::Game)
    return get_origin(g.input_tape,g.task)
end

function get_origin(grid,task)
    if in(task, ("copy","reverse"))
        return [1,1]
    elseif in(task, ("add","mul","radd"))
        return [1,size(grid,2)]
    elseif task == "walk"
        lastcol = grid[:,end]
        i = findfirst(x->x<0, lastcol)
        return [i,1]
    end
    error("invalid task: $task")
end

function move_timestep!(g::Game, write_symbol, move_action)
    write_action = write_symbol == NO_SYMBOL ? NOT_WRITE : WRITE
    write_action == WRITE && write!(g, write_symbol)
    move!(g, move_action)
    g.timestep += 1
    push!(g.prev_actions, (move_action, write_action))
    g.is_done = is_done(g)
end

function move!(g::Game, action)
    if action == "mr"
        g.head[2] += 1
    elseif action == "ml"
        g.head[2] -= 1
    elseif action == "<s>"
        g.is_done = true
    elseif action == "up"
        g.head[1] -= 1
    elseif action == "down"
        g.head[1] += 1
    else
        error("invalid action: $action")
    end
end

function write!(g::Game, symbol)
    (g.task in ("copy","reverse","walk")?push!:unshift!)(g.output_tape, symbol)
end

# FIXME: when it is done? right thing is to check the last action is <s> or not
function is_done(g::Game)
    len = length(g.output_tape)
    len >= length(g.gold_tape) && return true
    gold = nothing
    if g.task in ("copy","reverse","walk")
        gold = g.gold_tape[1:len]
    else
        gold = g.gold_tape[end-length(len):end]
    end
    return !(g.output_tape == gold)
end

function make_grid(x)
    x = collect(x)
    x = map(string, x)
    longest = length(x[1]); nelements = length(x)
    grid = -ones(Int64, nelements, longest)
    for (i,xi) in enumerate(x)
        for k = 1:length(xi)
            grid[i,end-length(xi)+k] = parse(Int64, xi[k])
        end
    end
    return grid
end

function reset!(g::Game)
    g.prev_actions = [STOP_ACTION]
    g.head = init_head(g.input_tape, g.task)
    g.timestep = 1
    empty!(g.output_tape)
    g.is_done = false
end

function read_symbol(grid, head)
    if 0 < head[1] <= size(grid,1) && 0 < head[2] <= size(grid,2)
        return grid[head...]
    end
    return -1
end

function read_symbol(g::Game)
    return read_symbol(g.input_tape, g.head)
end

# transitions
type Transition
    reward # +1 for true symbol output, 0 for otherwise
    input_symbol # read symbol - controller input
    prev_action # previous action - controller input
    nsteps # remaining steps
    h; c # controller states
    action # action has been taken
    output_symbol # symbol written to output tape
    is_done # last step or not
end

function take_action(w, b, s, steps_done; o=Dict())
    ei = get(o, :epsinit,  EPS_INIT)
    ef = get(o, :epsfinal, EPS_FINAL)
    ed = get(o, :epsdecay, EPS_DECAY)
    et = ef + (ei-ef) * exp(-steps_done/ed)

    # I think they did not use GLIE, they just unroll
    if rand() > et
        # @show w,s
        s = ndims(s) == 1 ? reshape(s,length(s),1) : s
        y = predict(w,b,s)
        return mapslices(indmax, Array(y), 1)
    else
        return rand(1:size(w,1), size(s)...)
    end
end

function take_action(w,b,s; eps=0.05)
    if rand() > eps
        s = ndims(s) == 1 ? reshape(s, length(s), 1) : s
        y = predict(w,b,s)
        return mapslices(indmax, Array(y), 1)
    else
        return rand(1:size(w,1), size(s,2))
    end
end

function get_reward(g::Game)
    length(g.output_tape) > length(g.gold_tape) && return 0
    if in(g.task, ("copy","reverse","walk"))
        is_desired = g.output_tape == g.gold_tape[1:length(g.output_tape)]
    else
        is_desired = g.output_tape == g.gold_tape[end-length(g.output_tape):end]
    end
    return Int(is_desired && length(g.output_tape) != 0)
end

function get_remaining_steps(g::Game)
    return length(g.gold_tape) - length(g.output_tape)
end

# deprecated old input functions, let me keep them for a while

# x: input tape, y: output tape, a: actions
function get_symgold(x,y,a,task)
    if task == "copy"
        return [y..., -1]
    elseif task == "reverse"
        return [-1, map(yi->-1, y)..., y..., -1]
    elseif task == "walk"
        return [map(i->-1, 1:size(x,2))..., y..., -1]
    elseif task == "add"
        x1digits = size(x,2)
        y1digits = length(y)
        symgold = mapreduce(i->[-1,y[i]], vcat, 1:x1digits)
        return vcat(symgold, [-1, (y1digits>x1digits?y[end]:-1)])
    else
        error("$task is not implemented yet")
    end
end

function make_data(games,s2i,a2i,actions)
    all_done = false
    inputs, outputs, masks = [], [], []

    t = 1
    while !all_done
        input = zeros(Cuchar, length(s2i)+length(a2i), length(games))
        mask = falses(1, length(games))
        symgold = NO_SYMBOL*ones(Int64, length(games))
        actgold = length(a2i)*ones(Int64, length(games))

        for (i,game) in enumerate(games)
            # skip if game is finished
            game.is_done && continue

            input_symbol = read_symbol(game.input_tape, game.head)
            input[s2i[input_symbol],i] = 1
            input[length(s2i)+a2i[game.prev_actions[end]],i] = 1

            move_action, write_action = actions[i][game.timestep]
            y = NO_SYMBOL
            if write_action == WRITE
                mask[1,i] = 1
                if game.task in ("copy","reverse","walk")
                    y = game.gold_tape[length(game.output_tape)+1]
                else
                    y = game.gold_tape[end-length(game.output_tape)]
                end
                symgold[i] = s2i[y]
            end

            actgold[i] = a2i[(move_action,write_action)]
            move_timestep!(game,y,move_action)
        end

        push!(inputs, input)
        push!(outputs, (symgold,actgold))
        push!(masks, mask)

        all_done = true
        for k = 1:length(games)
            if !games[k].is_done
                all_done = false
            end
        end

        # @show t,all_done
        t+=1
    end

    for k = 1:length(games); reset!(games[k]); end
    return inputs,outputs,masks
end

function make_input(g::Game, s2i, a2i)
    x1  = zeros(Float32, length(s2i), g.ninstances)

    # x1 => onehots, x11 => values, x12 => decoded (actions)
    x11 = map(i->read_symbol(g.input_tapes[i],g.pointers[i]), 1:g.ninstances)
    x12 = map(v->s2i[v], x11)
    for k = 1:length(x12); x1[x12[k],k] = 1; end

    # x2 => onehots, x21 => values, x22 => decoded (actions)
    x2  = zeros(Float32, length(a2i), g.ninstances)
    x21 = map(i->g.prev_actions[i][g.timestep], 1:g.ninstances)
    x22 = map(v->a2i[v], x21)
    for k = 1:length(x22); x2[x22[k],k] = 1; end

    return x1,x2
end

function make_inputs(g::Game, s2i, a2i)
    reset!(g)
    inputs = []
    for k = 1:length(g.prev_actions[1])
        push!(inputs, make_input(g,s2i,a2i))
        move_timestep!(g)
    end
    reset!(g)
    return inputs
end

function make_output(g::Game, s2i, a2i)
    y10 = map(i->g.symgold[i][g.timestep], 1:g.ninstances)
    y11 = map(yi->s2i[yi], y10)

    y20 = map(i->g.next_actions[i][g.timestep], 1:g.ninstances)
    y21 = map(yi->a2i[yi], y20)

    return y11, y21
end

function make_outputs(g, s2i, a2i)
    reset!(g)
    outputs = []
    for k = 1:length(g.next_actions[1])
        push!(outputs, make_output(g,s2i,a2i))
        move_timestep!(g)
    end
    reset!(g)
    return outputs
end

max_stars_repo_name: ereday/algorithm-learning

max_issues_repo_name: ereday/algorithm-learning

max_forks_repo_name: ereday/algorithm-learning


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2017/avx2/Compiler/gcc7.3/julia/1.3.0/bin/julia
#SBATCH --array=1-200
#SBATCH --mail-user=mkschleg@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o reproduce_test.%A_%a.out # Standard output
#SBATCH -e reproduce_test.%A_%a.err # Standard error
#SBATCH --mem-per-cpu=512M # Memory request of 2 GB
#SBATCH --time=0:01:00 # Running time of 12 hours
#SBATCH --account=rrg-whitem

using Pkg
Pkg.activate(".")

using Reproduce
const config_file = "./examples/arg_iter_config.toml"

function test_experiment()

    as = ArgParseSettings()
    @add_arg_table as begin
        "task_id"
        arg_type=Int
    end
    parsed = parse_args(as)

    experiment = Experiment(config_file)
    create_experiment_dir(experiment; tldr="hello tldr")
    ret = job(experiment, parsed["task_id"])
end

test_experiment()


text:
export setup_constraints

"""
input:

# Arguments
- constraint: dictionary with information about which constraints to use and their specifications
- comp_grid:  structure with computational grid information; comp_grid.n = (nx,ny,nz) number of gridpoints in each direction; comp_grid.d = (dx,dy,dz) spacing between gridpoints in each direction.
- TF:         Floating point precision (either Float32 or Float64)

output:
- P_sub    - vector of projection functions onto C: P_sub = [P_C_1(.); P_C_2(.); ... ; P_C_p(.)]
- TD_OP    - vector of linear operators [A_1;A_2;...;A_p]
- set_Prop - set properties, structure where each property is vector which has a lenght of the number of sets
                             (see SetIntersectionProjection.jl and setup_constraints.jl)
"""
function setup_constraints(constraint,comp_grid,TF)

# if    TF == Float64
#   TI = Int64
# elseif TF == Float32
#   TI = Int32
# end
TI=Int64

[@spawnat pid comp_grid for pid in workers()] #send computational grid to all workers in parallel mode

nr_constraints = length(constraint)

#convert all integers and floats to the desired precision
for i=1:nr_constraints
  if length(constraint[i].min) == 1 #scalar
    if typeof(constraint[i].min) in [Int32,Int64]
      #do nothing
    elseif typeof(constraint[i].min) in [Float32,Float64]
      constraint[i].min = TF(constraint[i].min)
      constraint[i].max = TF(constraint[i].max)
    end
  elseif length(constraint[i].min) > 1 #vector
    constraint[i].min = convert(Vector{TF},constraint[i].min)
    constraint[i].max = convert(Vector{TF},constraint[i].max)
  end
end

#allocate
P_sub = Vector{Any}(undef,nr_constraints)
TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef,nr_constraints)
AtA   = Vector{SparseMatrixCSC{TF,TI}}(undef,nr_constraints)

#initialize storage for set properties
set_Prop=set_properties(zeros(nr_constraints),zeros(nr_constraints),zeros(nr_constraints),Vector{Tuple{TI,TI}}(undef,nr_constraints),Vector{Tuple{String,String,String,String}}(undef,nr_constraints),zeros(nr_constraints),Vector{Vector{TI}}(undef,nr_constraints))

#other initialization
special_operator_list = ["DFT", "DCT", "wavelet"] #complex valued operators that are orthogonal
N =  prod(comp_grid.n)

for i = 1:nr_constraints

      #verify input consistency
      if (constraint[i].set_type in ["nuclear", "rank"]) && (constraint[i].app_mode[1] in ["matrix","tensor"]) && (length(comp_grid.n)==3)
        error("requested rank or nuclear norm constraints on a tensor, use mode=(slice,x) e.t.c. to define constraints per slice")
      end

      # #catch a few things not currently implemented, yet..
      if constraint[i].set_type in ["l1", "l2"] && constraint[i].app_mode[1] in ["slice", "fiber"]
        error("l1 and l2 constraints only available for matrix or tensor mode, currently")
      end

      (A,AtA_diag,dense,TD_n,banded) = get_TD_operator(comp_grid,constraint[i].TD_OP,TF)

      P_sub[i] = get_projector(constraint[i],comp_grid,special_operator_list,A,TD_n,TF)

      if constraint[i].TD_OP in special_operator_list
        TD_OP[i] =  SparseMatrixCSC{TF}(LinearAlgebra.I,N,N)
      else
        TD_OP[i] = A
      end

      set_Prop.AtA_diag[i] = AtA_diag
      set_Prop.dense[i]    = dense
      set_Prop.TD_n[i]     = TD_n
      set_Prop.banded[i]   = banded
      set_Prop.tag[i]      = (constraint[i].set_type,constraint[i].TD_OP,constraint[i].app_mode[1],constraint[i].app_mode[2])

      #determine set convexity
      if constraint[i].set_type in ["rank","cardinality"]
        set_Prop.ncvx[i]     = true
      elseif constraint[i].set_type == "bounds" && constraint[i].TD_OP != "identity" && TF(maximum(constraint[i].min))>TF(0.0)
        set_Prop.ncvx[i]     = true
      elseif constraint[i].set_type == "histogram" && constraint[i].TD_OP != "identity" && TF(maximum(constraint[i].min))>TF(0.0)
        set_Prop.ncvx[i]     = true
      else
        set_Prop.ncvx[i]     = false
      end

end

return P_sub,TD_OP,set_Prop
end #end setup_constraints


text:
function demo_pulsed_deck(cur_T_bar::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = Dict(
    :deck => :demo_pulsed,
    :T_bar => cur_T_bar,
    :H => 1.1,
    :Q => 39.86,
    :epsilon => 0.3226,
    :kappa_95 => 1.590,
    :delta_95 => 0.333,
    :nu_n => 0.27,
    :nu_T => 1.094,
    :Z_eff => 2.584,
    :f_D => 0.7753,
    :A => 2.735,
    :l_i => 1.155,
    :N_G => 1.2,
    :eta_T => 0.3531,
    :eta_RF => 0.4,
    :tau_FT => 7273.0,
    :C_ejima => 0.3,
    :max_beta_N => 0.02591,
    :max_q_95 => 2.5,
    :max_P_E => 1500.0,
    :max_P_W => 8.0,
    :max_q_DV => 10.0,
    :B_CS => 12.77,
    :sigma_CS => 660,
    :sigma_TF => 660,
    :J_CS => 50,
    :J_TF => 230,
    :k => 0.333,
    :eta_CD => 0.2721,
    :wave_theta => 135.0,
    :eta_LH => 0.6
  )

  merge!(cur_dict, Dict(cur_kwargs))

  BaseReactor(cur_dict)
end


text:
module CsminWel

using Reexport
using ForwardDiff
@reexport using Optim
using Calculus
using LinearAlgebra

import Base.minimum
import Random: AbstractRNG, MersenneTwister
import Printf.@printf
import Optim: AbstractOptimizer,
              OptimizationResults,
              MultivariateOptimizationResults,
              optimize, initial_state, minimizer,
              minimum, iterations, converged, x_converged,
              f_converged, f_tol, g_tol, g_converged, iteration_limit_reached,
              f_calls, method

#=
This code is based on a routine originally copyright Chris Sims.
See http://sims.princeton.edu/yftp/optimize/
It was transcripted in julia copyright Giuseppe Ragusa
Copyright (c) 2015, Federal Reserve Bank of New York Copyright (c) 2015, Chris Sims Copyright (c) 2015, Giuseppe Ragusa
=#

struct Csminwel <: Optim.AbstractOptimizer end

## Changed respoct to default to store Hessian matrix
mutable struct MultivariateOptimizationResultsCs{O<:AbstractOptimizer,T,N} <: OptimizationResults
    method::String
    initial_x::Array{T,N}
    minimizer::Array{T,N}
    minimum::T
    iterations::Int
    iteration_converged::Bool
    x_converged::Bool
    x_tol::Float64
    f_converged::Bool
    f_tol::Float64
    g_converged::Bool
    g_tol::Float64
    trace::OptimizationTrace{T,O}
    f_calls::Int
    g_calls::Int
    h_calls::Int
    invH::Array{T, 2}
end

initial_state(r::MultivariateOptimizationResultsCs) = r.initial_x
minimum(r::MultivariateOptimizationResultsCs) = r.minimum
minimizer(r::MultivariateOptimizationResultsCs) = r.minimizer
converged(r::MultivariateOptimizationResultsCs) = r.x_converged || r.f_converged || r.g_converged

x_converged(r::MultivariateOptimizationResultsCs) = r.x_converged
f_converged(r::MultivariateOptimizationResultsCs) = r.f_converged
g_converged(r::MultivariateOptimizationResultsCs) = r.g_converged

x_tol(r::MultivariateOptimizationResultsCs) = r.x_tol
f_tol(r::MultivariateOptimizationResultsCs) = r.f_tol
g_tol(r::MultivariateOptimizationResultsCs) = r.g_tol

function Base.show(io::IO, r::MultivariateOptimizationResultsCs)
    @printf io "Results of Optimization Algorithm\n"
    @printf io " * Algorithm: %s\n" summary(r)
    if length(join(initial_state(r), ",")) < 40
        @printf io " * Starting Point: [%s]\n" join(initial_state(r), ",")
    else
        @printf io " * Starting Point: [%s, ...]\n" join(initial_state(r)[1:2], ",")
    end
    if length(join(minimizer(r), ",")) < 40
        @printf io " * Minimizer: [%s]\n" join(minimizer(r), ",")
    else
        @printf io " * Minimizer: [%s, ...]\n" join(minimizer(r)[1:2], ",")
    end
    @printf io " * Minimum: %e\n" minimum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: %s\n" converged(r)
    @printf io "   * |x - x'| < %.1e: %s\n" x_tol(r) x_converged(r)
    @printf io "   * |f(x) - f(x')| / |f(x)| < %.1e: %s\n" f_tol(r) f_converged(r)
    @printf io "   * |g(x)| < %.1e: %s\n" g_tol(r) g_converged(r)
    @printf io "   * Reached Maximum Number of Iterations: %s\n" iteration_limit_reached(r)
    @printf io " * Objective Function Calls: %d\n" f_calls(r)
    return
end

const VERBOSITY = Dict(:none => 0, :low => 1, :high => 2)
const rc_messages = Dict(0 => "Standard Iteration",
1 => "zero gradient",
2 => "back and forth on step length never finished",
3 => "smallest step still improving too slow",
4 => "back and forth on step length never finished",
5 => "largest step still improving too fast",
6 => "smallest step still improving too slow, reversed gradient",
7 => "warning: possible inaccuracy in H matrix")

function Optim.optimize(f::Function,
    initial_x::Array,
    method::Csminwel,
    options::Optim.Options = Optim.Options();
    H0::Matrix             = 1e-5 .* Matrix{eltype(initial_x)}(I, length(initial_x), length(initial_x)),
    rng::AbstractRNG       = MersenneTwister(0))

    opts = extract_options(options)
    d = Optim.OnceDifferentiable(f, initial_x; autodiff = :finite)
    function gradwrap(g!, x)
        stor = Array{Float64}(undef, length(x))
        g!(stor, x)
        bad_grads = abs.(stor) .>= 1e15
        stor[bad_grads] .= 0.0
        return stor, any(bad_grads)
    end
    csminwel(f, x -> gradwrap(d.df,x), initial_x; H0 = H0, rng = rng, opts...)
end

function Optim.optimize(d::Optim.OnceDifferentiable,
    initial_x::Array,
    method::Csminwel,
    options::Optim.Options = Optim.Options();
    H0::Matrix             = 1e-5 .* Matrix{eltype(d.x_f)}(I, length(d.x_f), length(d.x_f)),
    rng::AbstractRNG       = MersenneTwister(0))

    opts = extract_options(options)

    function gradwrap(g!, x)
        stor = Array{Float64}(undef, length(x))
        g!(stor, x)
        bad_grads = abs.(stor) .>= 1e15
        stor[bad_grads] .= 0.0
        return stor, any(bad_grads)
    end
    csminwel(d.f, x -> gradwrap(d.df,x), initial_x; H0 = H0, rng = rng, opts...)
end

function extract_options(opt::Optim.Options)
    opts = Array{Any}(undef, 0)
    for name in fieldnames(typeof(opt))
        push!(opts, (name, getfield(opt, name)))
    end
    opts
end

function Optim.optimize(f::Function,
    g!::Function,
    initial_x::Array{T,1},
    method::Csminwel,
    options::Optim.Options = Optim.Options();
    H0::Matrix             = 1e-5.* Matrix{eltype(initial_x)}(I, length(initial_x), length(initial_x)),
    rng::AbstractRNG       = MersenneTwister(0)) where T <: Real

    opts = extract_options(options)

    function gradwrap(g!, x)
        stor = Array{Float64}(undef, length(x))
        g!(stor, x)
        bad_grads = abs.(stor) .>= 1e15
        stor[bad_grads] .= 0.0
        return stor, any(bad_grads)
    end
    csminwel(f, x->gradwrap(g!, x), initial_x; H0 = H0, rng = rng, opts...)
end

function csminwel(fcn::Function,
    grad::Function,
    x0::Vector;
    H0::Matrix           = 1e-5 .* Matrix{eltype(x0)}(I, length(x0), length(x0)),
    rng::AbstractRNG     = MersenneTwister(0),
    verbose::Symbol      = :high,
    x_tol::Real          = 1e-32,  # default from Optim.jl
    f_tol::Float64       = 1e-5, # 1e-7 # 1e-14,  # Default from csminwel
    g_tol::Real          = 1e-8,   # default from Optim.jl
    iterations::Int      = 1000,
    store_trace::Bool    = false,
    show_trace::Bool     = false,
    extended_trace::Bool = false,
    kwargs...)

    xtol, ftol, grtol  = x_tol, f_tol, g_tol

    if show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end

    # unpack dimensions
    nx = size(x0, 1)

    # Count function and gradient calls
    f_calls, g_calls, h_calls = 0, 0, 0

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(x0), copy(x0)

    # start with Initial Hessian
    H = H0

    # start rc parameter at 0
    rc = 0

    f_x = fcn(x0)
    f_calls += 1

    if f_x > 1e50
        throw(ArgumentError("Bad initial guess. Try again"))
    end

    gr, badg = grad(x0)
    g_calls += 1

    # Count iterations
    iteration = 0

    # Maintain a trace
    tr = OptimizationTrace{typeof(f_x), Csminwel}()
    tracing = show_trace || store_trace || extended_trace

    if tracing
        dt = Dict()
        if extended_trace
            dt["x"] = copy(x)
            dt["g(x)"] = copy(gr)
            dt["H(x)"] = copy(H)
            dt["rc"] = copy(rc)
            dt["rc_message"] = rc_messages[rc]
        end
        grnorm = norm(gr, Inf)
        Optim.update!(tr, iteration, f_x, grnorm, dt, store_trace, show_trace)
    end

    # set objects to their starting values
    retcode3 = 101

    # set up return variables so they are available outside while loop
    fh = copy(f_x)
    xh = copy(x0)
    gh = copy(x0)
    retcodeh = 1000

    # Assess multiple types of convergence
    x_converged, f_converged, gr_converged = false, false, false

    # Iterate until convergence or exhaustion
    converged = false
    while !converged && iteration < iterations
        iteration += 1

        f1, x1, fc, retcode1 = csminit(fcn, x, f_x, gr, badg, H; verbose=verbose)
        f_calls += fc

        if retcode1 != 1
            if retcode1 == 2 || retcode1 == 4
                wall1, badg1 = true, true
            else
                g1, badg1 = grad(x1)
                g_calls += 1
                wall1 = badg1
            end

            # Bad gradient or back and forth on step length.  Possibly at cliff edge. Try
            # perturbing search direction if problem not 1D
            if wall1 && (length(H) > 1)

                Hcliff = H + Diagonal(diag(H)).*rand(rng, nx)

                if VERBOSITY[verbose] >= VERBOSITY[:low]
                    @printf "Cliff.  Perturbing search direction.\n"
                end

                f2, x2, fc, retcode2 = csminit(fcn, x, f_x, gr, badg, Hcliff; verbose=verbose)
                f_calls += fc

                if f2 < f_x
                    if retcode2==2 || retcode2==4
                        wall2 = true; badg2 = true
                    else
                        g2, badg2 = grad(x2)
                        g_calls += 1
                        wall2 = badg2
                        badg2
                    end

                    if wall2
                        if VERBOSITY[verbose] >= VERBOSITY[:low]
                            @printf "Cliff again.  Try traversing\n"
                        end

                        if norm(x2-x1) < 1e-13
                            f3 = f_x
                            x3 = x
                            badg3 = true
                            retcode3 = 101
                        else
                            gcliff = ( (f2-f1) / ((norm(x2-x1))^2) )*(x2-x1)
                            if (size(x0 , 2)>1)
                                gcliff = gcliff'
                            end
                            f3, x3, fc, retcode3 = csminit(fcn, x, f_x, gcliff,
                            false, Matrix{eltype(nx)}(I, nx, nx); verbose=verbose)
                            f_calls += fc

                            if retcode3==2 || retcode3==4
                                wall3 = true
                                badg3 = true
                            else
                                g3, badg3 = grad(x3)
                                g_calls += 1
                                wall3 = badg3
                            end
                        end
                    else
                        f3 = f_x
                        x3 = x
                        badg3 = true
                        retcode3 = 101
                    end
                else
                    f3 = f_x
                    x3 = x
                    badg3 = true
                    retcode3 = 101
                end
            else
                # normal iteration, no walls, or else 1D, or else we're finished here.
                f2, f3 = f_x, f_x
                badg2, badg3 = true, true
                retcode2, retcode3 = 101, 101
            end
        else
            f1, f2, f3 = f_x, f_x, f_x
            retcode2, retcode3 = retcode1, retcode1
        end

        # how to pick gh and xh
        if f3 < f_x - ftol && badg3==0
            ih = 3
            fh = f3
            xh = x3
            gh = g3
            badgh = badg3
            retcodeh = retcode3
        elseif f2 < f_x - ftol && badg2==0
            ih = 2
            fh = f2
            xh = x2
            gh = g2
            badgh = badg2
            retcodeh = retcode2
        elseif f1 < f_x - ftol && badg1==0
            ih = 1
            fh = f1
            xh = x1
            gh = g1
            badgh = badg1
            retcodeh = retcode1
        else
            fh, ih = findmin([f1 , f2 , f3])

            if ih == 1
                xh = x1
                retcodeh = retcode1
            elseif ih == 2
                xh = x2
                retcodeh = retcode2
            elseif ih == 3
                xh = x3
                retcodeh = retcode3
            end

            if @isdefined gh
                nogh = isempty(gh)
            else
                nogh = true
            end

            if nogh
                gh, badgh = grad(xh)
                g_calls += 1
            end

            badgh = true
        end

        stuck = (abs(fh-f_x) < ftol)
        if !badg && !badgh && !stuck
            H = bfgsi(H , gh-gr , xh-x; verbose=verbose)
        end

        if VERBOSITY[verbose] >= VERBOSITY[:high]
            @printf "Improvement on iteration %d = %18.9f\n" iteration fh-f_x
        end

        if stuck
            if VERBOSITY[verbose] >= VERBOSITY[:low]
                @printf "improvement < ftol -- terminating\n"
            end
        end

        # record# retcodeh of previous x
        copyto!(x_previous, x)

        # update before next iteration
        f_x_previous, f_x = f_x, fh
        x = xh
        gr = gh
        badg = badgh

        # Check convergence
        x_converged, f_converged, gr_converged, converged =
            assess_convergence(x, x_previous, f_x, f_x_previous, gr, xtol, ftol, grtol)

        if tracing
            dt = Dict()
            if extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(gr)
                dt["H(x)"] = copy(H)
                dt["rc"] = copy(rc)
                dt["rc_message"] = rc_messages[rc]
            end
            grnorm = norm(gr, Inf)
            Optim.update!(tr, iteration, f_x, grnorm, dt, store_trace, show_trace)
        end

    end

    return MultivariateOptimizationResultsCs("Csminwel", x0, x, convert(Float64, f_x),
    iteration, iteration==iterations, x_converged, xtol, f_converged, ftol,
    gr_converged, grtol, tr, f_calls, g_calls, h_calls, H)   # also return H
end

#
#=
Version of `csminwel` that will use finite differencing methods to
approximate the gradient numerically. This is convenient for cases where
you cannot supply an analytical derivative, but it is not as robust as
using the true derivative.
=#

# """
# ```
# csminwel(fcn::Function, grad::Function, x0::Vector, H0::Matrix=1e-5.*eye(length(x0)), args...;
# xtol::Real=1e-32, ftol::Float64=1e-14, grtol::Real=1e-5, iterations::Int=1000,
# store_trace::Bool = false, show_trace::Bool = false, extended_trace::Bool = false,
# verbose::Symbol = :none, rng::AbstractRNG = MersenneTwister(), kwargs...)
# ```

# Minimizes `fcn` using the csminwel algorithm.

# ### Arguments

# * `fcn::Function` : The objective function
# * `grad::Function` : The gradient of the objective function. This argument can be omitted if
# an analytical gradient is not available, which will cause a numerical gradient to be
# calculated.
# * `x0::Vector`: The starting guess for the optimizer

# ### Optional Arguments

# * `H0::Matrix`: An initial guess for the Hessian matrix -- must be
# positive definite. If none is given, then a scaled down identity
# matrix is used.
# * `args...`:  Other positional arguments to be passed to `f` on each
# function call

#     ### Keyword Arguments

#     * `ftol::{T<:Real}=1e-14`: Threshold for convergence in terms of change
#     in function value across iterations.
#         * `iterations::Int=100`: Maximum number of iterations
#         * `kwargs...`: Other keyword arguments to be passed to `f` on each
#         function call
#             """

function csminwel(fcn::Function, x0::Vector{T}; kwargs...) where T <: Real
    x0 = float(x0)
    autodiff = false
    for fn in kwargs
        if fn[1] == :autodiff
            autodiff = fn[2]
        end
    end
    grad(x::Array{T}) where T<:Number = csminwell_grad(fcn, x, Val{autodiff})
    csminwel(fcn, grad, x0; kwargs...)
end

function csminwell_grad(fcn, x, ::Type{Val{false}})
    gr = Calculus.gradient(fcn, x)
    bad_grads = abs.(gr) .>= 1e15
    gr[bad_grads] .= 0.0
    return gr, any(bad_grads)
end

function csminwell_grad(fcn, x, ::Type{Val{true}})
    gr = ForwardDiff.gradient(fcn, x)
    bad_grads = abs.(gr) .>= 1e15
    gr[bad_grads] .= 0.0
    return gr, any(bad_grads)
end

function csminit(fcn, x0::Vector{T}, f0, g0, badg, H0, args...; verbose::Symbol=:high, kwargs...) where T <: Real
    x0 = float(x0)
    angle = .005

    #(0<THETA<.5) THETA near .5 makes long line searches, possibly fewer iterations.
    theta = .3
    fchange = 1000
    minlamb = 1e-9
    mindfac = .01
    f_calls = 0
    lambda = 1.0
    xhat = x0
    f = f0
    fhat = f0
    gr = g0
    gnorm = norm(gr)

    if gnorm < 1e-12 && !badg
        # gradient convergence
        retcode = 1
        dxnorm = 0.0
    else
        # with badg true, we don't try to match rate of improvement to directional
        # derivative.  We're satisfied just to get some improvement in f.
        dx = vec(-H0*gr)
        dxnorm = norm(dx)

        if dxnorm > 1e12

            if VERBOSITY[verbose] >= VERBOSITY[:low]
                @printf "Near singular H problem.\n"
            end

            dx = dx * fchange / dxnorm
        end

        dfhat = dot(dx, g0)

        if !badg
            # test for alignment of dx with gradient and fix if necessary
            a = -dfhat / (gnorm*dxnorm)

            if a < angle
                dx -= (angle*dxnorm/gnorm + dfhat/(gnorm*gnorm)) * gr
                dx *= dxnorm/norm(dx)
                dfhat = dot(dx, gr)

                if VERBOSITY[verbose] >= VERBOSITY[:low]
                    @printf "Correct for low angle %f\n" a
                end
            end
        end

        if VERBOSITY[verbose] >= VERBOSITY[:high]
            @printf "Predicted Improvement: %18.9f\n" (-dfhat/2)
        end
        # Have OK dx, now adjust length of step (lambda) until min and max improvement rate
        # criteria are met.
        done = false
        fact = 3.0
        shrink = true
        lambda_min = 0.0
        lambda_max = Inf
        lambda_peak = 0.0
        f_peak = f0
        lambda_hat = 0.0

        while !done
            if size(x0, 2) > 1
                dxtest = x0 + dx' * lambda
            else
                dxtest = x0 + dx * lambda
            end

            f = fcn(dxtest, args...; kwargs...)

            if VERBOSITY[verbose] >= VERBOSITY[:high]
                @printf "lambda = %10.5f; f = %20.7f\n" lambda f
            end

            if f < fhat
                fhat = f
                xhat = dxtest
                lambdahat = lambda
            end

            f_calls += 1
            shrink_signal = (!badg & (f0-f < maximum([-theta*dfhat*lambda 0]))) || (badg & ((f0-f) < 0))

            grow_signal = !badg && (lambda > 0)  &&
            (f0-f > -(1-theta)*dfhat*lambda)

            if shrink_signal && ((lambda > lambda_peak) || lambda < 0 )
                if (lambda > 0) && ((!shrink) || (lambda/fact <= lambda_peak))
                    shrink = true
                    fact = fact^.6
                    while lambda / fact <= lambda_peak
                        fact = fact^.6
                    end

                    if abs(fact - 1.0) < mindfac
                        if abs(lambda) < 4
                            retcode = 2
                        else
                            retcode = 7
                        end

                        done = true
                    end
                end

                if lambda < lambda_max && lambda > lambda_peak
                    lambda_max = lambda
                end

                lambda /= fact
                if abs(lambda) < minlamb
                    if (lambda > 0) && (f0 <= fhat)
                        lambda = -lambda*fact^6
                    else
                        if lambda < 0
                            retcode = 6
                        else
                            retcode = 3
                        end
                        done = true
                    end
                end

            elseif (grow_signal && lambda > 0) || (shrink_signal &&
                ((lambda <= lambda_peak) && (lambda > 0)))
                if shrink
                    shrink = false
                    fact = fact^.6
                    if abs(fact - 1) < mindfac
                        if abs(lambda) < 4
                            retcode = 4
                        else
                            retcode = 7
                        end
                        done = true
                    end
                end

                if f < f_peak && lambda > 0
                    f_peak = f
                    lambda_peak = lambda
                    if lambda_max <= lambda_peak
                        lambda_max = lambda_peak * fact^2
                    end
                end

                lambda *= fact
                if abs(lambda) > 1e20
                    retcode = 5
                    done = true
                end
            else
                done = true
                if fact < 1.2
                    retcode = 7
                else
                    retcode = 0
                end
            end
        end
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        @printf "Norm of dx %10.5f\n" dxnorm
    end

    return fhat, xhat, f_calls, retcode
end

"""
```
bfgsi(H0, dg, dx)
```
### Arguments
- `H0`: hessian matrix
- `dg`: previous change in gradient
- `dx`: previous change in x
"""
function bfgsi(H0, dg, dx; verbose::Symbol = :high)
    if size(dg, 2) > 1
        dg = dg'
    end

    if size(dx, 2) > 1
        dx = dx'
    end

    Hdg = H0*dg
    dgdx = dot(dx, dg)

    H = H0
    if abs(dgdx) > 1e-12
        H += (dgdx.+(dg'*Hdg)).*(dx*dx')/(dgdx^2) - (Hdg*dx'.+dx*Hdg')/dgdx
    elseif norm(dg) < 1e-7
        # gradient is super small so don't worry updating right now
        # do nothing
    else
        if VERBOSITY[verbose] >= VERBOSITY[:high]
            @warn("bfgs update failed")
            @printf "|dg| = %f, |dx| = %f\n" (norm(dg)) (norm(dx))
            @printf "dg'dx = %f\n" dgdx
            @printf "|H*dg| = %f\n" (norm(Hdg))
        end
    end

    return H
end

function assess_convergence(x::Array,
    x_previous::Array,
    f_x::Real,
    f_x_previous::Real,
    gr::Array,
    xtol::Real,
    ftol::Real,
    grtol::Real)
    x_converged, f_converged, gr_converged = false, false, false

    if Optim.maxdiff(x, x_previous) < xtol
        x_converged = true
    end

    # Relative Tolerance
    # if abs(f_x - f_x_previous) / (abs(f_x) + ftol) < ftol || nextfloat(f_x) >= f_x_previous
    # Absolute Tolerance
    if abs(f_x - f_x_previous) < ftol
        f_converged = true
    end

    if norm(vec(gr), Inf) < grtol
        gr_converged = true
    end

    converged = x_converged || f_converged || gr_converged

    return x_converged, f_converged, gr_converged, converged
end

## CSMINWEL
function getgradient(fcn, grad_f, s::Csminwel, ::Type{Val{true}}, ::Type{Val{false}})
    function gradient(x)
        gr = similar(x)
        grad_f(gr, x)
        bad_grads = abs.(gr) .>= 1e15
        gr[bad_grads] .= 0.0
        return gr, any(bad_grads)
    end
end

function getgradient(fcn, grad_f, s::Csminwel, ::Type{Val{false}}, ::Type{Val{true}})
    function gradient(x)
        gr = similar(x)
        ForwardDiff.gradient!(gr, fcn, x)
        bad_grads = abs.(gr) .>= 1e15
        gr[bad_grads] .= 0.0
        return gr, any(bad_grads)
    end
end

function getgradient(fcn, grad_f, s::Csminwel, ::Type{Val{false}}, ::Type{Val{false}})
    function gradient(x)
        gr = Calculus.gradient(fcn, x)
        bad_grads = abs.(gr) .>= 1e15
        gr[bad_grads] .= 0.0
        return gr, any(bad_grads)
    end
end

export Csminwel

end


text:
onednn_url = "https://github.com/oneapi-src/OneDNN"
onednn_tag = "v2.4.2"
install_dir = joinpath(@__DIR__, "dnnl")

onednn_src = joinpath(@__DIR__, "oneDNN")
onednn_build = joinpath(onednn_src, "build")

# Clone and checkout.
if !isdir(onednn_src)
    run(`git clone $onednn_url $onednn_src`)
    cd(onednn_src)
    run(`git checkout $onednn_tag`)
end

# Apply performance patch.
patchdir = joinpath(@__DIR__, "patches")
patchfile = joinpath(patchdir, "cnn.patch")
cd(onednn_src)

run(`git apply $patchfile`)
# try
#     run(`git apply $patchfile`)
# catch e
#     isa(e, ProcessFailedException) || rethrow(e)
# end

# Run cmake
cc = "clang"
cxx = "clang++"
cmake_args = [
    "-DCMAKE_BUILD_TYPE=Release",
    "-DCMAKE_INSTALL_PREFIX=$(install_dir)",
    "-DDNNL_LIBRARY_TYPE=SHARED",
    "-DCMAKE_CXX_COMPILER=$cc",
    "-DCMAKE_C_COMPILER=$cc",
    "-DDNNL_BUILD_TESTS=OFF",
    "-DDNNL_BUILD_EXAMPLES=OFF",
    "-DDNNL_CPU_RUNTIME=THREADPOOL",
]

isdir(onednn_build) || mkdir(onednn_build)
cd(onednn_build)
run(`cmake .. $cmake_args`)
nproc = parse(Int, read(`nproc`, String))
run(`make -j$nproc`)
run(`make install`)

# Cleanup
cd(@__DIR__)

#####
##### threadpool
#####

import CxxWrap

# Paths for linking
cxxhome = CxxWrap.prefix_path()
juliahome = dirname(Base.Sys.BINDIR)
dnnlhome = joinpath(@__DIR__, "dnnl")

# Use Clang since it seems to get along better with Julia
cxx = "clang++"

cxxflags = [
    "-g",
    "-O3",
    "-Wall",
    "-fPIC",
    "-std=c++17",
    "-DPCM_SILENT",
    "-DJULIA_ENABLE_THREADING",
    "-Dexcept_EXPORTS",
    # Surpress some warnings from Cxx
    "-Wno-unused-variable",
    "-Wno-unused-lambda-capture",
]

includes = [
    "-I$(joinpath(cxxhome, "include"))",
    "-I$(joinpath(juliahome, "include", "julia"))",
    "-I$(joinpath(dnnlhome, "include"))",
]

loadflags = [
    # Linking flags for Julia
    "-L$(joinpath(juliahome, "lib"))",
    "-Wl,--export-dynamic",
    "-Wl,-rpath,$(joinpath(juliahome, "lib"))",
    "-ljulia",
    # Linking Flags for CxxWrap
    "-L$(joinpath(cxxhome, "lib"))",
    "-Wl,-rpath,$(joinpath(cxxhome, "lib"))",
    "-lcxxwrap_julia",
]

src = joinpath(@__DIR__, "threadpool.cpp")
so = joinpath(@__DIR__, "libthreadpool.so")

cmd = `$cxx $cxxflags $includes -shared $src -lpthread -o $so $loadflags`
@show cmd
run(cmd)


text:
abstract type Estimator end
abstract type LGBMEstimator <: Estimator end

mutable struct LGBMRegression <: LGBMEstimator
    booster::Booster
    model::String
    application::String
    boosting::String

    num_iterations::Int
    learning_rate::Float64
    num_leaves::Int
    max_depth::Int
    tree_learner::String
    num_threads::Int
    histogram_pool_size::Float64

    min_data_in_leaf::Int
    min_sum_hessian_in_leaf::Float64
    lambda_l1::Float64
    lambda_l2::Float64
    min_gain_to_split::Float64
    feature_fraction::Float64
    feature_fraction_seed::Int
    bagging_fraction::Float64
    bagging_freq::Int
    bagging_seed::Int
    early_stopping_round::Int

    max_bin::Int
    data_random_seed::Int
    init_score::String
    is_sparse::Bool
    save_binary::Bool
    categorical_feature::Vector{Int}

    is_unbalance::Bool

    drop_rate::Float64
    max_drop::Int
    skip_drop:: Float64
    xgboost_dart_mode::Bool
    uniform_drop::Bool
    drop_seed::Int
    top_rate::Float64
    other_rate::Float64

    metric::Vector{String}
    metric_freq::Int
    is_training_metric::Bool
    ndcg_at::Vector{Int}

    num_machines::Int
    local_listen_port::Int
    time_out::Int
    machine_list_file::String

    num_class::Int
    device_type::String
end

"""
    LGBMRegression(; [
        objective = "regression",
        boosting = "gbdt",
        num_iterations = 10,
        learning_rate = .1,
        num_leaves = 127,
        max_depth = -1,
        tree_learner = \"serial\",
        num_threads = Sys.CPU_THREADS,
        histogram_pool_size = -1.,
        min_data_in_leaf = 100,
        min_sum_hessian_in_leaf = 10.,
        lambda_l1 = 0.,
        lambda_l2 = 0.,
        min_gain_to_split = 0.,
        feature_fraction = 1.,
        feature_fraction_seed = 2,
        bagging_fraction = 1.,
        bagging_freq = 0,
        bagging_seed = 3,
        early_stopping_round = 0,
        max_bin = 255,
        data_random_seed = 1,
        init_score = \"\",
        is_sparse = true,
        save_binary = false,
        categorical_feature = Int[],
        is_unbalance = false,
        drop_rate = 0.1,
        max_drop = 50,
        skip_drop = 0.5,
        xgboost_dart_mode = false,
        uniform_drop = false,
        drop_seed = 4,
        top_rate = 0.2,
        other_rate = 0.1,
        metric = [\"l2\"],
        metric_freq = 1,
        is_training_metric = false,
        ndcg_at = Int[],
        num_machines = 1,
        local_listen_port = 12400,
        time_out = 120,
        machine_list_file = \"\",
        device_type=\"cpu\",
    ])

Return a LGBMRegression estimator.
"""
function LGBMRegression(;
    objective = "regression",
    boosting = "gbdt",
    num_iterations = 10,
    learning_rate = .1,
    num_leaves = 127,
    max_depth = -1,
    tree_learner = "serial",
    num_threads = Sys.CPU_THREADS,
    histogram_pool_size = -1.,
    min_data_in_leaf = 100,
    min_sum_hessian_in_leaf = 10.,
    lambda_l1 = 0.,
    lambda_l2 = 0.,
    min_gain_to_split = 0.,
    feature_fraction = 1.,
    feature_fraction_seed = 2,
    bagging_fraction = 1.,
    bagging_freq = 0,
    bagging_seed = 3,
    early_stopping_round = 0,
    max_bin = 255,
    data_random_seed = 1,
    init_score = "",
    is_sparse = true,
    save_binary = false,
    categorical_feature = Int[],
    is_unbalance = false,
    drop_rate = 0.1,
    max_drop = 50,
    skip_drop = 0.5,
    xgboost_dart_mode = false,
    uniform_drop = false,
    drop_seed = 4,
    top_rate = 0.2,
    other_rate = 0.1,
    metric = ["l2"],
    metric_freq = 1,
    is_training_metric = false,
    ndcg_at = Int[],
    num_machines = 1,
    local_listen_port = 12400,
    time_out = 120,
    machine_list_file = "",
    device_type="cpu",
)

    return LGBMRegression(
        Booster(), "", objective, boosting, num_iterations, learning_rate, num_leaves,
        max_depth, tree_learner, num_threads, histogram_pool_size,
        min_data_in_leaf, min_sum_hessian_in_leaf, lambda_l1, lambda_l2,
        min_gain_to_split, feature_fraction, feature_fraction_seed,
        bagging_fraction, bagging_freq, bagging_seed, early_stopping_round,
        max_bin, data_random_seed, init_score,
        is_sparse, save_binary, categorical_feature,
        is_unbalance, drop_rate, max_drop, skip_drop, xgboost_dart_mode,
        uniform_drop, drop_seed, top_rate, other_rate,
        metric, metric_freq,
        is_training_metric, ndcg_at, num_machines, local_listen_port, time_out,
        machine_list_file, 1, device_type
    )
end

mutable struct LGBMClassification <: LGBMEstimator
    booster::Booster
    model::String
    application::String
    boosting :: String

    num_iterations::Int
    learning_rate::Float64
    num_leaves::Int
    max_depth::Int
    tree_learner::String
    num_threads::Int
    histogram_pool_size::Float64

    min_data_in_leaf::Int
    min_sum_hessian_in_leaf::Float64
    lambda_l1::Float64
    lambda_l2::Float64
    min_gain_to_split::Float64
    feature_fraction::Float64
    feature_fraction_seed::Int
    bagging_fraction::Float64
    bagging_freq::Int
    bagging_seed::Int
    early_stopping_round::Int

    max_bin::Int
    data_random_seed::Int
    init_score::String
    is_sparse::Bool
    save_binary::Bool
    categorical_feature::Vector{Int}

    is_unbalance::Bool

    drop_rate::Float64
    max_drop::Int
    skip_drop:: Float64
    xgboost_dart_mode::Bool
    uniform_drop::Bool
    drop_seed::Int
    top_rate::Float64
    other_rate::Float64

    metric::Vector{String}
    metric_freq::Int
    is_training_metric::Bool
    ndcg_at::Vector{Int}

    num_machines::Int
    local_listen_port::Int
    time_out::Int
    machine_list_file::String

    num_class::Int

    device_type::String
end

"""
    LGBMClassification(;[
        objective = "multiclass",
        boosting = "gbdt",
        num_iterations = 10,
        learning_rate = .1,
        num_leaves = 127,
        max_depth = -1,
        tree_learner = \"serial\",
        num_threads = Sys.CPU_THREADS,
        histogram_pool_size = -1.,
        min_data_in_leaf = 100,
        min_sum_hessian_in_leaf = 10.,
        lambda_l1 = 0.,
        lambda_l2 = 0.,
        min_gain_to_split = 0.,
        feature_fraction = 1.,
        feature_fraction_seed = 2,
        bagging_fraction = 1.,
        bagging_freq = 0,
        bagging_seed = 3,
        early_stopping_round = 0,
        max_bin = 255,
        data_random_seed = 1,
        init_score = \"\",
        is_sparse = true,
        save_binary = false,
        categorical_feature = Int[],
        is_unbalance = false,
        drop_rate = 0.1,
        max_drop = 50,
        skip_drop = 0.5,
        xgboost_dart_mode = false,
        uniform_drop = false,
        drop_seed = 4,
        top_rate = 0.2,
        other_rate = 0.1,
        metric = [\"multi_logloss\"],
        metric_freq = 1,
        is_training_metric = false,
        ndcg_at = Int[],
        num_machines = 1,
        local_listen_port = 12400,
        time_out = 120,
        machine_list_file = \"\",
        num_class = 1,
        device_type=\"cpu\",
    ])

Return a LGBMClassification estimator.
"""
function LGBMClassification(;
    objective = "multiclass",
    boosting = "gbdt",
    num_iterations = 10,
    learning_rate = .1,
    num_leaves = 127,
    max_depth = -1,
    tree_learner = "serial",
    num_threads = Sys.CPU_THREADS,
    histogram_pool_size = -1.,
    min_data_in_leaf = 100,
    min_sum_hessian_in_leaf = 10.,
    lambda_l1 = 0.,
    lambda_l2 = 0.,
    min_gain_to_split = 0.,
    feature_fraction = 1.,
    feature_fraction_seed = 2,
    bagging_fraction = 1.,
    bagging_freq = 0,
    bagging_seed = 3,
    early_stopping_round = 0,
    max_bin = 255,
    data_random_seed = 1,
    init_score = "",
    is_sparse = true,
    save_binary = false,
    categorical_feature = Int[],
    is_unbalance = false,
    drop_rate = 0.1,
    max_drop = 50,
    skip_drop = 0.5,
    xgboost_dart_mode = false,
    uniform_drop = false,
    drop_seed = 4,
    top_rate = 0.2,
    other_rate = 0.1,
    metric = ["None"],
    metric_freq = 1,
    is_training_metric = false,
    ndcg_at = Int[],
    num_machines = 1,
    local_listen_port = 12400,
    time_out = 120,
    machine_list_file = "",
    num_class = 2,
    device_type="cpu",
)

    return LGBMClassification(
        Booster(), "", objective, boosting, num_iterations, learning_rate,
        num_leaves, max_depth, tree_learner, num_threads, histogram_pool_size,
        min_data_in_leaf, min_sum_hessian_in_leaf, lambda_l1, lambda_l2,
        min_gain_to_split, feature_fraction, feature_fraction_seed,
        bagging_fraction, bagging_freq, bagging_seed, early_stopping_round,
        max_bin, data_random_seed, init_score,
        is_sparse, save_binary, categorical_feature,
        is_unbalance, drop_rate, max_drop, skip_drop, xgboost_dart_mode,
        uniform_drop, drop_seed, top_rate, other_rate, metric, metric_freq,
        is_training_metric, ndcg_at, num_machines, local_listen_port, time_out,
        machine_list_file, num_class, device_type,
    )
end


text:
export exitwithteststatus, asfloat16, asfloat32, asfloat64, asint, histdict, asdict
export tryasfloat32, tryasfloat64, tryasint, fromimage
export serve

function exitwithteststatus()
    s = FactCheck.getstats()
    if s["nNonSuccessful"] == 0
        print("   PASSED!")
    else
        print("   FAILED: $(s["nFailures"]) failures and $(s["nErrors"]) errors")
    end
    println("   ( $(s["nNonSuccessful"]+s["nSuccesses"]) tests for runtests.jl $(join(ARGS, " ")))")
    exit(s["nNonSuccessful"])
end

asint(a) = round(Int, a)
asint(a::AbstractString) = parse(Int, a)
asfloat16(a) = map(Float16, a)
asfloat32(a) = map(Float32, a)
asfloat32(a::AbstractString) = parse(Float32, a)
if isinstalled("Images")
    function asfloat32{T<:AbstractImage}(a::T)
        a = raw(a)
        if ndims(a) == 3
            a = permutedims(a,[3,2,1])
            a = a[:,:,1:min(sizeo(a),3)]
        else
            a = a'
        end
        @p asfloat32 a | divby 255 | clamp 0 1
    end
end

asfloat64(a) = map(Float64, a)
asfloat64(a::AbstractString) = parse(Float64, a)

tryasint(a, d = a) = try asint(a) catch return d end
tryasfloat32(a, d = a) = try asfloat32(a) catch return d end
tryasfloat64(a, d = a) = try asfloat64(a) catch return d end

histdict(a, field) = @p extract a field | histdict
function histdict(a)
    d = Dict()
    @p map a x->d[x] = get(d,x,0) + 1
    d
end

asdict(a) = @p fieldnames a | map (x->Pair(x,getfield(a,x))) | Dict
asdict{T<:Union{Dict,AbstractString,Array}}(a::T) = a

function serve(basepath::AbstractString, args...;kargs...)
    basepath = abspath(basepath)
    serve(args...; kargs...) do req, res
        filename = @p concat basepath req.resource | normpath
        @show filename
        if startswith(filename, basepath) && existsfile(filename) && !isdir(filename)
            HttpServer.FileResponse(filename)
        else
            Response(403)
        end
    end
end

function serve(f::Function, portrange = 8080:8199; ngrok = false)
    eval(:(using HttpServer))
    server = Server(HttpHandler(f))
    port, tcp = listenany(fst(portrange))
    close(tcp)
    @async run(server, port)
    if ngrok
        (stdout,stdin,process) = readandwrite(`ngrok http $port -log stdout --log-level "debug"`);
        a = eachline(stdout)
        b = []
        @async for x in a
            if contains(x,"Hostname:")
                push!(b,x)
                break
            end
        end
        for x in 1:50
            isempty(b) || break
            sleep(0.1)
        end
        return @p fst b | split "Hostname:" | last | split ".ngrok.io" | fst | concat ".ngrok.io"
    end
    Int(port)
end


text:
# This checks that the library has been successfully linked against an executable.
@testset "test-libgrasp" begin
    test_libgrasp_bin = joinpath(@__DIR__, "test-libgrasp")
    @test isfile(test_libgrasp_bin)

    test_libgrasp_arg = joinpath(@__DIR__, "../grasp/mixing/rmix.out")
    @test isfile(test_libgrasp_arg)
    @test success(`$(test_libgrasp_bin) $(test_libgrasp_arg)`)
end

@testset "test-readrwfn" begin
    test_readrwfn_bin = joinpath(@__DIR__, "test-readrwfn")
    @test isfile(test_readrwfn_bin)

    test_readrwfn_arg = joinpath(@__DIR__, "../grasp/mixing/rwfn.inp")
    @test isfile(test_readrwfn_arg)
    @test success(`$(test_readrwfn_bin) $(test_readrwfn_arg)`)

    test_readrwfn_arg = joinpath(@__DIR__, "../grasp/mixing/rwfn.out")
    @test isfile(test_readrwfn_arg)
    @test success(`$(test_readrwfn_bin) $(test_readrwfn_arg)`)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angles to quaternion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_quat.jl
# ========================================

# Functions: angle_to_quat
# ------------------------

@testset "Euler angles => Quaternion (Float64)" begin
    T = Float64

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to quaternion.
        q = angle_to_quat(ea)
        @test eltype(q) === T

        # Create the quaternion by composing the rotations.
        rot_seq_str = string(rot_seq)

        q₁ = if rot_seq_str[1] == 'X'
            Quaternion(cos(ea.a1 / 2), T[1, 0, 0] * sin(ea.a1 / 2))
        elseif rot_seq_str[1] == 'Y'
            Quaternion(cos(ea.a1 / 2), T[0, 1, 0] * sin(ea.a1 / 2))
        elseif rot_seq_str[1] == 'Z'
            Quaternion(cos(ea.a1 / 2), T[0, 0, 1] * sin(ea.a1 / 2))
        end

        q₂ = if rot_seq_str[2] == 'X'
            Quaternion(cos(ea.a2 / 2), T[1, 0, 0] * sin(ea.a2 / 2))
        elseif rot_seq_str[2] == 'Y'
            Quaternion(cos(ea.a2 / 2), T[0, 1, 0] * sin(ea.a2 / 2))
        elseif rot_seq_str[2] == 'Z'
            Quaternion(cos(ea.a2 / 2), T[0, 0, 1] * sin(ea.a2 / 2))
        end

        q₃ = if rot_seq_str[3] == 'X'
            Quaternion(cos(ea.a3 / 2), T[1, 0, 0] * sin(ea.a3 / 2))
        elseif rot_seq_str[3] == 'Y'
            Quaternion(cos(ea.a3 / 2), T[0, 1, 0] * sin(ea.a3 / 2))
        elseif rot_seq_str[3] == 'Z'
            Quaternion(cos(ea.a3 / 2), T[0, 0, 1] * sin(ea.a3 / 2))
        end
        qe = q₁ * q₂ * q₃

        # Make sure q0 is positive.
        (qe.q0 < 0) && (qe = -qe)

        # Compare.
        @test q ≈ qe
    end
end

@testset "Euler angles => Quaternion (Float32)" begin
    T = Float32

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to quaternion.
        q = angle_to_quat(ea)
        @test eltype(q) === T

        # Create the quaternion by composing the rotations.
        rot_seq_str = string(rot_seq)

        q₁ = if rot_seq_str[1] == 'X'
            Quaternion(cos(ea.a1 / 2), T[1, 0, 0] * sin(ea.a1 / 2))
        elseif rot_seq_str[1] == 'Y'
            Quaternion(cos(ea.a1 / 2), T[0, 1, 0] * sin(ea.a1 / 2))
        elseif rot_seq_str[1] == 'Z'
            Quaternion(cos(ea.a1 / 2), T[0, 0, 1] * sin(ea.a1 / 2))
        end

        q₂ = if rot_seq_str[2] == 'X'
            Quaternion(cos(ea.a2 / 2), T[1, 0, 0] * sin(ea.a2 / 2))
        elseif rot_seq_str[2] == 'Y'
            Quaternion(cos(ea.a2 / 2), T[0, 1, 0] * sin(ea.a2 / 2))
        elseif rot_seq_str[2] == 'Z'
            Quaternion(cos(ea.a2 / 2), T[0, 0, 1] * sin(ea.a2 / 2))
        end

        q₃ = if rot_seq_str[3] == 'X'
            Quaternion(cos(ea.a3 / 2), T[1, 0, 0] * sin(ea.a3 / 2))
        elseif rot_seq_str[3] == 'Y'
            Quaternion(cos(ea.a3 / 2), T[0, 1, 0] * sin(ea.a3 / 2))
        elseif rot_seq_str[3] == 'Z'
            Quaternion(cos(ea.a3 / 2), T[0, 0, 1] * sin(ea.a3 / 2))
        end
        qe = q₁ * q₂ * q₃

        # Make sure q0 is positive.
        (qe.q0 < 0) && (qe = -qe)

        # Compare.
        @test q ≈ qe
    end
end

@testset "Euler angles => Quaternion (Errors)" begin
    @test_throws ArgumentError angle_to_quat(1, 2, 3, :ZZY)
end

# Functions: smallangle_to_quat
# -----------------------------

@testset "Small Euler angles => Quaternion (Float64)" begin
    T = Float64

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Create the quaternion.
    q = smallangle_to_quat(θx, θy, θz)

    # Expected result.
    qe = Quaternion(1, θx / 2, θy / 2, θz / 2)
    qe = qe / norm(qe)

    @test qe ≈ q
end

@testset "Small Euler angles => Quaternion (Float32)" begin
    T = Float32

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Create the quaternion.
    q = smallangle_to_quat(θx, θy, θz)

    # Expected result.
    qe = Quaternion(1, θx / 2, θy / 2, θz / 2)
    qe = qe / norm(qe)

    @test qe ≈ q
end


text:
using Revise

using Tiger
using JuMP

# =================== TESTING ===================

cobra = read_cobra("/Users/jensen/Dropbox/repos/COBRA_models/iSMU.mat", "iSMU")
cobra.lb[cobra.lb .> 0.0] .= 0.0  # remove NGAM
model = build_base_model(cobra)
optimize!(model)

set_media_bounds!(model, "./test/CDM.toml")
optimize!(model)

add_gprs_cnf!(model, cobra)
optimize!(model)

#set_silent(model)
#@time single_deletions(model, variable_by_name.(model, cobra.genes))

ext = extend_cobra_cnf(cobra)
m_ext = build_base_model(ext)
set_media_bounds!(m_ext, "./test/CDM.toml")
optimize!(m_ext)


text:
using StableRNGs

"""
Function to initialise the parameters of the desired network.
"""
function initialise_model(layer_dims::Vector{Int}, seed::Any)
    params = Dict()

    for l = 2:length(layer_dims)
        params[string("W_", (l-1))] = ( rand( StableRNG(seed), layer_dims[l], layer_dims[l-1] ) .- 0.5) * 200#* sqrt(layer_dims[l-1])
        params[string("b_", (l-1))] = zeros(layer_dims[l], 1)
    end # for

    return params
end

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
# 5 cubes example --- Example 2.4.4 (Solid algebra 3D)
# Original generation data were random and no seed was previously set
V = [
    1.01181 0.215639 0.91979 0.123616 1.02252 0.226347 0.930498 0.134324 0.0458309 -0.301827 0.348275 0.0006172 0.579367 0.23171 0.881811 0.534154 -0.0521776 0.627953 -0.190635 0.489496 -0.0233884 0.656742 -0.161846 0.518285 0.27652 -0.0875132 0.52527 0.161237 0.509324 0.145291 0.758074 0.394041 0.27631 0.564484 0.0058279 0.294002 1.01221 1.30039 0.741732 1.02991;
    0.160033 0.0680099 0.956278 0.864255 0.160649 0.0686266 0.956895 0.864872 -0.200245 0.102199 0.417839 0.720283 -0.35354 -0.0510965 0.264543 0.566987 0.682359 0.543901 0.0592036 -0.0792537 0.956374 0.817917 0.333219 0.194762 -0.102028 0.146722 0.324834 0.573584 -0.16916 0.0795901 0.257702 0.506452 0.452356 0.181874 1.15396 0.883481 0.816151 0.545669 1.51776 1.24728;
    0.196256 0.206963 0.196872 0.20758 0.997729 1.00844 0.998346 1.00905 0.0677451 0.601282 -0.0855504 0.447986 0.502301 1.03584 0.349005 0.882542 0.159301 0.18809 0.433316 0.462105 0.797002 0.825792 1.07102 1.09981 0.1446 0.377404 0.0774682 0.310272 0.580364 0.813168 0.513232 0.746036 0.403805 1.13971 0.767599 1.5035 0.249344 0.985248 0.613139 1.34904
];

EV = Array{Int64,1}[
    [ 1, 2], [ 3, 4], [ 5, 6], [ 7, 8], [ 1, 3], [ 2, 4],
    [ 5, 7], [ 6, 8], [ 1, 5], [ 2, 6], [ 3, 7], [ 4, 8],
    [ 9,10], [11,12], [13,14], [15,16], [ 9,11], [10,12],
    [13,15], [14,16], [ 9,13], [10,14], [11,15], [12,16],
    [17,18], [19,20], [21,22], [23,24], [17,19], [18,20],
    [21,23], [22,24], [17,21], [18,22], [19,23], [20,24],
    [25,26], [27,28], [29,30], [31,32], [25,27], [26,28],
    [29,31], [30,32], [25,29], [26,30], [27,31], [28,32],
    [33,34], [35,36], [37,38], [39,40], [33,35], [34,36],
    [37,39], [38,40], [33,37], [34,38], [35,39], [36,40]
];

FV = Array{Int64,1}[
    [ 1, 2, 3, 4], [ 5, 6, 7, 8], [ 1, 2, 5, 6], [ 3, 4, 7, 8], [ 1, 3, 5, 7], [ 2, 4, 6, 8],
    [ 9,10,11,12], [13,14,15,16], [ 9,10,13,14], [11,12,15,16], [ 9,11,13,15], [10,12,14,16],
    [17,18,19,20], [21,22,23,24], [17,18,21,22], [19,20,23,24], [17,19,21,23], [18,20,22,24],
    [25,26,27,28], [29,30,31,32], [25,26,29,30], [27,28,31,32], [25,27,29,31], [26,28,30,32],
    [33,34,35,36], [37,38,39,40], [33,34,37,38], [35,36,39,40], [33,35,37,39], [34,36,38,40]
];

cop_EV = Lar.coboundary_0(EV);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV, EV); W = convert(Lar.Points, V');
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(W, cop_EW, cop_FE)

V,CVs,FVs,EVs = Lar.pols2tria(convert(Lar.Points, V'), copEV, copFE, copCF);

if todisplay
    GL.VIEW(GL.GLExplode(V, FVs,  5,5,5,  99,1));
    GL.VIEW(GL.GLExplode(V, CVs[2:end],  5,5,25,  99,0.5));
end;


text:
using BinaryBuilder

name = "Triangle"
version = v"1.6.1"

# Collection of sources required to build Triangle by J. Shewchuk
#
# More info:
# https://www.cs.cmu.edu/~quake/triangle.html
# https://github.com/JuliaGeometry/Triangulate.jl
 #
# Please be aware that triunsuitable.c is not part of the original Triangle distribution.
# It provides the possibility to pass a cfunction created in Julia as  user refinement callback.
# For this reason at least triunsiutable.c must be downloaded from Triangulate.jl repo, and
# for simplicity, we do this for the whole code.
sources = [
    GitSource("https://github.com/JuliaGeometry/Triangulate.jl.git","b2ffb23ca7d89c567fd31367882bd216757cdb9c")
]

script = raw"""
cd $WORKSPACE/srcdir/Triangulate.jl/deps/src
mkdir -p "${libdir}"

# Patch the code in order to replace exit() calls in the C code  so that we can
# throw an error instead.
sed -e "s/  exit/extern void error_exit(int); error_exit/g" triangle/triangle.c > triangle_patched.c

# Concerning the suppression of int - pointer cast warnings,
# see the following comment in triangle.c:
#
# encode() compresses an oriented triangle into a single pointer.  It
# relies on the assumption that all triangles are aligned to four-byte
# boundaries, so the two least significant bits of (otri).tri are zero.
#
# So there is an inherent hack in the code which we have to allow and hope things keep working.
$CC -Itriangle  -Wno-int-to-pointer-cast  -Wno-pointer-to-int-cast -DREAL=double -DTRILIBRARY -O3 -fPIC -DNDEBUG -DNO_TIMER -DEXTERNAL_TEST $LDFLAGS --shared -o "${libdir}/libtriangle.${dlext}" triangle_patched.c triwrap.c

install_license triangle/README
"""

platforms = supported_platforms(; experimental=true)
products = [LibraryProduct("libtriangle", :libtriangle)]
dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Base.isless
import PyPlot
#Allow sorting between Numbers and Strings
isless(A::Number,B::AbstractString) = false
isless(A::AbstractString,B::Number) = true

StringOrNumArray = Union{AbstractString,Array,Number}

#Creation of Assoc require StrUnique to split Single-Character-separated String Sequence.
include("StrUnique.jl")

#=
Type Assoc (Associative Array)
Support a
=#

type Assoc
    row::Array{Union{AbstractString,Number}}
    col::Array{Union{AbstractString,Number}}
    val::Array{Union{AbstractString,Number}}
    A::AbstractSparseMatrix

    Assoc(rowIn::StringOrNumArray,colIn::StringOrNumArray,valIn::StringOrNumArray) = Assoc(rowIn,colIn,valIn,min)
    Assoc(row::Array{Int64}, col::Array{Int64},val::Array{Int64},A::AbstractSparseMatrix) = new(row,col,val,A)
    #Setting Default Function

    function Assoc(rowIn::Array{Union{AbstractString,Number}}, colIn::Array{Union{AbstractString,Number}}, valIn::Array{Union{AbstractString,Number}}, AIn::AbstractSparseMatrix)
        return new(rowIn,colIn,valIn,AIn)
        end

    function Assoc(rowIn::StringOrNumArray,colIn::StringOrNumArray,valIn::StringOrNumArray,funcIn::Function)
        if isempty(rowIn) || isempty(colIn) || isempty(valIn)  #testing needed for isemtpy, for Matlab isemtpy is always possible TODO  Seems to work okay with String or NumArray type hard defined, Union type untested.  Should keep an eye.
            x = Array{Union{AbstractString,Number}}()
            return Assoc(x,x,x,spzeros(1,1));
            end
        if isa(rowIn,Number)
            rowIn = Array{Union{AbstractString,Number},1}([rowIn])
        end

        if isa(colIn,Number)

            colIn = Array{Union{AbstractString,Number},1}([colIn])
        end

        if isa(valIn,Number)
            valIn = Array{Union{AbstractString,Number},1}([valIn])
        end
        i = rowIn;
        j = colIn;
        v = valIn;
        row = rowIn;
        col = colIn;
        val = valIn;

        if isa(rowIn,AbstractString)
            row, i_out2in, i = StrUnique(rowIn);
            else
            row = unique(i)
            sort!(row)
            i = AbstractArray([Int64(searchsortedfirst(row,x)) for x in i])
            i = convert(AbstractArray{Int64},i)

            end

        if isa(colIn,AbstractString)
            col, j_out2in, j = StrUnique(colIn);
            else
            col = unique(j)
            sort!(col)
            j = [searchsortedfirst(col,x) for x in j]
            j = convert(AbstractArray{Int64},j)

            end

        if isa(valIn,AbstractString)
            val, v_out2in, v = StrUnique(valIn);
            else
            val = unique(v)
            sort!(val)
            if (isa(valIn[1],AbstractString))
                k = 0
                if val[1] == ""
                    k =-1
                end
                v = [searchsortedfirst(val,x)+ k  for x in v]
                v = convert(AbstractArray{Int64},v)
                if val[1] == ""
                    val = val[2:end]
                end
            end
           end

        NMax = maximum([length(i) length(j) length(v)]);
        if length(i) == 1
            x = i[1]
            i = Array{typeof(x)}(NMax)
            for n = 1:NMax
                i[n] = x
            end
            end

        if length(j) == 1
            x = j[1]
            j = Array{typeof(x)}(NMax)
            for n = 1:NMax
                j[n] = x
            end
            end

        if length(v) == 1
            x = v[1]
            v = Array{typeof(x)}(NMax)
            for n = 1:NMax
                v[n] = x
            end
            end

        #i = convert(AbstractArray{Int64},i)
        if isa(val[1],AbstractString) #If the values are string, assume that there are duplicates and take the earliest one ( the numbers should be the same)
            A = sparse(i,j,v,length(row),length(col),min);
        else
            A = sparse(i,j,v,length(row),length(col),(+));
        end
        #Accumarray isn't in Julia, use "push" for a more rapid array generation (acccumarray is too slow and cumbersome for Julia)  Sparse matrix generation condition with summation combine seem would do the trick.

        #End bit with val string.  Unknown purpose.
        return new(row,col,val,A)
        end
    end

#=
Adding related operations for Assoc_orig
=#
include("./Assoc_orig/getindex.jl")
include("./Assoc_orig/condense.jl")
include("./Assoc_orig/no.jl")
include("./Assoc_orig/sum.jl")
include("./Assoc_orig/isempty.jl")
include("./Assoc_orig/logical.jl")
include("./Assoc_orig/and.jl")
include("./Assoc_orig/print.jl")
include("./Assoc_orig/transpose.jl")
include("./Assoc_orig/multiply.jl")
include("./Assoc_orig/rdivide.jl")
include("./Assoc_orig/sqIn.jl")
include("./Assoc_orig/sqOut.jl")
include("./Assoc_orig/accessor.jl")
include("./Assoc_orig/find.jl")
include("./Assoc_orig/diag.jl")
include("./Assoc_orig/plus.jl")
include("./Assoc_orig/deepcondense.jl")
include("./Assoc_orig/lt.jl")
include("./Assoc_orig/gt.jl")
include("./Assoc_orig/spy.jl")
include("./Assoc_orig/jld.jl")
include("./Assoc_orig/put.jl")
include("./Assoc_orig/full.jl")
########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using Documenter
using RandomVariates

makedocs(
    sitename = "RandomVariates",
    format = Documenter.HTML(),
    modules = [RandomVariates],
    pages = [
        "Home" => "index.md",
        "prn.md",
        "Uniform" => "uniform.md",
        "Normal" => "normal.md",
        "Discrete" => "discrete.md",
        "Continuous" => "continuous.md"
    ]
)

deploydocs(
    repo = "github.com/chris-santiago/RandomVariates.jl.git",
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

# DocumenterTools.genkeys(user="chris-santiago",
#        repo="git@github.com:chris-santiago/RandomVariates.jl.git")


text:
using DataDeps

register(DataDep(
  "tv-sales",
  "http://www.stat.ufl.edu/~winner/data/tvsales.txt",
  "http://www.stat.ufl.edu/~winner/data/tvsales.dat",
  "6d51cfd2b0eaed03f0bcb578f8a2b972227b83c33b6334ea47dbf2d3965aa87f",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=2,
        feature_indices=3:10,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
@testset "Multiroute flow" begin
    #### Graphs for testing
    graphs = [
    # Graph with 8 vertices
    (8,
    [
    (1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
    (2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
    (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10),
    (8, 1, 8)                        # Reverse edge to test the slope in EMRF
    ],
    1, 8,                              # source/target
    [28, 28, 15, 0],                   # answer for 1 to 4 route(s) flows
    [(0., 0., 3), (5., 15., 2),        # breaking points
    (10., 25., 1), (13., 28., 0)],
    28.                                # value for 1.5 routes
    ),

    # Graph with 6 vertices
    (6,
    [
    (1, 2, 9), (1, 3, 9), (2, 3, 10), (2, 4, 8), (3, 4, 1),
    (3, 5, 3), (5, 4, 8), (4, 6, 10), (5, 6, 7)
    ],
    1, 6,                                      # source/target
    [12, 6, 0],                                # answer for 1 to 3 route(s) flows
    [(0., 0., 2), (3., 6., 1), (9., 12., 0)],  # breaking points
    9.                                         # value for 1.5 routes
    ),

    # Graph with 7 vertices
    (7,
    [
    (1, 2, 1), (1, 3, 2), (1, 4, 3), (1, 5, 4), (1, 6, 5),
    (2, 7, 1), (3, 7, 2), (4, 7, 3), (5, 7, 4), (6, 7, 5)
    ],
    1, 7,                                      # source/target
    [15, 15, 15, 12, 5, 0],                    # answer for 1 to 6 route(s) flows
    [(0., 0., 5), (1., 5., 4), (2., 9., 3),    # breaking points
    (3., 12., 2), (4., 14., 1), (5., 15., 0)],
    15.                                        # value for 1.5 routes
    ),

    # Graph with 6 vertices
    (6,
    [
    (1, 2, 1), (1, 3, 1), (1, 4, 2), (1, 5, 2),
    (2, 6, 1), (3, 6, 1), (4, 6, 2), (5, 6, 2),
    ],
    1, 6,                                      # source/target
    [6, 6, 6, 4, 0],                           # answer for 1 to 5 route(s) flows
    [(0., 0., 4), (1., 4., 2), (2., 6., 0)],   # breaking points
    6.                                         # value for 1.5 routes
    )
    ]

    for (nvertices, flow_edges, s, t, froutes, breakpts, ffloat) in graphs
        flow_graph = Graphs.DiGraph(nvertices)
        for g in testdigraphs(flow_graph)
            capacity_matrix = zeros(Int, nvertices, nvertices)
            for e in flow_edges
                u, v, f = e
                Graphs.add_edge!(g, u, v)
                capacity_matrix[u, v] = f
            end
            # Test ExtendedMultirouteFlowAlgorithm when the number of routes is either
            # Noninteger or 0 (the algorithm returns the breaking points)
            @test multiroute_flow(g, s, t, capacity_matrix) == breakpts
            @test multiroute_flow(g, s, t, capacity_matrix, routes=1.5)[1] ≈ ffloat
            @test multiroute_flow(breakpts, 1.5)[2] ≈ ffloat

            # Test all other algorithms - PR is unstable - see #553
            for AlgoFlow in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
                # When the input are breaking points and routes number
                @test multiroute_flow(breakpts, 1.5, g, s, t, capacity_matrix)[1] ≈ ffloat
                for AlgoMrf in [ExtendedMultirouteFlowAlgorithm, KishimotoAlgorithm]
                    for (k, val) in enumerate(froutes)
                        @test multiroute_flow(g, s, t, capacity_matrix,
                            flow_algorithm = AlgoFlow(), mrf_algorithm = AlgoMrf(),
                            routes = k)[1] ≈ val
                    end
                end
            end
        end
    end
end