text:
# operators imported and overridden from Base
import Base: +, -, *, /, ^, |, %, ==, !=, <, >, <<, >>, ~

# functions imported and overridden from Base
# convention: prefix `Base.` should always be present when overriding
# so comment the following will still work
import Base: show,convert, getproperty, propertynames, setproperty!
if isdefined(Base, :hasproperty) # Julia 1.2
    import Base: hasproperty
end
import Base: abs, inv, adjoint, exp, getindex

export Mv

export norm, rev, dual, involute, proj, refl, rot, exp_with_hint, scalar, even, odd
# export \cdot, \wedge, \intprod, \intprodr, \odot, \boxtimes, \circledast, \times
export +, -, *, /, ^, |, %, ==, !=, <, >, <<, >>, ~, ⋅, ∧, ⨼, ⨽, ⊙, ⊠, ⊛, ×
# Operator precedence: they have the same precedence, unlike in math
# julia> for op ∈ [:* :⋅ :∧ :⨼ :⨽ :⊙ :⊠ :⊛ :×]; println(String(op), "  ", Base.operator_precedence(op)) end
# *  13
# ⋅  13
# ∧  13
# ⨼  13
# ⨽  13
# ⊙  13
# ⊠  13
# ⊛  13
# ×  13

# experimental export \bar\times
# export ×̄

# \^-\^1
@define_postfix_symbol(⁻¹)
# # \^T
# @define_postfix_symbol(ᵀ)
# \^x
@define_postfix_symbol(ˣ)
# \doublepipe
@define_postfix_symbol(ǂ)
# \_+
@define_postfix_symbol(₊)
# \_-
@define_postfix_symbol(₋)
# \bot
# @define_postfix_symbol(⊥)

@doc raw"""
A wrapper class for `galgebra.mv.Mv`:

- all methods of `galgebra.mv.Mv` are delegated and can be called like in Python.
- enhanced with operator overriding and extra methods.
- automatically supports pretty printing like in Python.
"""
mutable struct Mv
    o::PyCall.PyObject
end

@define_show(Mv)
@delegate_properties(Mv, :o)
@delegate_doc(Mv)

"Addition."
@define_op(Mv, +, __add__)

"Subtraction."
@define_op(Mv, -, __sub__)

@doc raw"""
Geometric product.

``A * B \equiv A B``.
"""
@define_op(Mv, *, __mul__)

@doc raw"""
Division.

``A / B \equiv A B^{-1}``. Only valid when ``B`` has inverse.
"""
@define_op(Mv, /, __truediv__)

@doc raw"""
Comparisons of equality.

``A = B \equiv \mathrm{simplify}(A - B) = 0``
"""
@define_op(Mv, ==, __eq__)

@doc raw"""
Comparisons of inequality.

Hint: Type ≠ with `\neq`. Alternatively, use `!=`.

``A \neq B \equiv \mathrm{simplify}(A - B) \neq 0``
"""
@define_op(Mv, ≠, __ne__)
@define_op(Mv, !=, __ne__)

@doc raw"""
Wedge product.

Hint: type ∧ with `\wedge`.
"""
@define_op(Mv, ∧, __xor__)

@doc raw"""
Hestenes' inner product.

Hint: type ⋅ with `\cdot`. Alternatively, use `|`.

``\boldsymbol{A}_{r} \cdot \boldsymbol{B}_{s} \equiv \left\{\begin{array}{lr}{r \text { and } s \neq 0 :} & {\left\langle\boldsymbol{A}_{r} \boldsymbol{B}_{s}\right\rangle_{|r-s|}} \\ {r \text { or } s=0 :} & {0}\end{array}\right.``
"""
@define_op(Mv, ⋅, __or__)
@define_op(Mv, |, __or__)

@doc raw"""Left contraction, i.e. "contraction onto".

Hint: type ⨼ with `\intprod`. Alternatively, use `<`.

``A \rfloor B \equiv \sum\limits_{r, s}\left\langle\langle A\rangle_{r}\langle B\rangle_{s}\right\rangle_{s-r}``

In literature the notation is usually ``A \rfloor B``, but `\rfloor` is reserved by Julia.
"""
@define_op(Mv, ⨼, __lt__)
@define_op(Mv, <, __lt__)

@doc raw"""Right contraction, i.e. "contraction by".

Hint: type ⨽ with `\intprodr`. Alternatively, use `>`.

``A \lfloor B \equiv \sum\limits_{r, s}\left\langle\langle A\rangle_{r}\langle B\rangle_{s}\right\rangle_{r-s}``

In literature the notation is usually ``A \lfloor B``, but `\lfloor` is reserved by Julia.
"""
@define_op(Mv, ⨽, __gt__)
@define_op(Mv, >, __gt__)

@doc raw"""Commutator product.

Hint: type ⊠ with `\boxtimes`. Alternatively, use `>>`.

``A \underline{\times} B \equiv \dfrac{1}{2}(AB-BA)``.
"""
@define_op(Mv, ⊠, __rshift__)
@define_op(Mv, >>, __rshift__)

@doc raw"""Anti-commutator product.

Hint: type ⊙ with `\odot`. Alternatively, use `<<`.

``A \bar{\times} B \equiv \dfrac{1}{2}(AB+BA)``.
"""
@define_op(Mv, ⊙, __lshift__)
@define_op(Mv, <<, __lshift__)

# # experimental symbol for anti-comutator product: \bar\times
# A×̄B = (AB+BA)/2
# @define_op(Mv, ×̄, __lshift__)

@doc raw"""
Cross product for vectors in 3D.
"""
@define_op_with_impl(Mv, ×, mv.cross(x, y))

@doc raw"""Scalar product.

Hint: type ⊛ with `\circledast`. Alternatively, use `%`.

``A \circledast B \equiv \langle A B^{\dagger} \rangle``.

In literature the notation is usually ``\ast`` , but it's visually indistinguishable from `*`.
"""
@define_op_with_impl(Mv, ⊛, (x * ~y).scalar())
@define_op_with_impl(Mv, %, x ⊛ y)

