text:
# API definition file for data layer abstraction.
# The idea is for a one line call to change from internal Graphs.jl to CloudGraphs.jl
# switch for as data layer.

type DataLayerAPI
  addvertex!::Function
  getvertex::Function
  # setupvertgraph!::Function
  # setupfncvertgraph!::Function
  makeaddedge!::Function
  getedge::Function
  outneighbors::Function
  updatevertex!::Function
  updateedge!::Function
  deletevertex!::Function
  deleteedge!::Function
  cgEnabled::Bool
end

# global dlapi

dlapi = DataLayerAPI(addGraphsVert!,          # addvertex
                     getVertNode,             # getvertex
                     makeAddEdge!,            # makeaddedge
                     graphsGetEdge,           # getedge
                     graphsOutNeighbors,      # outneighbors
                     updateFullVertData!,     # updatevertex!
                     +,                       # updateedge!
                     graphsDeleteVertex!,     # deletevertex!
                     +,                       # deleteedge!
                     false )

localapi = DataLayerAPI(addGraphsVert!,          # addvertex
                        getVertNode,             # getvertex
                        makeAddEdge!,            # makeaddedge
                        graphsGetEdge,           # getedge
                        graphsOutNeighbors,      # outneighbors
                        updateFullVertData!,     # updatevertex!
                        +,                       # updateedge!
                        graphsDeleteVertex!,     # deletevertex!
                        +,                       # deleteedge!
                        false )
#

# setCloudDataLayerAPI!
function setdatalayerAPI!(;addvertex::Function = addGraphsVert!,
      getvertex::Function = getVertNode,
      makeaddedge::Function = makeAddEdge!,
      getedge::Function = graphsGetEdge,
      outneighbors::Function = graphsOutNeighbors,
      updatevertex::Function = updateFullVertData!,
      updateedge::Function = +,
      deletevertex::Function = graphsDeleteVertex!,
      deleteedge::Function = +,
      cgEnabled::Bool = false  )

  dlapi.addvertex! = addvertex
  dlapi.getvertex = getvertex
  dlapi.makeaddedge! = makeaddedge
  dlapi.getedge = getedge
  dlapi.outneighbors = outneighbors
  dlapi.updatevertex! = updatevertex
  dlapi.updateedge! = updateedge
  dlapi.deletevertex! = deletevertex
  dlapi.deleteedge! = deleteedge
  dlapi.cgEnabled = cgEnabled

  # dlapi.addvertex! = addCloudVert!
  # dlapi.getvertex = getExVertFromCloud
  # dlapi.makeaddedge! = makeAddCloudEdge!
  # dlapi.getedge = getEdgeFromCloud
  # dlapi.updatevertex! = updateFullCloudVertData!
  # dlapi.outneighbors = getCloudOutNeighbors
  # dlapi.deletevertex! = deleteCloudVertex!
  # dlapi.deleteedge! = deleteCloudEdge!
  # dlapi.cgEnabled = true

  println("Changed internal API calls to use outside calls.")
  nothing
end

function showcurrentdlapi()
  @show dlapi
end

# Remember 3rd party users interact with
# addNode!
# addFactor!
# prepBatchTree!
# inferOverTree!
#
## Visualization functions -- Visualization may well be separated or abstracted out
# writeGraphPdf
# drawHorBeliefsList


text:
# # Logistic regression
using DataFrames
using Plots
using RDatasets
using Convex
using SCS

# This is an example logistic regression using `RDatasets`'s iris data.
# Our goal is to gredict whether the iris species is versicolor
# using the sepal length and width and petal length and width.
iris = dataset("datasets", "iris");
iris[1:10, :]

# We'll define `Y` as the outcome variable: +1 for versicolor, -1 otherwise.
Y = [species == "versicolor" ? 1.0 : -1.0 for species in iris.Species]

# We'll create our data matrix with one column for each feature
# (first column corresponds to offset).
X = hcat(
    ones(size(iris, 1)),
    iris.SepalLength,
    iris.SepalWidth,
    iris.PetalLength,
    iris.PetalWidth,
);

# Now to solve the logistic regression problem.
n, p = size(X)
beta = Variable(p)
problem = minimize(logisticloss(-Y .* (X * beta)))
solve!(problem, SCS.Optimizer; silent_solver = true)

# Let's see how well the model fits.
using Plots
logistic(x::Real) = inv(exp(-x) + one(x))
perm = sortperm(vec(X * evaluate(beta)))
plot(1:n, (Y[perm] .+ 1) / 2, st = :scatter)
plot!(1:n, logistic.(X * evaluate(beta))[perm])

ext: jl


text:
# To-do: make user-defined objective function and constraints work (no metter what)
#  -> try declaring x1, x2, ... separatelly. Ugly, but may work.
using DwsimOpt
using JuMP
using Ipopt
using FiniteDifferences
using LinearAlgebra
using PyCall

path2sim = "c:/Users/lfsfr/Desktop/DwsimOpt.jl/examples/SMR/SMR.dwxmz"
path2dwsim = "C:/Users/lfsfr/AppData/Local/DWSIM7/"

include("c:/Users/lfsfr/Desktop/DwsimOpt.jl/src/SimOpt.jl")
sim_jl = py"sim_smr"

# just for testing: connect to another simulation

py"""
sim2 = SimulationOptimization(dof=np.array([]), path2sim= "c:/Users/lfsfr/Desktop/DwsimOpt.jl/examples/SMR/SMR.dwxmz",
                    path2dwsim = $path2dwsim)
sim2.connect(interf)
"""

py"""
sim_smr = sim2

# Import dwsim-python data exchange interface (has to be after Automation2)
from dwsimopt.py2dwsim import create_pddx, assign_pddx

# Assign DoF:
create_pddx( ["MR-1", "CompoundMassFlow", "Nitrogen", "kg/s"],    sim_smr, element="dof" )
create_pddx( ["MR-1", "CompoundMassFlow", "Methane", "kg/s"],     sim_smr, element="dof" )
create_pddx( ["MR-1", "CompoundMassFlow", "Ethane", "kg/s"],      sim_smr, element="dof" )
create_pddx( ["MR-1", "CompoundMassFlow", "Propane", "kg/s"],     sim_smr, element="dof" )
create_pddx( ["MR-1", "CompoundMassFlow", "Isopentane", "Pa"],    sim_smr, element="dof" )
create_pddx( ["VALV-01", "OutletPressure", "Mixture", "Pa"],      sim_smr, element="dof" )
create_pddx( ["COMP-1", "OutletPressure", "Mixture", "Pa"],       sim_smr, element="dof" )
create_pddx( ["COOL-08", "OutletTemperature", "Mixture", "K"],    sim_smr, element="dof" )

# Assign F
create_pddx( ["Sum_W", "EnergyFlow", "Mixture", "kW"], sim_smr, element="fobj" )

# adding constraints (g_i <= 0):
g1 = create_pddx( ["MITA1-Calc", "OutputVariable", "mita", "°C"], sim_smr, element="constraint", assign=False )
assign_pddx( lambda: 3-g1[0]() , ["MITA1-Calc", "OutputVariable", "mita", "°C"], sim_smr, element="constraint" )
g2 = create_pddx( ["MITA2-Calc", "OutputVariable", "mita", "°C"], sim_smr, element="constraint", assign=False )
assign_pddx( lambda: 3-g2[0]() , ["MITA2-Calc", "OutputVariable", "mita", "°C"], sim_smr, element="constraint" )

# decision variables bounds
x0 = np.array( [0.25/3600, 0.70/3600, 1.0/3600, 1.10/3600, 1.80/3600, 2.50e5, 50.00e5, -60+273.15] )
bounds_raw = np.array( [0.5*np.asarray(x0), 1.5*np.asarray(x0)] )   # 50 % around base case
bounds_raw[0][-1] = 153     # precool temperature low limit manually
bounds_raw[1][-1] = 253     # precool temperature upper limit manually

# regularizer calculation
regularizer = np.zeros(x0.size)
for i in range(len(regularizer)):
    regularizer[i] = 10**(-1*math.floor(math.log(x0[i],10))) # regularizer for magnitude order of 1e0

# bounds regularized
bounds_reg = regularizer*bounds_raw

# objective and constraints lambda definitions
f = lambda x: sim_smr.calculate_optProblem(np.asarray(x)/regularizer)[0:sim_smr.n_f]
g = lambda x: sim_smr.calculate_optProblem(np.asarray(x)/regularizer)[sim_smr.n_f:(sim_smr.n_f+sim_smr.n_g)]
"""

f = py"f"
g = py"g"
x0 = py"x0*regularizer"
searchSpace = py"bounds_reg"
dim = sim_jl.n_dof
op1 = optProblem(f, g, x0, searchSpace, dim, sim_jl, py"sim.n_dof", py"sim.n_f", py"sim.n_g")
op1.sim_jl.verbose = false;
py"sim_smr.verbose=False";
save_sim() = py"""sim.interface.SaveFlowsheet(sim.flowsheet,$pwd()+"/examples/PRICO/PRICO_composite2.dwxmz",True)"""

# model = Model(optimizer_with_attributes(Ipopt.Optimizer))
model = Model(Ipopt.Optimizer)
set_optimizer_attribute(model, "max_cpu_time", 60.0)
set_optimizer_attribute(model, "print_level", 0)

