text:
using Test, LinearAlgebra, SparseArrays
using JuMP
using Polyhedra

mutable struct MockOptimizer{T} <: Polyhedra.AbstractPolyhedraOptimizer{T}
    lphrep::LPHRep{T}
    # Feasible set defined by user
    rep::Union{Rep{T}, Nothing}
    # Either `rep` or `polyhedron(lphrep, library)`.
    # It is kept between consecutive solve if not modified,
    # e.g. if only the objective is changed.
    feasible_set::Union{Rep{T}, Nothing}

    objective_sense::MOI.OptimizationSense
    objective_func::Union{SparseVector{T, Int64}, Nothing}
    objective_constant::T

    optimize!::Function
    status::MOI.TerminationStatusCode
    solution::Union{AbstractVector{T}, Nothing}

    function MockOptimizer{T}(optimize!::Function) where T
        return new(LPHRep(Polyhedra._MOIModel{T}()), nothing, nothing,
                   MOI.FEASIBILITY_SENSE, nothing, zero(T),
                   optimize!, MOI.OPTIMIZE_NOT_CALLED, nothing)
    end
end
coefficient_type(::MockOptimizer{T}) where {T} = T

function MOI.empty!(mock::MockOptimizer{T}) where T
    mock.lphrep = LPHRep(Polyhedra._MOIModel{T}())
    mock.rep = nothing
    mock.feasible_set = nothing
    mock.objective_sense = MOI.FEASIBILITY_SENSE
    mock.objective_func = nothing
    mock.objective_constant = zero(T)
    mock.status = MOI.OPTIMIZE_NOT_CALLED
    mock.solution = nothing
end
function MOI.is_empty(mock::MockOptimizer{T}) where T
    MOI.is_empty(mock.lphrep.model) &&
    mock.rep === nothing &&
    mock.feasible_set === nothing &&
    mock.objective_sense == MOI.FEASIBILITY_SENSE &&
    mock.objective_func === nothing &&
    iszero(mock.objective_constant) &&
    mock.status == MOI.OPTIMIZE_NOT_CALLED &&
    mock.solution === nothing
end

function MOI.optimize!(mock::MockOptimizer)
    mock.optimize!(mock)
end

MOI.get(mock::MockOptimizer, ::MOI.TerminationStatus) = mock.status
function MOI.get(mock::MockOptimizer, ::MOI.ResultCount)
    if mock.status == MOI.OPTIMAL || mock.status == MOI.DUAL_INFEASIBLE
        return 1
    else
        return 0
    end
end
function MOI.get(mock::MockOptimizer{T}, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                         <:Union{MOI.EqualTo{T},
                                                 MOI.LessThan{T}}}) where T
    return MOIU.get_fallback(mock, attr, ci)
end
function MOI.get(mock::MockOptimizer, ::MOI.ObjectiveValue)
    if mock.status == MOI.OPTIMAL
        return mock.objective_func ⋅ mock.solution + mock.objective_constant
    elseif mock.status == MOI.DUAL_INFEASIBLE
        return mock.objective_func ⋅ mock.solution
    else
        error("No objective value available when termination status is $(mock.status).")
    end
end
function MOI.get(mock::MockOptimizer, ::MOI.PrimalStatus)
    if mock.status == MOI.OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif mock.status == MOI.DUAL_INFEASIBLE
        return MOI.INFEASIBILITY_CERTIFICATE
    else
        return MOI.NO_SOLUTION
    end
end
function MOI.get(mock::MockOptimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if mock.status != MOI.OPTIMAL && mock.status != MOI.DUAL_INFEASIBLE
        error("No primal value available when termination status is $(mock.status).")
    end
    return mock.solution[vi.value]
end

function polyhedra_model_test(factory::JuMP.OptimizerFactory)
    model = Model(factory)
    n = 3
    @variable(model, 0 <= x[1:n] <= 1)
    @variable(model, y)
    @constraint(model, 2y == 2)
    csum = @constraint(model, sum(x[i] for i=1:n) == 1)
    c12 = @constraint(model, x[1] - x[2] <= 1)
    c23 = @constraint(model, x[2] - x[3] >= 1)
    @objective(model, Max, sum(x) - 2y)

    optimize!(model)
    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) == -1.0
    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value.(x) ≈ [0.0, 1.0, 0.0]
    @test value(y) ≈ 1.0
    @test value(csum) ≈ 1.0
    @test value(c12) ≈ -1.0
    @test value(c23) ≈ 1.0

    # TODO dual
end

@testset "AbstractPolyhedraOptimizer" begin
    polyhedra_model_test(with_optimizer(
        MockOptimizer{Float64},
        mock -> begin
            mock.status = MOI.OPTIMAL
            mock.solution = [0.0, 1.0, 0.0, 1.0]
        end))
end


text:
function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    ODESystem(eqs, t, [], [], systems=[g], name=name)
end

function Resistor(;name, R = 1.0)
    val = R

    @named p = Pin()
    @named n = Pin()
    @parameters R
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ p.i * R
          ]
    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
end

function Capacitor(; name, C = 1.0)
    val = C

    @named p = Pin()
    @named n = Pin()
    @parameters C
    @variables v(t)

    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(v) ~ p.i / C
          ]
    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
end

function Inductor(; name, L = 1.0)
    val = L

    @named p = Pin()
    @named n = Pin()
    @parameters L
    @variables v(t)

    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(p.i) ~ v / L
          ]
    ODESystem(eqs, t, [v], [L], systems=[p, n], defaults=Dict(L => val), name=name)
end

function IdealOpAmp(; name)
    @named p1 = Pin()
    @named p2 = Pin()
    @named n1 = Pin()
    @named n2 = Pin()
    @variables v1(t) v2(t) # u"v"
    @variables i1(t) i2(t) # u"A"

    eqs = [
        v1 ~ p1.v - n1.v
        v2 ~ p2.v - n2.v
        0 ~ p1.i + n1.i
        0 ~ p2.i + n2.i
        i1 ~ p1.i
        i2 ~ p2.i
        v1 ~ 0
        i1 ~ 0
    ]
    ODESystem(eqs, t, [i1, i2, v1, v2], [], systems=[p1, p2, n1, n2], name=name)
end

ext: jl


text:
using Hyperscript, Markdown
using JSServe, Observables
using JSServe: Application, Session, evaljs, linkjs, update_dom!, div, active_sessions
using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session
using JSServe.DOM

with_session() do session, req
    s1 = Slider(1:100)
    s2 = Slider(1:100)
    b = Button("hi")
    t = TextField("Write!")
    linkjs(session, s1.value, s2.value)
    onjs(session, s1.value, js"(v)=> console.log(v)")
    on(t) do text
        println(text)
    end
    # scene = scatter(
    #     1:100, rand(100) .* 100,
    #     markersize = s1, axis = (names = (title = t,),)
    # )
    return md = md"""
    # IS THIS REAL?

    My first slider: $(s1)

    My second slider: $(s2)

    Test: $(s1.value)

    The BUTTON: $(b)

    Type something for the list: $(t)

    some list $(t.value)
    """
end


text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setup.jl");

# periodic boundary conditions
pbc = [1; 1; 1]

# path to the database
datapath = "/Users/ngoccuongnguyen/Dropbox (MIT)/Research/Software/MDP/fitting/snap/Ta/JSON/";
dataformat = "json"
fileextension = "json"
atomspecies = ["Ta"];

# folders in the datapath
folders = ["Displaced_A15", "Displaced_BCC", "Displaced_FCC", "Elastic_BCC",
            "Elastic_FCC", "GSF_110", "GSF_112", "Liquid", "Surface",
            "Volume_A15", "Volume_BCC", "Volume_FCC"];

# weights for energies, forces, and stresses in the linear fit
weightinner=[100            1               1.00E-08
            100             1               1.00E-08
            100             1               1.00E-08
            1.00E-08        1.00E-08        0.0001
            1.00E-09        1.00E-09        1.00E-09
            100             1               1.00E-08
            100             1               1.00E-08
            4.67E+02        1               1.00E-08
            100             1               1.00E-08
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09]

# randomly selecting the configurations in the database
randomize = false;

# use all the data
percentage = 100.0;

# translate atom positions
translationvector = nothing

# rotate atom positions
rotationmatrix = nothing

# transpose lattice vectors because the lattice ordering is python-based
transposelattice = true

# training data
for i = 1:length(folders)
    traindata[i] = adddata(datapath * folders[i], dataformat, fileextension,
                percentage, randomize, atomspecies, weightinner[i,:], translationvector,
                rotationmatrix, transposelattice)
end

rcut = 4.67637   # cut-off radiuss

# single descriptors to catpure energy of isolated pure elements
descriptors[1] = PotentialStruct("single", "nonbonded", [1], "single", [0.0], rcut);
include(descriptors[1].potentialfunction * ".jl");  # include the descriptors file

besselscale = 2.5;
eta = [rcut; besselscale];
kappa = [6];      # number of Bessel descriptors

# Bessel descriptors
descriptors[2] = PotentialStruct("pair", "nonbonded", [1, 1], "scaledbessel", [0.0], rcut);
include(descriptors[2].potentialfunction * ".jl");  # include the descriptors file

# SNAP parameters
twojmax = 6
ntypes = 1
rcutfac = rcut;
rfac0 = 0.99363;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)

# SNAP descriptors
descriptors[3] = initsna(snaparam, elemradius, elemweight)

# loss function style must be: "energy", "force", "energyforce", "energyforcestress"
lossfunc = "energyforcestress"
# use least-square method
method = "lsq"
# divide energies by the number of atoms
normalizeenergy = true
# divide virials by the box volume
normalizestress = true
# optimization parameters
optim = setoptim(lossfunc, method, normalizeenergy, normalizestress, eta, kappa)

# linear fit to compute SNAP coefficients
coeff,~ = linearfit2(traindata, descriptors, potentials, optim, pbc)

print("SNAP Coeffficients: "), show(stdout, "text/plain", coeff)

# compute unweighted MAE, RMSE, RSQ errors
energyerrors, forceerrors, stresserrors = validate2(traindata, descriptors, potentials, optim, coeff, pbc)

printerrors(folders, energyerrors, "Energy Errors")
printerrors(folders, forceerrors, "Force Errors")
printerrors(folders, stresserrors, "Stress Errors")

# ----------------------------------------------------------------------------------------
# Energy Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.07752835884217    |    0.22705475009736    |    0.99596448956663   |
# Displaced_A15 |   0.00049192909655    |    0.00058897195332    |    0.95326317157895   |
# Displaced_BCC |   0.00179938127019    |    0.00240364187430    |    0.99758372333633   |
# Displaced_FCC |   0.00229495671435    |    0.00240573323954    |    0.63958365406662   |
# Elastic_BCC   |   0.01227624775096    |    0.01229967922620    |    -735.09749881215   |
# Elastic_FCC   |   0.00281118185418    |    0.00369493891024    |    -2.3502044635248   |
# GSF_110       |   0.00473898165805    |    0.00555933602248    |    0.77894211169715   |
# GSF_112       |   0.00556624505684    |    0.00617549243995    |    0.89842350462463   |
# Liquid        |   0.00160658042785    |    0.00180561530815    |    0.99960899813060   |
# Surface       |   0.01298185518923    |    0.01588038702008    |    0.93534466141737   |
# Volume_A15    |   0.20707055289531    |    0.25553937978906    |    0.99375825428829   |
# Volume_BCC    |   0.29043994595705    |    0.39373648391637    |    0.99855177887171   |
# Volume_FCC    |   0.45029070090874    |    0.65941196368284    |    0.98171731349904   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Force Errors  |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.07435126021971    |    0.12848376278885    |    0.96636534845803   |
# Displaced_A15 |   0.09159588698128    |    0.11543137635400    |    0.96709623259060   |
# Displaced_BCC |   0.11631098454386    |    0.14961150381635    |    0.98750172427778   |
# Displaced_FCC |   0.08406970708774    |    0.10543003421090    |    0.95619944089030   |
# Elastic_BCC   |   0.05860645433728    |    0.07186259158171    |    -47819.762694298   |
# Elastic_FCC   |   0.04317009096991    |    0.05538270877199    |    -18726.563076726   |
# GSF_110       |   0.02727025920704    |    0.04553237860642    |    0.87898335977552   |
# GSF_112       |   0.04987378103932    |    0.07922494107237    |    0.82408980275735   |
# Liquid        |   0.27366116024599    |    0.34985442714443    |    0.95494518581465   |
# Surface       |   0.04259395186436    |    0.09563703931568    |    0.82445969671280   |
# Volume_A15    |   3.37240113741997    |    6.15659213169514    |    -Inf000000000000   |
# Volume_BCC    |   3.32863251469136    |    6.98673598115050    |    -Inf000000000000   |
# Volume_FCC    |   1.97839064823913    |    3.85646088762313    |    -Inf000000000000   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Stress Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   46622.9349273975    |    235182.209214730    |    0.98144926194258   |
# Displaced_A15 |   1221.04850726689    |    1553.55605717264    |    0.98563925123290   |
# Displaced_BCC |   7341.61299136141    |    9909.83900242742    |    0.86467137841218   |
# Displaced_FCC |   29864.5485879079    |    41957.0091551672    |    -2.6414119710704   |
# Elastic_BCC   |   425.430825216202    |    504.480804809002    |    0.99619587954770   |
# Elastic_FCC   |   39292.2213685134    |    54711.6947589831    |    -11.003934696752   |
# GSF_110       |   2084.38759763698    |    3183.70255860428    |    0.50643293174301   |
# GSF_112       |   1707.54624727980    |    2859.07610929099    |    0.71666105292148   |
# Liquid        |   42390.0156468890    |    58778.5888928048    |    -6.0401903250365   |
# Surface       |   1662.20520087705    |    2917.99248883842    |    0.90799433870165   |
# Volume_A15    |   70617.2614162872    |    140043.228853478    |    0.99284697772326   |
# Volume_BCC    |   239272.502211254    |    672446.741169593    |    0.98672888320201   |
# Volume_FCC    |   169066.027472050    |    558413.504176687    |    0.94963712438035   |
# ----------------------------------------------------------------------------------------

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
struct Symplectic2ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  b1::T
  b2::T