"Unary negation."
@define_unary_op(Mv, -, __neg__)

@doc raw"""
Norm.

`norm(A)` = `A.norm()` ``\equiv \left\lVert A \right\rVert \equiv \sqrt{A \tilde{A}}``

Alternatively:

- `A.norm(hint="+")` ``\equiv \sqrt{A \tilde{A}}``
- `A.norm(hint="-")` ``\equiv \sqrt{- A \tilde{A}}``
- `A.norm(hint="0")` ``\equiv \sqrt{\left| A \tilde{A} \right|}``

Only valid when the result is a scalar.
"""
@define_unary_op(Mv, norm, norm)
@define_unary_op(Mv, Base.abs, norm)

@doc raw"""
Inverse.

`(A)⁻¹ = A^-1 = inv(A) = A.inv()` ``\equiv A^{-1}``

Hint: type ⁻¹ with `\^-\^1`.
"""
@define_unary_op(Mv, Base.inv, inv)
@define_postfix_op(Mv, ⁻¹, Base.inv)

@doc raw"""
Reversion.

`~A = A[:~] = rev(A) = A.rev()` ``\equiv \tilde{A} \equiv A^{\dagger}``

In literature the notation is usually ``\tilde{A}`` or ``A^{\dagger}``, the former is illegal syntax and `\dagger` in the latter is is reserved by Julia.
"""
@define_unary_op(Mv, ~, rev)
@define_unary_op(Mv, rev, rev)

# @deprecated
# @define_postfix_op(Mv, ᵀ, rev)

@doc raw"""
Dual, i.e. orthogonal complement, ``\Lambda^p \to \Lambda^{n-p}``.

`A'` ``\equiv A^{\bot} \equiv A I``

Note: call `Ga.dual_mode(mode)` to globally specify a different dual mode (`I+` is the default):

| dual_mode | ``A^{\bot}`` |
|-----------|--------------|
|   `+I`    |     ``IA``   |
|   `-I`    |    ``-IA``   |
|   `I+`    |     ``AI``   |
|   `I-`    |    ``-AI``   |
|  `+Iinv`  |  ``I^{-1}A`` |
|  `-Iinv`  | ``-I^{-1}A`` |
|  `Iinv+`  |  ``AI^{-1}`` |
|  `Iinv-`  | ``-AI^{-1}`` |
"""
@define_unary_op(Mv, Base.adjoint, dual)
@define_unary_op(Mv, dual, dual)
# @define_postfix_op(Mv, ⊥, dual)

@doc raw"""
Grade involution.

`(A)ˣ = A[:*] = involute(A)` ``\equiv A_+ - A_- \equiv`` `A.even() - A.odd()`

Hint: type ˣ with `\^x`.

In literature the notation is usually ``A^{*}``.
"""
@define_unary_op_with_impl(Mv, involute, x.even() - x.odd())
@define_postfix_op(Mv, ˣ, involute)

@doc raw"""
Clifford conjugate.

`(A)ǂ = A[:ǂ]` ``\equiv A^{*\dagger}``

Hint: type ǂ with `\doublepipe`.

In literature the notation is usually ``A^{\ddagger}``, but `\ddagger` is reserved by Julia.
"""
@define_unary_op_with_impl(Mv, Base.conj, involute(x).rev())
@define_postfix_op(Mv, ǂ, Base.conj)

@doc raw"""
Projection.

`proj(B, A)` `` \equiv P_{B}(A) \equiv`` `A.project_in_blade(B)`

Only valid if B is a blade.
"""
@define_op_with_impl(Mv, proj, mv.proj(x, y))

@doc raw"""
Reflection.

`refl(B, A)` `` \equiv \mathrm{Refl}_{B}(A) \equiv`` `A.reflect_in_blade(B)`

Only valid if B is a blade.
"""
@define_op_with_impl(Mv, refl, mv.refl(x, y))

@doc raw"""
Rotation.

Rotate the multivector `A` by the 2-blade `itheta`.

`rot(itheta, A)` ``\equiv A e^{I \theta} \equiv`` `A.rotate_multivector(itheta)`
"""
@define_op_with_impl(Mv, rot, rot_with_hint(x, y))
@pure rot_with_hint(itheta::Mv, A::Mv, hint::AbstractString="-") = mv.rot(itheta, A, hint)

@doc raw"Natural base exponential of X: ``e^X``"
@define_unary_op_with_impl(Mv, Base.exp, exp_with_hint(x))
@pure exp_with_hint(x::Mv, hint::AbstractString="-") = mv.exp(x, hint)

@doc raw"""
The `i`-th grade part.

`A[i] = A.grade(i)` ``\equiv \langle A B^{\dagger} \rangle_i``
"""
@pure Base.getindex(x::Mv, i::Integer) = x.grade(i)
@pure function Base.getindex(x::Mv, sym::Symbol)
    if sym == :+
        even(x)
    elseif sym == :-
        odd(x)
    elseif sym == :~
        rev(x)
    elseif sym == :⁻¹
        inv(x)
    elseif sym == :*
        involute(x)
    elseif sym == :ǂ
        conj(x)
    else
        throw(DomainError(x, "argument can only be one of :+, :-, :~, :⁻¹, :*, :ǂ"))
    end
end

@doc raw"""
Scalar (grade-0) part.

`scalar(A) = A.scalar()` ``\equiv \langle A B^{\dagger} \rangle \equiv \langle A B^{\dagger} \rangle_0``

Note: it returns a SymPy expression unlike A[0] which returns a Mv object
"""
@define_unary_op(Mv, scalar, scalar)

@doc raw"""
Even-grade part.

`A[:+] = (A)₊ = even(A) = A.even()` ``\equiv A_+``
"""
@define_unary_op(Mv, even, even)
@define_postfix_op(Mv, ₊, even)