# VARIABLES
# @variable(model, op1.searchSpace[1, i]' <= x[i = 1:op1.dim] <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 1;
@variable(model, op1.searchSpace[1, i]' <= x1 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 2;
@variable(model, op1.searchSpace[1, i]' <= x2 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 3;
@variable(model, op1.searchSpace[1, i]' <= x3 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 4;
@variable(model, op1.searchSpace[1, i]' <= x4 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 5;
@variable(model, op1.searchSpace[1, i]' <= x5 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 6;
@variable(model, op1.searchSpace[1, i]' <= x6 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 7;
@variable(model, op1.searchSpace[1, i]' <= x7 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
i = 8;
@variable(model, op1.searchSpace[1, i]' <= x8 <= op1.searchSpace[2, i]', start = (op1.searchSpace[1, i] + op1.searchSpace[2, i]) ./ 2)
x0 = 0.95 * (op1.searchSpace[1, :] + op1.searchSpace[2, :]) ./ 2
function f_block(x)
    # println("hi")
    # global x0, ∇f_bkp, ∇g1_bkp, ∇g2_bkp
    # if norm(x0 - x) > 1e-5
    #     println("bye")
    #     ∇f_bkp, ∇g1_bkp, ∇g2_bkp = op1.f(x)[1], op1.g(x)[1], op1.g(x)[2]
    #     x0 = x
    # end
    return [op1.f(x)[1], op1.g(x)[1], op1.g(x)[2]]
end
function jac_block(x)
    println("hi")
    global x0, ∇f_bkp, ∇g1_bkp, ∇g2_bkp
    if norm(x0 - x) > 1e-5
        println("bye")
        tmp = jacobian(forward_fdm(2, 1), f_block, x0)[1]
        ∇f_bkp = tmp[1, :]
        ∇g1_bkp = tmp[2, :]
        ∇g2_bkp = tmp[3, :]
        x0 = x
    end
end

# OBJECTIVE FUNCTION
# fobj_bb(x) = op1.f([x])
fobj_bb(x1, x2, x3, x4, x5, x6, x7, x8) = op1.f([x1, x2, x3, x4, x5, x6, x7, x8])[1]
function ∇fobj_bb(g::AbstractVector{T}, x1::T, x2::T, x3::T, x4::T, x5::T, x6::T, x7::T, x8::T) where {T}
    println("I am in ∇fobj_bb with x=", x1)
    global x0, ∇f_bkp, ∇g1_bkp, ∇g2_bkp
    jac_block([x1, x2, x3, x4, x5, x6, x7, x8])
    ∇f_tmp = ∇f_bkp
    # ∇f_tmp = grad(forward_fdm(2, 1), op1.f, [x1, x2, x3, x4, x5, x6, x7, x8])[1]
    for i = 1:length([x1, x2, x3, x4, x5, x6, x7, x8])
        g[i] = ∇f_tmp[i]
    end
end
# ∇fobj_bb(x) = grad(central_fdm(5, 1), op1.f, x)
register(model, :fobj_bb, op1.dim, fobj_bb, ∇fobj_bb; autodiff=false)
# @NLobjective(model, Min, fobj_bb(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8]))
@NLobjective(model, Min, fobj_bb(x1, x2, x3, x4, x5, x6, x7, x8))

# CONSTRAINTS
g1(x) = op1.g(x)[1]
g2(x) = op1.g(x)[2]
g1_bb(x1, x2, x3, x4, x5, x6, x7, x8) = op1.g([x1, x2, x3, x4, x5, x6, x7, x8])[1]
g2_bb(x1, x2, x3, x4, x5, x6, x7, x8) = op1.g([x1, x2, x3, x4, x5, x6, x7, x8])[2]
function ∇g1_bb(g::AbstractVector{T}, x1::T, x2::T, x3::T, x4::T, x5::T, x6::T, x7::T, x8::T) where {T}
    println("I am in ∇g1_bb with x=", x1)
    global x0, ∇f_bkp, ∇g1_bkp, ∇g2_bkp
    jac_block([x1, x2, x3, x4, x5, x6, x7, x8])
    ∇g1_tmp = ∇g1_bkp
    # ∇g1_tmp = grad(forward_fdm(2, 1), g1, [x1, x2, x3, x4, x5, x6, x7, x8])[1]
    for i = 1:length([x1, x2, x3, x4, x5, x6, x7, x8])
        g[i] = ∇g1_tmp[i]
    end
end
function ∇g2_bb(g::AbstractVector{T}, x1::T, x2::T, x3::T, x4::T, x5::T, x6::T, x7::T, x8::T) where {T}
    println("I am in ∇g2_bb with x=", x1)
    global x0, ∇f_bkp, ∇g1_bkp, ∇g2_bkp
    jac_block([x1, x2, x3, x4, x5, x6, x7, x8])
    ∇g2_tmp = ∇g2_bkp
    # ∇g2_tmp = grad(forward_fdm(2, 1), g2, [x1, x2, x3, x4, x5, x6, x7, x8])[1]
    for i = 1:length([x1, x2, x3, x4, x5, x6, x7, x8])
        g[i] = ∇g2_tmp[i]
    end
end
register(model, :g1_bb, op1.dim, g1_bb, ∇g1_bb; autodiff=false)

@NLconstraint(model, g1_bb(x1, x2, x3, x4, x5, x6, x7, x8) <= 0)
register(model, :g2_bb, op1.dim, g2_bb, ∇g2_bb; autodiff=false)
@NLconstraint(model, g2_bb(x1, x2, x3, x4, x5, x6, x7, x8) <= 0)

JuMP.optimize!(model)

@show JuMP.termination_status(model)
@show JuMP.primal_feasibility_report(model)
@show JuMP.dual_status(model)

# retrieve the objective value, corresponding x values and the status
println(JuMP.value.([x1, x2, x3, x4, x5, x6, x7, x8]))
println(JuMP.objective_value(model))
println(JuMP.termination_status(model))


text:
using ReverseDiff: InstructionTape, GradientConfig, JacobianConfig, HessianConfig,
                   value, deriv, tape, valtype,
                   derivtype, track, track!

using Random

const COMPILED_TAPE_LIMIT = 5000

# These functions correctly emit NaNs for certain arguments, but ReverseDiff's test
# machinery is currently too dumb to handle them properly.
const SKIPPED_BINARY_SCALAR_TESTS = Symbol[:hankelh1, :hankelh1x, :hankelh2, :hankelh2x,
                                           :pow, :besselj, :besseli, :bessely, :besselk,
                                           :polygamma]

# make RNG deterministic, and thus make result inaccuracies
# deterministic so we don't have to retune EPS for arbitrary inputs
Random.seed!(1)

test_println(kind, f, pad = "  ") = println(pad, "testing $(kind): `$(f)`...")

@inline test_approx(A, B, _atol = 1e-5) = @test isapprox(A, B, atol = _atol)

tracked_is(a, b) = value(a) === value(b) && deriv(a) === deriv(b) && tape(a) === tape(b)
tracked_is(a::AbstractArray, b::AbstractArray) = all(map(tracked_is, a, b))
tracked_is(a::Tuple, b::Tuple) = all(map(tracked_is, a, b))


text:
module Timestamps

using Dates
import Dates: month, year, day, value
import Base: show, size, getindex, push!, isempty
import Tables
# import Base: getindex, show, convert, sum, prod, mean, var, std, maximum, minimum, merge

include("timestamp.jl")
include("timestamparrays.jl")
include("tables.jl")
include("conversion.jl")
include("operators.jl")
include("utils.jl")

export Timestamp, TimestampArray

end


text:
const PROJ_ROOT = dirname(dirname(@__DIR__))
const DEPS_DIR = joinpath(PROJ_ROOT, "deps")
const DATA_DIR = joinpath(PROJ_ROOT, "data")
    const RAW_DATA_DIR = joinpath(DATA_DIR, "raw")
        const HUMAN1_PUBLICATION_RAW_DATA_DIR = joinpath(RAW_DATA_DIR, HUMAN1_PUBLICATION_PROJ_NAME)
        const HUMAN1_PUBLICATION_ZIP_FILE = joinpath(RAW_DATA_DIR, "$HUMAN1_PUBLICATION_PROJ_NAME.zip")
    const PROCESSED_DATA_DIR = joinpath(DATA_DIR, "processed")
        const HUMAN1_MODELS_DIR = joinpath(PROCESSED_DATA_DIR, "human1_models")
        const FVA_PP_MODELS_DIR = joinpath(PROCESSED_DATA_DIR, "fva_pp_models")
        const HUMAN1_RESULTS_DIR = joinpath(PROCESSED_DATA_DIR, "human1_results")
        const HUMAN1_NCI60_RESULTS_DIR = joinpath(HUMAN1_RESULTS_DIR, "11_cellLines_NCI60")
        const CACHE_DIR = joinpath(PROCESSED_DATA_DIR, "cache")
    const FIGURES_DATA_DIR = joinpath(DATA_DIR, "figures")

function _make_dirs()
    for dir in [DATA_DIR, RAW_DATA_DIR, DEPS_DIR, HUMAN1_MODELS_DIR,
                HUMAN1_RESULTS_DIR, PROCESSED_DATA_DIR,
                FIGURES_DATA_DIR, FVA_PP_MODELS_DIR, CACHE_DIR]
        if !isdir(dir)
            mkpath(dir)
        end
    end
end


text:
module Tst

include("src/transport/transport.jl")

end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
const DDMVector{T} = DDMArray{T,1}

"""
    CartesianDDMVector{T,C,A}

Objects of `CartesianDDMVector` behave like a normal vector, except
that the storage is not stored continuously but by DDM subdomain.

For now, a Boolean partition of unity is assumed but this could be
generalized at a later time.

# Note on `DenseArray`s

1. Would `CartesianDDMVector <: DenseVector` if contiguous storage?
1. Optional field to signal/ensure coherence?

"""
struct CartesianDDMVector{T,C,A} <: DDMVector{T}
    context::C
    parent::A
end

parent(x::CartesianDDMVector) = getproperty(x, :parent)
getcontext(x::CartesianDDMVector) = getproperty(x, :context)

function size(x::CartesianDDMVector)
    (; context, parent) = x
    (; dims) = context
    (mapreduce(getdof, *, dims),)
end

"""
    getindex(x::CartesianDDMVector, i)

Assumes partition of unity is Boolean for now. This is a very
inefficient way to read data. It is included only for debugging
purposes.

"""
function getindex(x::CartesianDDMVector, i::Int)
    @boundscheck checkbounds(x, i)

    (; context, parent) = x

    (; dims) = context
    ci = getindex(CartesianIndices(getdof.(dims)), i)

    nocontext = removeoverlap(context)
    noglb, _ = ranges(nocontext)

    index = findfirst(noglb) do el
        in(ci, CartesianIndices(el))
    end

    glb, _, lcl = ranges(context)
    domain = CartesianIndices(glb[index])

    cj = findfirst(domain) do cj
        cj == ci
    end

    j = LinearIndices(CartesianIndices(lcl[index]))[cj]

    parent[index][j]
end

"""
    setindex!(x::CartesianDDMVector, val, i::Int)

Writes data in every subdomain that contains element i. This is
a very inefficient way to write data. It is included only for
debugging purposes.

"""
function setindex!(x::CartesianDDMVector, val, i::Int)
    @boundscheck checkbounds(x, i)

    (; context, parent) = x

    (; dims) = context
    ci = getindex(CartesianIndices(getdof.(dims)), i)

    glb, _, lcl = ranges(context)

    indices = findall(glb) do el
        in(ci, CartesianIndices(el))
    end

    for index in indices
        domain = CartesianIndices(glb[index])

        cj = findfirst(domain) do cj
            cj == ci
        end

        domain = CartesianIndices(lcl[index])

        j = LinearIndices(domain)[cj]

        parent[index][j] = val
    end

    val
end

#=
"""
Do not use until `makecoherent` implemented.

"""
function CartesianDDMVector{T}(::UndefInitializer, context) where {T}
    _, lcl = ranges(context)

    parent = map(lcl) do el
        Array{T}(undef, prod(length.(el)))
    end

    CartesianDDMVector{T,typeof(context),typeof(parent)}(context, parent)
end
=#

"""
Do not use until `makecoherent` implemented.

"""
function CartesianDDMVector(init::Function, context)
    _, lcl = ranges(context)

    parent = map(lcl) do el
        init(prod(length.(el)))
    end

    T = eltype(eltype(parent))

    x = CartesianDDMVector{T,typeof(context),typeof(parent)}(context, parent)

    makecoherent!(x)
end

function decompose(context::CartesianDDMContext, x::AbstractVector)
    (; dims) = context

    y = reshape(x, getdof.(dims))

    glb, _ = ranges(context)

    parent = map(glb) do el
        reshape(y[CartesianIndices(el)], :)
    end

    CartesianDDMVector{eltype(x),typeof(context),typeof(parent)}(context, parent)
end

function similar(x::CartesianDDMVector)
    (; context, parent) = x

    out = map(parent) do el
        similar(el)
    end

    CartesianDDMVector{eltype(x),typeof(context),typeof(parent)}(context, out)
end


text:
# ## Atoms and Molecules

export AtoBohr!, Atom, Molecule, AtoZ

mutable struct Atom
    atno::Int64
    x::Float64
    y::Float64
    z::Float64
end

mutable struct Molecule
    atomlist::Array{Atom,1}
end

function push!(mol::Molecule,at::Atom)
    Base.push!(atomlist,at)
end

nuclear_repulsion(a::Atom,b::Atom)= a.atno*b.atno/sqrt(dist2(a.x-b.x,a.y-b.y,a.z-b.z))
function nuclear_repulsion(mol::Molecule)
    nr = 0
    for (i,j) in spairs(nat(mol))
        nr += nuclear_repulsion(mol.atomlist[i],mol.atomlist[j])
    end
    return nr
end

nel(mol::Molecule) = sum([at.atno for at in mol.atomlist])
nat(mol::Molecule) = length(mol.atomlist)

# Other molecule methods to implement
# nocc, nclosed, nopen, nup, ndown, stoich, mass,
# center_of_mass, center!

# Array of symbols, masses
# First we need to dereference atom name to atomic charge (Z)
AtoZ=Dict(j=>i for (i,j) in enumerate(["H","He",
    "Li","Be","B","C","N","O","F","Ne",
    "Na","Mg","Al","Si","P","S","Cl","Ar",
    "K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
    "Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe",
    "Cs","Ba","La",
    "Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu",
    "Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn",
    "Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr",
    "Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"]))
# Magic comprehension ends up like:
# ZfromAtom = Dict("H" => 1, "He" => 2, "Li" => 3)

# Sample molecules for tests - all from PyQuante
# All units currently Angstrom; converted into Bohr below
h2 = Molecule([Atom(1,  0.00000000,     0.00000000,     0.36628549),
               Atom(1,  0.00000000,     0.00000000,    -0.36628549)])

h2o = Molecule([Atom(8,   0.00000000,     0.00000000,     0.04851804),
                Atom(1,   0.75300223,     0.00000000,    -0.51923377),
                Atom(1,  -0.75300223,     0.00000000,    -0.51923377)])

ch4 = Molecule([Atom(6,   0.00000000,     0.00000000,     0.00000000),
                Atom(1,   0.62558332,    -0.62558332,     0.62558332),
                Atom(1,  -0.62558332,     0.62558332,     0.62558332),
                Atom(1,   0.62558332,     0.62558332,    -0.62558332),
                Atom(1,  -0.62558332,    -0.62558332,    -0.62558332)])

c6h6 = Molecule([ Atom(6,  0.98735329,     0.98735329,     0.00000000),
                  Atom(6,  1.34874967,    -0.36139639,     0.00000000),
                  Atom(6,  0.36139639,    -1.34874967,     0.00000000),
                  Atom(6, -0.98735329,    -0.98735329,     0.00000000),
                  Atom(6, -1.34874967,     0.36139639,     0.00000000),
                  Atom(6, -0.36139639,     1.34874967,     0.00000000),
                  Atom(1,  1.75551741,     1.75551741,     0.00000000),
                  Atom(1,  2.39808138,    -0.64256397,     0.00000000),
                  Atom(1,  0.64256397,    -2.39808138,     0.00000000),
                  Atom(1, -1.75551741,    -1.75551741,     0.00000000),
                  Atom(1, -2.39808138,     0.64256397,     0.00000000),
                  Atom(1, -0.64256397,     2.39808138,     0.00000000)])

lih = Molecule([Atom(3,    0.00000000,     0.00000000,    -0.53999756),
                Atom(1,    0.00000000,     0.00000000,     1.08999756)])

# Convert to atomic units (bohr)
AtoBohr(x::Float64) = x*BohrInA # These assume units start as Angstrom.
function AtoBohr!(at::Atom)
    at.x *= BohrInA
    at.y *= BohrInA
    at.z *= BohrInA
end
function AtoBohr!(mol::Molecule)
    for at in mol.atomlist
        AtoBohr!(at)
    end
end

AtoBohr!(h2)
AtoBohr!(h2o)
AtoBohr!(ch4)
AtoBohr!(c6h6)
AtoBohr!(lih)


text:
using Combinatorics, StaticArrays, LinearAlgebra

# miscallaneous
Base.ntuple(n::Integer, ::Val{M}) where {M} = ntuple(_->n, M)

# --------------- helpful functions to connect ---------------
#                 permutations of corners and edges

# 4-body = 4-simplex has 4 corners (atom positions) and 6 edges  rᵢⱼ
#
# edge index : A[1]  A[2]  A[3]  A[4]  A[5]  A[6]
# edge length: r12   r13   r14   r23   r24   r34
# (where rij = |xᵢ - xⱼ| with xᵢ the corner positions)

const b4_e_inds = [0 1 2 3
                   1 0 4 5
                   2 4 0 6
                   3 5 6 0]

# 5-body = 5-simplex has 5 corners (atom positions) and 10 edges  rᵢⱼ
#
# edge index : A[1]  A[2]  A[3]  A[4]  A[5]  A[6] A[7]  A[8]  A[9]  A[10]
# edge length: r12   r13   r14   r15   r23   r24  r25   r34   r35   r45
# (where rij = |xᵢ - xⱼ| with xᵢ the corner positions)

const b5_e_inds = [0 1 2 3 4
                   1 0 5 6 7
                   2 5 0 8 9
                   3 6 8 0 10
                   4 7 9 10 0]

# preparation for true n-body terms only
 #for 3-body
 const b3_e_proj = [1 1 0
                    1 0 1
                    0 1 1]
#for 4-body
 const b4_e_proj = [1 1 1 0 0 0
                    1 0 0 1 1 0
                    0 1 0 1 0 1
                    0 0 1 0 1 1]
#for 5-body
const b5_e_proj = [1 1 1 1 0 0 0 0 0 0
                   1 0 0 0 1 1 1 0 0 0
                   0 1 0 0 1 0 0 1 1 0
                   0 0 1 0 0 1 0 1 0 1
                   0 0 0 1 0 0 1 0 1 1]

const πb3 = collect(permutations(1:3))

"""
convert a permutation of simplex corners into a permutation of
simplex edges (for 4 body)
"""
S4_to_S6(π::Vector{Int}) = Int[
   b4_e_inds[π[1], π[2]], b4_e_inds[π[1], π[3]], b4_e_inds[π[1], π[4]],
   b4_e_inds[π[2], π[3]], b4_e_inds[π[2], π[4]], b4_e_inds[π[3], π[4]] ]

"""
convert a permutation of simplex corners into a permutation of
simplex edges (for 5 body)
"""
S5_to_S10(π::Vector{Int}) = Int[
   b5_e_inds[π[1], π[2]], b5_e_inds[π[1], π[3]], b5_e_inds[π[1], π[4]],
   b5_e_inds[π[1], π[5]], b5_e_inds[π[2], π[3]], b5_e_inds[π[2], π[4]],
   b5_e_inds[π[2], π[5]], b5_e_inds[π[3], π[4]], b5_e_inds[π[3], π[5]],
   b5_e_inds[π[4], π[5]] ]

"""
generate all permutations of A that correspond to permutations of corners,
then keep only the unique ones so that we don't double-count.
"""
simplex_permutations(::Val{5}, A) =
   (  [ A[S5_to_S10(πX)]
              for πX in permutations(1:5) ]  )

simplex_permutations(::Val{4}, A) =
   (  [ A[S4_to_S6(πX)]
              for πX in permutations(1:4) ]  )

simplex_permutations(::Val{3}, A) = ( [ A[π] for π in πb3 ] )

simplex_permutations(::Val{2}, A) = [ A ]

simplex_permutations(x::SVector{3}) =
   simplex_permutations(Val(3), x)

simplex_permutations(x::SVector{6}) =
   simplex_permutations(Val(4), x)

simplex_permutations(x::SVector{10}) =
   simplex_permutations(Val(5), x)

# # ==== TESTS =======
# using NBodyIPs, JuLIP, StaticArrays
#
# R = rand(JVecF, 10)
# J = @SVector [3,5,7]
# s = Array(norm.(R[J]))
# for i = 1:2, j = i+1:3
#    push!(s, norm(R[J[i]]-R[J[j]]))
# end
# s2, _ = NBodyIPs.simplex_edges(R, J)
# s2 == s

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# abstracttensor.jl
#
# Abstract Tensor type
#----------------------
"""
    abstract type AbstractTensorMap{T, S<:IndexSpace, N₁, N₂} end

Abstract supertype of all tensor maps, i.e. linear maps between tensor products
of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from
an input space of type `ProductSpace{S,N₂}` to an output space of type
`ProductSpace{S,N₁}` and has entries (with respect to the canonical basis) of
type `T`.
"""
abstract type AbstractTensorMap{T<:Number, S<:IndexSpace, N₁, N₂} end
"""
    AbstractTensor{T<:Number, S<:IndexSpace, N} = AbstractTensorMap{T,S,N,0}

Abstract supertype of all tensors, i.e. elements in the tensor product space
of type `ProductSpace{S,N}`, built from elementary spaces of type `S<:IndexSpace`.
The entries of the tensor are of type `T<:Number` (with respect to the canonical basis).

An `AbstractTensor{T,S,N}` is actually a special case `AbstractTensorMap{T,S,N,0}`,
i.e. a tensor map with only a non-trivial output space.
"""
const AbstractTensor{T<:Number, S<:IndexSpace, N} = AbstractTensorMap{T, S, N, 0}

# tensor characteristics
Base.eltype(t::AbstractTensorMap) = eltype(typeof(t))
spacetype(t::AbstractTensorMap) = spacetype(typeof(t))
fieldtype(t::AbstractTensorMap) = fieldtype(typeof(t))
numind(t::AbstractTensorMap) = numind(typeof(t))

Base.eltype(::Type{<:AbstractTensorMap{T}}) where {T<:Number} = T
spacetype(::Type{<:AbstractTensorMap{T,S}}) where {T<:Number, S<:IndexSpace} = S
fieldtype(::Type{<:AbstractTensorMap{T,S}}) where {T<:Number, S<:IndexSpace} = fieldtype(S)
numind(::Type{<:AbstractTensorMap{S,T,N₁,N₂}}) where {T<:Number, S<:IndexSpace, N₁, N₂} = N₁ + N₂

const order = numind

space(t::AbstractTensorMap) = codomain(t) ⊗ dual(domain(t))

# Defining vector spaces:
#------------------------
const TensorSpace{S,N} = ProductSpace{S,N}
const TensorMapSpace{S,N₁,N₂} = Pair{ProductSpace{S,N₂},ProductSpace{S,N₁}}

# Little unicode hack to define TensorMapSpace
→(dom::ProductSpace{S}, codom::ProductSpace{S}) where {S<:IndexSpace} = dom => codom
→(dom::S, codom::ProductSpace{S}) where {S<:IndexSpace} = ProductSpace(dom) => codom
→(dom::ProductSpace{S}, codom::S) where {S<:IndexSpace} = dom => ProductSpace(codom)
→(dom::S, codom::S) where {S<:IndexSpace} = ProductSpace(dom) => ProductSpace(codom)

←(codom::ProductSpace{S}, dom::ProductSpace{S}) where {S<:IndexSpace} = dom => codom
←(codom::S, dom::ProductSpace{S}) where {S<:IndexSpace} = ProductSpace(dom) => codom
←(codom::ProductSpace{S}, dom::S) where {S<:IndexSpace} = dom => ProductSpace(codom)
←(codom::S, dom::S) where {S<:IndexSpace} = ProductSpace(dom) => ProductSpace(codom)

# dims
dims(P::TensorSpace) = map(dim, P.spaces)
dims(P::TensorMapSpace) = (dims(P[2])..., dims(P[1])...)

# Basic algebra
#---------------
Base.copy(t::AbstractTensorMap) = Base.copy!(similar(t), t)

Base.:-(t::AbstractTensorMap) = scale!(t, -one(eltype(t)))
function Base.:+(t1::AbstractTensorMap{T1}, t2::AbstractTensorMap{T2}) where {T1,T2}
    T = promote_type(T1,T2)
    return add!(copy!(similar(t1,T), t1), t2, one(T))
end
function Base.:-(t1::AbstractTensorMap{T1}, t2::AbstractTensorMap{T2}) where {T1,T2}
    T = promote_type(T1,T2)
    return add!(copy!(similar(t1,T), t1), t2, -one(T))
end

Base.:*(t::AbstractTensorMap, a::Number) = scale!(similar(t, promote_type(eltype(t),typeof(a))), t, a)
Base.:*(a::Number, t::AbstractTensorMap) = *(t,a)
Base.:/(t::AbstractTensorMap, a::Number) = *(t, one(a)/a)
Base.:\(a::Number, t::AbstractTensorMap) = *(t, one(a)/a)

Base.scale!(t::AbstractTensorMap, a::Number) = scale!(t, t, a)
Base.scale!(a::Number, t::AbstractTensorMap) = scale!(t, t, a)
Base.scale!(tdest::AbstractTensorMap, a::Number, tsrc::AbstractTensorMap) = scale!(tdest, tsrc, a)

Base.LinAlg.axpy!(a::Number, tx::AbstractTensorMap, ty::AbstractTensorMap) = add!(ty, tx, a)

# Base.conj(t::AbstractTensor) = Base.conj!(similar(t, conj(space(t))), t)
# Base.transpose(t::AbstractTensor) = Base.transpose!(similar(t, space(t).'), t)
# Base.ctranspose(t::AbstractTensor) = Base.ctranspose!(similar(t, space(t)'), t)

# Tensor operations
#-------------------
# convenience definition which works for vectors and matrices but also sometimes useful in general case
# *{S,T1,T2,N1,N2}(t1::AbstractTensor{S,T1,N1},t2::AbstractTensor{S,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[1:N1-1] ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(1:N1-1,0),'N',t2,vcat(0,numind(t1)-1+(1:N2-1)),'N',0,t3,1:(N1+N2-2)))
# Base.At_mul_B{S,T1,T2,N1,N2}(t1::AbstractTensor{S,T1,N1},t2::AbstractTensor{S,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[2:N1].' ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(0,reverse(1:N1-1)),'N',t2,vcat(0,N1-1+(1:N2-1)),'N',0,t3,1:(numind(t1)+numind(t2)-2)))
# Base.Ac_mul_B{S,T1,T2,N1,N2}(t1::AbstractTensor{S,T1,N1},t2::AbstractTensor{S,T2,N2})=(t3=similar(t1,promote_type(T1,T2),space(t1)[2:N1]' ⊗ space(t2)[2:N2]);tensorcontract!(1,t1,vcat(0,reverse(1:N1-1)),'C',t2,vcat(0,N1-1+(1:N2-1)),'N',0,t3,1:(numind(t1)+numind(t2)-2)))
#
# ⊗{S}(t1::AbstractTensor{S},t2::AbstractTensor{S})=tensorproduct(t1,1:numind(t1),t2,numind(t1)+(1:numind(t2)))
# Base.trace{S,T}(t::AbstractTensor{S,T,2})=scalar(tensortrace(t,[1,1],[]))
#
# # general tensor operations: no error checking, pass to mutating methods
# function tensorcopy(A::AbstractTensor,labelsA,outputlabels=labelsA)
#     spaceA=space(A)
#     spaceC=spaceA[indexin(outputlabels,labelsA)]
#     C=similar(A,spaceC)
#     tensorcopy!(A,labelsA,C,outputlabels)
#     return C
# end
# function tensoradd{S,TA,TB,N}(A::AbstractTensor{S,TA,N},labelsA,B::AbstractTensor{S,TB,N},labelsB,outputlabels=labelsA)
#     spaceA=space(A)
#     spaceC=spaceA[indexin(outputlabels,labelsA)]
#     T=promote_type(TA,TB)
#     C=similar(A,T,spaceC)
#     tensorcopy!(A,labelsA,C,outputlabels)
#     tensoradd!(1,B,labelsB,1,C,outputlabels)
#     return C
# end
# function tensortrace(A::AbstractTensor,labelsA,outputlabels)
#     T=eltype(A)
#     spaceA=space(A)
#     spaceC=spaceA[indexin(outputlabels,labelsA)]
#     C=similar(A,spaceC)
#     tensortrace!(1,A,labelsA,0,C,outputlabels)
#     return C
# end
# function tensortrace(A::AbstractTensor,labelsA) # there is no one-line method to compute the default outputlabels
#     ulabelsA=unique(labelsA)
#     labelsC=similar(labelsA,0)
#     sizehint(labelsC,length(ulabelsA))
#     for j=1:length(ulabelsA)
#         ind=findfirst(labelsA,ulabelsA[j])
#         if findnext(labelsA,ulabelsA[j],ind+1)==0
#             push!(labelsC,ulabelsA[j])
#         end
#     end
#     return tensortrace(A,labelsA,labelsC)
# end
# function tensorcontract{S}(A::AbstractTensor{S},labelsA,conjA::Char,B::AbstractTensor{S},labelsB,conjB::Char,outputlabels=symdiff(labelsA,labelsB);method::Symbol=:BLAS,buffer::TCBuffer=defaultcontractbuffer)
#     spaceA=conjA=='C' ? conj(space(A)) : space(A)
#     spaceB=conjB=='C' ? conj(space(B)) : space(B)
#     spaceC=(spaceA ⊗ spaceB)[indexin(outputlabels,vcat(labelsA,labelsB))]
#     T=promote_type(eltype(A),eltype(B))
#     C=similar(A,T,spaceC)
#     tensorcontract!(1,A,labelsA,conjA,B,labelsB,conjB,0,C,outputlabels;method=method,buffer=buffer)
#     return C
# end
# tensorcontract{S}(A::AbstractTensor{S},labelsA,B::AbstractTensor{S},labelsB,outputlabels=symdiff(labelsA,labelsB);
#     method::Symbol=:BLAS,buffer::TCBuffer=defaultcontractbuffer)=tensorcontract(A,labelsA,'N',B,labelsB,'N',outputlabels;method=method,buffer=buffer)
#
# function tensorproduct{S}(A::AbstractTensor{S},labelsA,B::AbstractTensor{S},labelsB,outputlabels=vcat(labelsA,labelsB))
#     spaceA=space(A)
#     spaceB=space(B)
#
#     spaceC=(spaceA ⊗ spaceB)[indexin(outputlabels,vcat(labelsA,labelsB))]
#     T=promote_type(eltype(A),eltype(B))
#     C=similar(A,T,spaceC)
#     tensorproduct!(1,A,labelsA,B,labelsB,0,C,outputlabels)
#     return C
# end

# Factorization
#---------------
# Base.svd{S,T}(t::AbstractTensor{S,T,2}, truncation::TruncationScheme = notrunc()) = svd(t, 1, 2, truncation)
# leftorth{S,T}(t::AbstractTensor{S,T,2}, truncation::TruncationScheme = notrunc()) = leftorth(t, 1, 2, truncation)
# rightorth{S,T}(t::AbstractTensor{S,T,2}, truncation::TruncationScheme = notrunc()) = rightorth(t, 1, 2, truncation)
#
# # general tensor factorizations: permute to correct order and pass to in place methods
# """
#     svd(t::AbstractTensor, leftind, rightind = setdiff(1:numind(t),leftind), truncation::TruncationScheme = notrunc()) -> U,S,V'
#
# Create orthonormal basis `U` for indices in `leftind` and orthonormal basis `V'`for indices in
# `rightind`, and a diagonal tensor with singular values `S`, such that tensor `t`
# (permuted into index order `vcat(leftind,rightind)`) can be written as `U*S*V'`.
#
# A truncation parameter can be specified for the new internal dimension, in which case
# a singular value decomposition will be performed. See `svd(!)` for further information.
# """
# Base.svd(t::AbstractTensor,leftind,truncation::TruncationScheme=notrunc())=svd(t,leftind,setdiff(1:numind(t),leftind),truncation)
# function Base.svd(t::AbstractTensor,leftind,rightind,truncation::TruncationScheme=notrunc())
#     N = numind(t)
#     p = vcat(leftind,rightind)
#     (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
#     newt = tensorcopy(t, 1:N, p)
#     return svd!(newt, length(leftind), truncation)
# end
#
# """
#     leftorth(t::AbstractTensor, leftind, rightind = setdiff(1:numind(t),leftind), truncation::TruncationScheme = notrunc()) -> Q,R
#
# Create orthonormal basis `Q` for indices in `leftind`, and remainder `R` such that tensor `t`
# (permuted into index order `vcat(leftind,rightind)`) can be written as `Q*R`.
#
# This decomposition should be unique, such that it always returns the same result for the
# same input tensor `t`. The QR decomposition is fastest but only unique after correcting for
# phases. A truncation parameter can be specified for the new internal dimension, in which case
# a singular value decomposition will be performed. See `svd(!)` for further information.
# """
# leftorth(t::AbstractTensor, leftind, truncation::TruncationScheme) = leftorth(t, leftind, setdiff(1:numind(t),leftind), truncation)
# function leftorth(t::AbstractTensor, leftind, rightind=setdiff(1:numind(t),leftind), truncation::TruncationScheme = notrunc())
#     N = numind(t)
#     p = vcat(leftind,rightind)
#     (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
#     newt = tensorcopy(t, 1:N, p)
#     return leftorth!(newt, length(leftind), truncation)
# end
#
# """
#     rightorth(t::AbstractTensor, leftind, rightind = setdiff(1:numind(t),leftind), truncation::TruncationScheme = notrunc()) -> L,Q
#
# Create orthonormal basis `Q'` for indices in `rightind`, and remainder `L` such that tensor `t`
# (permuted into index order `vcat(leftind,rightind)`) can be written as `L*Q`.
#
# This decomposition should be unique, such that it always returns the same result for the
# same input tensor `t`. The LQ decomposition is fastest but only unique after correcting for
# phases. A truncation parameter can be specified for the new internal dimension, in which case
# a singular value decomposition will be performed. See `svd(!)` for further information.
# """
# rightorth(t::AbstractTensor, leftind, truncation::TruncationScheme) = rightorth(t, leftind, setdiff(1:numind(t), leftind), truncation)
# function rightorth(t::AbstractTensor, leftind, rightind = setdiff(1:numind(t), leftind), truncation::TruncationScheme = notrunc())
#     N = numind(t)
#     p = vcat(leftind, rightind)
#     (isperm(p) && length(p)==N) || throw(IndexError("Not a valid bipartation of the tensor indices"))
#     newt = tensorcopy(t, 1:N, p)
#     return rightorth!(newt, length(leftind), truncation)
# end
#
# svd!(t::AbstractTensor, n::Int) = svd!(t, n, notrunc())
#
# """
#     leftorth!(t::AbstractTensor, n::Int, truncation::TruncationScheme = notrunc()) -> Q,R
#
# Create orthonormal basis `Q` for the first `n` indices and remainder `R` such that tensor `t` can be written as `Q*R`.
#
# This decomposition should be unique, such that it always returns the same result for the
# same input tensor `t`. The QR decomposition is fastest but only unique after correcting for
# phases. A truncation parameter can be specified for the new internal dimension, in which case
# a singular value decomposition will be performed. See `svd(!)` for further information.
#
# The data in input tensor `t` is overwritten for the computation of `Q` and `R`.
# """
# leftorth!(t::AbstractTensor, n::Int) = leftorth!(t, n, notrunc())
#
# """
#     rightorth!(t::AbstractTensor, n::Int, truncation::TruncationScheme = notrunc()) -> Q,R
#
# Create orthonormal basis `Q'` for the indices `n+1:numind(t)` and remainder `L` such that tensor `t` can be written as `L*Q`.
#
# This decomposition should be unique, such that it always returns the same result for the
# same input tensor `t`. The LQ decomposition is fastest but only unique after correcting for
# phases. A truncation parameter can be specified for the new internal dimension, in which case
# a singular value decomposition will be performed. See `svd(!)` for further information.
#
# The data in input tensor `t` is overwritten for the computation of `L` and `Q`.
# """
# rightorth!(t::AbstractTensor, n::Int) = rightorth!(t, n, notrunc())


text:
struct DTWWorkspace{T,N,AT<:AbstractArray,D}
    q::AT
    dist::D
    r::Int
    l::Vector{T}
    u::Vector{T}
    l_buff::Vector{T}
    u_buff::Vector{T}
    cb::Vector{T}
    c1::Vector{T}
    c2::Vector{T}
    function DTWWorkspace(q::AbstractArray{QT}, dist, r::Int, ::Type{N}=Nothing) where {QT, N}
        T      = floattype(QT)
        m      = lastlength(q)
        n      = 2r + 1
        l      = zeros(T, m)
        u      = zeros(T, m)
        l_buff = zeros(T, m)
        u_buff = zeros(T, m)
        cb     = zeros(T, m)
        c1     = zeros(T, n)
        c2     = zeros(T, n)
        new{T, N, typeof(q), typeof(dist)}(q, dist, r, l, u, l_buff, u_buff, cb, c1, c2)
    end
end

struct DTWSearchResult{QT,C,D} <: AbstractSearchResult{QT}
    q::QT
    cost::C
    loc::Int
    prunestats
    dists::D
end

SlidingDistancesBase.value(r::DTWSearchResult) = r.cost
SlidingDistancesBase.location(r::DTWSearchResult) = r.loc
SlidingDistancesBase.payload(r::DTWSearchResult) = r.dists
SlidingDistancesBase.target(r::DTWSearchResult) = r.q

Base.findmin(results::Vector{<:DTWSearchResult}) = (i=argmin(results); (results[i].cost,i))
Base.findmax(results::Vector{<:DTWSearchResult}) = _findres(results, >)
Base.minimum(results::Vector{<:DTWSearchResult}) = findmin(results)[1]
Base.maximum(results::Vector{<:DTWSearchResult}) = findmax(results)[1]

function _findres(results::Vector{<:DTWSearchResult}, comp)
    mapreduce((a,b)->comp(a[1], b[1]) ? a : b, enumerate(results)) do (i,r)
        minimum(r.dists), i
    end
end

function lower_upper_envs!(w::DTWWorkspace{T}, q, bsf, query = false) where {T}
    du, dl = Deque{Int}(), Deque{Int}()
    push!(du, 0)
    push!(dl, 0)
    r = w.r
    if query
        u,l = w.u, w.l
    else
        u,l = w.u_buff, w.l_buff
    end
    m = lastlength(w.q)
    for i = 1:m-1
        if i > r
            u[i-r] = q[first(du)+1]
            l[i-r] = q[first(dl)+1]
        end
        if q[i+1] > q[i]
            pop!(du)
            while (!isempty(du) && q[i+1] > q[last(du)+1])
                pop!(du)
            end
        else
            pop!(dl)
            while (!isempty(dl) && q[i+1] < q[last(dl)+1])
                pop!(dl)
            end
        end
        push!(du, i)
        push!(dl, i)
        if i == 2r + 1 + first(du)
            popfirst!(du)
        elseif (i == 2r + 1 + first(dl))
            popfirst!(dl)
        end
    end
    for i = m:m+r
        u[i-r] = q[first(du)+1]
        l[i-r] = q[first(dl)+1]
        if i - first(du) >= 2r + 1
            popfirst!(du)
        end
        if i - first(dl) >= 2r + 1
            popfirst!(dl)
        end
    end
end

function lb_endpoints(dist, q, buffer, best_so_far; kwargs...)
    m = lastlength(q)

    x1 = buffer[!,1]
    y1 = buffer[!,m]
    lb = dist(q[!,1], x1; kwargs...) + dist(q[!,m], y1; kwargs...)
    lb >= best_so_far && return lb

    x2 = buffer[!,2]
    d = min(dist(x2, q[!,1]; kwargs...), dist(x1,q[!,2]; kwargs...), dist(x2,q[!,2]); kwargs...)
    lb += d
    lb >= best_so_far && return lb

    y2 = buffer[!,m-1]
    d = min(dist(y2, q[!,m]; kwargs...), dist(y1,q[!,m-1]; kwargs...), dist(y2,q[!,m-1]); kwargs...)
    lb += d
    lb >= best_so_far && return lb

    return lb
    # TODO: can add more comparisons here
end

function lb_env!(w::DTWWorkspace{T}, buffer, best_so_far; kwargs...) where T
    lb = zero(T)
    q, dist, u, l = w.q, w.dist, w.u, w.l
    for i in 1:lastlength(q)
        x = buffer[!,i] # This function only supports data with natural ordering
        d = zero(T)
        if x > u[i]
            d = dist(x, u[i]; kwargs...)
        elseif x < l[i]
            d = dist(x, l[i]; kwargs...)
        end
        lb += d
        w.cb[i] = d
        lb > best_so_far && return lb
    end
    return lb
end

function rev_cumsum!(cb)
    @inbounds for k = length(cb)-1:-1:1
        cb[k] = cb[k+1] + cb[k]
    end
end

"""
    search_result = dtwnn(q, y, dist, rad, [normalizer::Type{Nothing}]; kwargs...)

Compute the nearest neighbor to `q` in `y`. An optinal normalizer type can be supplied, see, `ZNormalizer, DiagonalZNormalizer, NormNormalizer`.

# Arguments:
- `q`: query (the short time series)
- `y`: data ( the long time series)
- `dist`: distance
- `rad`: radius
- `showprogress`: Defaults to true
- `prune_endpoints = true`: use endpoint heuristic
- `prune_envelope  = true`: use envelope heuristic
- `bsf_multiplier  = 1`: If > 1, require lower bound to exceed `bsf_multiplier*best_so_far`.
- `saveall = false`: compute a dense result (takes longer, no early stopping methods used). If false, then a vector of lower bounds on the distance is stored in `search_result.dists`, if true, all distances are computed and stored.
- `avoid`: If an integer index (or set of indices) is provided, this index will be avoided in the search. This is useful in case `q` is a part of `y`.
"""
function dtwnn(q::AbstractArray{QT}, y, dist, rad, ::Type{N}=Nothing; kwargs...) where {QT,N}
    q, y = setup_normalizer(N, q, y)
    w = DTWWorkspace(q, dist, rad, N)::DTWWorkspace{floattype(QT),N,typeof(q),typeof(dist)}
    dtwnn(w, y; kwargs...)
end

function dtwnn(w::DTWWorkspace{T,normalizer}, y::AbstractArray;
    prune_endpoints = true,
    prune_envelope  = true,
    saveall         = false,
    bsf_multiplier  = 1,
    transportcost   = 1,
    showprogress    = true,
    avoid           = nothing,
    kwargs...) where {T,normalizer}

    normalizer !== Nothing && !isa(y, AbstractNormalizer) && @warn("Normalizer in use but `y` is not wrapped in a normalizer object. This will result in highly suboptimal performance", maxlog=10)
    bsf_multiplier >= 1 || throw(DomainError("It does not make sense to have the bsf_multiplier < 1"))
    best_so_far = typemax(T)
    best_loc    = 1
    q           = w.q
    m           = lastlength(q)
    my          = actuallastlength(y)
    my >= m || throw(ArgumentError("q must be shorter than y, swap inputs."))
    onedim      = ndims(q) == 1 && eltype(q) <: Real
    onedim && prune_envelope && lower_upper_envs!(w, q, best_so_far, true) # Result stored in w

    # Counters to keep track of how many times lb helps
    prune_end   = 0
    prune_env   = 0
    dists       = fill(typemax(T), my-m+1)

    prog = Progress((my-m)÷max(my÷100,1), dt=1, desc="DTW NN")
    # @inbounds @showprogress 1.5 "DTW NN" for it = 1:my-m
    @inbounds for it = 1:my-m+1
        showprogress && it % max(my÷100,1) == 0 && next!(prog)
        advance!(y)
        avoid !== nothing && it ∈ avoid && continue
        bsf = bsf_multiplier*best_so_far
        ym = getwindow(y, m, it)
        if prune_endpoints && !saveall
            lb_end = lb_endpoints(w.dist, w.q, ym, bsf; kwargs...)
            if lb_end > bsf
                prune_end += 1
                continue
            end
        end
        if onedim && prune_envelope && !saveall # This bound only works when there is a natural ordering
            # lower_upper_envs!(w, ym, bsf) # This step is only required for reverse bound
            lb_env = lb_env!(w, ym, bsf; kwargs...) # updates w.cb
            rev_cumsum!(w.cb)
            if lb_env > bsf
                prune_env += 1
                continue
            end
        end
        # If we get here, we must normalize the entire y
        buffern = normalize(normalizer, ym) # This only normalizes what's not already normalized

        newdist = dtw_cost(q, buffern, w.dist, w.r;
            cumulative_bound = w.cb,
            best_so_far      = saveall ? typemax(T) : bsf,
            s1               = w.c1,
            s2               = w.c2,
            transportcost    = transportcost,
            kwargs...
        )
        dists[it] = newdist
        if newdist < best_so_far
            best_so_far = newdist
            best_loc = it
        end
    end
    prunestats = (prune_end=prune_end, prune_env=prune_env)
    DTWSearchResult(q, best_so_far, best_loc, prunestats, dists)
end

struct Neighbor{T}
    i::Int
    d::T
end

Base.isless(n1::Neighbor,n2::Neighbor) = isless(n1.d, n2.d)
Base.isless(n1, n2::Neighbor) = isless(n1, n2.d)
Base.isless(n1::Neighbor,n2) = isless(n1.d, n2)

"""
    dists, inds = sparse_distmat(y::Vector{<:AbstractVector{T}}, k, dist, radius; kwargs...) where T

Compute the `k` nearest neighbors between signals in `y`, corresponding to the `k` smallest entries in each row of the pairwise distance matrix. The return values are vectors of length-k vectors with the calculated distances and neighbor indices.

#Arguments:
- `y`: Vector of vectors containing the signals
- `k`: number of neighbors
- `dist`: the inner metric, e.g., `SqEuclidean()`
- `kwargs`: these are sent to `dtw_cost`.
- `showprogress = true`
"""
function sparse_distmat(
    y::AbstractVector{<:AbstractArray{S}},
    k,
    dist,
    rad;
    showprogress::Bool = true,
    kwargs...,
) where {S}
    T     = floattype(S)
    N     = length(y)
    INDS  = [zeros(Int, k) for _ = 1:N]
    DISTS = [zeros(T, k) for _ = 1:N]
    showprogress && (p = Progress(N^2, 1, "sparse_distmat"))
    for i = 1:N
        bsf   = typemax(T)
        dists = BinaryMaxHeap{Neighbor{T}}()
        for j = 1:N
            j == i && (showprogress && next!(p); continue)
            d = lb_endpoints(dist, y[i], y[j], bsf; kwargs...)
            if d < bsf
                d = dtw_cost(y[i], y[j], dist, rad; best_so_far = bsf, kwargs...)
            end
            push!(dists, Neighbor(j, T(d)))
            if length(dists) > k
                bsf = pop!(dists)
            end
            showprogress && next!(p)
        end

        for j = k:-1:1
            n = pop!(dists)
            INDS[i][j] = n.i
            DISTS[i][j] = n.d
        end
    end
    DISTS, INDS
end


text:
using Primes

function solve()
    n = parse(Int, readline())
    fs = Primes.factor(n)
    ret = 0
    for x in values(fs)
        i = 1
        while x >= i
            ret += 1
            x -= i
            i += 1
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
export prime_plane

"""
`prime_plane(n)` returns a finite projective plane of order `n` provided
`n` is prime.
"""
function prime_plane(p::Int)
    @assert isprime(p) "$p is not prime"
    alist = [ [a,b,1] for a=0:p-1 for b=0:p-1 ]
    blist = [ [a,1,0] for a=0:p-1 ]
    clist = [ [1,0,0] ]

    list = [alist; blist; clist ]

    n = length(list)

    A = zeros(Int,n,n)
    for i=1:n
        for j=1:n
            d = list[i]' * list[j]
            A[i,j] = mod(d[1],p) == 0
        end
    end
    return A
end


text:
arrayshape(::Type{Array{T,N}}, sz...) where {T,N} =
  Partial{Array{T,N}}(convert(Array{Any}, fill(T, sz)))

arrayshape(T::Type, sz...) = arrayshape(Array{T,length(sz)}, sz...)

@abstract Basic getindex(xs::Const{<:Array}, i::Const...) =
  Const(xs.value[map(i -> i.value, i)...])

@abstract Basic length(xs::Const) = Const(length(xs.value))
@abstract Basic length(xs::Partial{<:Array}) = Const(length(xs.value))
@abstract Basic length(xs::Shape{<:Array}) = Const(prod(size(xs)))
@abstract Basic length(xs::Array) = Int

@abstract Basic size(xs::Const) = Const(size(xs.value))
@abstract Basic size(xs::Partial{<:Array}) = Const(size(xs.value))
@abstract Basic size(xs::Shape{Array{T,N}}) where {T,N} = Const(size(xs))
@abstract Basic size(xs::AType{Array{T,N}}) where {T,N} = NTuple{N,Int}

@abstract Basic eltype(xs::AbstractArray{T}) where T = T

@partial Basic getindex(xs::Partial{<:Array}, i::Const...) =
  xs.value[map(i -> i.value, i)...]

@abstract Basic Vector{Any}() = Partial{Vector{Any}}([])

@partial Basic push!(xs::Partial{Vector{T}}, x::T) where T = (push!(xs.value, x); xs)

@partial Basic pop!(xs::Partial{Vector{T}}) where T = pop!(xs.value)

@pure Basic Colon(), similar

@abstract Basic function Broadcast.broadcasted(::Broadcast.AbstractArrayStyle, f, args...)
  A = Core.Compiler.return_type(broadcast, Tuple{widen(f),widen.(args)...})
  if args isa Tuple{Vararg{Union{Const,Shape,AType{<:Number}}}} && !(args isa Tuple{Vararg{AType{<:Number}}})
    return Shape{A}(Broadcast.broadcast_shape(size.(args)...))
  else
    return A
  end
end

@abstract Basic mapreduce(f, op, A; dims = :) =
  Core.Compiler.return_type(mapreduce, Tuple{widen(f),widen(op),widen(A)})


text:
using Base.Test
using Mimi

@defcomp Foo begin
    input = Parameter()
    intermed = Variable(index=[time])
end

function run_timestep(c::Foo, tt)
    c.Variables.intermed[tt] = c.Parameters.input
end

@defcomp Bar begin
    intermed = Parameter(index=[time])
    output = Variable(index=[time])
end

function run_timestep(c::Bar, tt)
    c.Variables.output[tt] = c.Parameters.intermed[tt]
end

m = Model()
setindex(m, :time, 1)
foo = addcomponent(m, Foo)
bar = addcomponent(m, Bar)

foo[:input] = 3.14
bar[:intermed] = foo[:intermed]

run(m)

@test m[:Bar, :output][1] == 3.14


text:
"""
hatfun(x,t,k)

Evaluate a piecewise linear "hat" function at `x`, where `t` is a
vector of n+1 interpolation nodes and `k` is an integer in 0:n
giving the index of the node where the hat function equals one.
"""

function hatfun(x,t,k)
    n = length(t)-1
    # Return correct node given mathematical index k, including
    # fictitious choices.
    function node(k)
        if k < 0
            2t[1]-t[2]
        elseif k > n
            2t[n+1]-t[n]
        else
            t[k+1]
        end
    end

    H1 = (x-node(k-1))/(node(k)-node(k-1))   # upward slope
    H2 = (node(k+1)-x)/(node(k+1)-node(k))   # downward slope

    H = min(H1,H2)
    return max(0,H)
end

"""
plinterp(t,y)

Create a piecewise linear interpolating function for data values in
`y` given at nodes in `t`.
"""
function plinterp(t,y)
n = length(t)-1
return x -> sum( y[k+1]*hatfun(x,t,k) for k in 0:n )
end

"""
spinterp(t,y)

Create a cubic not-a-knot spline interpolating function for data
values in `y` given at nodes in `t`.
"""
function spinterp(t,y)

    n = length(t)-1
    h = diff(t)         # differences of all adjacent pairs

    # Preliminary definitions.
    Z = zeros(n,n);
    In = I(n);  E = In[1:n-1,:];
    J = diagm(0=>ones(n),1=>-ones(n-1))
    H = diagm(0=>h)

    # Left endpoint interpolation:
    AL = [ In Z Z Z ]
    vL = y[1:n]

    # Right endpoint interpolation:
    AR = [ In H H^2 H^3 ];
    vR = y[2:n+1]

    # Continuity of first derivative:
    A1 = E*[ Z J 2*H 3*H^2 ]
    v1 = zeros(n-1)

    # Continuity of second derivative:
    A2 = E*[ Z Z J 3*H ]
    v2 = zeros(n-1)

    # Not-a-knot conditions:
    nakL = [ zeros(1,3*n) [1 -1 zeros(1,n-2)] ]
    nakR = [ zeros(1,3*n) [zeros(1,n-2) 1 -1] ]

    # Assemble and solve the full system.
    A = [ AL; AR; A1; A2; nakL; nakR ]
    v = [ vL; vR; v1; v2; 0; 0 ]
    z = A\v

    # Break the coefficients into separate vectors.
    rows = 1:n
    a = z[rows]
    b = z[n.+rows];  c = z[2*n.+rows];  d = z[3*n.+rows]
    S = [ Polynomial([a[k],b[k],c[k],d[k]]) for k = 1:n ]
    # This function evaluates the spline when called with a value
    # for x.
    function evaluate(x)
        k = findfirst(@. x<t)   # one greater than interval x belongs to
        k==1 && return NaN
        if isnothing(k)
            return x==t[end] ? y[end] : NaN
        end
        return S[k-1](x-t[k-1])
    end
    return evaluate
end

"""
fdweights(t,m)

Return weights for the `m`th derivative of a function at zero using
values at the nodes in vector `t`.
"""
function fdweights(t,m)
    # This is a compact implementation, not an efficient one.

    function weight(t,m,r,k)
        # Recursion for one weight. Input: t   nodes (vector) m
        # order of derivative sought r   number of nodes to use from
        # t (<= length(t)) k   index of node whose weight is found

        if (m<0) || (m>r)        # undefined coeffs must be zero
            c = 0
        elseif (m==0) && (r==0)  # base case of one-point interpolation
            c = 1
        else                     # generic recursion
            if k<r
                c = (t[r+1]*weight(t,m,r-1,k) -
                    m*weight(t,m-1,r-1,k))/(t[r+1]-t[k+1])
            else
                numer = r > 1 ? prod(t[r]-x for x in t[1:r-1]) : 1.0
                denom = r > 0 ? prod(t[r+1]-x for x in t[1:r]) : 1.0
                β = numer/denom
                c = β*(m*weight(t,m-1,r-1,r-1) - t[r]*weight(t,m,r-1,r-1))
            end
        end
        return c
    end

    r = length(t)-1
    w = zeros(size(t))
    return [ weight(t,m,r,k) for k=0:r ]
end

"""
trapezoid(f,a,b,n)

Apply the trapezoid integration formula for integrand `f` over
interval [`a`,`b`], broken up into `n` equal pieces. Returns
estimate, vector of nodes, and vector of integrand values at the
nodes.
"""
function trapezoid(f,a,b,n)
    h = (b-a)/n
    t = LinRange(a,b,n+1)
    y = f.(t)
    T = h * ( sum(y[2:n]) + 0.5*(y[1] + y[n+1]) )

    return T,t,y
end

"""
intadapt(f,a,b,tol)

Do adaptive integration to estimate the integral of `f` over
[`a`,`b`] to desired error tolerance `tol`. Returns estimate and a
vector of evaluation nodes used.
"""
function intadapt(f,a,b,tol)
    # Use error estimation and recursive bisection.
    function do_integral(a,fa,b,fb,m,fm,tol)
        # These are the two new nodes and their f-values.
        xl = (a+m)/2;  fl = f(xl);
        xr = (m+b)/2;  fr = f(xr);
        t = [a,xl,m,xr,b]          # all 5 nodes at this level

        # Compute the trapezoid values iteratively.
        h = (b-a)
        T = [0.,0.,0.]
        T[1] = h*(fa+fb)/2
        T[2] = T[1]/2 + (h/2)*fm
        T[3] = T[2]/2 + (h/4)*(fl+fr)

        S = (4*T[2:3]-T[1:2]) / 3      # Simpson values
        E = (S[2]-S[1]) / 15           # error estimate

        if abs(E) < tol*(1+abs(S[2]))  # acceptable error?
            Q = S[2]                   # yes--done
        else
            # Error is too large--bisect and recurse.
            QL,tL = do_integral(a,fa,m,fm,xl,fl,tol)
            QR,tR = do_integral(m,fm,b,fb,xr,fr,tol)
            Q = QL + QR
            t = [tL;tR[2:end]]   # merge the nodes w/o duplicate
        end
        return Q,t
    end

    m = (b+a)/2
    Q,t = do_integral(a,f(a),b,f(b),m,f(m),tol)
    return Q,t
end


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "X-3-5 scenario" begin

@time for X in 6:7
    filepath = "./data/quick_adjacency_decomposition/"
    filename = "$X-3-5"

    println("$filename scenario")

    PM = LocalSignaling(X,5,3)
    BG_seed = BellGame(cat(
        [1 0 0 0 0;1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0],
        zeros(Int64,(5,X-5)),
        dims=2
    ), 3)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices=200)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
include("services/s3.jl")
using .aws_s3

function high_level_main()
    buckets = aws_s3.ListBuckets()
    println(buckets)
end

max_stars_repo_name: nicoleepp/JuliaCloud-AWS-State

max_issues_repo_name: nicoleepp/JuliaCloud-AWS-State

max_forks_repo_name: nicoleepp/JuliaCloud-AWS-State


text:
function solve_008(n::Int64)

    # Hardcode string
    str = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"

    max_prod = 0
    # Iterate over string and multiply digits
    for j = 1:(1000-n)
        sub_str = str[j:j+n]
        prod = 1
        for i = 1:n
            prod *= parse(Int64,sub_str[i])
        end
        if prod > max_prod
            max_prod = prod
        end
    end

    return max_prod
end


text:
const Point{T} = SArray{Tuple{d},T,1,d}       # point in Rd
const Unc{T} = SArray{Tuple{d,d},T,2,d*d}     # Matrix presenting uncertainty
const PointF = Point{Float64}
const UncF = Unc{Float64}

struct NState{P, M<:AbstractMatrix{P}}
    x::M
end
const State = NState

NState(x::Vector) = NState(reshape(x, (2, length(x)>>1)))

import Base: axes, #=iterate,=# eltype, copy, copyto!, zero, eachindex, getindex, setindex!, size, vec

q(x::NState, i) = x.x[1, i]
p(x::NState, i) = x.x[2, i]

eltype(x::NState) = eltype(x.x)
deepeltype(x::Union{NState,Vector}) = eltype(eltype(x))
q(x::NState) = @view x.x[1:2:end]
p(x::NState) = @view x.x[2:2:end]

import LinearAlgebra: norm
norm(x::NState) = norm(vec(x))

vecofpoints2state(x::Vector) = NState(x)

size(s::NState) = size(s.x)
axes(s::NState, i) = axes(s.x, i)
deepvec(x::NState) = vec(reinterpret(deepeltype(x), x.x))
deepvec(x::Vector) = reinterpret(deepeltype(x), x)

function deepvec2state(x::Vector)
    x = reinterpret(Point{eltype(x)}, x)
    #dump(length(x))
    NState(reshape(x, (2, length(x)>>1)))
end

function deepvec2state(x::Base.ReshapedArray{<:Any,1,<:Base.ReinterpretArray})
    NState(x.parent.parent)
end

function NState(q::AbstractVector, p::AbstractVector)
    length(p) != length(q) && throw(DimensionMismatch())
    NState([((q,p)[i])[j] for i in 1:2, j in 1:length(p)])
end

vec(x::NState) = vec(x.x)
Base.broadcastable(x::NState) = x
Broadcast.BroadcastStyle(::Type{<:NState}) = Broadcast.Style{NState}()

"""
	split_state(x::NState)

Split a state into its `q` and `p` parts, where the latter are reinterpreted as vectors.

Note that merge does the inverse operation

## Example
```
y = State(rand(PointF,3),rand(PointF,3))
q, p = split_state(y)
yy = merge_state(q,p)
yy-y
```
"""
function split_state(x::NState)
	q =  vec(reinterpret(deepeltype(x), x.x[1,:]))
	p =  vec(reinterpret(deepeltype(x), x.x[2,:]))
	q, p
end

"""
	merge_state(q,p)

Merge `q` and `p` parts of state to an element of type `NState`.

Note that merge does the inverse operation

## Example
```
y = State(rand(PointF,3),rand(PointF,3))
q, p = split_state(y)
yy = merge_state(q,p)
yy-y
```
"""
function merge_state(q,p)
	NState(reinterpret(Point{eltype(q)}, q),reinterpret(Point{eltype(p)}, p))
end

function Base.getproperty(u::NState, s::Symbol)
    x = getfield(u, :x)
    s == :x && return x

    if s == :q
        @view x[1:2:end]
    elseif s == :p
        @view x[2:2:end]
    else
        throw(ArgumentError("NState has properties `p`, `q` or `x`"))
    end
end

zero!(v::NState) = zero!(v.x)
zero!(v::AbstractArray) = v .= Ref(zero(eltype(v)))
zero(v::NState) = NState(zero(v.x))

copy(x::NState) = NState(copy(x.x))
function copyto!(x::NState, y::NState)
    copyto!(x.x, y.x)
    x
end

function copyto!(x::NState, y)
    for i in eachindex(x)
        x[i] = y[i]
    end
    x
end

getindex(x::NState, I) = getindex(x.x, I)
function setindex!(x::NState, val, I)
    x.x[I] = val
end

eachindex(x::NState) = CartesianIndices(axes(x.x))

import Base: *, +, /, -
import LinearAlgebra: dot

#    import Bridge: outer, inner

function outer(x::NState, y::NState)
    [outer(x[i],y[j]) for i in eachindex(x), j in eachindex(y)]
end

*(c::Number, x::NState) = NState(c*x.x)
*(x::NState, c::Number) = NState(x.x*c)
+(x::NState, y::NState) = NState(x.x + y.x)
-(x::NState, y::NState) = NState(x.x - y.x)
/(x::NState, y) = NState(x.x/y)

dot(x::NState, y::NState) = dot(vec(x),vec(y))

function Base.:*(A::AbstractArray{<:Unc,2}, x::State)
    vecofpoints2state(A*vec(x))
end

q(i::Int) = 2i - 1
p(i::Int) = 2i

flipmomenta(x::NState) = NState(x.q, -x.p)

onemask(x) = onemask.(x)

"""
    onemask(x::Number) = one(x)
"""
onemask(x::Number) = one(x)


text:
using Statistics, Printf

using Oceananigans
using Oceananigans.OutputWriters
using Oceananigans.Diagnostics
using Oceananigans.AbstractOperations

using Oceananigans: Cell, Face, cell_advection_timescale

# Float type and architecture
FT = Float64
arch = GPU()

# Resolution
Nx = 256
Ny = 256
Nz = 256

# Domain size
Lx = 100
Ly = 100
Lz = 100

# Initial stratification
∂θ∂z = 0.01

# Simulation time
days = 8
end_time = day * days

# Data directory
base_dir = "ocean_convection_data"
mkpath(base_dir)

# Physical constants.
const f₀ = 1e-4  # Coriolis parameter [s⁻¹]
const ρ₀ = 1027  # Density of seawater [kg/m³]
const cₚ = 4000  # Specific heat capacity of seawater at constant pressure [J/(kg·K)]

# Surface heat flux (negative values imply cooling)
@inline Q(x, y, t) = -100
Q_str = "-100"

# Surface wind stress along the x-direction
@inline τx(x, y, t) = 0
τx_str = "0"

# Convert heat flux and wind stress to boundary conditions.
@inline Fu(x, y, t) = τx(x, y, t) / ρ₀
@inline Fθ(x, y, t) = -Q(x, y, t) / (ρ₀ * cₚ)

Fu_str = "τx(x, y, t) / ρ₀"
Fθ_str = "Q(x, y, t) / (ρ₀ * cₚ)"

u_top_bc = BoundaryFunction{:z, Face, Cell}(Fu)
θ_top_bc = BoundaryFunction{:z, Cell, Cell}(Fθ)

# Define boundary conditions
ubcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, u_top_bc))
θbcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, θ_top_bc),
                               bottom = BoundaryCondition(Gradient, ∂θ∂z))
Sbcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Gradient, 0),
                               bottom = BoundaryCondition(Gradient, ∂θ∂z))