end

function PseudoVerletLeapfrogConstantCache(T,T2)
  a1 = convert(T,1)
  a2 = convert(T,0)
  b1 = convert(T,1//2)
  b2 = convert(T,1//2)
  Symplectic2ConstantCache{T,T2}(a1,a2,b1,b2)
end

function McAte2ConstantCache(T,T2)
  a2 = convert(T,1 - (1/2)*sqrt(convert(T,2)))
  a1 = convert(T,1 - a2)
  b2 = convert(T,1/(2*(1-a2)))
  b1 = convert(T,1 - b2)
  Symplectic2ConstantCache{T,T2}(a1,a2,b1,b2)
end

function VerletLeapfrogConstantCache(T,T2)
  a1 = convert(T,1//2)
  a2 = convert(T,1//2)
  b1 = convert(T,0)
  b2 = convert(T,1)
  Symplectic2ConstantCache{T,T2}(a1,a2,b1,b2)
end

struct Symplectic3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  b1::T
  b2::T
  b3::T
end

function Ruth3ConstantCache(T,T2)
  a1 = convert(T,2//3)
  a2 = convert(T,-2//3)
  a3 = convert(T,1)
  b1 = convert(T,7//24)
  b2 = convert(T,3//4)
  b3 = convert(T,-1//24)
  Symplectic3ConstantCache{T,T2}(a1,a2,a3,b1,b2,b3)
end

function McAte3ConstantCache(T,T2)
  a1 = convert(T,0.9196615230173999)
  a2 = convert(T,0.25/a1 - a1/2)
  a3 = convert(T,1 - a1 - a2)
  b1 = convert(T,a3)
  b2 = convert(T,a2)
  b3 = convert(T,a1)
  Symplectic3ConstantCache{T,T2}(a1,a2,a3,b1,b2,b3)
end

struct Symplectic4ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  b1::T
  b2::T
  b3::T
  b4::T
end

function CandyRoz4ConstantCache(T,T2)
  a1 = convert(T,(2 + T(2)^(1//3) + convert(T,2)^(-1//3))/6)
  a2 = convert(T,(1 - T(2)^(1//3) - convert(T,2)^(-1//3))/6)
  a3 = convert(T,a2)
  a4 = convert(T,a1)
  b1 = convert(T,0)
  b2 = convert(T,(2 - T(2)^(1//3))^-1)
  b3 = convert(T,(1-T(2)^(2//3))^-1)
  b4 = convert(T,b2)
  Symplectic4ConstantCache{T,T2}(a1,a2,a3,a4,b1,b2,b3,b4)
end

function McAte4ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1 =convert(T, 0.515352837431122936)
  a2 =convert(T,-0.085782019412973646)
  a3 =convert(T, 0.441583023616466524)
  a4 =convert(T, 0.128846158365384185)
  b1 =convert(T, 0.134496199277431089)
  b2 =convert(T,-0.224819803079420806)
  b3 =convert(T, 0.756320000515668291)
  b4 =convert(T, 0.334003603286321425)
  Symplectic4ConstantCache{T,T2}(a1,a2,a3,a4,b1,b2,b3,b4)
end

function McAte4ConstantCache(T::Type,T2::Type)
  a1 =convert(T,big"0.515352837431122936")
  a2 =convert(T,big"-0.085782019412973646")
  a3 =convert(T,big" 0.441583023616466524")
  a4 =convert(T,big" 0.128846158365384185")
  b1 =convert(T,big" 0.134496199277431089")
  b2 =convert(T,big"-0.224819803079420806")
  b3 =convert(T,big" 0.756320000515668291")
  b4 =convert(T,big" 0.334003603286321425")
  Symplectic4ConstantCache{T,T2}(a1,a2,a3,a4,b1,b2,b3,b4)
end

struct Symplectic45ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
end

function CalvoSanz4ConstantCache(T,T2)
  a1= convert(T,0.205177661542290)
  a2= convert(T,0.403021281604210)
  a3=-convert(T,0.12092087633891)
  a4= convert(T,0.512721933192410)
  a5= convert(T,0.0)
  b1= convert(T,0.061758858135626)
  b2= convert(T,0.33897802655364)
  b3= convert(T,0.61479130717558)
  b4=-convert(T,0.14054801465937)
  b5= convert(T,0.12501982279453)
  Symplectic45ConstantCache{T,T2}(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5)
end

# Broken
# http://epubs.siam.org/doi/pdf/10.1137/0916010
# On the numerical integration of ordinary differential equations by symmetric composition methods
function McAte42ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1 =convert(T,0.40518861839525227722)
  a2 =convert(T,-0.28714404081652408900)
  a3 = 1 - 2a1 - 2a2
  a4 = a2
  a5 = a1
  b1 = convert(T,-3//73)
  b2 = convert(T,17//59)
  b3 = 1 - 2b1 - 2b2
  b4 = b2
  b5 = b1
  Symplectic45ConstantCache{T,T2}(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5)
end

function McAte42ConstantCache(T::Type,T2::Type)
  a1 =convert(T,big"0.40518861839525227722")
  a2 =convert(T,big"-0.28714404081652408900")
  a3 = 1 - 2a1 - 2a2
  a4 = a2
  a5 = a1
  b1 = convert(T,-3//73)
  b2 = convert(T,17//59)
  b3 = 1 - 2b1 - 2b2
  b4 = b2
  b5 = b1
  Symplectic45ConstantCache{T,T2}(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5)
end

struct Symplectic5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
end

function McAte5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1=convert(T,0.339839625839110000)
  a2=convert(T,-0.088601336903027329)
  a3=convert(T,0.5858564768259621188)
  a4=convert(T,-0.603039356536491888)
  a5=convert(T,0.3235807965546976394)
  a6=convert(T,0.4423637942197494587)
  b1=convert(T,0.1193900292875672758)
  b2=convert(T,0.6989273703824752308)
  b3=convert(T,-0.1713123582716007754)
  b4=convert(T,0.4012695022513534480)
  b5=convert(T,0.0107050818482359840)
  b6=convert(T,-0.0589796254980311632)
  Symplectic5ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6)
end

function McAte5ConstantCache(T::Type,T2::Type)
  a1=convert(T,big"0.339839625839110000")
  a2=convert(T,big"-0.088601336903027329")
  a3=convert(T,big"0.5858564768259621188")
  a4=convert(T,big"-0.603039356536491888")
  a5=convert(T,big"0.3235807965546976394")
  a6=convert(T,big"0.4423637942197494587")
  b1=convert(T,big"0.1193900292875672758")
  b2=convert(T,big"0.6989273703824752308")
  b3=convert(T,big"-0.1713123582716007754")
  b4=convert(T,big"0.4012695022513534480")
  b5=convert(T,big"0.0107050818482359840")
  b6=convert(T,big"-0.0589796254980311632")
  Symplectic5ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6)
end

struct Symplectic6ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  a7::T
  a8::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
end

function Yoshida6ConstantCache(T,T2)
  a1=convert(T, 0.78451361047756)
  a2=convert(T, 0.23557321335936)
  a3=convert(T,-1.1776799841789)
  a4=convert(T, 1.3151863206839)
  a5=convert(T, a3)
  a6=convert(T, a2)
  a7=convert(T, a1)
  a8=convert(T, 0.0)
  b1= a1/2
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T, b4)
  b6=convert(T, b3)
  b7=convert(T, b2)
  b8=convert(T, b1)
  Symplectic6ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8)
end

struct Symplectic62ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  a7::T
  a8::T
  a9::T
  a10::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
end

function KahanLi6ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1=convert(T,0.39216144400731413927925056)
  a2=convert(T,0.33259913678935943859974864)
  a3=convert(T,-0.70624617255763935980996482)
  a4=convert(T,0.08221359629355080023149045)
  a5=convert(T,0.79854399093482996339895035)
  a6= a4
  a7= a3
  a8= a2
  a9= a1
  a10 = convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=b5
  b7=b4
  b8=b3
  b9=b2
  b10=b1
  Symplectic62ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10)
end

function KahanLi6ConstantCache(T::Type,T2::Type)
  a1=convert(T,big"0.39216144400731413927925056")
  a2=convert(T,big"0.33259913678935943859974864")
  a3=convert(T,big"-0.70624617255763935980996482")
  a4=convert(T,big"0.08221359629355080023149045")
  a5=convert(T,big"0.79854399093482996339895035")
  a6= a4
  a7= a3
  a8= a2
  a9= a1
  a10 = convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=b5
  b7=b4
  b8=b3
  b9=b2
  b10=b1
  Symplectic62ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10)
end

struct McAte8ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  a7::T
  a8::T
  a9::T
  a10::T
  a11::T
  a12::T
  a13::T
  a14::T
  a15::T
  a16::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
  b11::T
  b12::T
  b13::T
  b14::T
  b15::T
  b16::T
end

function McAte8ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1=convert(T,0.74167036435061295344822780)
  a2=convert(T,-0.40910082580003159399730010)
  a3=convert(T,0.19075471029623837995387626)
  a4=convert(T,-0.57386247111608226665638773)
  a5=convert(T,0.29906418130365592384446354)
  a6=convert(T,0.33462491824529818378495798)
  a7=convert(T,0.31529309239676659663205666)
  a8=convert(T,-0.79688793935291635401978884)
  a9=a7
  a10=a6
  a11=a5
  a12=a4
  a13=a3
  a14=a2
  a15=a1
  a16=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=b8
  b10=b7
  b11=b6
  b12=b5
  b13=b4
  b14=b3
  b15=b2
  b16=b1
  McAte8ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,
                            b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16)
end

function McAte8ConstantCache(T::Type,T2::Type)
  a1=convert(T,big"0.74167036435061295344822780")
  a2=convert(T,big"-0.40910082580003159399730010")
  a3=convert(T,big"0.19075471029623837995387626")
  a4=convert(T,big"-0.57386247111608226665638773")
  a5=convert(T,big"0.29906418130365592384446354")
  a6=convert(T,big"0.33462491824529818378495798")
  a7=convert(T,big"0.31529309239676659663205666")
  a8=convert(T,big"-0.79688793935291635401978884")
  a9=a7
  a10=a6
  a11=a5
  a12=a4
  a13=a3
  a14=a2
  a15=a1
  a16=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=b8
  b10=b7
  b11=b6
  b12=b5
  b13=b4
  b14=b3
  b15=b2
  b16=b1
  McAte8ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,
                            b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16)
end

struct KahanLi8ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  a7::T
  a8::T
  a9::T
  a10::T
  a11::T
  a12::T
  a13::T
  a14::T
  a15::T
  a16::T
  a17::T
  a18::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
  b11::T
  b12::T
  b13::T
  b14::T
  b15::T
  b16::T
  b17::T
  b18::T
end

function KahanLi8ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1=convert(T,0.13020248308889008087881763)
  a2=convert(T,0.56116298177510838456196441)
  a3=convert(T,-0.38947496264484728640807860)
  a4=convert(T,0.15884190655515560089621075)
  a5=convert(T,-0.39590389413323757733623154)
  a6=convert(T,0.18453964097831570709183254)
  a7=convert(T,0.25837438768632204729397911)
  a8=convert(T,0.29501172360931029887096624)
  a9=convert(T,-0.60550853383003451169892108)
  a10=a8
  a11=a7
  a12=a6
  a13=a5
  a14=a4
  a15=a3
  a16=a2
  a17=a1
  a18=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=convert(T,(a8+a9)/2)
  b10=b9
  b11=b8
  b12=b7
  b13=b6
  b14=b5
  b15=b4
  b16=b3
  b17=b2
  b18=b1
  KahanLi8ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
                              b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18)
end

function KahanLi8ConstantCache(T::Type,T2::Type)
  a1=convert(T,big"0.13020248308889008087881763")
  a2=convert(T,big"0.56116298177510838456196441")
  a3=convert(T,big"-0.38947496264484728640807860")
  a4=convert(T,big"0.15884190655515560089621075")
  a5=convert(T,big"-0.39590389413323757733623154")
  a6=convert(T,big"0.18453964097831570709183254")
  a7=convert(T,big"0.25837438768632204729397911")
  a8=convert(T,big"0.29501172360931029887096624")
  a9=convert(T,big"-0.60550853383003451169892108")
  a10=a8
  a11=a7
  a12=a6
  a13=a5
  a14=a4
  a15=a3
  a16=a2
  a17=a1
  a18=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=convert(T,(a8+a9)/2)
  b10=b9
  b11=b8
  b12=b7
  b13=b6
  b14=b5
  b15=b4
  b16=b3
  b17=b2
  b18=b1
  KahanLi8ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
                              b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18)
end

struct SofSpa10ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a1::T
  a2::T
  a3::T
  a4::T
  a5::T
  a6::T
  a7::T
  a8::T
  a9::T
  a10::T
  a11::T
  a12::T
  a13::T
  a14::T
  a15::T
  a16::T
  a17::T
  a18::T
  a19::T
  a20::T
  a21::T
  a22::T
  a23::T
  a24::T
  a25::T
  a26::T
  a27::T
  a28::T
  a29::T
  a30::T
  a31::T
  a32::T
  a33::T
  a34::T
  a35::T
  a36::T
  b1::T
  b2::T
  b3::T
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T
  b9::T
  b10::T
  b11::T
  b12::T
  b13::T
  b14::T
  b15::T
  b16::T
  b17::T
  b18::T
  b19::T
  b20::T
  b21::T
  b22::T
  b23::T
  b24::T
  b25::T
  b26::T
  b27::T
  b28::T
  b29::T
  b30::T
  b31::T
  b32::T
  b33::T
  b34::T
  b35::T
  b36::T
end

function SofSpa10ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a1=convert(T,0.07879572252168641926390768)
  a2=convert(T,0.31309610341510852776481247)
  a3=convert(T,0.02791838323507806610952027)
  a4=convert(T,-0.22959284159390709415121340)
  a5=convert(T,0.13096206107716486317465686)
  a6=convert(T,-0.26973340565451071434460973)
  a7=convert(T,0.07497334315589143566613711)
  a8=convert(T,0.11199342399981020488957508)
  a9=convert(T,0.36613344954622675119314812)
  a10=convert(T,-0.39910563013603589787862981)
  a11=convert(T,0.10308739852747107731580277)
  a12=convert(T,0.41143087395589023782070412)
  a13=convert(T,-0.00486636058313526176219566)
  a14=convert(T,-0.39203335370863990644808194)
  a15=convert(T,0.05194250296244964703718290)
  a16=convert(T,0.05066509075992449633587434)
  a17=convert(T,0.04967437063972987905456880)
  a18=convert(T,0.04931773575959453791768001)
  a19=a17
  a20=a16
  a21=a15
  a22=a14
  a23=a13
  a24=a12
  a25=a11
  a26=a10
  a27=a9
  a28=a8
  a29=a7
  a30=a6
  a31=a5
  a32=a4
  a33=a3
  a34=a2
  a35=a1
  a36=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=convert(T,(a8+a9)/2)
  b10=convert(T,(a9+a10)/2)
  b11=convert(T,(a10+a11)/2)
  b12=convert(T,(a11+a12)/2)
  b13=convert(T,(a12+a13)/2)
  b14=convert(T,(a13+a14)/2)
  b15=convert(T,(a14+a15)/2)
  b16=convert(T,(a15+a16)/2)
  b17=convert(T,(a16+a17)/2)
  b18=convert(T,(a17+a18)/2)
  b19=b18
  b20=b17
  b21=b16
  b22=b15
  b23=b14
  b24=b13
  b25=b12
  b26=b11
  b27=b10
  b28=b9
  b29=b8
  b30=b7
  b31=b6
  b32=b5
  b33=b4
  b34=b3
  b35=b2
  b36=b1
  SofSpa10ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
                              a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,
                              a35,a36,
                              b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,
                              b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,
                              b35,b36)
end

function SofSpa10ConstantCache(T::Type,T2::Type)
  a1=convert(T,big"0.07879572252168641926390768")
  a2=convert(T,big"0.31309610341510852776481247")
  a3=convert(T,big"0.02791838323507806610952027")
  a4=convert(T,big"-0.22959284159390709415121340")
  a5=convert(T,big"0.13096206107716486317465686")
  a6=convert(T,big"-0.26973340565451071434460973")
  a7=convert(T,big"0.07497334315589143566613711")
  a8=convert(T,big"0.11199342399981020488957508")
  a9=convert(T,big"0.36613344954622675119314812")
  a10=convert(T,big"-0.39910563013603589787862981")
  a11=convert(T,big"0.10308739852747107731580277")
  a12=convert(T,big"0.41143087395589023782070412")
  a13=convert(T,big"-0.00486636058313526176219566")
  a14=convert(T,big"-0.39203335370863990644808194")
  a15=convert(T,big"0.05194250296244964703718290")
  a16=convert(T,big"0.05066509075992449633587434")
  a17=convert(T,big"0.04967437063972987905456880")
  a18=convert(T,big"0.04931773575959453791768001")
  a19=a17
  a20=a16
  a21=a15
  a22=a14
  a23=a13
  a24=a12
  a25=a11
  a26=a10
  a27=a9
  a28=a8
  a29=a7
  a30=a6
  a31=a5
  a32=a4
  a33=a3
  a34=a2
  a35=a1
  a36=convert(T,0)
  b1= convert(T,a1/2)
  b2=convert(T,(a1+a2)/2)
  b3=convert(T,(a2+a3)/2)
  b4=convert(T,(a3+a4)/2)
  b5=convert(T,(a4+a5)/2)
  b6=convert(T,(a5+a6)/2)
  b7=convert(T,(a6+a7)/2)
  b8=convert(T,(a7+a8)/2)
  b9=convert(T,(a8+a9)/2)
  b10=convert(T,(a9+a10)/2)
  b11=convert(T,(a10+a11)/2)
  b12=convert(T,(a11+a12)/2)
  b13=convert(T,(a12+a13)/2)
  b14=convert(T,(a13+a14)/2)
  b15=convert(T,(a14+a15)/2)
  b16=convert(T,(a15+a16)/2)
  b17=convert(T,(a16+a17)/2)
  b18=convert(T,(a17+a18)/2)
  b19=b18
  b20=b17
  b21=b16
  b22=b15
  b23=b14
  b24=b13
  b25=b12
  b26=b11
  b27=b10
  b28=b9
  b29=b8
  b30=b7
  b31=b6
  b32=b5
  b33=b4
  b34=b3
  b35=b2
  b36=b1
  SofSpa10ConstantCache{T,T2}(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
                              a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,
                              a35,a36,
                              b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,
                              b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,
                              b35,b36)
end


text:
t = Tuple{UInt8, Bool}          # a tuple type
@assert typeof(t) == DataType   # tuple types aren't special
@assert isa(t, Type)            # tuple types are types

@assert isa((0x08, false), t)   # tuples of the right type belong
@assert isa((0x7A, true), t)
@assert !isa((3, "wrong"), t)   # tuples of the wrong type do not
@assert !isa(0x25, t)           # non-tuples do not belong

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
function to_magma(f::IOStream, a::Vector{T}; name::String="R") where T
  print(f, name, " := [\n")
  for i=1:(length(a)-1)
    try
      to_magma(f, a[i]);
    catch b
      print(f, a[i])
    end
    print(f, ",\n")
  end
  try
    to_magma(f, a[end])
  catch b
    print(f, a[end])
  end
  print(f, "];\n")
end

function to_magma(f::IOStream, t::Tuple)
  print(f, "<")
  for i=1:length(t)
    try
      to_magma(f, t[i])
    catch e
      print(f, t[i])
    end
    if i<length(t)
      print(f, ", ")
    else
      print(f, ">\n")
    end
  end
end

function to_magma(s::String, a::Vector{T}; name::String="R", mode::String ="w") where T
  f = open(s, mode)
  to_magma(f, a, name = name)
  close(f)
end

################################################################################
# fmpz_mat -> magma file
# use as read(...)
################################################################################
function to_magma(io::IOStream, A::fmpz_mat; name = "A")
  println(io, name, " := Matrix(Integers(), ", nrows(A), ", ", ncols(A), ", [")
  for i = 1:nrows(A)
    for j = 1:ncols(A)
      print(io, A[i,j])
      if j < ncols(A)
        print(io, ", ")
      end
    end
    if i<nrows(A)
      println(io, ",")
    end
  end
  println(io, "]);")
  println(io, "\"Loaded ", name, "\";")
end

function to_magma(s::String, A::fmpz_mat)
  f = open(s, "w")
  to_magma(f, A)
  close(f)
end

################################################################################
# SMat -> magma file
# use as read(...)
################################################################################
function to_magma(io::IOStream, A::SMat; name = "A")
  println(io, name, " := SparseMatrix(Integers(), ", nrows(A), ", ", ncols(A), ", [")
  for i = 1:nrows(A)
    for xx = 1:length(A.rows[i].pos)
      print(io, "<", i, ", ", A.rows[i].pos[xx], ", ", A.rows[i].values[xx], ">")
      if xx < length(A.rows[i].pos) || i<nrows(A)
        print(io, ", ")
      end
    end
    println(io, "")
  end
  println(io, "]);")
  println(io, "\"Loaded ", name, "\";")
end

function to_magma(s::String, A::SMat)
  f = open(s, "w")
  to_magma(f, A)
  close(f)
end

################################################################################
# MPoly -> magma file
# use as read(...)
################################################################################
function to_magma(io::IOStream, R::AbstractAlgebra.MPolyRing; base_name::String = "S", name::String = "R")
  print(io, "$name<")
  S = symbols(R)
  for i = 1:length(S)-1
    print(io, "$(S[i]),")
  end
  print(io, "$(S[end])> := PolynomialRing($base_name, $(length(S)));\n")
end

function to_magma(p::String, R::AbstractAlgebra.MPolyRing; base_name::String = "S", name::String = "R", make::String = "w")
  f = open(p, mode)
  Hecke.to_magma(f, R, base_name = base_name, name = name)
  close(f)
end

function to_magma(io::IOStream, f::Generic.MPolyElem)
  S = symbols(parent(f))
  for i=1:length(f)
    if i>1
      print(io, "+")
    end
    s = "$(coeff(f, i))"
    s = replace(s, "//" => "/")
    print(io, "($s)")
    e = exponent_vector(f, i)
    if iszero(e)
      continue
    end
    print(io, "*")
    fi = true
    for j=1:length(S)
      if e[j] > 0
        if !fi
          print(io, "*")
        else
          fi = false
        end
        print(io, "$(S[j])^$(e[j])")
      end
    end
  end
end

function to_magma(io::IOStream, k::AnticNumberField; name::String = "S", gen_name::String="_a")
  print(io, "$name<$gen_name> := NumberField($(k.pol));\n")
end

function to_magma(io::IOStream, s::Symbol, v::Any)
  print(io, "$s := ")
  Hecke.to_magma(io, v)
  print(io, ";\n")
end


text:
using Latexify
using Test

@test latexify(((1.0, 2), (3, 4)); env=:align) ==
raw"\begin{align}
1.0 =& 3 \\
2 =& 4
\end{align}
"

@test latexify(((1.0, 2), (3, 4)); separator = [" =& ", " ∈& "], env = :align) ==
raw"\begin{align}
1.0 =& 3 \\
2 ∈& 4
\end{align}
"

# @test_throws MethodError latexify(rn; bad_kwarg="should error")


text:
isabspath(ARGS[1]) || error("Please use an absolute path to the config file.")
include(ARGS[1])

using Distributed

@info "Checking packages for $([h for (h,_) in worker_list])."

addprocs([(w, 1) for (w,cpu) in worker_list], sshflags=sshflags, exeflags=exeflags)
@everywhere using Pkg

@everywhere global CUSTOM_PACKAGES = [PackageSpec(url="https://github.com/JuliaOpt/JuMP.jl", rev="v0.19-beta"),
                                      PackageSpec(url="https://github.com/englhardt/SVDD.jl", rev=("subsvdd-update")),
                                      PackageSpec(url="https://github.com/englhardt/OneClassActiveLearning.jl", rev="subsvdd-update"),
                                      PackageSpec(url="https://github.com/holtri/OCALPlots.jl", rev="master")]

@everywhere function update_or_install_custom_packages()
    for spec in CUSTOM_PACKAGES
        @info "Installing $(spec.repo.url) on revision $(spec.repo.rev) on $(gethostname())."
        Pkg.add(spec)
    end
end

@info "Bringing packages up to date..."
@info "instantiate"
# @everywhere update_or_install_custom_packages()
# @everywhere Pkg.instantiate()
for wid in workers()
    # remotecall_fetch(update_or_install_custom_packages, wid)
    remotecall_fetch(Pkg.instantiate, wid)
end
@info "Trigger precompile."
@everywhere using SVDD, OneClassActiveLearning, JuMP
if foldl(==, remotecall_fetch.(Pkg.installed, workers()))
    @error "Package update did not work."
else
    @info "All hosts are up to date."
end

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
# This file is based on automatically generated code from gen/wrap.jl

const NC_NAT = 0
const NC_BYTE = 1
const NC_CHAR = 2
const NC_SHORT = 3
const NC_INT = 4
const NC_LONG = NC_INT
const NC_FLOAT = 5
const NC_DOUBLE = 6
const NC_UBYTE = 7
const NC_USHORT = 8
const NC_UINT = 9
const NC_INT64 = 10
const NC_UINT64 = 11
const NC_STRING = 12
const NC_MAX_ATOMIC_TYPE = NC_STRING
const NC_VLEN = 13
const NC_OPAQUE = 14
const NC_ENUM = 15
const NC_COMPOUND = 16
const NC_FIRSTUSERTYPEID = 32

const NC_FILL = 0
const NC_NOFILL = 0x0100
const NC_NOWRITE = 0x0000
const NC_WRITE = 0x0001
const NC_CLOBBER = 0x0000
const NC_NOCLOBBER = 0x0004
const NC_DISKLESS = 0x0008
const NC_MMAP = 0x0010
const NC_64BIT_DATA = 0x0020
const NC_CDF5 = NC_64BIT_DATA
const NC_CLASSIC_MODEL = 0x0100
const NC_64BIT_OFFSET = 0x0200
const NC_LOCK = 0x0400
const NC_SHARE = 0x0800
const NC_NETCDF4 = 0x1000
const NC_MPIIO = 0x2000
const NC_MPIPOSIX = 0x4000
const NC_INMEMORY = 0x8000
const NC_PNETCDF = NC_MPIIO
const NC_FORMAT_CLASSIC = 1
const NC_FORMAT_64BIT_OFFSET = 2
const NC_FORMAT_64BIT = NC_FORMAT_64BIT_OFFSET
const NC_FORMAT_NETCDF4 = 3
const NC_FORMAT_NETCDF4_CLASSIC = 4
const NC_FORMAT_64BIT_DATA = 5
const NC_FORMAT_CDF5 = NC_FORMAT_64BIT_DATA
const NC_FORMATX_NC3 = 1
const NC_FORMATX_NC_HDF5 = 2
const NC_FORMATX_NC4 = NC_FORMATX_NC_HDF5
const NC_FORMATX_NC_HDF4 = 3
const NC_FORMATX_PNETCDF = 4
const NC_FORMATX_DAP2 = 5
const NC_FORMATX_DAP4 = 6
const NC_FORMATX_UNDEFINED = 0
const NC_FORMAT_NC3 = NC_FORMATX_NC3
const NC_FORMAT_NC_HDF5 = NC_FORMATX_NC_HDF5
const NC_FORMAT_NC4 = NC_FORMATX_NC4
const NC_FORMAT_NC_HDF4 = NC_FORMATX_NC_HDF4
const NC_FORMAT_PNETCDF = NC_FORMATX_PNETCDF
const NC_FORMAT_DAP2 = NC_FORMATX_DAP2
const NC_FORMAT_DAP4 = NC_FORMATX_DAP4
const NC_FORMAT_UNDEFINED = NC_FORMATX_UNDEFINED
const NC_SIZEHINT_DEFAULT = 0

const NC_UNLIMITED = Int32(0)
const NC_GLOBAL = -1
const NC_MAX_DIMS = 1024
const NC_MAX_ATTRS = 8192
const NC_MAX_VARS = 8192
const NC_MAX_NAME = 256
const NC_MAX_VAR_DIMS = 1024
const NC_MAX_HDF4_NAME = 64
const NC_ENDIAN_NATIVE = 0
const NC_ENDIAN_LITTLE = 1
const NC_ENDIAN_BIG = 2
const NC_CHUNKED = 0
const NC_CONTIGUOUS = 1
const NC_NOCHECKSUM = 0
const NC_FLETCHER32 = 1
const NC_NOSHUFFLE = 0
const NC_SHUFFLE = 1

const NC_NOERR = 0
const NC2_ERR = -1
const NC_EBADID = -33
const NC_ENFILE = -34
const NC_EEXIST = -35
const NC_EINVAL = -36
const NC_EPERM = -37
const NC_ENOTINDEFINE = -38
const NC_EINDEFINE = -39
const NC_EINVALCOORDS = -40
const NC_EMAXDIMS = -41
const NC_ENAMEINUSE = -42
const NC_ENOTATT = -43
const NC_EMAXATTS = -44
const NC_EBADTYPE = -45
const NC_EBADDIM = -46
const NC_EUNLIMPOS = -47
const NC_EMAXVARS = -48
const NC_ENOTVAR = -49
const NC_EGLOBAL = -50
const NC_ENOTNC = -51
const NC_ESTS = -52
const NC_EMAXNAME = -53
const NC_EUNLIMIT = -54
const NC_ENORECVARS = -55
const NC_ECHAR = -56
const NC_EEDGE = -57
const NC_ESTRIDE = -58
const NC_EBADNAME = -59
const NC_ERANGE = -60
const NC_ENOMEM = -61
const NC_EVARSIZE = -62
const NC_EDIMSIZE = -63
const NC_ETRUNC = -64
const NC_EAXISTYPE = -65
const NC_EDAP = -66
const NC_ECURL = -67
const NC_EIO = -68
const NC_ENODATA = -69
const NC_EDAPSVC = -70
const NC_EDAS = -71
const NC_EDDS = -72
const NC_EDATADDS = -73
const NC_EDAPURL = -74
const NC_EDAPCONSTRAINT = -75
const NC_ETRANSLATION = -76
const NC_EACCESS = -77
const NC_EAUTH = -78
const NC_ENOTFOUND = -90
const NC_ECANTREMOVE = -91
const NC4_FIRST_ERROR = -100
const NC_EHDFERR = -101
const NC_ECANTREAD = -102
const NC_ECANTWRITE = -103
const NC_ECANTCREATE = -104
const NC_EFILEMETA = -105
const NC_EDIMMETA = -106
const NC_EATTMETA = -107
const NC_EVARMETA = -108
const NC_ENOCOMPOUND = -109
const NC_EATTEXISTS = -110
const NC_ENOTNC4 = -111
const NC_ESTRICTNC3 = -112
const NC_ENOTNC3 = -113
const NC_ENOPAR = -114
const NC_EPARINIT = -115
const NC_EBADGRPID = -116
const NC_EBADTYPID = -117
const NC_ETYPDEFINED = -118
const NC_EBADFIELD = -119
const NC_EBADCLASS = -120
const NC_EMAPTYPE = -121
const NC_ELATEFILL = -122
const NC_ELATEDEF = -123
const NC_EDIMSCALE = -124
const NC_ENOGRP = -125
const NC_ESTORAGE = -126
const NC_EBADCHUNK = -127
const NC_ENOTBUILT = -128
const NC_EDISKLESS = -129
const NC_ECANTEXTEND = -130
const NC_EMPI = -131
const NC4_LAST_ERROR = -131
const DIM_WITHOUT_VARIABLE = "This is a netCDF dimension but not a netCDF variable."
const NC_HAVE_NEW_CHUNKING_API = 1
const NC_EURL = NC_EDAPURL
const NC_ECONSTRAINT = NC_EDAPCONSTRAINT

const NC_ENTOOL = NC_EMAXNAME
const NC_EXDR = -32
const NC_SYSERR = -31
const NC_FATAL = 1

const nc_type = Cint

mutable struct nc_vlen_t
    len::Integer
    p::Ptr{Cvoid}
end

const nclong = Cint

function nc_inq_libvers()
    ccall((:nc_inq_libvers, libnetcdf), Ptr{UInt8}, ())
end

function nc_strerror(ncerr::Integer)
    ccall((:nc_strerror, libnetcdf), Ptr{UInt8}, (Cint,), ncerr)
end


text:
"""
$(TYPEDSIGNATURES)

Return [`ToolInfo`](@ref) with tool & version attributes and content.

The content can be one or more well-formed xml elements.
Each are wrapped in a [`PnmlLabel`](@ref).
"""
function parse_toolspecific(node, pntd; kw...)
    nn = nodename(node)
    nn == "toolspecific" || error("element name wrong: $nn")
    EzXML.haskey(node, "tool") || throw(MalformedException("$(nn) missing tool attribute", node))
    EzXML.haskey(node, "version") || throw(MalformedException("$(nn) missing version attribute", node))

    d = PnmlDict(
                :tool    => node["tool"],
                :version => node["version"])

    d[:content] = AnyElement[] #

    foreach(elements(node)) do child
        #TODO: Specialize/verify on tool, version. User supplied?
        #TODO: Register additional tool specific parsers?
        push!(d[:content], anyelement(child, pntd; kw...))
    end
    ToolInfo(d[:tool], d[:version], d[:content], node)
end


text:
import Base: show

const DFBC = DefaultPrognosticFieldBoundaryCondition

bc_str(::FBC)     = "Flux"
bc_str(::PBC)     = "Periodic"
bc_str(::OBC)     = "Open"
bc_str(::VBC)     = "Value"
bc_str(::GBC)     = "Gradient"
bc_str(::ZFBC)    = "ZeroFlux"
bc_str(::DFBC)    = "Default"
bc_str(::Nothing) = "Nothing"

#####
##### BoundaryCondition
#####

print_condition(n::Union{Nothing, Number}) = "$n"
print_condition(A::AbstractArray) = "$(Base.dims2string(size(A))) $(typeof(A))"
print_condition(bf::Union{DiscreteBoundaryFunction, ContinuousBoundaryFunction}) = print_condition(bf.func)

function print_condition(f::Function)
    ms = methods(f).ms
    length(ms) == 1 && return "$(ms[1])"
    return "$(ms)"
end

show(io::IO, bc::BoundaryCondition) =
    print(io, "BoundaryCondition: classification=$(bc_str(bc)), condition=$(print_condition(bc.condition))")

#####
##### FieldBoundaryConditions
#####

Base.summary(fbcs::FieldBoundaryConditions) =
    string("west=$(bc_str(fbcs.west)), ",
           "east=$(bc_str(fbcs.east)), ",
           "south=$(bc_str(fbcs.south)), ",
           "north=$(bc_str(fbcs.north)), ",
           "bottom=$(bc_str(fbcs.bottom)), ",
           "top=$(bc_str(fbcs.top)), ",
           "immersed=$(bc_str(fbcs.immersed))")

show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") =
    string("Oceananigans.FieldBoundaryConditions, with boundary conditions", '\n',
           padding, "├── west: ",     typeof(bcs.west), '\n',
           padding, "├── east: ",     typeof(bcs.east), '\n',
           padding, "├── south: ",    typeof(bcs.south), '\n',
           padding, "├── north: ",    typeof(bcs.north), '\n',
           padding, "├── bottom: ",   typeof(bcs.bottom), '\n',
           padding, "├── top: ",      typeof(bcs.top), '\n',
           padding, "└── immersed: ", typeof(bcs.immersed))

Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs))


text:
import ForwardDiff

import ReverseDiffPrototype
const _RD_ = ReverseDiffPrototype

using JuLIP: vecs

export FDPotential, FDPotential_r, RDPotential_r

# Implementation of a ForwardDiff site potential
# ================================================

function ad_evaluate end

"""
`abstract FDPotential <: SitePotential`

To implement a site potential that uses ForwardDiff.jl for AD, create
a concrete type and overload `JuLIP.Potentials.ad_evaluate`.

Example
```julia
@pot type P1 <: FDPotential
end
JuLIP.Potentials.ad_evaluate{T<:Real}(pot::P1, R::Matrix{T}) =
               sum( exp(-norm(R[:,i])) for i = 1:size(R,2) )
```

Notes:

* For an `FDPotential`, the argument to `ad_evaluate` will be
`R::Matrix` rather than `R::JVecsF`; this is an unfortunate current limitation of
`ForwardDiff`. Hopefully it can be fixed.
* TODO: allow arguments `(r, R)` then use chain-rule to combine them.
"""
abstract FDPotential <: SitePotential

evaluate(pot::FDPotential, r, R) = ad_evaluate(pot, mat(collect(R)))

evaluate_d(pot::FDPotential, r, R) =
   ForwardDiff.gradient( R_ -> ad_evaluate(pot, R_), mat(collect(R)) ) |> vecs

abstract FDPotential_r <: SitePotential

evaluate(pot::FDPotential_r, r, R) = ad_evaluate(pot, collect(r))

function evaluate_d(pot::FDPotential_r, r, R)
   d = ForwardDiff.gradient( r_ -> ad_evaluate(pot, r_), collect(r) )
   return [ (d[i]/r[i]) * R[i] for i = 1:length(r) ]
end

# Implementation of a ReverseDiffPrototype site potential
# ========================================================

abstract RDPotential_r <: SitePotential

evaluate(pot::RDPotential_r, r, R) = ad_evaluate(pot, collect(r))

function evaluate_d(pot::RDPotential_r, r, R)
   d = _RD_.gradient( r_ -> ad_evaluate(pot, r_), collect(r) )
   return [ (d[i]/r[i]) * R[i] for i = 1:length(r) ]
end

# # Implementation of a ReverseDiffSource site potential
# # ========================================================
#
# abstract RDSPotential_r <: SitePotential
#
# init(pot::RDSPotential_r,
#
# evaluate(pot::RDPotential_r, r, R) = ad_evaluate(pot, collect(r))
#
# function evaluate_d(pot::RDPotential_r, r, R)
#    d = _RD_.gradient( r_ -> ad_evaluate(pot, r_), collect(r) )
#    return [ (d[i]/r[i]) * R[i] for i = 1:length(r) ]
# end


text:
"""
	regress(cmd0::String="", arg1=nothing, kwargs...)

Linear regression of 1-D data sets.

Full option list at [`regress`]($(GMTdoc)gmtregress.html)

Parameters
----------

- **A** | **all_slopes** :: [Type => Str | List]        ``Arg = min/max/inc``

    Instead of determining a best-fit regression we explore the full range of regressions.
    ($(GMTdoc)gmtregress.html#a)
- **C** | **confidence_level** :: [Type => Int]      ``Arg = level``

    Set the confidence level (in %) to use for the optional calculation of confidence
    bands on the regression [95].
    ($(GMTdoc)gmtregress.html#c)
- **E** | **regression_type** :: [Type => Str]   ``Arg = x|y|o|r``

    Type of linear regression, i.e., select the type of misfit we should calculate.
    ($(GMTdoc)gmtregress.html#e)
- **F** | **column_combination** :: [Type => Str]   ``Arg = x|y|m|l|c``

    Append a combination of the columns you wish returned;
    ($(GMTdoc)gmtregress.html#f)
- **N** | **norm** :: [Type => Str | Int]          ``Arg = 1|2|r|w``

    Selects the norm to use for the misfit calculation.
    ($(GMTdoc)gmtregress.html#n)
- **S** | **restrict** :: [Type => Str | []]        ``Arg = [r]``

    Restricts which records will be output.
    ($(GMTdoc)gmtregress.html#s)
- **T** | **equi_space** :: [Type => Str | List]     ``Arg = [min/max/]inc[+a|n]] or file|list``

    Evaluate the best-fit regression model at the equidistant points implied by the arguments.
    ($(GMTdoc)gmtregress.html#t)
- **W** | **weighted** :: [Type => Str | []]     ``Arg = [w][x][y][r]``

    Specifies weighted regression and which weights will be provided.
    ($(GMTdoc)gmtregress.html#w)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function regress(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtregress", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:V_params :b :d :e :g :h :i :o :yx])
    cmd = parse_these_opts(cmd, d, [[:A :all_slopes], [:C :confidence_level], [:E :regression_type], [:N :norm],
                [:F :column_combination], [:S :restrict], [:T :equi_space], [:W :weighted]])

	common_grd(d, cmd0, cmd, "gmtregress ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
regress(arg1, cmd0::String=""; kw...) = regress(cmd0, arg1; kw...)

gmtregress = regress 		# Alias


text:
function na_to_zero(input_value)
	if isna(input_value)
		return 0
	else
		return input_value
	end
end

function dict_to_df(input_dict, ID)
	new_df = DataFrame(PFAM_Model=collect(keys(input_dict)))
	new_df[symbol(ID)] = convert(Array{Int64},collect(values(input_dict)))
	return new_df
end

function read_hmm_tabular_output(filename, ID,seqinfo)
	hmm_output = open(filename)
	matched_sequences_dict = Dict()
	for line in eachline(hmm_output)
		if startswith(line, ['#'])
			continue
		else
			split_line = split(line)
			family = split_line[3]
			best_domain_evalue = float(split_line[8])
			target_sequence_name = join(split(split_line[1], '.')[1:end-1],'.')
			translation_frame = split(split_line[1], '.')[end]
			if haskey(matched_sequences_dict, target_sequence_name)
				if best_domain_evalue < matched_sequences_dict[target_sequence_name][1]
					matched_sequences_dict[target_sequence_name] = [best_domain_evalue,family,translation_frame]
				end
			else
				matched_sequences_dict[target_sequence_name] = [best_domain_evalue,family,translation_frame]
			end
		end
	end
	Pfam_dict = Dict()
	seq_df = DataFrame()
	for sequence in matched_sequences_dict
		if seqinfo
			newline_seq_df = DataFrame()
		end
		family = sequence[2][2]
		if haskey(Pfam_dict, family)
			Pfam_dict[family] = Pfam_dict[family] + 1
			if seqinfo
				newline_seq_df = DataFrame(seqID = sequence[1], translation = sequence[2][3], model = family, evalue = sequence[2][1])
			end
		else
			Pfam_dict[family] = 1
			if seqinfo
				newline_seq_df = DataFrame(seqID = sequence[1], translation = sequence[2][3], model = family, evalue = sequence[2][1])
			end
		end
		if seqinfo
			seq_df = vcat(seq_df, newline_seq_df)
		end
	end

	output_dict = Dict()
	output_dict["pfam_df"] = dict_to_df(Pfam_dict, ID)
	output_dict["seq_df"] = seq_df

	return output_dict
end

function add_newsample_to_pfam_table(filenames, IDs, pfam_table,seqinfo)
	new_pfam_table = pfam_table
	full_seq_df = DataFrame()
	for (i, filename) in enumerate(filenames)
		print("Reading "*filename*"\n")
		column_to_add = read_hmm_tabular_output(filename, IDs[i],seqinfo)
		new_pfam_table = join(column_to_add["pfam_df"], new_pfam_table, on = :PFAM_Model, kind = :outer)
		if seqinfo
			seq_df = column_to_add["seq_df"]
			seq_df[:Sample] = fill(IDs[i], size(seq_df,1))
			full_seq_df = vcat(full_seq_df, seq_df)
		end
	end
#	for i in 1:size(new_pfam_table,1)
#		for j in 1:size(new_pfam_table,2)
#			if isna(new_pfam_table[i,j])
#				new_pfam_table[i,j] = 0
#			end
#		end
#	end
	for column in names(new_pfam_table)
		for (i, value) in enumerate(new_pfam_table[column])
			new_pfam_table[column][i] = na_to_zero(value)
		end
	end
	output_dict = Dict()
	output_dict["new_pfam_table"] = new_pfam_table
	output_dict["seq_df"] = full_seq_df
	return output_dict
end

function tabdict(filename)
	desc_table = Dict()
	for line in eachline(open(filename))
		desc_table[split(line, '\t')[1]] = strip(split(line, '\t')[2],'\n')
	end

	return desc_table
end

function reads_to_rpob_equiv(pfam_df, pfam_db_stats)
	#Get the names of the PFAM models for RpoB
	RpoB_models = ["RNA_pol_Rpb2_1","RNA_pol_Rpb2_2","RNA_pol_Rpb2_3","RNA_pol_Rpb2_4","RNA_pol_Rpb2_45","RNA_pol_Rpb2_5","RNA_pol_Rpb2_6","RNA_pol_Rpb2_7"]

	#Get their lengths
	RpoB_model_lengths = Array(Float64,0)
	for model in RpoB_models
		append!(RpoB_model_lengths, [float(pfam_db_stats[pfam_db_stats[:NAME] .== model, :LENG][1])])
	end

	#Calculate RpoB reads per aa for each sample
	RpoB_reads_per_aa = Array(Float64,0)
	for sample in names(pfam_df)
		RpoB_reads = Array(Float64,0)
		if sample != :PFAM_Model
			for RpoB_model in RpoB_models
				append!(RpoB_reads,float(pfam_df[pfam_df[:PFAM_Model] .== RpoB_model, sample]))
			end
			if length(RpoB_reads) != length(RpoB_model_lengths)
				print("It looks like there's some RpoB models missing from your dataset.")
			end

			data_for_linear_regression = DataFrame(reads=RpoB_reads,aa=RpoB_model_lengths)
			append!(RpoB_reads_per_aa,[coef(lm(reads~aa,data_for_linear_regression))[2]])
		end
	end

	#Calculate RpoB equivalents for every PFAM model in your dataset
	rpoB_table_output = DataFrame(PFAM_Model = [])
	for colnum in 2:size(pfam_df)[2]
		rpoB_table_output[names(pfam_df)[colnum]] = []
	end
	for rownum in 1:size(pfam_df)[1]
		model_name = pfam_df[rownum,:PFAM_Model]
		readnums = Array(pfam_df[rownum,2:size(pfam_df)[2]])
		if 	model_name in pfam_db_stats[:NAME]
			pfam_length = float(pfam_db_stats[pfam_db_stats[:NAME] .== model_name, :LENG][1])
		else
			error("At least one of the models in your hmmsearch output wasn't in the specified HMM database. Are you sure that the same HMM database was specified for both the run_hmmsearch() step and the new_project() step?\n")
		end
		reads_per_aa = []
		for readnum in readnums
			push!(reads_per_aa, readnum / pfam_length)
		end
		RpoB_equivalents = []
		for (i, reads) in enumerate(reads_per_aa)
			push!(RpoB_equivalents, reads / RpoB_reads_per_aa[i])
		end
		newrow = DataFrame(PFAM_Model = model_name)
		for colnum in 2:size(pfam_df)[2]
			newrow[names(pfam_df)[colnum]] = RpoB_equivalents[colnum-1]
		end
		rpoB_table_output = vcat(rpoB_table_output,newrow)
	end
	return rpoB_table_output
end

#function import_hmm_output(project_table)
#	blank_table = DataFrame(PFAM_Model=[])
#	pfam_table = add_newsample_to_pfam_table(project_table[:Filename],project_table[:SampleID],blank_table)
#	output_dict = Dict()
#	output_dict["pfam_table"] = pfam_table["new_pfam_table"]
#	output_dict["seq_df"] = pfam_table[seq_df]
#	return output_dict
#end

function new_project(;project_table_file="apfamet_project_table.txt", hmm_database="", seqinfo=true)
	project_table = readtable(project_table_file, separator='\t')

	blank_table = DataFrame(PFAM_Model=[])
	pfam_table = add_newsample_to_pfam_table(project_table[:Filename],project_table[:SampleID],blank_table,seqinfo)
	read_counts_table = pfam_table["new_pfam_table"]
	seq_info_table = pfam_table["seq_df"]

	hmm_database_info = prepare_pfam_database(hmm_database)
	rpoB_equiv_table = reads_to_rpob_equiv(read_counts_table, hmm_database_info)
	return Dict("project_table"=>project_table,"read_counts_table"=>read_counts_table,"hmm_database_info"=>hmm_database_info,"rpoB_equiv_table"=>rpoB_equiv_table, "seq_info_table"=>seq_info_table)
end

function save_project(project; base_filename="apfamet_project", overwrite=false)
	if isfile(base_filename)
		if overwrite == false
			error("Project "*base_filename*" already exists. Please select a new filename or set overwrite=true")
		else
			print("Overwriting previous project "*base_filename*" because overwrite=true.")
		end
	end
	writetable(base_filename*".apfamet_project_table", project["project_table"], separator='\t')
	writetable(base_filename*".apfamet_read_counts_table", project["read_counts_table"], separator='\t')
	writetable(base_filename*".apfamet_hmm_database_info", project["hmm_database_info"], separator='\t')
	writetable(base_filename*".apfamet_rpoB_equiv_table", project["rpoB_equiv_table"], separator='\t')
	writetable(base_filename*".apfamet_seq_info_table", project["seq_info_table"], separator='\t')
end

function load_project(base_filename)
	project_table = readtable(base_filename*".apfamet_project_table", separator='\t')
	read_counts_table = readtable(base_filename*".apfamet_read_counts_table", separator='\t')
	hmm_database_info = readtable(base_filename*".apfamet_hmm_database_info", separator='\t')
	rpoB_equiv_table = readtable(base_filename*".apfamet_rpoB_equiv_table", separator='\t')
	seq_info_table = readtable(base_filename*".apfamet_seq_info_table", separator='\t')
	return Dict("project_table"=>project_table,"read_counts_table"=>read_counts_table,"hmm_database_info"=>hmm_database_info,"rpoB_equiv_table"=>rpoB_equiv_table, "seq_info_table"=>seq_info_table)
end

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
mutable struct BackoffSchedulerEntry
    capacity::Int
    fuel::Int
    bantime::Int
    banremaining::Int
end

isbanned(d::BackoffSchedulerEntry) = d.banremaining > 0

"""
A Rewrite Scheduler that implements exponential rule backoff.
For each rewrite, there exists a configurable initial match limit.
If a rewrite search yield more than this limit, then we ban this rule
for number of iterations, double its limit, and double the time it
will be banned next time.

This seems effective at preventing explosive rules like
associativity from taking an unfair amount of resources.
"""
struct BackoffScheduler <: AbstractScheduler
    data::Dict{Rule, BackoffSchedulerEntry}
    G::EGraph
    theory::Vector{Rule}
end

shouldskip(s::BackoffScheduler, r::Rule) = s.data[r].banremaining > 0

function BackoffScheduler(g::EGraph, theory::Vector{Rule})
    BackoffScheduler(g, theory, 8, 2)
end

function BackoffScheduler(G::EGraph, theory::Vector{Rule}, fuel::Int, bantime::Int)
    gsize = length(G.U)
    data = Dict{Rule, BackoffSchedulerEntry}()

    # These numbers seem to fit
    for rule ∈ theory
        data[rule] = BackoffSchedulerEntry(fuel, fuel, bantime, 0)
    end

    return BackoffScheduler(data, G, theory)
end

# can saturate if there's no banned rule
cansaturate(s::BackoffScheduler) = all(kv -> !isbanned(last(kv)), s.data)

function readstep!(s::BackoffScheduler)
    for rule ∈ s.theory
        rd = s.data[rule]
        if rd.banremaining > 0
            rd.banremaining -= 1

            if rd.banremaining == 0
                # @info "unbanning rule" rule
                rd.bantime *= 2
                rd.capacity *= 2
                rd.fuel = rd.capacity
            end
        end
    end
end

function writestep!(s::BackoffScheduler, rule::Rule)
    rd = s.data[rule]

    # decrement fuel, ban rule if fuel is empty
    rd.fuel -= 1
    if rd.fuel == 0
        # @info "banning rule $rule for $(rd.bantime)!"
        rd.banremaining = rd.bantime
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ 67ce08a4-9a04-11eb-2800-d5325a874022
begin
	using Revise
	using Pkg
	Pkg.activate("../..")
	Pkg.instantiate()
	using ACCESS
	using CSV
	using DelimitedFiles
	using StatsPlots
	using Measurements
	using CCDReduction
	using PlutoUI
	using BenchmarkTools
	using Glob, FITSIO, DataFramesMeta
	using Statistics
	using Unitful, UnitfulAstro
	using MultivariateStats, LinearAlgebra, StatsBase
end

# ╔═╡ 94595695-b439-4032-924d-fba32a481199
TableOfContents()

# ╔═╡ 67366046-900f-43c9-9ede-359bca71ef55
base_dir = "/home/mango/data/data_detrending/WASP50b/out_l/WASP50/w50_131219/white-light"

# ╔═╡ 670ad85c-6ab7-4652-ab78-ea83b2993927
lc_path = "$base_dir/lc.dat"

# ╔═╡ 4a438e18-edb9-4992-b859-70a5bf31834b
BMA_results_path = "$base_dir/results.dat"

# ╔═╡ eef8fa8c-4989-4721-bcc9-42167bbfd325
comps_path = "$base_dir/comps.dat"

# ╔═╡ b33bdd2a-d4c6-44ca-ab14-fc67552e3f2c
eparams_path = "$base_dir/../eparams.dat"

# ╔═╡ 099a0f42-70c6-499d-a7dd-fb6d1d8832f4
md"""
## Raw data
"""

# ╔═╡ 763e1621-7922-4bc7-ad7d-6d6888ed778d
lc_data = CSV.File(lc_path, header=["t", "f", "f_index"])

# ╔═╡ ea2d99ac-feb6-4b3a-bee0-72c60c070755
begin
	idx = lc_data.f_index .== 0;
	t, f = lc_data.t[idx], lc_data.f[idx]
end

# ╔═╡ 92bb0deb-d2b9-4826-ad50-6a4c11a86c59
md"""
## External params
"""

# ╔═╡ 13cb7356-6e36-4565-b586-f05e3fab4b2f
eparams = CSV.read(eparams_path, DataFrame, normalizenames=true)

# ╔═╡ e85b316c-f0a3-4892-964b-53e7577cac25
X = norm_data(eparams, dims=1)

# ╔═╡ a119f4b6-ab36-493b-b732-a055f09cdc54
md"""
## Comp stars
"""

# ╔═╡ c04b2b40-d3dd-47ed-bee7-a4e42ece3f96
comps = CSV.read(comps_path, DataFrame, header=false)

# ╔═╡ c4f1f656-ce41-4319-a394-cae2094fdd14
Xc = norm_data(comps, dims=1)

# ╔═╡ 51aebd47-156c-4bf7-b991-f1c4ff0ee8b5
plot_params(A, A_names) = plot(A, label=reshape(names(A_names), (1, :)))

# ╔═╡ 29a716f5-740a-4e5e-abea-f902b6d195f0
plot_params(X, eparams)

# ╔═╡ d621bf6d-c6ab-4626-89fd-6f5d5ac68a47
plot_params(Xc, comps)

# ╔═╡ 519cd15e-e6a0-4d7f-b9c7-40be545fa373
md"""
## PCA
"""

# ╔═╡ c52b38fb-bd77-4176-a570-0923aa99cc07
V, S, A = classic_PCA(Xc)

# ╔═╡ 83760c7c-fe94-4b4b-a7c9-ca20c051a973
# 	# Load pickle
# 	gpts_pkl_to_dict = load_pickle(

# )
# 	# Show keys
# 	println(keys(gpts_pkl_to_dict))

# 	# Load specified params
# 	const gpts_params = ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"]
# 	gpts_params_dict = sub_dict(
# 		gpts_pkl_to_dict, gpts_params
# 	)

# 	# Add desired offsets
# 	gpts_params_dict[:t0] .-= 2.45485e6; gpts_params_dict

# 	# Load into DataFrame
# 	df_gpts = DataFrame(gpts_params_dict)

# ╔═╡ 030a5a05-8009-4a3d-8e83-a1979fadf9f4
import DarkMode

# ╔═╡ 23987d7b-187d-47d8-a3d3-99c1a9fa8b8c
#DarkMode.enable()

# ╔═╡ be317a00-2c87-4b38-8ab5-fb4cf7c32bf2
begin
	const pal = palette(
		[
			"#f5d300", # yellow
			"#d55e00", # orange
			"#08f7fe", # cyan
			"#0173b2", # blue
			"#b2df8a", # light green
			"#029e73", # green
		]
	)
	theme(:dark)
	default(
		titlefont = "Lato",
		guidefont = "Lato",
		markerstrokewidth = 0,
		palette = pal,
		dpi = 250,
		fmt = :png,
	)
end

# ╔═╡ Cell order:
# ╟─94595695-b439-4032-924d-fba32a481199
# ╠═67366046-900f-43c9-9ede-359bca71ef55
# ╠═670ad85c-6ab7-4652-ab78-ea83b2993927
# ╠═4a438e18-edb9-4992-b859-70a5bf31834b
# ╠═eef8fa8c-4989-4721-bcc9-42167bbfd325
# ╠═b33bdd2a-d4c6-44ca-ab14-fc67552e3f2c
# ╟─099a0f42-70c6-499d-a7dd-fb6d1d8832f4
# ╠═763e1621-7922-4bc7-ad7d-6d6888ed778d
# ╠═ea2d99ac-feb6-4b3a-bee0-72c60c070755
# ╟─92bb0deb-d2b9-4826-ad50-6a4c11a86c59
# ╠═13cb7356-6e36-4565-b586-f05e3fab4b2f
# ╠═e85b316c-f0a3-4892-964b-53e7577cac25
# ╠═29a716f5-740a-4e5e-abea-f902b6d195f0
# ╟─a119f4b6-ab36-493b-b732-a055f09cdc54
# ╠═c04b2b40-d3dd-47ed-bee7-a4e42ece3f96
# ╠═c4f1f656-ce41-4319-a394-cae2094fdd14
# ╠═d621bf6d-c6ab-4626-89fd-6f5d5ac68a47
# ╠═51aebd47-156c-4bf7-b991-f1c4ff0ee8b5
# ╟─519cd15e-e6a0-4d7f-b9c7-40be545fa373
# ╠═c52b38fb-bd77-4176-a570-0923aa99cc07
# ╠═83760c7c-fe94-4b4b-a7c9-ca20c051a973
# ╠═67ce08a4-9a04-11eb-2800-d5325a874022
# ╠═030a5a05-8009-4a3d-8e83-a1979fadf9f4
# ╠═23987d7b-187d-47d8-a3d3-99c1a9fa8b8c
# ╠═be317a00-2c87-4b38-8ab5-fb4cf7c32bf2


text:
using Base.Iterators: product

"""
    RectilinearGrid{Tl, Tr} <: AbstractVector{Tuple{Tl, Tr}}

A `RectilinearGrid` is parametrised by `AbstractVector`s of points `xl` and `xr`, whose
element types are `Tl` and `Tr` respectively, comprising `length(xl) * length(xr)`
elements. Linear indexing is the same as `product(eachindex(xl), eachindex(xr))` - `xl`
iterates more quickly than `xr`.
"""
struct RectilinearGrid{
    Tl, Tr, Txl<:AbstractVector{Tl}, Txr<:AbstractVector{Tr},
} <: AbstractVector{Tuple{Tl, Tr}}
    xl::Txl
    xr::Txr
end

get_space(x::RectilinearGrid) = Zygote.literal_getfield(x, Val(:xl))

get_times(x::RectilinearGrid) = Zygote.literal_getfield(x, Val(:xr))

Base.size(X::RectilinearGrid) = (length(X.xl) * length(X.xr),)

function Base.collect(X::RectilinearGrid{Tl, Tr}) where{Tl, Tr}
    return vec(
        map(
            ((p, q),) -> (X.xl[p], X.xr[q]),
            product(eachindex(X.xl), eachindex(X.xr)),
        )
    )
end

function Base.getindex(X::RectilinearGrid, n::Integer)
    return (X.xl[mod(n - 1, length(X.xl)) + 1], X.xr[div(n - 1, length(X.xl)) + 1])
end

Base.show(io::IO, x::RectilinearGrid) = Base.show(io::IO, collect(x))

"""
    SpaceTimeGrid{Tr, Tt<:Real}

A `SpaceTimeGrid` is a `RectilinearGrid` in which the left vector corresponds to space, and
the right `time`. The left eltype is arbitrary, but the right must be `Real`.
"""
const SpaceTimeGrid{Tr, Tt<:Real} = RectilinearGrid{
    Tr, Tt, <:AbstractVector{Tr}, <:AbstractVector{Tt},
}

#
# Implement internal API for transforming between "flat" representation, which is useful for
# GPs, an a time-centric representation, which is useful for state-space models.
#

# See docstring elsewhere for context.
function inputs_to_time_form(x::SpaceTimeGrid)
    return Fill(get_space(x), length(get_times(x)))
end

# See docstring elsewhere for context.
function merge_inputs(x1::SpaceTimeGrid, x2::SpaceTimeGrid)
    if get_space(x1) != get_space(x2)
        throw(error("Space coords of inputs not compatible, cannot merge."))
    end
    return RectilinearGrid(get_space(x1), vcat(get_times(x1), get_times(x2)))
end

# See docstring elsewhere for context.
function sort_in_time(x::SpaceTimeGrid)
    idx = sortperm(get_times(x))
    return idx, RectilinearGrid(get_space(x), get_times(x)[idx])
end

# See docstring elsewhere for context.
function observations_to_time_form(x::SpaceTimeGrid, y::AbstractVector{<:Union{Real, Missing}})
    return restructure(y, Fill(length(get_space(x)), length(get_times(x))))
end

function observations_to_time_form(x::SpaceTimeGrid, ::AbstractVector{Missing})
    return fill(missing, length(get_times(x)))
end

function get_zeros(x::SpaceTimeGrid{T}) where {T<:Real}
    return fill(Diagonal(zeros(T, length(get_space(x)))), length(get_times(x)))
end

# See docstring elsewhere for context.
function noise_var_to_time_form(x::RectilinearGrid, S::Diagonal{<:Real})
    vs = restructure(
        diag(S),
        Zygote.ignore() do
            Fill(length(get_space(x)), length(get_times(x)))
        end,
    )
    return zygote_friendly_map(v -> Diagonal(collect(v)), vs)
end

destructure(::RectilinearGrid, y::AbstractVector) = reduce(vcat, y)


text:
import .LibCImPlot: IsPlotHovered, GetPlotMousePos, GetPlotLimits, IsPlotQueried, GetPlotQuery
import .LibCImPlot: SetNextPlotLimits, SetNextPlotLimitsX, SetNextPlotLimitsY, SetPlotYAxis
import .LibCImPlot: ShowColormapScale, PushPlotClipRect, PopPlotClipRect

function SetNextPlotTicksX(values::Vector{<:Real}, n_ticks::Integer;
                           labels::Vector{String} = [""], show_default::Bool = false)

    eltype(values) !== Float64 && (values = Float64.(values))

    LibCImPlot.SetNextPlotTicksXdoublePtr(values, Cint(n_ticks), labels, show_default)
end

function SetNextPlotTicksX(x_min, x_max, n_ticks::Integer;
                           labels::Vector{String} = [""], show_default::Bool = false)

    typeof(x_min) !== Float64 && (x_min = Float64(x_min))
    typeof(x_max) !== Float64 && (x_max = Float64(x_max))

    LibCImPlot.SetNextPlotTicksXdouble(x_min, x_max, Cint(n_ticks), labels, show_default)
end

function SetNextPlotTicksY(values::Vector{<:Real}, n_ticks::Integer;
                           labels::Vector{String} = [""], show_default::Bool = false,
                           y_axis::Integer = 0)

    eltype(values) !== Float64 && (values = Float64.(values))

    LibCImPlot.SetNextPlotTicksYdoublePtr(values, Cint(n_ticks), labels, show_default,
                                          Cint(y_axis))
end

function SetNextPlotTicksY(y_min, y_max, n_ticks::Integer;
                           labels::Vector{String} = [""], show_default::Bool = false,
                           y_axis::Integer = 0)

    typeof(y_min) !== Float64 && (y_min = Float64(y_min))
    typeof(y_max) !== Float64 && (y_max = Float64(y_max))

    LibCImPlot.SetNextPlotTicksYdouble(y_min, y_max, Cint(n_ticks), labels, show_default,
                                       Cint(y_axis))
end

function GetPlotPos()
    out = Ref{ImVec2}()
    LibCImPlot.GetPlotPos(out)
    return out[]
end

function GetPlotSize()
    out = Ref{ImVec2}()
    LibCImPlot.GetPlotSize(out)
    return out[]
end

function PlotToPixels(plt::ImPlotPoint, y_axis::Integer = -1)
    out = Ref{ImVec2}()
    LibCImPlot.PlotToPixels(out, plt, Cint(y_axis))
    return out[]
end


text:
## full test to run a GASS optimization with the output
##

rootdir = "/home/stephane/Science/ALMA/ArrayConfig/GASS"

push!(LOAD_PATH,"$rootdir/master/src")
using GASS

## directory
wdir    = "$rootdir/products"

cd(wdir)

function main(inpfile)
    println("# GASS subarray optimization")
    println("## Parameter file: $inpfile")
    cfg = read_cfg(inpfile)

    res= full_run("testing" , cfg , true)
end

pop= @time main("../master/data/GA_Parameters_testing.txt.julia")

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
using Statistics

function pooling!(
    cLayer::OneD,
    Ai::AbstractArray{T,3},
    Ao::AbstractArray{T,3},
) where {OneD<:Union{MaxPool1D,AveragePool1D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, c, m = size(Aip)
    s_H = cLayer.s
    f_H = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1

    if cLayer isa MaxPoolLayer
        pool = maximum
    else
        pool = mean
    end #if cLayer isa MaxPoolLayer

    # @simd
    for mi = 1:m
        # @simd for
        Threads.@threads for ci = 1:c
            # @simd
            for hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                @inbounds ai = view(Aip,h_start:h_end, ci, mi)
                @inbounds Ao[hi, ci, mi] = pool(ai)
                # ai = nothing
                # Base.GC.gc()
            end #for
        end #for ci=1:c
    end #for mi=1:m, ci=1:c

    return nothing
end #function pooling!(cLayer::OneD

function pooling!(
    cLayer::TwoD,
    Ai::AbstractArray{T,4},
    Ao::AbstractArray{T,4},
) where {TwoD<:Union{MaxPool2D,AveragePool2D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, c, m = size(Aip)
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    if cLayer isa MaxPoolLayer
        pool = maximum
    else
        pool = mean
    end #if cLayer isa MaxPoolLayer
    # @simd
    for mi = 1:m
        # @simd for
        Threads.@threads for ci = 1:c
            # @simd for

                # @simd
                for wi = 1:n_W, hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                    w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                    @inbounds ai =
                        view(Aip, h_start:h_end, w_start:w_end, ci, mi)
                    @inbounds Ao[hi, wi, ci, mi] = pool(ai)
                    # ai = nothing
                    # Base.GC.gc()
                end #for hi=1:n_H
        #     end #for wi=1:n_W
        end #for ci=1:c
    end #for mi=1:m,
    return nothing
end #function pooling!(cLayer::TwoD,

function pooling!(
    cLayer::ThreeD,
    Ai::AbstractArray{T,5},
    Ao::AbstractArray{T,5},
) where {ThreeD<:Union{MaxPool3D,AveragePool3D},T}

    Aip = padding(cLayer, Ai)
    n_Hi, n_Wi, n_Di, c, m = size(Aip)
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    if cLayer isa MaxPoolLayer
        pool = maximum
    else
        pool = mean
    end #if cLayer isa MaxPoolLayer

    # @simd
    for mi = 1:m
        # @simd for
        Threads.@threads for ci = 1:c
            # @simd for

                # @simd for

                    # @simd
                    for di = 1:n_D, wi = 1:n_W, hi = 1:n_H
                        h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                        h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                        w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                        w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                        d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                        d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                        @inbounds ai = view(Aip,
                            h_start:h_end,
                            w_start:w_end,
                            d_start:d_end,
                            ci,
                            mi,
                        )
                        @inbounds Ao[hi, wi, di, ci, mi] = pool(ai)
                        # ai = nothing
                        # Base.GC.gc()
                    end #for hi=1:n_H
        #         end #for wi=1:n_W
        #     end #for di=1:n_D
        end #for ci=1:c
    end #for mi=1:m, ci=1:c

    return nothing
end #function pooling!(cLayer::ThreeD,

export pooling!

###dpooling

function dpooling!(
    cLayer::OneD,
    Ai::AbstractArray{T,3},
    dAi::AbstractArray{T,3},
    Ao::AbstractArray{T,3},
    dAo::AbstractArray{T,3},
) where {OneD<:Union{MaxPool1D,AveragePool1D},T}

    padS = paddingSize(cLayer, Ai)
    paddedS = paddedSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    dAip = zeros(eltype(Ai), paddedS)
    n_Hi, c, m = paddedS
    s_H = sS = cLayer.s
    f_H = fS = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1

    # if sS == fS #to use the @simd the parallele the operation
    #     # @simd
    #     for mi = 1:m,
    #         # @simd for
    #         ci = 1:c
    #             # @simd
    #             Threads.@threads for hi = 1:n_H
    #                 h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
    #                 h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
    #                 @inbounds ai = Aip[h_start:h_end, ci, mi]
    #                 if cLayer isa MaxPoolLayer
    #                     mask = (Ao[hi,ci,mi] .== ai)
    #                     @inbounds dAip[h_start:h_end, ci, mi] .+=
    #                         (dAo[hi, ci, mi] .* mask)
    #                 else
    #                     mask = similar(ai)
    #                     mask .= 1 / prod(size(ai))
    #                     @inbounds dAip[h_start:h_end, ci, mi] .+=
    #                         (dAo[hi, ci, mi] .* mask)
    #                 end #if cLayer isa MaxPoolLayer
    #             end #for hi=1:n_H
    #         # end #for ci=1:c
    #     end #for mi=1:m
    # else
        # @simd
        for mi = 1:m
            # @simd
            Threads.@threads for ci = 1:c
                for hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    @inbounds ai = Aip[h_start:h_end, ci, mi]
                    if cLayer isa MaxPoolLayer
                        mask = (Ao[hi,ci,mi] .== ai)
                        @inbounds dAip[h_start:h_end, ci, mi] .+=
                            (dAo[hi, ci, mi] .* mask)
                    else
                        mask = similar(ai)
                        mask .= 1 / prod(size(ai))
                        @inbounds dAip[h_start:h_end, ci, mi] .+=
                            (dAo[hi, ci, mi] .* mask)
                    end #if cLayer isa MaxPoolLayer
                end #for hi=1:n_H
            end #for ci=1:c
        end #for mi=1:m
    # end #if sS == fS

    dAi .= dAip[1+padS[1]:end-padS[2], :, :]

    return nothing
end #function dpooling!(cLayer::OneD

function dpooling!(
    cLayer::TwoD,
    Ai::AbstractArray{T,4},
    dAi::AbstractArray{T,4},
    Ao::AbstractArray{T,4},
    dAo::AbstractArray{T,4},
) where {TwoD<:Union{MaxPool2D,AveragePool2D},T}

    padS = paddingSize(cLayer, Ai)
    paddedS = paddedSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    dAip = zeros(eltype(Ai), paddedS)
    n_Hi, n_Wi, c, m = paddedS
    s_H, s_W = sS = cLayer.s
    f_H, f_W = fS = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    # if sS == fS
    #     # @simd
    #     for mi = 1:m,
    #         # @simd for
    #         ci = 1:c,
    #             # @simd for
    #             wi = 1:n_W
    #                 # @simd
    #                 Threads.@threads for hi = 1:n_H
    #                     h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
    #                     h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
    #                     w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
    #                     w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
    #                     @inbounds ai = Aip[h_start:h_end, w_start:w_end, ci, mi]
    #                     if cLayer isa MaxPoolLayer
    #                         mask = (Ao[hi,wi,ci,mi] .== ai)
    #                         @inbounds dAip[
    #                             h_start:h_end,
    #                             w_start:w_end,
    #                             ci,
    #                             mi,
    #                         ] .+= (dAo[hi, wi, ci, mi] .* mask)
    #                     else
    #                         mask = similar(ai)
    #                         mask .= 1 / prod(size(ai))
    #                         @inbounds dAip[
    #                             h_start:h_end,
    #                             w_start:w_end,
    #                             ci,
    #                             mi,
    #                         ] .+= (dAo[hi, wi, ci, mi] .* mask)
    #                     end #if cLayer isa MaxPoolLayer
    #                 end #for hi=1:n_H
    #         #     end #for wi=1:n_W
    #         # end #for ci=1:c
    #     end #for mi=1:m
    # else
        # @simd
        for mi = 1:m
            # @simd
            Threads.@threads for ci = 1:c
                for wi = 1:n_W, hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                    w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                    @inbounds ai = Aip[h_start:h_end, w_start:w_end, ci, mi]
                    if cLayer isa MaxPoolLayer
                        mask = (Ao[hi,wi,ci,mi] .== ai)
                        @inbounds dAip[h_start:h_end, w_start:w_end, ci, mi] .+=
                            (dAo[hi, wi, ci, mi] .* mask)
                    else
                        mask = similar(ai)
                        mask .= 1 / prod(size(ai))
                        @inbounds dAip[h_start:h_end, w_start:w_end, ci, mi] .+=
                            (dAo[hi, wi, ci, mi] .* mask)
                    end #if cLayer isa MaxPoolLayer
                end #for wi=1:n_W, hi=1:n_H
            end #for ci=1:c
        end #for mi=1:m
    # end #if sS == fS

    dAi .= dAip[1+padS[1]:end-padS[2], 1+padS[3]:end-padS[4], :, :]

    return nothing
end #function dpooling!(cLayer::TwoD,

function dpooling!(
    cLayer::ThreeD,
    Ai::AbstractArray{T,5},
    dAi::AbstractArray{T,5},
    Ao::AbstractArray{T,5},
    dAo::AbstractArray{T,5},
) where {ThreeD<:Union{MaxPool3D,AveragePool3D},T}

    padS = paddingSize(cLayer, Ai)
    paddedS = paddedSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    dAip = zeros(eltype(Ai), paddedS)
    n_Hi, n_Wi, n_Di, c, m = paddedS
    s_H, s_W, s_D = sS = cLayer.s
    f_H, f_W, f_D = fS = cLayer.f

    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    # if sS == fS
    #     # @simd
    #     for mi = 1:m,
    #         # @simd for
    #         ci = 1:c,
    #             # @simd for
    #             di = 1:n_D,
    #                 # @simd for
    #                 wi = 1:n_W
    #                     # @simd
    #                     Threads.@threads for hi = 1:n_H
    #                         h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
    #                         h_end =
    #                             hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
    #                         w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
    #                         w_end =
    #                             wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
    #                         d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
    #                         d_end =
    #                             di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
    #                         @inbounds ai = Aip[
    #                             h_start:h_end,
    #                             w_start:w_end,
    #                             d_start:d_end,
    #                             ci,
    #                             mi,
    #                         ]
    #                         if cLayer isa MaxPoolLayer
    #                             mask = (Ao[hi,wi,di,ci,mi] .== ai)
    #                             @inbounds dAip[
    #                                 h_start:h_end,
    #                                 w_start:w_end,
    #                                 d_start:d_end,
    #                                 ci,
    #                                 mi,
    #                             ] .+= (dAo[hi, wi, di, ci, mi] .* mask)
    #                         else
    #                             mask = similar(ai)
    #                             mask .= 1 / prod(size(ai))
    #                             @inbounds dAip[
    #                                 h_start:h_end,
    #                                 w_start:w_end,
    #                                 d_start:d_end,
    #                                 ci,
    #                                 mi,
    #                             ] .+= (dAo[hi, wi, di, ci, mi] .* mask)
    #                         end #if cLayer isa MaxPoolLayer
    #                     end #for hi=1:n_H
    #         #         end #for wi=1:n_W
    #         #     end #for di=1:n_D
    #         # end #for ci=1:c
    #     end #for mi=1:m
    #
    # else
        # @simd
        for mi = 1:m
            # @simd
            Threads.@threads for ci = 1:c
                for di = 1:n_D, wi = 1:n_W, hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                    w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                    d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                    d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                    @inbounds ai =
                        Aip[h_start:h_end, w_start:w_end, d_start:d_end, ci, mi]
                    if cLayer isa MaxPoolLayer
                        mask = (Ao[hi,wi,di,ci,mi] .== ai)
                        @inbounds dAip[
                            h_start:h_end,
                            w_start:w_end,
                            d_start:d_end,
                            ci,
                            mi,
                        ] .+= (dAo[hi, wi, di, ci, mi] .* mask)
                    else
                        mask = similar(ai)
                        mask .= 1 / prod(size(ai))
                        @inbounds dAip[
                            h_start:h_end,
                            w_start:w_end,
                            d_start:d_end,
                            ci,
                            mi,
                        ] .+= (dAo[hi, wi, di, ci, mi] .* mask)
                    end #if cLayer isa MaxPoolLayer
                end #for wi=1:n_W, hi=1:n_H, di=1:n_D
            end #for ci=1:c
        end #for mi=1:m
    # end #if sS == fS

    dAi .= dAip[1+padS[1]:end-padS[2],
        1+padS[3]:end-padS[4],
        1+padS[5]:end-padS[6],
        :,
        :,
    ]

    return nothing
end #function dpooling!(cLayer::ThreeD,

export dpooling!


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "breakpad"
version = v"2020.02.21"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/google/breakpad.git", "815497495eee92e028320ef96dc5b7ec13d85216"),
    GitSource("https://chromium.googlesource.com/linux-syscall-support.git", "f70e2f1641e280e777edfdad7f73a2cfa38139c7"),
    DirectorySource("./bundled"),

    # These are listed in the `.gclient` generated by `depot_tools`, but they dont' appear to actually be needed.
    #GitSource("https://github.com/google/googletest.git", "5ec7f0c4a113e2f18ac2c6cc7df51ad6afc24081"),
    #GitSource("https://github.com/google/protobuf.git", "cb6dd4ef5f82e41e06179dcd57d3b1d9246ad6ac"),
    #GitSource("https://chromium.googlesource.com/external/gyp.git", "324dd166b7c0b39d513026fa52d6280ac6d56770"; unpack_target="gyp"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/breakpad

# Apply old glibc patch for aarch64
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/glibc_aarch64.patch"

# Link necessary subprojects into their appropriate locations
ln -s ${WORKSPACE}/srcdir/linux-syscall-support src/third_party/lss
./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE}
make -j${nproc}
make install
install_license LICENSE
"""

# musl is broken, same with ppc64le
platforms = filter(p -> arch(p) != "powerpc64le" && libc(p) != "musl", supported_platforms())
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("minidump_dump", :minidump_dump),
    # Unfortunately, this is not provided on every platform.  Sigh.
    #ExecutableProduct("minidump-2-core", :minidump_2_core),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DataFrames

"""
    padItem(item, nchar)

Pads an item with an appropriate amount of "." characters so that n-grams can be
created. Provides two methods:

    padItem(item::Array{String}, nchar::Integer)
    padItem(item::String, nchar::Integer)

# Arguments
* **item** The item that needs to be padded; may be an `Array` of `Strings` or
    simply a `String`
* **nchar** The number of characters contained in the n-gram that is being
    examined (e.g., 2 for diphones)

# Returns
The item with appropriate padding
"""
function padItem(item::Array{String}, nchar::Integer)
    padding = ["." for x in 1:(nchar-1)]
    item = vcat(padding, item, padding)
    return item
end

function padItem(item::String, nchar::Integer)
    padding = repeat(".", nchar-1)
    item = padding * item * padding
    return item
end

"""
    createGrams(item, nchar::Integer)

Creates the n-grams that will be analyzed. If looking at diphones, will produce
diphones from the input, for example. Here, n-grams is used in a general sense
as an extension of the n-gram in natural language processing.

# Arguments
* **item** The item for which to create the list of n-grams
* **nchar** The number of characters for each gram (e.g., 2 for diphones)

# Returns
A `Vector` of the created n-grams
"""
function createGrams(item, nchar::Integer)
    grams = Vector()
    for i=1:(length(item)-(nchar-1))
        push!(grams, item[i:i+(nchar-1)])
    end
    return grams
end

"""
    countSymbols(corpus::Array, frequencies::Array, positional::Bool,
        nchar::Integer, pad::Bool)

Counts the number of occurrences of each symbol in the passed in corpus;
produces, for example, the number of times "K" appears in the corpus.

# Arguments
* **corpus** The corpus from which to derive the counts
* **frequencies** The number of times each item in the corpus has been observed;
    this variable must correspond with `corpus` such that the n-th element of
    `corpus` is related to the n-th element of `frequencies`.
* **positional** Whether to consider where in a word a given phone appears
(e.g., should "K" as the first sound be considered a different category than "K"
    as the second sound?)
* **nchar** The number of characters for each n-gram that will be examined
    (e.g., 2 for diphones)
* **pad** Whether to add padding to each query or not

# Returns
A `Dict` where each key is a differnt observed "gram," and each paired value is
the number of times that n-gram was observed
"""
function countSymbols(corpus::Array, frequencies::Array, positional::Bool,
    nchar::Integer, pad::Bool)

    counts = Dict()
    total = 0

    for item in corpus

        if pad
            item = padItem(item, nchar)
        end
        item = createGrams(item, nchar)

        for (i, ch) in enumerate(item)

            key = positional ? (i, ch) : ch

            if haskey(counts, key)
                counts[key] += frequencies[i]
            else
                counts[key] = frequencies[i]
            end

            total += 1
        end
    end

    return counts
end

"""
    calcProb(counts, query, positional, total)

Calculates the phonotactic probability for a query by determining the count of
each n-gram contained in the query, dividing each count by the total number of
n-grams observed in the corpus, and taking the product of those proportions

# Arguments
* **counts** A `Dict` that can be indexed with n-grams and produce the number
    of occurrences for that item
* **query** The item for which the phonotactic probability should be calculated
* **positional** Whether to consider where in the query a given phone appears
    (e.g., should "K" as the first sound be considered a different category
    than "K" as the second sound?)

# Returns
The calculated phonotactic probability for the given query item
"""
function calcProb(counts::Dict, query, positional::Bool, total)
    unitCounts = Int64[]

    for (i, ch) in enumerate(query)

        key = positional ? (i, ch) : ch

        if haskey(counts, key)
            push!(unitCounts, counts[key])
        else
            push!(unitCounts, 0)
        end
    end

    return prod(unitCounts / total)
end

"""
    phnprb(corpus::Array, frequencies::Array, queries::Array; positional=false,
        nchar=1, pad=true)

Calculates the phonotactic probability for each item in a list of queries based
on a corpus

# Arguments
* **corpus** The corpus on which to base the probability calculations
* **frequencies** The frequencies associated with each element in `corpus`
* **queries** The items for which the probability should be calculated

## Keyword arguments
* **positional**  Whether to consider where in the query a given phone appears
(e.g., should "K" as the first sound be considered a different category than "K"
    as the second sound?)
* **nchar** The number of characters for each n-gram that will be examined
    (e.g., 2 for diphones)
* **pad** Whether to add padding to each query or not

# Returns
A `DataFrame` with the queries in the first column and the probability values
    in the second
"""
function phnprb(corpus::Array, frequencies::Array{Int64}, queries::Array;
    positional=false, nchar=1, pad=true)

    counts = countSymbols(corpus, frequencies, positional, nchar, pad)
	total = sum(values(counts))

    results = Vector()

    for query in queries
        local paddedQuery
        if pad
            paddedQuery = padItem(query, nchar)
        else
            paddedQuery = query
        end
        grams = createGrams(paddedQuery, nchar)
        prob = calcProb(counts, grams, positional, total)
        push!(results, prob)
    end

    return DataFrame(Query=queries, Probability=results)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/LibUV_jll.jl

module LibUV_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libuv

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libuv_handle = C_NULL
libuv_path = ""

if Sys.iswindows()
    const libuv = "libuv-2.dll"
elseif Sys.isapple()
    const libuv = "@rpath/libuv.2.dylib"
else
    const libuv = "libuv.so.2"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libuv_handle = dlopen(libuv)
    global libuv_path = dlpath(libuv_handle)
end

is_available() = true

end  # module LibUV_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
module TestOTSIN
    using Test
    using EcologicalNetworks

    M = [1 2 1 1 1;
         0 5 0 0 0;
         3 1 1 0 0;
         1 1 1 2 1]

    a = [0.25, 0.25, 0.25, 0.25]
    b = [0.3, 0.1, 0.1, 0.25, 0.25]

    Q = optimaltransportation(M; a, b)
    @test Q isa AbstractMatrix
    @test sum(Q) ≈ 1.0
    @test sum(Q, dims=2)[:] ≈ a
    @test sum(Q, dims=1)[:] ≈ b
    @test optimaltransportation(M; a, b, λ=0.0) ≈ a * b'

    Q = optimaltransportation(M; a)
    @test sum(Q) ≈ 1.0
    @test sum(Q, dims=2)[:] ≈ a

    Q = optimaltransportation(M; b)
    @test sum(Q) ≈ 1.0
    @test sum(Q, dims=1)[:] ≈ b

    Q = optimaltransportation(M)
    @test sum(Q) ≈ 1.0

    N = BipartiteNetwork(M.>0)
    Q = optimaltransportation(N; a, b)

    @test Q isa BipartiteQuantitativeNetwork

end


text:
using GeoData, Test, ArchGDAL, GeoFormatTypes
using GeoData: resample

include(joinpath(dirname(pathof(GeoData)), "../test/test_utils.jl"))

@testset "resample" begin
    raster_path = maybedownload("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")

    output_res = 0.0027
    output_crs = EPSG(4326)
    resample_method = "near"

    ## Resample cea.tif manually with ArchGDAL
    wkt = convert(String, convert(WellKnownText, output_crs))
    AG_output = ArchGDAL.read(raster_path) do dataset
        ArchGDAL.gdalwarp([dataset], ["-t_srs", "$(wkt)",
                                "-tr", "$(output_res)", "$(output_res)",
                                "-r", "$(resample_method)"]) do warped
            ArchGDAL.read(ArchGDAL.getband(warped, 1))
        end
    end

    ## Resample cea.tif using resample
    cea = GeoArray(GDALarray(raster_path))
    crs(cea)
    GD_output = resample(cea, output_res, crs=output_crs, method=resample_method)

    ## Compare the two
    @test AG_output == GD_output.data[:, :, 1]
    @test abs(step(dims(GD_output, Y))) ≈ abs(step(dims(GD_output, X))) ≈ output_res

    @testset "snapped size and dim index match" begin
        snaptarget = GD_output
        snapped = resample(cea, snaptarget)
        @test size(snapped) == size(snaptarget)
        @test isapprox(index(snaptarget, Y), index(snapped, Y))
        @test isapprox(index(snaptarget, X), index(snapped, X))
    end
end