@doc raw"""
Odd-grade part.

`A[:-] = (A)₋ = odd(A) = A.odd()` ``\equiv A_-``
"""
@define_unary_op(Mv, odd, odd)
@define_postfix_op(Mv, ₋, odd)

@define_lop(Mv, Number, +, __add__)
@define_rop(Mv, Number, +, __radd__)
@define_lop(Mv, Number, -, __sub__)
@define_rop(Mv, Number, -, __rsub__)
@define_lop(Mv, Number, *, __mul__)
@define_rop(Mv, Number, *, __rmul__)
@define_lop(Mv, Number, /, __truediv__)
@define_rop(Mv, Number, /, __rdiv__)
@define_lop(Mv, Number, ==, __eq__)
@define_lop(Mv, Number, !=, __ne__)

@pure function ^(x::Mv, y::Integer)
    if y < 0
        x.__pow__(abs(y)).inv()
    elseif y == 0
        1
    else
        x.__pow__(y)
    end
end


text:
using CoupledODETools
using Test

import CoupledODETools
import CoupledODETools: flaginsert
import CoupledODETools: flagreplace

@testset "CoupledODETools.jl" begin
    #test metaprogramming functions
    @test flaginsert(:f, :(f+a+t), :t) == :(t+f+a+t)
    @test flaginsert(:f, :(f-a-t), :t) == :(t+f-a-t)
    @test flaginsert(:f, :(√(f-a*t)), :q) == :(√(q+f-a*t))
    @test flaginsert(:f, :(a*f*b), :t) == :(t*a*f*b)
    @test flaginsert(:f, :(x^f), :a) == :(x^(a+f))
    @test flaginsert(:f, :(t/f), :a) == :(t/(a*f))
    @test flaginsert(:f, :(t/(1+f)), :a) == :(t/(a+1+f))

end
flagreplace(:s, :((f+a+t)), :a) == :(f+a+t)
flagreplace(:a, :((f+a+t)), :s) == :(f+s+t)

@testset "free parameters -> SLArray" begin

    @Component FitzughCell begin
        Dv = v - v^3 - w + Iext + gsyn, .0
        Dw = v + a -b*w, .1
        to_alpha = vpre -> ~v
        elec = gsyn -> gelec*(~v - v)
    end begin
        Iext = .5
        a = :free
        b = .7
        gelec = 3
    end
    @Component AlphaSynapse begin
        DS = α*(1-S)/(1+exp(♉-vpre)*k) - β*S, 0
        synout = gsyn -> g*~S*(v-E)
    end begin
        α = :free
        β = .2
        ♉ = -20
        k = 20
        E = -5
        g = .1
    end
    #test defaults
    c = FitzughCell(:c)
    @test c.couplings[2].targets == []
    @test c.pars[1].val == .5
    @test c.eqns[1].sym == :v
    @test c.eqns[1].val == .0
    #test coupling
    c1 = FitzughCell(:c1, elec = :c2)
    c2 = FitzughCell(:c2, v0 = .1, w0 = .2, a = 3, elec = :c1, to_alpha = :s1)
    s1 = AlphaSynapse(:s1, synout = [:c2])

    @test c1.eqns[1].val == .0
    @test c2.eqns[1].val == .1
    @test c2.pars[2].val == 3

    net = Network([c1,c2,s1])()
end

ext: jl


text:
@testset "test_pid_design" begin
@test stabregionPID("exp(-sqrt(s))", exp10.(range(-3, stop=3, length=50)), kd=1)[2][1] ≈ -1.022356911142034

P = tf(1,[1,1])^4
gangoffourplot(P,tf(1))
gangoffourplot(ss(P),ss(1))
ωp = 0.8
kp,ki,C = loopshapingPI(P,ωp,phasemargin=60, doplot=true)
@test kp ≈ 0.82274734724854
@test ki ≈ -0.45472580052281375
end


text:
# code for taking two tables of results and doing comparisions
using CSV, DataFrames,  PyPlot

function df_append_columns!(df::DataFrame,append::String;ignore=Array{Symbol,1}())
    for old_name in names(df)
        new_name = String(old_name) * append
        if !(old_name in ignore)
            rename!(df,old_name => Symbol(new_name))
        end
    end
end

function status_table(df_combine)
    # create a CSV with the statuses of each solver
end

function ratio_to_fastest(mat::Matrix)
    # loop through all columns of the dataframe and pick
    #@assert(mat .> 0.0)
    ratios = zeros(size(mat))
    for i = 1:size(mat,1)
        fastest_val = 10.0^10.0
        for j = 1:size(mat,2)
            if isnan(fastest_val)
                warn("NaN")
            else
                #@show mat[i,j]
                fastest_val = min(fastest_val,mat[i,j])
            end
        end

        for j = 1:size(mat,2)
            if mat[i,j] > 0.0 && fastest_val > 0.0
                ratios[i,j] = mat[i,j] / fastest_val
            end
        end
    end

    return ratios
end

function sort_mat!(mat::Matrix)
    for j = 1:size(mat,2)
        mat[:,j] = sort(mat[:,j])
    end
end

function plot_ratios!(ratios::Matrix;max_ratio=128.0,labels=Array{String,1}(),colors=Array{String,1}(),lines=Array{String,1}())
    sort_mat!(ratios)
    n = size(ratios,1)

    for j = 1:size(ratios,2)
        vals = ratios[:,j]
        only_plot = vals[vals .> 1.0];
        theta = (length(vals) - length(only_plot)) / length(vals)
        y = linspace(theta,1,length(only_plot))
        if length(labels) > 0
            label = labels[j]
        else
            label = ""
        end
        semilogx(only_plot,y,basex=2,label=label)
    end

    ylim([0.0,1.0])
    xlim([1.0,max_ratio])
    if length(labels) > 0
        legend()
    end
end

function geomean(vals::Vector)
    n = length(vals)
    return exp((1/n) * sum(log.(vals)))
end