# Create model
model = Model(float_type = FT,
            architecture = arch,
                    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                coriolis = FPlane(FT; f=f₀),
                buoyancy = SeawaterBuoyancy(FT; equation_of_state=LinearEquationOfState(β=0)),
                 closure = AnisotropicMinimumDissipation(FT),
     boundary_conditions = BoundaryConditions(u=ubcs, T=θbcs, S=Sbcs))

# Set initial condition.
ε(μ) = μ * randn() # noise
θ₀(x, y, z) = 20 + ∂θ∂z * z + ε(1e-10) * exp(4z/Lz)
S₀(x, y, z) = θ₀(x, y, z)

set_ic!(model, T=θ₀, S=S₀)

# Function that saves metadata with every output file.
function init_save_parameters_and_bcs(file, model)
    file["parameters/coriolis_parameter"] = f₀
    file["parameters/density"] = ρ₀
    file["parameters/specific_heat_capacity"] = cₚ
    file["parameters/viscosity"] = model.closure.ν
    file["parameters/diffusivity"] = model.closure.κ.T
    file["parameters/diffusivity_T"] = model.closure.κ.T
    file["parameters/diffusivity_S"] = model.closure.κ.S
    file["parameters/surface_cooling"] = Q_str
    file["parameters/temperature_gradient"] = ∂θ∂z
    file["parameters/wind_stress_x"] = τx_str
    file["boundary_conditions/top/FT"] = Fθ_str
    file["boundary_conditions/top/Fu"] = Fu_str
    file["boundary_conditions/bottom/dTdz"] = ∂θ∂z