function summary(df::DataFrame)
    data = Dict{Symbol,Array{Union{String,Float64}}}()
    data[:stat] = ["median","mean","geomean"];

    for h in names(df)
        data[h]=[median(df[:,h]),mean(df[:,h]),geomean(df[:,h])]
    end

    df_new = DataFrame(data)

    return df_new
end

function count(vec::Vector)
    count_dic = Dict()
    for el in vec
        if !haskey(count_dic,el)
            count_dic[el] = 1
        else
            count_dic[el] += 1
        end
    end
    return count_dic
end

function success(vec::Vector)
    return (vec .== "Optimal") .| (vec .== "dual_infeasible") .| (vec .== "primal_infeasible") .| (vec .== "Optimal_SS")
end

function compare_status(vec1::Vector,vec2::Vector)

end

test_name = "CUTEst/test-July-1"
RESULTS_DIR = "../benchmark-2/results/$test_name"

df_Ipopt = readtable("$RESULTS_DIR/Ipopt_summary.csv");
df_append_columns!(df_Ipopt,"_Ipopt",ignore=[:name])
df_OnePhase = readtable("$RESULTS_DIR/OnePhase_summary.csv");
df_append_columns!(df_OnePhase,"_OnePhase",ignore=[:name])

df_combine = join(df_Ipopt, df_OnePhase, on = :name) #, makeunique=true)
CSV.write("$RESULTS_DIR/combined.csv",df_combine)
df_combine = df_combine[df_combine[:status_Ipopt] .!= "TOO_FEW_DEGREES_OF_FREEDOM",:]
df_combine[df_combine[:status_Ipopt] .== "Optimal_SS",:status_Ipopt] = :Optimal;

count(df_combine[:status_Ipopt])
count(df_combine[:status_OnePhase])
count(df_combine[:status_Ipopt] .* df_combine[:status_OnePhase])

both_optimal_bools = (df_combine[:status_Ipopt].=="Optimal") .& (df_combine[:status_OnePhase].=="Optimal")
df_both_optimal = df_combine[both_optimal_bools,:];

fac_ratios = ratio_to_fastest([df_both_optimal[:n_fac_Ipopt] df_both_optimal[:n_fac_OnePhase]]);
plot_ratios!(fac_ratios,labels=["#fac Ipopt","#fac One Phase"])

fval_ratios = ratio_to_fastest([df_both_optimal[:n_fval_Ipopt] df_both_optimal[:n_fval_OnePhase]]);
plot_ratios!(fval_ratios,labels=["#fval Ipopt","#fval One Phase"])

n_hess_ratios = ratio_to_fastest([df_both_optimal[:n_Hess_Ipopt] df_both_optimal[:n_Hess_OnePhase]]);
plot_ratios!(n_hess_ratios,labels=["#hess Ipopt","#hess One Phase"])

df_n = df_both_optimal[:,[:iter_Ipopt,:iter_OnePhase,:n_fac_Ipopt,:n_fac_OnePhase,:n_fval_Ipopt,:n_fval_OnePhase,:n_Hess_Ipopt,:n_Hess_OnePhase]];
summary(df_n)

##################

function test_ratio_to_fastest()
    mat = rand(10,3)
    ratios = ratio_to_fastest(mat)
end

ratio_to_fastest([df_combine[:iter_Ipopt] df_combine[:iter_OnePhase]])

df_combine[:iter] ./ df_combine[:iter_1]

df_opt = df_combine[df_combine[:status] .== "Optimal",:]

it_ratio = log.(df_opt[:iter]) - log.(df_opt[:iter_1])
fac_ratio = log.(df_opt[:n_fac]) - log.(df_opt[:n_fac_1])
f_ratio = log.(df_opt[:n_fval]) - log.(df_opt[:n_fval_1])
c_ratio = log.(df_opt[:n_cons]) - log.(df_opt[:n_cons_1])

median(exp(it_ratio))
median(exp(fac_ratio))
median(exp(f_ratio))
median(exp(c_ratio))

sum(df_combine[:status] .== "Optimal")
sum(df_combine[:status_1] .== "Optimal")
sum(df_combine[:status] .== "primal_infeasible")
sum(df_combine[:status_1] .== "primal_infeasible")
sum(df_combine[:status] .== "dual_infeasible")
sum(df_combine[:status_1] .== "dual_infeasible")


text:
export ZabbixHost, get_hosts

struct ZabbixHost
  id::Int
  name::String
end
ZabbixHost(dict) =
  ZabbixHost(
    parse(Int, dict["hostid"]),
    dict["name"]
  )

showres(h::ZabbixHost; offset=0) =
  println(repeat(" ", offset), "ZabbixHost: ", h.id, " => ", h.name)

function get_hosts(z)
  _hosts = @checkerr make_request(z, "host.get")
  hosts = ZabbixHost[]
  for h in _hosts
    push!(hosts, ZabbixHost(h))
  end
  hosts
end