end

# Saving 3D fields to JLD2 output files.
fields = Dict(
     :u => model -> Array(model.velocities.u.data.parent),
     :v => model -> Array(model.velocities.v.data.parent),
     :w => model -> Array(model.velocities.w.data.parent),
     :T => model -> Array(model.tracers.T.data.parent),
     :S => model -> Array(model.tracers.S.data.parent),
    :nu => model -> Array(model.diffusivities.νₑ.data.parent),
:kappaT => model -> Array(model.diffusivities.κₑ.T.data.parent),
:kappaS => model -> Array(model.diffusivities.κₑ.S.data.parent)
)

field_writer = JLD2OutputWriter(model, fields; dir=base_dir, prefix="ocean_convection_fields",
                                init=init_save_parameters_and_bcs,
                                max_filesize=100GiB, interval=6hour, force=true, verbose=true)
push!(model.output_writers, field_writer)

####
#### Set up diagnostics.
####

# NaN checker will abort simulation if NaNs are produced.
push!(model.diagnostics, NaNChecker(model; frequency=1000, fields=Dict(:w => model.velocities.w)))

# Time interval for computing and saving profiles.
Δtₚ = 10minute

# Define horizontal average diagnostics.
 Up = HorizontalAverage(model.velocities.u;       return_type=Array)
 Vp = HorizontalAverage(model.velocities.v;       return_type=Array)
 Wp = HorizontalAverage(model.velocities.w;       return_type=Array)
 Tp = HorizontalAverage(model.tracers.T;          return_type=Array)
 Sp = HorizontalAverage(model.tracers.S;          return_type=Array)
 νp = HorizontalAverage(model.diffusivities.νₑ;   return_type=Array)
κTp = HorizontalAverage(model.diffusivities.κₑ.T; return_type=Array)
κSp = HorizontalAverage(model.diffusivities.κₑ.S; return_type=Array)
dTp = HorizontalAverage(model.timestepper.Gⁿ.T;   return_type=Array)

u = model.velocities.u
v = model.velocities.v
w = model.velocities.w
T = model.tracers.T
S = model.tracers.S

uu = HorizontalAverage(u*u, model; return_type=Array)
vv = HorizontalAverage(v*v, model; return_type=Array)
ww = HorizontalAverage(w*w, model; return_type=Array)
uv = HorizontalAverage(u*v, model; return_type=Array)
uw = HorizontalAverage(u*w, model; return_type=Array)
vw = HorizontalAverage(v*w, model; return_type=Array)
wT = HorizontalAverage(w*T, model; return_type=Array)
wS = HorizontalAverage(w*S, model; return_type=Array)

# Create output writer that writes vertical profiles to JLD2 output files.
profiles = Dict(
     :u => model -> Up(model),
     :v => model -> Vp(model),
     :w => model -> Wp(model),
     :T => model -> Tp(model),
     :S => model -> Sp(model),
    :nu => model -> νp(model),
:kappaT => model -> κTp(model),
:kappaS => model -> κSp(model),
  :dTdt => model -> dTp(model),
    :uu => model -> uu(model),
    :vv => model -> vv(model),
    :ww => model -> ww(model),
    :uv => model -> uv(model),
    :uw => model -> uw(model),
    :vw => model -> vw(model),
    :wT => model -> wT(model),
    :wS => model -> wS(model)
)