text:
function gsgp(cfg, X_train, y_train, X_test, y_test; rec=Record(), seed=rand(UInt), verbose=true)
    Random.seed!(seed)

    # these should be executed before main process of GSGP
    init_id()
    init_ref_counter(cfg)
    init_entry_pool()

    # for simplicity
    data_train = Dataset(X_train, y_train)
    data_test  = Dataset(X_test, y_test)

    # initial population (individuals are treated as "entry")
    init_entries = cfg.create_population(cfg, data_train, data_test)
    add_entry!.(init_entries)
    pop = getfield.(init_entries, :index)
    new_pop = similar(pop)
    old_pop = similar(pop)

    avoid_overwrite.(pop)

    fitness = getfield.(init_entries, :train_fitness)

    register!(rec, "train_fitness", train_fitness, Float64)
    register!(rec, "test_fitness", test_fitness, Float64)
    add_record!(rec, 1, pop, cfg)

    verbose && print_latest_fitness(rec, 1)

    # control whether entries that are no longer referenced can be overwritten
    overwrite = false

    for g in 2:cfg.num_generations
        overwrite = g > 2

        elite_idx = argmin(fitness)
        new_pop[1] = pop[elite_idx]
        avoid_overwrite(pop[elite_idx])

        for i in 2:cfg.population_size
            rand_num = rand()

            if rand_num < cfg.p_crossover
                p1 = tournament_selection(pop, fitness, cfg)
                p2 = tournament_selection(pop, fitness, cfg)
                offspring = gscrossover(p1, p2, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            else
                p = tournament_selection(pop, fitness, cfg)
                offspring = gsmutation(p, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            end
        end

        update_fitness!(fitness, new_pop)

        g > 2 && update_entry_pool!(pop)

        add_record!(rec, g, pop, cfg)

        old_pop = copy(pop)
        pop = copy(new_pop)

        verbose && print_latest_fitness(rec, g)
    end

    return pop[argmin(fitness)], rec
end

function gscrossover(t1, t2, data_train, data_test, cfg, overwrite)
    increment_ref(entry_pool[t1])
    increment_ref(entry_pool[t2])

    local r, offspring
    if overwrite
        r = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg)
        add_entry!(r)
        set_ref(r)
        increment_ref(r)

        offspring = Entry(entry_pool[t1], entry_pool[t2], r, crossover, data_train, data_test, cfg)
    else
        r = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg, index=get_counter_size()+1)
        add_entry!(r)
        set_ref(r)
        increment_ref(r)

        offspring = Entry(entry_pool[t1], entry_pool[t2], r, crossover, data_train, data_test, cfg, index=get_counter_size()+1)
    end

    add_entry!(offspring)
    set_ref(offspring)

    # to avoid overwritten
    increment_ref(offspring)

    return offspring
end

function gsmutation(t, data_train, data_test, cfg, overwrite)
    increment_ref(entry_pool[t])

    local offspring, r1, r2
    if overwrite
        r1 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg)
        add_entry!(r1)
        set_ref(r1)
        increment_ref(r1)

        r2 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg)
        add_entry!(r2)
        set_ref(r2)
        increment_ref(r2)

        offspring = Entry(entry_pool[t], r1, r2, mutation, data_train, data_test, cfg)
    else
        r1 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg, index=get_counter_size()+1)
        add_entry!(r1)
        set_ref(r1)
        increment_ref(r1)

        r2 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg, index=get_counter_size()+1)
        add_entry!(r2)
        set_ref(r2)
        increment_ref(r2)

        offspring = Entry(entry_pool[t], r1, r2, mutation, data_train, data_test, cfg, index=get_counter_size()+1)
    end

    add_entry!(offspring)
    set_ref(offspring)

    # to avoid overwritten
    increment_ref(offspring)

    return offspring
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
TOLERANCE = 0.0001

x = Variable()
y = Variable()
solve!(x == 3)
@test evaluate(x) == 3
@test_throws ErrorException solve!(x + 3 * y == 2)

y = Variable(3, 1)
A = rand(3, 3)
solve!(x == 3, A*y == 4)
@test evaluate(x) == 3
@test all(abs(A*evaluate(y) - 4) .<= TOLERANCE)

solve!([x == 5, y == 7])
@test evaluate(x) == 5
@test all(evaluate(y) .== 7)

x = Variable(3)
y = sum_squares(x)
optval = minimize!(y)
@test abs(optval - 0) < TOLERANCE
@test all(abs(evaluate(x)) .<= TOLERANCE)

x = Variable(3)
y = sum_squares(x)
optval = minimize!(y, x == 3)
@test abs(optval - 27) < TOLERANCE
@test all(abs(evaluate(x) - 3) .<= TOLERANCE)

x = Variable(3)
optval = minimize!(sum_squares(3), x == 3)
@test abs(optval - 9) < TOLERANCE
@test all(abs(evaluate(x) - 3) .<= TOLERANCE)

x = Variable(2)
A = [1 0; 1 0]
y = sum_squares(x)
@test_throws ErrorException optval = minimize!(y, A * x == [1; 2])

A = randn(5, 5)
x_real = randn(5, 5)
x = Variable(5, 5)
y = A * x
b = A * x_real
constraints = []
for i = 1:5
  constraints += y[:, i] == b[:, i]
end
optval = minimize!(sum_squares(x), constraints)
@test abs(sum(x_real.^2) - optval) <= TOLERANCE
@test all(abs(evaluate(x) .- x_real) .<= TOLERANCE)

info("All solve/minimize tests passed")


text:
import Base.getindex

"""
    getindex(x::Var, inds...)

```julia
x = Var(rand(Float32,10,5))
y = x[1:3]
y = x[2:2]
```
Note that `y = x[i]` throws an error since `y` is not a vector but a scholar.
Instead, use `y = x[i:i]`.
"""
function getindex(x::Var, I::Tuple)
    Var(x.data[I...], ∇getindex!, (x,I))
end
getindex(x::Var, inds...) = getindex(x, inds)
getindex(x::Node, inds...) = Node(getindex, (x,inds))

function ∇getindex!(y::Var, x::Var, I::Tuple)
    isnothing(x.grad) && return
    addto!(view(x.grad,I...), y.grad)
end


text:
export BC4BP

#---------------------------#
# BC4BP EQUATIONS OF MOTION #
#---------------------------#
struct _BC4BP_ODEFunctions{S,F,F2} <: Abstract_ModelODEFunctions
    ode_system :: S
    ode_f      :: F
    ode_stm_f  :: F2
end