profile_writer = JLD2OutputWriter(model, profiles; dir=base_dir, prefix="ocean_convection_profiles",
                                  init=init_save_parameters_and_bcs,
                                  interval=Δtₚ, max_filesize=25GiB, force=true, verbose=true)

push!(model.output_writers, profile_writer)

# Wizard utility that calculates safe adaptive time steps.
wizard = TimeStepWizard(cfl=0.25, Δt=3.0, max_change=1.2, max_Δt=30.0)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 50

while model.clock.time < end_time
    walltime = @elapsed time_step!(model; Nt=Ni, Δt=wizard.Δt)

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Calculate advective CFL number.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)
    CFL = wizard.Δt / cell_advection_timescale(model)

    # Calculate diffusive CFL number.
    νmax = maximum(model.diffusivities.νₑ.data.parent)
    κmax = maximum(model.diffusivities.κₑ.T.data.parent)

    Δ = min(model.grid.Δx, model.grid.Δy, model.grid.Δz)
    νCFL = wizard.Δt / (Δ^2 / νmax)
    κCFL = wizard.Δt / (Δ^2 / κmax)

    # Calculate a new adaptive time step.
    update_Δt!(wizard, model)

    # Print progress statement.
    @printf("[%06.2f%%] i: %d, t: %5.2f days, umax: (%6.3g, %6.3g, %6.3g) m/s, CFL: %6.4g, νκmax: (%6.3g, %6.3g), νκCFL: (%6.4g, %6.4g), next Δt: %8.5g, ⟨wall time⟩: %s\n",
            progress, model.clock.iteration, model.clock.time / day,
            umax, vmax, wmax, CFL, νmax, κmax, νCFL, κCFL,
            wizard.Δt, prettytime(walltime / Ni))
end


text: ~f=0∉3:length(f).|>i->f[i]==f[i-2]+f[i-1]

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus <: SwaggerModel
    certificate::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: certificate
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCertificatesV1beta1CertificateSigningRequestCondition} } # spec name: conditions

    function IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus(;certificate=nothing, conditions=nothing)
        o = new()
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus, Symbol("certificate"), certificate)
        setfield!(o, Symbol("certificate"), certificate)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        o
    end
end # type IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus

const _property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus = Dict{Symbol,Symbol}(Symbol("certificate")=>Symbol("certificate"), Symbol("conditions")=>Symbol("conditions"))
const _property_types_IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus = Dict{Symbol,String}(Symbol("certificate")=>"Vector{UInt8}", Symbol("conditions")=>"Vector{IoK8sApiCertificatesV1beta1CertificateSigningRequestCondition}")
Base.propertynames(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus }) = collect(keys(_property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus))
Swagger.property_type(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus }, property_name::Symbol) =  _property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus[property_name]

function check_required(o::IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus)
    true
end

function validate_property(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestStatus }, name::Symbol, val)
    if name === Symbol("certificate")
    end
end


text:
n = rand(1:38)
@test mdinfo(builtin(n)) != nothing
@test mdinfo(builtin(1:n)) != nothing


text:
## modeling of TCR-MHC-CD4 kinetics from adhesion frequency experiments
##
using Revise
src_dir = "C:/Users/stravaglino3/Downloads/julia_code/src/"
# includet(src_dir*"model_fit_global.jl")
includet(src_dir*"model_fit_global_v2.jl")
includet(src_dir*"models.jl")
includet(src_dir*"models_dissociation.jl")
includet(src_dir*"helper_functions.jl")
includet(src_dir*"binning_helper_fs.jl")
includet(src_dir*"global_plotting.jl")
includet(src_dir*"mle_fitting.jl")
using Plots
plot(0:10)
using StatsBase, KernelDensity, Distributions
using BlackBoxOptim
using DataFrames
# load data
data_dir = "C:/Users/stravaglino3/Downloads/julia_code/TCR_CD4_project/dissociation data/"
tm_j = CSV.read(data_dir*"jurkat_TCR_MHC.csv", DataFrame)
cm_j = CSV.read(data_dir*"jurkat_CD4_MHC.csv", DataFrame)
tmc_j = CSV.read(data_dir*"jurkat_TCR_MHC_CD4.csv", DataFrame)

tm_b = CSV.read(data_dir*"Bimolecular polylink interactions.txtbinned bond lifetime.csv", DataFrame)
cm_b = CSV.read(data_dir*"CD4HLA only.txtbinned bond lifetime.csv", DataFrame)
tmc_b = CSV.read(data_dir*"All from 2018-05-1_10_15.txtbinned bond lifetime.csv", DataFrame)

# clean data from 0s
tm_j_binned = clean_dissociation_data(tm_j)
cm_j_binned = clean_dissociation_data(cm_j)
tmc_j_binned = clean_dissociation_data(tmc_j)

tm_b = clean_dissociation_data(tm_b)
cm_b = clean_dissociation_data(cm_b)
tmc_b = clean_dissociation_data(tmc_b)
# cobine bins to complete dataset
F_tm_j, t_tm_j = combine_bins(tm_j_binned)
F_cm_j, t_cm_j = combine_bins(cm_j_binned)
F_tmc_j, t_tmc_j = combine_bins(tmc_j_binned)

F_tm_b, t_tm_b = combine_bins(tm_b)
F_cm_b, t_cm_b = combine_bins(cm_b)
F_tmc_b, t_tmc_b = combine_bins(tmc_b)

gr()
s1 = scatter(F_tm_j, t_tm_j,label="TCR-MHC");
s2 = scatter(F_cm_j, t_cm_j,label="CD4-MHC");
s3 = scatter(F_tmc_j, t_tmc_j,label="TCR-MHC-CD4");
plot(s1,s2,s3); title!("Jurkats")
s1 = scatter(F_tm_b, t_tm_b,label="TCR-MHC");
s2 = scatter(F_cm_b, t_cm_b,label="CD4-MHC");
s3 = scatter(F_tmc_b, t_tmc_b,label="TCR-MHC-CD4");
plot(s1,s2,s3); title!("Beads")

plotly()
plot()
Color=:black; plot_F_vs_t(tm_j)
Color=:blue; plot_F_vs_t(cm_j)
Color=:red; plot_F_vs_t(tmc_j)
title!("Jurkats")
plot()
Color=:black; plot_F_vs_t(tm_b)
Color=:blue; plot_F_vs_t(cm_b)
Color=:red; plot_F_vs_t(tmc_b)
title!("Beads")

##
# load new data from Kaitao and compare to Muaz data
cm_b_new = CSV.read(data_dir*"CD4_LT_bead_KL.csv", DataFrame)
F_cm_b_new = cm_b_new[:,1]; t_cm_b_new = cm_b_new[:,2]
scatter(F_cm_b, log.(t_cm_b),label="Muaz")
scatter!(F_cm_b_new, log.(t_cm_b_new), label="Kaitao")
xlabel!("F"); ylabel!("ln(t)"); title!("CD4-MHC (beads)")