function ModelingToolkit.ODESystem(::Type{_BC4BP_ODEFunctions})
    @parameters  μ   # Mass fraction (smaller 2 bodies)
    @parameters  μ2  # Mass fraction (note: inverse, as per [DeiTos2018])
    @parameters  a3  # Distance to 3rd body
    @parameters  α0  # Initial phase angle of 3rd body
    @parameters  ω3  # Angular velocity of 3rd body
    @parameters  t   # Time
    @variables   x(t) y(t) z(t)
    D, D2 = Differential(t), Differential(t)^2

    # [Dei Tos 2018]
    p   = [x, y, z]                  # Spacecraft position
    α   = ω3 * t
    p_S = a3 .* [cos(α + α0), sin(α + α0), 0]  # Position of 3rd body
    x̂   = [1, 0, 0]                  # X unit vector

    vectors = @. [
        # (Coefficient,     vector to be normalised)
        (-(μ2 - 1),         p - p_S),
        (-(μ2 - 1)*(1 - μ), μ*x̂ + p_S),
        (+(μ2 - 1)*μ,       (1 - μ)*x̂ - p_S),
        (-(1 - μ),          p + μ*x̂),
        (-μ,                p - (1 - μ)*x̂)
    ]

    # Add up the equations [DeiTos2018 Eq. 21]
    forces = sum([coeff .* (v / norm(v)^3) for (coeff, v) in vectors])
    eqs    = @. D2(p) ~ [+2D(y) + x, -2D(x) + y, 0] + forces

    # Simplify and return the system
    return ODESystem(simplify.(eqs), t, [x, y, z], [μ, μ2, a3, α0, a3, ω3]), []
end

# Build the equations at pre-compile time
const BC4BP_ODEFunctions = _BC4BP_ODEFunctions()

#-------------#
# BC4BP MODEL #
#-------------#
struct BC4BP{O<:_BC4BP_ODEFunctions,P<:R4BPSystemProperties} <: Abstract_R3BPModel
    ode   :: O
    props :: P
end
BC4BP(args...; kwargs...) = BC4BP(BC4BP_ODEFunctions, R4BPSystemProperties(args...; kwargs...))


text:
# const getf_master = Libdl.dlsym(polo_lib, :getf_master)
# const getx_master = Libdl.dlsym(polo_lib, :getx_master)

struct Master <: ParameterServer
    fval::Float64
    x::Vector{Float64}
    psopts::ParameterServerOptions

    function (::Type{Master})(; kw...)
        master =  new(0., Vector{Float64}(), ParameterServerOptions(; kw...))
    end
end

function initialize!(master::Master, x₀::AbstractVector)
    resize!(master.x,length(x₀))
    master.x .= x₀
end
getx(master::Master) = master.x

paramserver_handle(::Master) = POLO.proxgradient_master
execution_handle(::Master) = POLO.run_master
delete_handle(::Master) = POLO.delete_proxgradient_master
# getf_handle(::Master) = getf_master
# getx_handle(::Master) = getx_master


text:
# holds all functions related to the type MetaTuple

#all information about each field
struct MetaTuple{T<:Union{String, Bool, Symbol, Missing, Int, Float64}}
  fieldName::Symbol
  columnType::Type
  columnDict::S where S<:Union{Dict, Nothing}
  shortDescription::U where U<:MString
  longDescription::V where V<:MString
end

#for atomic type fields
function MetaTuple(fieldName::Symbol, shortDescription::U where U<:MString,
  longDescription::V where V<:MString, ::Type{T})::MetaTuple where T<:Any

  if !(T <: Union{Bool, String, Missing, Float64, Int})
    error("Invalid atomic metatuple type $T")
  end
  return MetaTuple{T}(fieldName, T, nothing, shortDescription, longDescription)
end

#for conformed fields
function MetaTuple(fieldName::Symbol, valueDict::Dict, shortDescription::MString,
  longDescription::MString, ::Type{T})::MetaTuple where T<:Any

  if !(T==Symbol)
    error("Invalid conformed metatuple type $T")
  end

  return MetaTuple{T}(fieldName::Symbol, CategoricalValue{Union{T, Missing}}#=Union{CategoricalValue{T}, Missing}=#,
    valueDict, shortDescription, longDescription)
end

#convenience method which creates the dictionary of valid values from the dataframe
function MetaTuple(fieldName::Symbol, columnDF::DataFrame,
    shortDescription::MString, longDescription::MString, ::Type{T})::MetaTuple where T<:Any

  NValues::Int = size(columnDF, 1)

  #convert the value column to symbols
  values::Vector{Union{Symbol, Missing}}=
    ((s::U where U<:Union{Int, Float64, String, Missing})->ismissing(s) ? missing : Symbol(s)).(columnDF[:value])

  columnDict::Dict = Dict(values[i] => columnDF[i,:description] for i::Int ∈ 1:NValues)

  return MetaTuple(fieldName, columnDict, shortDescription, longDescription, T)

end

#creates a MetaTuple purely using the fields from the metadata index
function MetaTuple(fileName::String, fieldName::Symbol, fieldType::Symbol, subTableInd::Bool,
    shortDescription::MString, longDescription::MString; subTablesPath::String = SUB_TABLES_PATH)::MetaTuple

    #longDescription = ifelse(ismissing(longDescription),missing,
    #  "\"$(replace(longDescription, "\n"=>"\r"))\"")

    if !subTableInd
      fieldType::Type = Union{VAR_TYPES[fieldType], Missing}
      return MetaTuple(fieldName, shortDescription, longDescription, fieldType)
    else
      fieldType = Symbol

      path::String =  "$subTablesPath\\$fileName\\$fileName - $fieldName.csv"
      #get the subtable info
      #iStream::IOStream = open()
      columnDF::DataFrame = CSV.read(path, strings=:raw, categorical=false)
      #columnDF::DataFrame = load(path) |> DataFrame
      #close(iStream)

      return MetaTuple(fieldName, columnDF, shortDescription, longDescription, fieldType)
    end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Add a calculate_γ from Bertsekas.
# This is when an initial α = 1 far oversteps the first point that hits a boundary
# on the piece-wise linear projected search path. This could be done always, never
# or if last search required a lot of line search reductions.

"""
# ActiveBox
## Constructor
```julia
    ActiveBox(; epsilon=1e-8)
```

`epsilon` determines the threshold for whether a bound is approximately active or not, see eqn. (32) in [1].

## Description
ActiveBox second order for bound constrained convex optimization. It's an active set and allows for rapid exploration of the constraint face. It employs a modified Armijo-line search that takes the active set into account. Details can be found in [1].

## References
- 1) http://www.mit.edu/~dimitrib/ProjectedNewton.pdf
- 2) Iterative Methods for Optimization https://archive.siam.org/books/textbooks/fr18_book.pdf
"""
struct ActiveBox{F,T}
    factorize::F
    ϵ::T
end
ActiveBox(; factorize = cholesky, epsilon = nothing) = ActiveBox(factorize, epsilon)
summary(::ActiveBox) = "ActiveBox"
modelscheme(::ActiveBox) = Newton()
"""
    diagrestrict(x, c, i)

Returns the correct element of the Hessian according to the active set and the diagonal matrix described in [1].

[1] http://www.mit.edu/~dimitrib/ProjectedNewton.pdf
"""
function diagrestrict(x::T, ci, cj, i) where {T}
    if !(ci | cj)
        # If not binding, then return the value
        return x
    else
        # If binding, then return 1 if the diagonal or 0 otherwise
        T(i)
    end
end

function is_ϵ_active(x, lower, upper, ∇fx, ϵ∇f = eltype(x)(0))
    # it is requied that ϵ ⩽ min(U_i - L_i)/2 to uniquely choose
    # an underestimate of the inactive set or else there would be
    # two ways of defining 𝓐^ϵ.
    lowerbinding = x <= lower + ϵ∇f
    upperbinding = x >= upper - ϵ∇f

    pointing_down = ∇fx >= 0
    pointing_up = ∇fx <= 0

    lower_active = lowerbinding && pointing_down
    upper_active = upperbinding && pointing_up

    lower_active || upper_active
end
isbinding(i, j) = i & j

factorize(ab::ActiveBox, M) = ab.factorize(M)
function solve(
    problem::OptimizationProblem,
    x0,
    approach::ActiveBox,
    options::OptimizationOptions,
)
    B0 = false * x0 * x0' + I
    s0 = (x0, B0)
    _solve(problem, s0, approach, options)
end
function solve(
    problem::OptimizationProblem,
    s0::Tuple,
    approach::ActiveBox,
    options::OptimizationOptions,
)
    _solve(problem, s0, approach, options)
end