## fit bimolecular TCR-MHC (on Jurkats)

tm_j  = F_t_Data(F_tm_b, t_tm_b, InterpKDE(kde(F_tm_b)))
cm_j   = F_t_Data(F_cm_b, t_cm_b, InterpKDE(kde(F_cm_b)))
tmc_j = F_t_Data(F_tmc_b, t_tmc_b, InterpKDE(kde(F_tmc_b)))

model = bimolecular_diss_bell_activated_state_w_cons!
u₀ = [] # don't provide initial states
u₀_type = "equilibrium" # calculate initial states based on equilibrium
tspan = (0.0,maximum(tm_j.t))
tspan = (0.0,10)
k₋ᵒf_tm = 0.74
M_tm_j = M(model,u₀,u₀_type,tspan,[],[k₋ᵒf_tm],0.0,Rodas5)

pinit = [   0.01,
        0.1,0.6,
        3.,0.2,
        150.0,4.0]
opt =  [        (0.,10.),
        (0.,1.),(0.,10.),
        (0.,10.),(0.,10.),
        (0.,500.),(0.,10.)]

# pinit = rand(8)*10
@time mle_loss(pinit,tm_j,M_tm_j)
using BlackBoxOptim

@time rates_tm_j, l_tm_j, hes_tm_j = optimization(M_tm_j,tm_j,opt)
rates_tm_j =[0.13772078920646594, 0.9999999882398696, 0.2216191952783784, 6.454133898001537e-5, 1.6550629048008287, 499.08524377458565, 9.935988576599074]
ks_tm_j = rates_tm_j
σ_tm_j = hessian2σ(hes_tm_j)
@show k_tm_j = rates_tm_j .± σ_tm_j

includet(src_dir*"mle_plotting.jl")
plotly()
plot_ft_fit(M_tm_j,ks_tm_j,tm_j)

# fitting is nonsensical

# try fitting again, this time allowing k₋ᵒf_tm to change as well
##
model = bimolecular_diss_bell!
u₀ = [] # don't provide initial states
u₀_type = "equilibrium" # calculate initial states based on equilibrium
tspan = (0.0,10)
M_tm_j  = M(model,u₀,u₀_type,tspan,[],[],0.0,Rodas5)

pinit = [14.1,0.01,
        0.1,0.6,
        3.,0.2,
        150.0,4.0]

opt =  [(0.,10.),(0.,10.),
        (0.,1.),(0.,10.),
        (0.,10.),(0.,10.),
        (0.,500.),(0.,10.)]

@time mle_loss(pinit,tm_j,M_tm_j )

@time rates_tm_j, l_tm_j, hes_tm_j = optimization(M_tm_j,tm_j,opt)
ks_tm_j = rates_tm_j
σ_tm_j = hessian2σ(hes_tm_j)
@show k_tm_j = rates_tm_j .± σ_tm_j

includet(src_dir*"mle_plotting.jl")
plotly()
plot_ft_fit(M_tm_j,ks_tm_j,tm_j)
plot(); Color = :blue
plot_mean_ft_fit(M_tm_j,ks_tm_j,tm_j,tm_j_binned;species="T-M")

## try fitting kaitao cd4 data


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

using Pkg # hideall
Pkg.activate("_literate/A-ensembles-3/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

using MLJ
using PyPlot
ioff() # hide
import Statistics

Xs = source()
ys = source()

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
atom = DecisionTreeRegressor()

machines = (machine(atom, Xs, ys) for i in 1:100)

Statistics.mean(v...) = mean(v)
Statistics.mean(v::AbstractVector{<:AbstractNode}) = node(mean, v...)

yhat = mean([predict(m, Xs) for  m in machines]);

surrogate = Deterministic()
mach = machine(surrogate, Xs, ys; predict=yhat)

@from_network mach begin
    mutable struct OneHundredModels
        atom=atom
    end
end

one_hundred_models = OneHundredModels()

X, y = @load_boston;

r = range(atom,
          :min_samples_split,
          lower=2,
          upper=100, scale=:log)

mach = machine(atom, X, y)

figure()
curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e1.svg")) # hide

r = range(one_hundred_models,
          :(atom.min_samples_split),
          lower=2,
          upper=100, scale=:log)

mach = machine(one_hundred_models, X, y)

figure()
curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)
plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e2.svg")) # hide

PyPlot.close_figs() # hide

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using BinaryProvider

# BinaryProvider support
const prefix = Prefix(joinpath(dirname(@__FILE__),"usr"))
const platform = platform_key()

# These are the two binary objects we care about
libfoo = LibraryProduct(prefix, "libfoo")
fooifier = ExecutableProduct(prefix, "fooifier")

# This is where we download things from, for different platforms
const bin_prefix = "https://github.com/staticfloat/small_bin/raw/74b7fd81e3fbc8963b14b0ebbe5421e270d8bdcf"
const download_info = Dict(
    Linux(:i686) =>     ("$bin_prefix/libfoo.i686-linux-gnu.tar.gz", "1398353bcbbd88338189ece9c1d6e7c508df120bc4f93afbaed362a9f91358ff"),
    Linux(:x86_64) =>   ("$bin_prefix/libfoo.x86_64-linux-gnu.tar.gz", "b9d57a6e032a56b1f8641771fa707523caa72f1a2e322ab99eeeb011f13ad9f3"),
    Linux(:aarch64) =>  ("$bin_prefix/libfoo.aarch64-linux-gnu.tar.gz", "19d9da0e6e7fb506bf4889eb91e936fda43493a39cd4fd7bd5d65506cede6f95"),
    Linux(:armv7l) =>   ("$bin_prefix/libfoo.arm-linux-gnueabihf.tar.gz", "8e33c1a0e091e6e5b8fcb902e5d45329791bb57763ee9cbcde49c1ec9bd8532a"),
    Linux(:ppc64le) =>  ("$bin_prefix/libfoo.powerpc64le-linux-gnu.tar.gz", "b48a64d48be994ec99b1a9fb60e0af7f4415a57596518cb90a340987b79fad81"),
    MacOS() =>          ("$bin_prefix/libfoo.x86_64-apple-darwin14.tar.gz", "661b71edb433ab334b0fef70db3b5c45d35f2b3bee0d244f54875f1ec899c10f"),
    Windows(:i686) =>   ("$bin_prefix/libfoo.i686-w64-mingw32.tar.gz", "3d4a8d4bf0169007a42d809a1d560083635b1540a1bc4a42108841dcb6d2aaea"),
    Windows(:x86_64) => ("$bin_prefix/libfoo.x86_64-w64-mingw32.tar.gz", "2d08fbc9a534cd021f36b6bbe86ddabb2dafbedeb589581240aa4a8c5b896055"),
)
if platform in keys(download_info)
    # Grab the url and tarball hash for this particular platform
    url, tarball_hash = download_info[platform]

    # Build a BinaryPackage from the metadata, and install it
    binpkg = BinaryPackage(url, tarball_hash, platform, [libfoo, fooifier])
    install(binpkg; prefix=prefix, force=true, verbose=true)

    # Finaly, write out a deps file containing paths to libfoo and fooifier
    @write_deps_file libfoo fooifier
else
    error("Your platform $(Sys.MACHINE) is not recognized, we cannot install Libfoo!")
end


text:
"""
    train_val_test_inds(indices, ratios=(0.6,0.2,0.2); seed=nothing)

Split indices.
"""
function train_val_test_inds(indices, ratios=(0.6,0.2,0.2); seed=nothing)
    (sum(ratios) ≈ 1 && length(ratios) == 3) ? nothing :
    	error("ratios must be a vector of length 3 that sums up to 1")

    # set seed
    (seed == nothing) ? nothing : Random.seed!(seed)

    # set number of samples in individual subsets
    n = length(indices)
    ns = cumsum([x for x in floor.(Int, n .* ratios)])

    # scramble indices
    _indices = sample(indices, n, replace=false)

    # restart seed
    (seed == nothing) ? nothing : Random.seed!()

    # return the sets of indices
    _indices[1:ns[1]], _indices[ns[1]+1:ns[2]], _indices[ns[2]+1:ns[3]]
end

"""
	train_val_test_split(data_normal, data_anomalous, ratios=(0.6,0.2,0.2); seed=nothing,
	    contamination::Real=0.0)

Split data.
"""
function train_val_test_split(data_normal, data_anomalous, ratios=(0.6,0.2,0.2);
	seed=nothing, contamination::Real=0.0)

	# split the normal data, add some anomalies to the train set and divide
	# the rest between validation and test
	(0 <= contamination <= 1) ? nothing : error("contamination must be in the interval [0,1]")
	nd = ndims(data_normal) # differentiate between 2D tabular and 4D image data

	# split normal indices
	indices = 1:size(data_normal, nd)
	split_inds = train_val_test_inds(indices, ratios; seed=seed)

	# select anomalous indices
	na = size(data_anomalous, nd)
	indices_anomalous = 1:na
	na_tr = floor(Int, length(split_inds[1])*contamination/(1-contamination))
	(na_tr > na) ? error("selected contamination rate $contamination is too high, not enough anomalies available") : nothing
	tr = na_tr/length(indices_anomalous) # training ratio
	vtr = (1 - tr)/2 # validation/test ratio
	split_inds_anomalous = train_val_test_inds(indices_anomalous, (tr, vtr, vtr); seed=seed)

	# this can be done universally - how?
	if nd == 2
		tr_n, val_n, tst_n = map(is -> data_normal[:,is], split_inds)
		tr_a, val_a, tst_a = map(is -> data_anomalous[:,is], split_inds_anomalous)
	elseif nd == 4
		tr_n, val_n, tst_n = map(is -> data_normal[:,:,:,is], split_inds)
		tr_a, val_a, tst_a = map(is -> data_anomalous[:,:,:,is], split_inds_anomalous)
	end

	# cat it together
	tr_x = cat(tr_n, tr_a, dims = nd)
	val_x = cat(val_n, val_a, dims = nd)
	tst_x = cat(tst_n, tst_a, dims = nd)

	# now create labels
	tr_y = vcat(zeros(Float32, size(tr_n, nd)), ones(Float32, size(tr_a,nd)))
	val_y = vcat(zeros(Float32, size(val_n, nd)), ones(Float32, size(val_a,nd)))
	tst_y = vcat(zeros(Float32, size(tst_n, nd)), ones(Float32, size(tst_a,nd)))

	(tr_x, tr_y), (val_x, val_y), (tst_x, tst_y)
end

"""
	load_data(dataset::String, ratios=(0.6,0.2,0.2); seed=nothing,
	method="leave-one-out", contamination::Real=0.0)

Returns 3 tuples of (data, labels) representing train/validation/test part. Arguments are the splitting
ratios for normal data, seed and training data contamination.

For a list of available datasets, check `GenerativeAD.Datasets.uci_datasets`, `GenerativeAD.Datasets.other_datasets`
and `GenerativeAD.Datasets.mldatasets`.
"""
function load_data(dataset::String, ratios=(0.6,0.2,0.2); seed=nothing,
	method="leave-one-out", contamination::Real=0.0, kwargs...)
	any(method .== ["leave-one-out","leave-one-in"]) ? nothing : error("unknown method, choose one of `leave-one-in`, `leave-one-out`")
	(method ==  "leave-one-in" && !(dataset in mldatasets)) ? error("`leave-one-in` only implemented for image datasets") : nothing

	# extract data and labels
	if dataset in uci_datasets # UCI Loda data, standardized
		data_normal, data_anomalous = load_uci_data(dataset; kwargs...)
	elseif dataset in mldatasets # MNIST,FMNIST, SVHN2, CIFAR10
		data_normal, data_anomalous = load_mldatasets_data(dataset; kwargs...)
	elseif dataset in other_datasets
		data_normal, data_anomalous = load_other_data(dataset; standardize=true, kwargs...)
	else
		error("Dataset $(dataset) not known, either not implemented or misspeled.")
		# TODO add the rest
	end

	# now do the train/validation/test split
	if method=="leave-one-in" # in this case, we swap the anomalous nad normal data
		return train_val_test_split(data_anomalous, data_normal, ratios; seed=seed, contamination=contamination)
	else
		return train_val_test_split(data_normal, data_anomalous, ratios; seed=seed, contamination=contamination)
	end
end

"""
	vectorize(data)

Vectorizes the 4D image data returned from `load_data`.
"""
vectorize(data) = map(d->(reshape(d[1], :, size(d[1],4)), d[2]),data)