function _solve(
    prob::OptimizationProblem,
    s0::Tuple,
    scheme::ActiveBox,
    options::OptimizationOptions,
)
    t0 = time()
    x0, B0 = s0

    lower, upper = bounds(prob)
    if isnothing(scheme.ϵ)
        ϵbounds = mapreduce(b -> (b[2] - b[1]) / 2, min, zip(lower, upper)) # [1, pp. 100], [2, 5.41]
    else
        ϵbounds = scheme.ϵ
    end
    !any(clamp.(x0, lower, upper) .!= x0) ||
        error("Initial guess not in the feasible region")

    linesearch = ArmijoBertsekas()
    mstyle = OutOfPlace()

    objvars = prepare_variables(prob, scheme, x0, copy(x0), B0)
    f0, ∇f0 = objvars.fz, norm(objvars.∇fz, Inf) # use user norm
    fz, ∇fz = objvars.fz, objvars.∇fz # use user norm
    fx, ∇fx = fz, copy(∇fz)
    B = B0
    x, z = copy(x0), copy(x0)
    Tf = typeof(fz)
    is_first = false
    Ix = Diagonal(z .* 0 .+ 1)
    for iter = 1:options.maxiter
        x = copy(z)
        fx = copy(fz)
        ∇fx = copy(∇fz)
        ϵ = min(norm(clamp.(x .- ∇fx, lower, upper) .- x), ϵbounds) # Kelley 5.41 and just after (83) in [1]
        activeset = is_ϵ_active.(x, lower, upper, ∇fx, ϵ)
        Hhat = diagrestrict.(B, activeset, activeset', Ix)
        # Update current gradient and calculate the search direction
        HhatFact = factorize(scheme, Hhat)
        d = -(HhatFact \ ∇fx) # use find_direction here

        φ = (; prob, ∇fz, z, x, p = d, φ0 = fz, dφ0 = dot(∇fz, d))

        # Perform line search along d
        # Also returns final step vector and update the state
        α, f_α, ls_success, z = find_steplength(
            mstyle,
            linesearch,
            φ,
            Tf(1),
            ∇fz,
            activeset,
            lower,
            upper,
            x,
            d,
            ∇fx,
            activeset,
        )
        # # Calculate final step vector and update the state
        s = @. x - z

        # Update approximation
        fz, ∇fz, B, s, y = update_obj(prob.objective, s, ∇fx, z, ∇fz, B, Newton(), is_first)
        if norm(x .- clamp.(x .- ∇fz, lower, upper), Inf) < options.g_abstol
            return ConvergenceInfo(
                scheme,
                (
                    prob = prob,
                    B = B,
                    ρs = norm(x .- z),
                    ρx = norm(x),
                    solution = z,
                    fx = fx,
                    minimum = fz,
                    ∇fz = ∇fz,
                    f0 = f0,
                    ∇f0 = ∇f0,
                    iter = iter,
                    time = time() - t0,
                ),
                options,
            )
        end
    end
    iter = options.maxiter
    return ConvergenceInfo(
        scheme,
        (
            prob = prob,
            B = B,
            ρs = norm(x .- z),
            ρx = norm(x),
            solution = z,
            fx = fx,
            minimum = fz,
            ∇fz = ∇fz,
            f0 = f0,
            ∇f0 = ∇f0,
            iter = iter,
            time = time() - t0,
        ),
        options,
    )
end

"""
# ArmijoBertsekas
## Constructor
```julia
    ArmijoBertsekas()
```
## Description
ArmijoBertsekas is the modified Armijo backtracking line search described in [1]. It takes into account whether an element of the gradient is active or not.

## References
[1] http://www.mit.edu/~dimitrib/ActiveBox.pdf
"""
struct ArmijoBertsekas{T1,T2,T3,TR} <: LineSearcher
    ratio::T1
    decrease::T1
    maxiter::T2
    interp::T3
    steprange::TR
    verbose::Bool
end
ArmijoBertsekas(;
    ratio = 0.5,
    decrease = 1e-4,
    maxiter = 50,
    steprange = (0.0, Inf),
    interp = FixedInterp(),
    verbose = false,
) = ArmijoBertsekas(ratio, decrease, maxiter, interp, steprange, verbose)

function find_steplength(
    mstyle,
    ls::ArmijoBertsekas,
    φ::T,
    λ,
    ∇fx,
    Ibool,
    lower,
    upper,
    x,
    p,
    g,
    activeset,
) where {T}
    #== unpack ==#
    φ0, dφ0 = φ.φ0, φ.dφ0
    Tf = typeof(φ0)
    ratio, decrease, maxiter, verbose =
        Tf(ls.ratio), Tf(ls.decrease), ls.maxiter, ls.verbose

    #== factor in Armijo condition ==#
    t0 = decrease * dφ0 # dphi0 should take into account the active set
    iter, α, β = 0, λ, λ # iteration variables
    x⁺ = box_retract.(lower, upper, x, p, α)
    f_α = (; ϕ = φ.prob.objective.f(x⁺))  # initial function value

    if verbose
        println("Entering line search with step size: ", λ)
        println("Initial value: ", φ0)
        println("Value at first step: ", f_α)
    end

    is_solved =
        isfinite(f_α.ϕ) &&
        f_α.ϕ <= φ0 - decrease * sum(bertsekas_R.(x, x⁺, g, p, α, activeset))
    while !is_solved && iter <= maxiter
        iter += 1
        β, α = α, α / 2
        x⁺ = box_retract.(lower, upper, x, p, α)
        f_α = (; ϕ = φ.prob.objective.f(x⁺))  # initial function value
        #        β, α, f_α = interpolate(ls.interp, x->φ, φ0, dφ0, α, f_α.ϕ, ratio)
        is_solved =
            isfinite(f_α.ϕ) &&
            f_α.ϕ <= φ0 - decrease * sum(bertsekas_R.(x, x⁺, g, p, α, activeset))
    end

    ls_success = iter >= maxiter ? false : true

    if verbose
        !ls_success && println("maxiter exceeded in backtracking")
        println("Exiting line search with step size: ", α)
        println("Exiting line search with value: ", f_α)
    end
    return α, f_α, ls_success, x⁺
end

bertsekas_R(x, x⁺, g, p, α, i) = i ? g * (x - x⁺) : α * p * g
# defined univariately
# should be a "manifodl"
box_retract(lower, upper, x, p, α) = min(upper, max(lower, x + α * p))

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
@time @everywhere include("../BioClim/src/required.jl")

# Import CSV file
sp_table = CSV.read("Canis_lupus_baileyi.csv")

# select only scientificName, decimalLatitude, decimalLongitude, and year
sp_table = sp_table[:, [:species, :infraspecificEpithet, :taxonRank,
                        :decimalLatitude, :decimalLongitude, :year]]
# create new column copying all info from species col
sp_table.new_name = copy(sp_table.species)
# replace spaces by "_"
sp_table.new_name .= replace.(sp_table.new_name, " " .=> "_")
# if subspecies, specify in new_name column
for i in 1:length(sp_table)
    if sp_table.taxonRank[i] == "SUBSPECIES"
        sp_table.new_name[i] = string(sp_table.species[i], "_", sp_table.infraspecificEpithet[i])
    end
end

# Subset with specific columns
df = sp_table[:, [:new_name, :year, :decimalLatitude, :decimalLongitude]]
# Rename coordinate names
rename!(df, :decimalLatitude => :latitude)
rename!(df, :decimalLongitude => :longitude)
rename!(df, :new_name => :species)
# Replace spaces by underscores in species names
df.species .= replace.(df.species, " " .=> "_")
# Remove entries with missing year, latitude and lontitude
dropmissing!(df, [:year, :latitude, :longitude])

# Get number of observatons per species/subspecies
newdf = by(df, :species, n = :species => length)
show(sort(newdf, order(:species)), allrows=true)

# Show result for subspecies observations
show(sp_table[sp_table.taxonRank .== "SUBSPECIES", [:species, :infraspecificEpithet, :new_name]], allrows=true)

# # select only occurences with included year date
# sp_table = sp_table[sp_table.year .!= "NA", :]
# # convert year "string" to "integer"
# sp_table.year = parse.(Int64, sp_table.year)

# Prepare csv data using fnc
# df = prepare_csvdata(sp_table)

## From Tim's BioClim/main.jl
# Get the worldclim data by their layer number
cd("/Users/kiristern/Documents/GitHub/BioClim/")

@info "Extract and crop bioclim variables"
@time wc_vars = [clip(worldclim(i), df) for i in 1:19];
# Make the prediction for each layer
@info "Predictions for each layer"
@time predictions = [bioclim(wc_vars[i], df) for i in 1:length(wc_vars)];
# Make the final prediction by taking the minimum
@info "Minimum-consensus aggregation"
@time prediction = reduce(minimum, predictions);
# Get the threshold for NaN given a percentile
@info "Threshold estimation"
@time threshold = first(quantile(prediction[df], [0.05]))
@info "5% threshold:\t$(round(threshold; digits=3))"
# Filter the predictions based on the threshold
@info "Final prediction filtering"
@time for i in eachindex(prediction.grid)
    prediction.grid[i] < threshold && (prediction.grid[i] = NaN)
end

worldmap = clip(worldshape(50), prediction)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = first(unique(df.species)))
xaxis!(sdm_plot, (prediction.left, prediction.right), "Longitude")
yaxis!(sdm_plot, (prediction.bottom, prediction.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(prediction), latitudes(prediction), prediction.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, prediction.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end

sdm_plot

savefig("../SDM/figures/sdm.png")
savefig("../SDM/figures/sdm_map-$(first(unique(sp_table.species))).pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM