text:
# # Iris data

# In this example, we create a logistic regression model that classifies iris flowers.
# It consists of a [single-layer neural network](https://sebastianraschka.com/faq/docs/logisticregr-neuralnet.html)
# that outputs **three** probabilities (one for each species of iris flowers).
# We use Fisher's classic dataset to train the model. This dataset is retrieved from
# the [UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/iris).

# In Machine Learning, a classification task can be performed by a logistic regression model.
# However, we also can create a logistic regression model as a single-layer neural network.
# This neural network has the following characteristics:

# * Uses the [logitcrossentropy](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.logitcrossentropy) loss function.
# * Expects the class labels of the iris flowers encoded using [One-Hot encoding](https://fluxml.ai/Flux.jl/stable/data/onehot/#One-Hot-Encoding).
# * Outputs the index in the output vector with the highest value as the class label using
# [onecold](https://fluxml.ai/Flux.jl/stable/data/onehot/#Flux.onecold) which is the inverse operation of One-Hot encoding.

# To run this example, we need the following packages:

using Flux
using Flux: logitcrossentropy, normalise, onecold, onehotbatch
using Statistics: mean
using Parameters: @with_kw

# We set default values for the learning rate *lr* (for the training routine) and the number of
# times that we repeat the train data (more information below):

@with_kw mutable struct Args
    lr::Float64 = 0.5
    repeat::Int = 110
end

# ## Data

# We create the function `get_processed_data` to load the iris data, preprocess
# it (normalize and One-Hot encode the class labels), and split it into train and test datasets.

function get_processed_data(args)
    labels = Flux.Data.Iris.labels()
    features = Flux.Data.Iris.features()

    ## Subract mean, divide by std dev for normed mean of 0 and std dev of 1.
    normed_features = normalise(features, dims=2)

    klasses = sort(unique(labels))
    onehot_labels = onehotbatch(labels, klasses)

    ## Split into training and test sets, 2/3 for training, 1/3 for test.
    train_indices = [1:3:150 ; 2:3:150]

    X_train = normed_features[:, train_indices]
    y_train = onehot_labels[:, train_indices]

    X_test = normed_features[:, 3:3:150]
    y_test = onehot_labels[:, 3:3:150]

    ## Repeat the data `args.repeat` times
    train_data = Iterators.repeated((X_train, y_train), args.repeat)
    test_data = (X_test,y_test)

    return train_data, test_data
end

# The iris data is a 4×150 matrix. This means that the iris data has 150 examples,
# and each example has four features as well as a class label.
# After normalizing and encoding the data, the `get_processed_data` function divides it into train and test data.
# Also, it repeats the examples in the train data so that we have more data to train the neural network.

# ## Metrics

# We use two functions to assess the output of the model: `accuracy` and `confusion matrix`.
# The [accuracy function](https://developers.google.com/machine-learning/crash-course/classification/accuracy)
# measures the percentage of the labels that the model classified correctly.
# On the other hand, the [confusion matrix](https://machinelearningmastery.com/confusion-matrix-machine-learning/)
# is a table that summarises how good the model is for predicting data.

accuracy(x, y, model) = mean(onecold(model(x)) .== onecold(y))

function confusion_matrix(X, y, model)
    ŷ = onehotbatch(onecold(model(X)), 1:3)
    y * transpose(ŷ)
end

# ## Train function

# We define the `train` function that defines the model and trains it:

function train(; kws...)
    ## Initialize hyperparameter arguments
    args = Args(; kws...)

    ## Load processed data
    train_data, test_data = get_processed_data(args)

    ## #Declare model taking 4 features as inputs and outputting 3 probabiltiies,
    ## one for each species of iris.
    model = Chain(Dense(4, 3))

    ## Define loss function to be used in training
    ## For numerical stability, we use here logitcrossentropy
    loss(x, y) = logitcrossentropy(model(x), y)

    ## Training
    ## Gradient descent optimiser with learning rate `args.lr`
    optimiser = Descent(args.lr)

    println("Starting training.")
    Flux.train!(loss, params(model), train_data, optimiser)

    return model, test_data
end

# The function above loads the train and test data.
# Then, it creates the model as a single-layer network that expects as an input
# a four-element vector (features) and outputs a three-element vector
# (the number of classes of species of iris flowers).
# Also, it sets [logitcrossentropy](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.logitcrossentropy)
# as the loss function and the Gradient descent optimiser
# [Descent](https://fluxml.ai/Flux.jl/stable/training/optimisers/#Flux.Optimise.Descent).
# Finally, it runs a training step with the
# [train! function](https://fluxml.ai/Flux.jl/stable/training/training/#Flux.Optimise.train!).

# ## Test function

# After training the model, we define the `test` function that
# computes the model performance on the test data.
# It calls the `accuracy` function and displays the confusion matrix.
# This function validates that the model should achieve at least a 0.8 accuracy score.

function test(model, test)
    ## Testing model performance on test data
    X_test, y_test = test
    accuracy_score = accuracy(X_test, y_test, model)

    println("\nAccuracy: $accuracy_score")

    ## Sanity check.
    @assert accuracy_score > 0.8

    ## To avoid confusion, here is the definition of a
    ## Confusion Matrix: https://en.wikipedia.org/wiki/Confusion_matrix
    println("\nConfusion Matrix:\n")
    display(confusion_matrix(X_test, y_test, model))
end

# ## Run the example

# We call the `train` function to run the iris data example and compute the model performance:

cd(@__DIR__)
model, test_data = train()
test(model, test_data)


text:
"""
linear_aa()

Test aasol.jl for a two-dimensional linear problem
"""
function linear_aa()
maxit = 10
maxm = 2
vdim = 3 * maxm + 3
Vstore = zeros(2, vdim)
vdimS = 2 * maxm + 4
VstoreS = zeros(2, vdimS)
eigs=[.1, .5]
xstar=ones(2,)
pdata=makeLinpdata(eigs)
m=0
#
# Test for termination on entry
#
x0=[1.0, 1.0]; m=0;
aout = aasol(GLin!, x0, m, Vstore;
     rtol = 1.e-10, pdata=pdata, maxit=maxit, keepsolhist=true);
# Same results with low storage method?
aoutS = aasol(GLin!, x0, m, VstoreS;
     rtol = 1.e-10, pdata=pdata, maxit=maxit, keepsolhist=true);
linerrS = (
        norm(aoutS.history-aout.history,Inf) +
        norm(aoutS.stats.condhist-aout.stats.condhist,Inf) +
        norm(aoutS.stats.alphanorm-aout.stats.alphanorm,Inf)
    )
tflag = (aout.errcode===-1) && aout.idid && (linerrS < 1.e-15)
tflag || println("Failure in aasol terminate on entry test.")
#
# Test for failure to converge
#
x0=[2.0, 10.0]
m=0
maxit=10
aout = aasol(GLin!, x0, m, Vstore;
     rtol = 1.e-10, pdata=pdata, maxit=maxit, keepsolhist=true);
failflag = ~aout.idid && (aout.errcode == 10)
failflag || println("Linear iteration failure in aasol test fails.")
#
# Test for convergence is two iterations.
#
m=2
aout = aasol(GLin!, x0, m, Vstore;
     rtol = 1.e-10, pdata=pdata, maxit=maxit, keepsolhist=true);
termflag=(length(aout.history) == 4) &&
         (norm(aout.solution - xstar,Inf) < 1.e-14)
termflag || println("Terminate in two aa iterations test fails.")
#
# Now set the eigenvalues to [2.0, 10.0] and beta=-1/9
#
eigs=[2.0, 10.0]
pdata=makeLinpdata(eigs)
beta=-1.0/9.0
maxit=10
m=1
aout1 = aasol(GLin!, x0, m, Vstore;
     rtol = 1.e-10, pdata=pdata, maxit=maxit)
aout2 = aasol(GLin!, x0, m, Vstore;
     rtol = 1.e-10, pdata=pdata, maxit=maxit, beta=beta)
bflag= ~aout1.idid && aout2.idid && (length(aout2.history)==8)
return tflag && failflag && termflag && bflag

end

function GLin!(gout,xin,pdata)
M=pdata.M
b=pdata.b
gout=M*xin+b
return gout
end

function GLinBeta!(gout,xin,pdata)
M=pdata.M
b=pdata.b
beta=pdata.beta
gout = Glin!(gout,xin,pdata)
gout .*=beta
gout .+= (1.0-beta)*xin
return gout
end

function makeLinpdata(eigs, beta=1.0)
U=[1 -1; 1 1]./sqrt(2.0)
V=[3 -4; 4 3]./5.0
S=diagm(eigs)
M=U*S*V'
b=(I-M)*ones(2,);
return (M=M, b=b, beta=beta)
end

ext: jl


text:
# * Dispersed fields
#   We calculate the effect of dispersion described by the Sellmeier
#   equations in the frequency domain, to which we transform via an
#   FFT. For this to be possible, we only allow the evaluation of the
#   field for a specified sampling frequency, i.e. we don't provide an
#   implementation for evaluating a =DispersedField= at an arbitrary
#   time point.

struct DispersedField <: AbstractField
    a::AbstractField
    m::Medium
    d::Unitful.Length
end

function show(io::IO, f::DispersedField)
    show(io, f.a)
    write(io, "\n  – dispersed through $(f.d) of $(f.m)")
end

disperse(a::AbstractField, m::Medium, d::Unitful.Length) =
    DispersedField(a, m, d)

(f::DispersedField)(t::Union{Unitful.Time,AbstractVector{<:Unitful.Time}}) =
    error("Dispersed fields can only be evaluated by specifying a sampling frequency")

function (f::DispersedField)(fs::Unitful.Frequency = default_sampling_frequency(f.a))
    F̂ = fft(f.a, fs)
    ifft(F̂.*dispersion(f.m, f.d, fftfreq(f.a, fs), frequency(f.a)))*(u"V"*u"m")
end

max_frequency(f::DispersedField) = max_frequency(f.a)

# This is, strictly speaking, not true, since dispersing a pulse will
# in general stretch it in the time domain. It is up to the user to
# ensure that the time window is large enough to contain even the
# stretched pulse, by setting the σmax (or friends) to a large enough
# value.
span(f::DispersedField) = span(f.a)

export disperse


text:
reload("SSH")
sock = connect("localhost", 10000)
session = SSH.connect(SSH.Session, sock)
SSH.negotiate_algorithms!(session)
SSH.client_dh_kex!(session)
SSH.enter_userauth!(session)
SSH.clientauth_list(session, "kfischer", "ssh-connection")
@assert SSH.clientauth_pubkey(session, "kfischer",
    joinpath(ENV["HOME"],".ssh","id_rsa.pub"),
    joinpath(ENV["HOME"],".ssh","id_rsa"))


text:
using Sobol

function sampling(lb::AbstractVector{T}, ub::AbstractVector{T}, ns::Int, gs::Function...) where T
    @assert length(lb) == length(ub)
    @assert all(lb .< ub)
    ndim = length(lb)
    s = SobolSeq(ndim)
    x = Matrix{T}(undef, ndim, ns)
    scalar = MinMaxScaler(lb, ub)
    for i in 1:ns
        while true
            xx = KrigingModel.inverse!(scalar,  next!(s))
            if all(map(g->g(xx)<=0, gs))
                x[:, i] .= xx
                break
            end
        end
    end
    return x
end


text:
using NNlib: conv, DenseConvDims, PoolDims
using Flux: maxpool

const to_NCHW = :([0, 3, 1, 2])
const to_NHWC = :([0, 2, 3, 1])

# Library of mathematical functions we consider primitive.
# TODO: store Julia functions and types, to avoid tensorflow.js-specific functions

matVecMul(args...) = *(args...)

@primitive Trace x::AbstractMatrix * y::AbstractMatrix =
  Staged(*, x, y)

jscall(::typeof(*), a, b) = jscall(:(math.matMul), b, a)

@primitive Trace x::AbstractMatrix * y::AbstractVector =
  Staged(matVecMul, x, y)

jscall(::typeof(matVecMul), a, b) = jscall(:(math.vectorTimesMatrix), b, a)

@primitive Trace x::AbstractArray + y::AbstractArray =
  Staged(+, x, y)

jscall(::typeof(+), a, b) = jscall(:(math.add), b, a)

# cat

concat1D(a, b) = vcat(a, b)
concat(a, b) = vcat(a, b)

@primitive Trace vcat(a::AbstractVector, b::AbstractVector) =
  Staged(concat1D, a, b)

@primitive Trace vcat(a::AbstractMatrix, b::AbstractMatrix) =
  Staged(concat, a, b)

jscall(::typeof(concat1D), a, b) = jscall(:(math.concat1d), jscall(tuple,a, b))
jscall(::typeof(concat), a, b) = jscall(:(math.concat), jscall(tuple, a, b), -1)

# softmax

@primitive Trace softmax(x::AbstractVecOrMat) =
  Staged(softmax, x)

jscall(::typeof(softmax), x) = jscall(:(math.softmax), x)

# Dense uses invoke which JuliaInterpreter does not enter
@primitive Trace function (a::Dense)(x::AbstractArray)
  W, b, σ = a.W, a.b, a.σ
  trace(() -> σ.(W*x .+ b))
end

# conv2d

@primitive Trace function (c::Conv)(x)
  cdims = DenseConvDims(val(x), c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  out = conv(val(x), c.weight, cdims)

  all(x-> x == c.pad[1], c.pad) ||
    error("Assymetric padding is unsupported by tf.conv2d")

  y = Staged(conv, stagedinputs(x)..., c.weight, c.stride, c.pad[1], c.dilation, v=out)
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  out = overdub(Trace(), (x) -> (σ).(x .+ b), y)
  wrap(out, vcall(vertex(DataFlow.Lambda(1, unwrap(out))), x))
end

Base.permutedims(x::Union{Staged,IVertex}, p) = jscall(:(math.transpose), x, p)
Base.reverse(x::Staged) = jscall(:(math.transpose), x)

# tf-js uses NHWC while js default is NCHW
function jscall(::typeof(conv), x, w, s, p, d)
  _x = permutedims(x, to_NHWC)
  _w = jscall(:(math.reverse), jscall(:(math.transpose), w, :([2, 3, 1,0]), x), :([0,1]))
  _s = reverse(s)
  _d = reverse(d)
  _out = jscall(:(math.conv2d), _x, _w, _s, p, "NHWC", _d, "floor")
  permutedims(_out, to_NCHW)
end

# maxpool

@primitive Trace function (m::MaxPool)(x)
  pdims = PoolDims(val(x), m.k; padding=m.pad, stride=m.stride)
  out = maxpool(val(x), pdims)

  all(x-> x == m.pad[1], m.pad) ||
    throw(error("Assymetric padding is unsupported by deeplearn-js"))

  Staged(maxpool, x, m.k, m.pad[1], m.stride, v=out)
end

function jscall(::typeof(maxpool), x, k, pad, stride)
  _x = permutedims(x, to_NHWC)
  _k = reverse(k)
  _s = reverse(stride)
  _out = jscall(:(math.maxPool), _x, _k, _s, pad, "floor")
  permutedims(_out, to_NCHW)
end

# broadcasted ops
bcastable(+, *, /, ^, tanh, σ, relu, leakyrelu, abs, exp, log, -, copy)
# copy for residual blocks

jscall(::typeof(broadcast), ::typeof(+), a, b) = jscall(:(math.add), a, b)
jscall(::typeof(broadcast), ::typeof(*), a, b) = jscall(:(math.mul), a, b)
jscall(::typeof(broadcast), ::typeof(/), a, b) = jscall(:(math.div), a, b)
jscall(::typeof(broadcast), ::typeof(^), a, b) = jscall(:(math.pow), a, b)
jscall(::typeof(broadcast), ::typeof(tanh), x) = jscall(:(math.tanh), x)
jscall(::typeof(broadcast), ::typeof(σ), x) = jscall(:(math.sigmoid), x)
jscall(::typeof(broadcast), ::typeof(relu), x) = jscall(:(math.relu), x)
jscall(::typeof(broadcast), ::typeof(leakyrelu), x) = jscall(:(math.leakyRelu), x, 0.01)
jscall(::typeof(broadcast), ::typeof(leakyrelu), x, a) = jscall(:(math.leakyRelu), x, a)
jscall(::typeof(broadcast), ::typeof(abs), x) = jscall(:(math.abs), x)
jscall(::typeof(broadcast), ::typeof(exp), x) = jscall(:(math.exp), x)
jscall(::typeof(broadcast), ::typeof(log), x) = jscall(:(math.log), x)
jscall(::typeof(broadcast), ::typeof(-), x, y) = jscall(:(math.sub), x, y)
jscall(::typeof(broadcast), ::typeof(copy), x) = jscall(:(flux.slice), x)

# reshape

Base.reshape(A::AbstractArray, i...) = nothing

@primitive Trace Base.reshape(parent::AbstractArray, dims::Union{Colon,StagedReal}...) =
  ! any(x -> x isa StagedType, (parent, dims...)) ?
  trace(reshape, parent, dims...) :
  begin
    dims = any(x -> val(x) isa Colon, dims) ? begin
      # resolve dimensions
      p = 1
      pos = 0
      for i=1:length(dims)
        !( val(dims[i]) isa Colon ) ?
         p = tracecall((p, v)-> p*v, p, dims[i]) : (pos = i)
      end
      c = tracecall((x, p) -> (count(x)/p), parent, p)
      (dims[1:pos-1]..., c, dims[pos+1:end]...)
    end : dims
    Staged(reshape, parent, dims..., v=reshape(val(parent), Int.(val.(dims))...))
  end

Base.count(x::AbstractArray) = prod(size(x))
@primitive Trace Base.count(x::AbstractArray) = Staged(getindex, x, :(String("size")), v=count(val(x)))

jscall(::typeof(reshape), p, dims...) =
  jscall(:(math.reshape), p, jscall(tuple, reverse(dims)...))

# size
@primitive Trace Base.size(x::StagedArray) =
  Staged(getindex, x, :("shape"), v=size(val(x)))

@primitive ctx::Trace Base.size(x, i) =
  !any(x -> x isa StagedType, (x, i)) ?
  trace(size, x, i) :
  begin
    s, n = invertedindex(x, i)
    wrap(s[i], unwrap(s[n]))
  end

# # gate ( for LSTM and GRU )
@primitive ctx::Trace function Flux.gate(x::AbstractArray, h, n)
  out = Flux.gate(val(x), val(h), val(n))
  _start =  trace((h, n) -> h * (n-1), h, n)
  Staged(view, x, _start, h, v=out)
end

jscall(::typeof(view), x, start, length) =
  jscall(:(math.slice), x, start, length)

add(x, y) = x + y
sub(x, y) = x - y
mul(x, y) = x * y
div(x, y) = x / y

# for Staged{Int}
function binary_op(op, sub)
  @eval @primitive Trace ($op)(x::Number, y::Number) =
    any(e-> e isa StagedReal, (x, y)) ? Staged($sub, x, y) : $(op)(x, y)
end

binary_op(+, add)
binary_op(*, mul)
binary_op(-, sub)
binary_op(/, div)

jscall(::typeof(add), x, y) = jscall(:(+), x, y)
jscall(::typeof(sub), x, y) = jscall(:(-), x, y)
jscall(::typeof(mul), x, y) = jscall(:(*), x, y)
jscall(::typeof(div), x, y) = jscall(:(/), x, y)

bn(args...) = nothing

@primitive Trace function (BN::BatchNorm)(x)
  μ, σ², ϵ, γ, β, λ = BN.μ, BN.σ², BN.ϵ, BN.γ, BN.β, BN.λ
  s = Staged(bn, stagedinputs(x)..., μ, σ², ϵ, γ, β, v=BN(val(x)))
  y = trace((l) -> λ.(l), s)
end

jscall(::typeof(bn), x, μ, σ², ϵ, γ, β) =
  jscall(:(tf.batchNormalization), x, μ, σ², ϵ, γ, β)

@primitive Trace copy(A::Staged{T}) where {T <: AbstractArray} =
  Staged(copy,A)

jscall(::typeof(copy), A) = jscall(:(flux.slice), A)

function invertedindex(x::AbstractArray, i)
  _size = trace((x) -> size(x), x)
  l = length(val(_size))
  index = trace((l, i)-> l - i + 1, l, i)
  return _size, index
end


text:
module EhrenfestMethods

using LinearAlgebra: lmul!, eigvecs, diag, dot

using NQCDynamics:
    NQCDynamics,
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    Calculators,
    DynamicsMethods,
    DynamicsUtils,
    Estimators

using NQCModels: NQCModels, Model
using NQCBase: Atoms

"""
Abstract type for Ehrenfest method.
"""
abstract type AbstractEhrenfest <: DynamicsMethods.Method end

function DynamicsMethods.motion!(du, u, sim::Simulation{<:AbstractEhrenfest}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    dσ = DynamicsUtils.get_quantum_subsystem(du)

    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)
    σ = DynamicsUtils.get_quantum_subsystem(u)

    DynamicsUtils.velocity!(dr, v, r, sim, t)
    Calculators.update_electronics!(sim.calculator, r)
    acceleration!(dv, v, r, sim, t, σ)
    DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim)
end

function DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim::AbstractSimulation{<:AbstractEhrenfest})
    V = DynamicsUtils.calculate_density_matrix_propagator!(sim, v)
    DynamicsUtils.commutator!(dσ, V, σ, sim.method.tmp_complex_matrix)
    lmul!(-im, dσ)
end

include("ehrenfest.jl")
include("ehrenfest_rpmd.jl")

end # module


text:
# synchronization

export synchronize_threads

"""
     synchronize_threads(ctx::AbstractKernelContext)

in CUDA terms `__synchronize`
in OpenCL terms: `barrier(CLK_LOCAL_MEM_FENCE)`
"""
function synchronize_threads(ctx::AbstractKernelContext)
    error("Not implemented") # COV_EXCL_LINE
end


text:
@testset "Debug Options" begin
    # helper to get debug as string
    io = IOBuffer()
    M = ManifoldsBase.DefaultManifold(2)
    x = [4.0, 2.0]
    o = GradientDescentOptions(x, StopAfterIteration(20), ConstantStepsize(1.0))
    f = y -> distance(M, y, x) .^ 2
    ∇f = y -> -2 * log(M, y, x)
    p = GradientProblem(M, f, ∇f)
    a1 = DebugDivider("|", io)
    @test Manopt.dispatch_options_decorator(DebugOptions(o, a1)) === Val{true}()
    # constructors
    @test DebugOptions(o, a1).debugDictionary[:All] == a1
    @test DebugOptions(o, [a1]).debugDictionary[:All].group[1] == a1
    @test DebugOptions(o, Dict(:A => a1)).debugDictionary[:A] == a1
    @test DebugOptions(o, ["|"]).debugDictionary[:All].group[1].divider == a1.divider
    # single AbstractOptionsActions
    # DebugDivider
    a1(p, o, 0)
    s = @test String(take!(io)) == "|"
    #DebugGroup
    DebugGroup([a1, a1])(p, o, 0)
    @test String(take!(io)) == "||"
    # Debug Every
    DebugEvery(a1, 10, false)(p, o, 9)
    @test String(take!(io)) == ""
    DebugEvery(a1, 10, true)(p, o, 10)
    @test String(take!(io)) == "|"
    @test DebugEvery(a1, 10, true)(p, o, -1) == nothing
    # Debug Cost
    @test DebugCost("A").prefix == "A"
    DebugCost(false, io)(p, o, 0)
    @test String(take!(io)) == "F(x): 0.0"
    DebugCost(false, io)(p, o, -1)
    @test String(take!(io)) == ""
    # entry
    DebugEntry(:x, "x:", io)(p, o, 0)
    @test String(take!(io)) == "x: $x"
    DebugEntry(:x, "x:", io)(p, o, -1)
    @test String(take!(io)) == ""
    # Change
    a2 = DebugChange(StoreOptionsAction((:x,)), "Last: ", io)
    a2(p, o, 0) # init
    o.x = [3.0, 2.0]
    a2(p, o, 1)
    @test String(take!(io)) == "Last: 1.0"
    # Iterate
    DebugIterate(io)(p, o, 0)
    @test String(take!(io)) == "x:$(o.x)"
    DebugIterate(io)(p, o, 1)
    @test String(take!(io)) == "x:$(o.x)"
    # Iteration
    DebugIteration(io)(p, o, 0)
    @test String(take!(io)) == "Initial"
    DebugIteration(io)(p, o, 23)
    @test String(take!(io)) == "# 23"
    # DEbugEntryChange - reset
    o.x = x
    a3 = DebugEntryChange(
        :x, (p, o, x, y) -> distance(p.M, x, y), StoreOptionsAction((:x,)), "Last: ", io
    )
    a4 = DebugEntryChange(
        x, :x, (p, o, x, y) -> distance(p.M, x, y), StoreOptionsAction((:x,)), "Last: ", io
    )
    a3(p, o, 0) # init
    @test String(take!(io)) == ""
    a4(p, o, 0) # init
    @test String(take!(io)) == ""
    #change
    o.x = [3.0, 2.0]
    a3(p, o, 1)
    @test String(take!(io)) == "Last: 1.0"
    a4(p, o, 1)
    @test String(take!(io)) == "Last: 1.0"
    # StoppingCriterion
    DebugStoppingCriterion(io)(p, o, 1)
    @test String(take!(io)) == ""
    o.stop(p, o, 20)
    DebugStoppingCriterion(io)(p, o, 20)
    @test String(take!(io)) == ""
    o.stop(p, o, 21)
    DebugStoppingCriterion(io)(p, o, 21)
    @test String(take!(io)) ==
          "The algorithm reached its maximal number of iterations (20).\n"
    #DebugFactory
    df = DebugFactory([:Stop, "|"])
    @test isa(df[:Stop], DebugStoppingCriterion)
    @test isa(df[:All], DebugGroup)
    @test isa(df[:All].group[1], DebugDivider)
    @test length(df[:All].group) == 1
    df = DebugFactory([:Stop, "|", 20])
    @test isa(df[:All], DebugEvery)
    @test all(isa.(
        DebugFactory([:Change, :Iteration, :Iterate, :Cost, :x])[:All].group,
        [DebugChange, DebugIteration, DebugIterate, DebugCost, DebugEntry],
    ))
    @test DebugActionFactory(a3) == a3
end


text:
#
# The type of the query output.
#

# Output cardinality.

immutable OutputMode
    optional::Bool
    plural::Bool
end

OutputMode() = OutputMode(false, false)

setoptional(omode::OutputMode, optional::Bool=true) =
    OutputMode(optional, omode.plural)
setplural(omode::OutputMode, plural::Bool=true) =
    OutputMode(omode.optional, plural)

setoptional(optional::Bool=true) =
    omode -> setoptional(omode, optional)

setplural(plural::Bool=true) =
    omode -> setplural(omode, plural)

isplain(omode::OutputMode) = !omode.optional && !omode.plural
isoptional(omode::OutputMode) = omode.optional
isplural(omode::OutputMode) = omode.plural

# The output type.

immutable Output
    dom::Domain
    mode::OutputMode
end

Output(dom) = Output(convert(Domain, dom), OutputMode())

decorate(otype::Output, d::Decoration) =
    Output(decorate(otype.dom, d), otype.mode)

setoptional(otype::Output, optional::Bool=true) =
    Output(otype.dom, setoptional(otype.mode, optional))

setplural(otype::Output, plural::Bool=true) =
    Output(otype.dom, setplural(otype.mode, plural))

convert(::Type{Output}, dom::Union{Type, Symbol, Tuple, Vector, Domain}) =
    Output(convert(Domain, dom))

function show(io::IO, otype::Output)
    print(io, otype.dom)
    if isplural(otype) && isoptional(otype)
        print(io, "*")
    elseif isplural(otype)
        print(io, "+")
    elseif isoptional(otype)
        print(io, "?")
    end
end

# Predicates and properties.

domain(otype::Output) = otype.dom
mode(otype::Output) = otype.mode

isdata(otype::Output) = isdata(otype.dom)
isany(otype::Output) = isany(otype.dom)
isvoid(otype::Output) = isvoid(otype.dom)
isnone(otype::Output) = isnone(otype.dom)
isentity(otype::Output) = isentity(otype.dom)
isrecord(otype::Output) = isrecord(otype.dom)

isplain(otype::Output) = isplain(otype.mode)
isoptional(otype::Output) = isoptional(otype.mode)
isplural(otype::Output) = isplural(otype.mode)

decorations(otype::Output) = decorations(otype.dom)
decoration(otype::Output, name, T, default) =
    decoration(otype.dom, name, T, default)

classname(otype::Output) = classname(otype.dom)
fields(otype::Output) = fields(otype.dom)

# The native Julia type that can represent the output.

datatype(otype::Output) =
    let T = datatype(otype.dom)
        isplural(otype) ?
            Vector{T} :
        isoptional(otype) ?
            Nullable{T} :
            T
    end

# Domain definitions that depend on Output.

convert(::Type{Domain}, desc::Union{Tuple, Vector}) =
    Domain(Output[convert(Output, otype) for otype in desc])

isdata(::Vector{Output}) = false
isany(::Vector{Output}) = false
isvoid(::Vector{Output}) = false
isnone(::Vector{Output}) = false
isentity(::Vector{Output}) = false
isrecord(::Vector{Output}) = true

datatype(desc::Vector{Output}) =
    Tuple{(datatype(otype) for otype in desc)...}

classname(::Vector{Output}) = NO_CLASSNAME

const NO_FIELDS = Output[]

fields(desc::Type) = NO_FIELDS
fields(desc::Symbol) = NO_FIELDS
fields(desc::Vector{Output}) = desc

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using SemanticModels
using SemanticModels.ModelTools

expr =quote
    function f(a::Int64, b::Int64)
    c = a + b
    d = g(float(c))

    return d
end
function g(res::Float64)
    res = abs(res)
    res -= 1

    while res > 0
        res = g(res)
    end

    return res
end
@show f(1,2)
end

# +
funcs = Expr[]
for ex in expr.args
    if isa(ex, LineNumberNode)
        continue
    end
    if ex.head == :function
        push!(funcs, ex)
    end
end

funcs

# +
nametype(ex::Expr) = begin
    ex.head == :(::) || error("$ex is not a type assertion")
    avar = ex.args[1]
    atyp = ex.args[2]
    return avar, atyp
end

nametype(s::Symbol) = (s, :Any)

function describeargs(fu::Expr)
    f = fu
    @show argl = f.args[1].args[2:end]
    @show argl
    @show body = f.args[2].args
    fname = string(f.args[1])
    for a in reverse(argl)
        avar, atyp = nametype(a)
        varname = string(avar)
        pushfirst!(body,
            :(println("F: ", $fname,";",
                $varname,"=", $avar,"::",typeof($avar), "<:", $atyp))
        )
    end
    return f
end

for f in funcs
    describeargs(f)
    body = f.args[2].args
    fname = string(f.args[1])
    for ex in body
        if isa(ex,LineNumberNode)
            continue
        end
        if ex.head == :(=)
            @show a = ex.args[1]
            varname = string(a)
            insert!(body, length(body)-1,:(println("A: ", $fname,";",$varname, "=", $a)))
        end
    end
end
expr
# -

eval(expr)


text:
module TestOperators
    using Base.Test
    using DataArrays
    using DataFrames
    using Stats

    macro test_da_pda(da, code)
        esc(quote
            let $da = copy($da)
                $code
            end
            let $da = PooledDataArray($da)
                $code
            end
        end)
    end

    unary_operators = [(+), (-), (!)]

    numeric_unary_operators = [(+), (-)]

    logical_unary_operators = [(!)]

    elementary_functions = [abs, sign, acos, acosh, asin,
                            asinh, atan, atanh, sin, sinh,
                            cos, cosh, tan, tanh, ceil, floor,
                            round, trunc, exp, exp2, expm1, log, log10, log1p,
                            log2, exponent, sqrt, gamma, lgamma, digamma,
                            erf, erfc]

    special_comparison_operators = [isless]

    comparison_operators = [(==), (.==), (!=), (.!=),
                            (>), (.>), (>=), (.>=), (<), (.<),
                            (<=), (.<=)]

    scalar_comparison_operators = [(==), (!=), (>), (>=),
                                   (<), (<=)]

    array_comparison_operators = [(.==), (.!=), (.>), (.>=), (.<), (.<=)]

    binary_operators = [(+), (.+), (-), (.-), (*), (.*), (/), (./),
                        (.^), (div), (mod), (fld), (rem),
                        (&), (|), ($)]

    induced_binary_operators = [(^)]

    arithmetic_operators = [(+), (.+), (-), (.-), (*), (.*), (/), (./),
                            (.^), (div), (mod), (fld), (rem)]

    induced_arithmetic_operators = [(^)]

    biscalar_operators = [(max), (min)]

    scalar_arithmetic_operators = [(+), (-), (*), (/),
                                   (div), (mod), (fld), (rem)]

    induced_scalar_arithmetic_operators = [(^)]

    array_arithmetic_operators = [(+), (.+), (-), (.-), (.*), (./), (.^)]

    bit_operators = [(&), (|), ($)]

    unary_vector_operators = [minimum, maximum, prod, sum, mean, median, std, var, norm]

    ex_tunary_vector_operators = [mad, skewness, kurtosis]

    pairwise_vector_operators = [diff]#, percent_change]

    cumulative_vector_operators = [cumprod, cumsum, cumsum_kbn, cummin, cummax]

    ffts = [fft]

    binary_vector_operators = [dot, cor, cov, corspearman]

    columnar_operators = [colmins, colmaxs, colprods, colsums,
                          colmeans, colmedians, colstds, colvars,
                          colnorms]

    boolean_operators = [any, all]

    # Unary operators on DataFrame's should be equivalent to elementwise
    # application of those same operators
    df = DataFrame(quote
                       A = [1, 2, 3, 4]
                       B = [1.0, pi, pi, e]
                   end)
    for f in numeric_unary_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end
    df = DataFrame(quote
                       A = [true, false, true, false]
                   end)
    for f in logical_unary_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end

    # Elementary functions on DataFrames's
    N = 5
    df = DataFrame(A = DataArray(ones(N)),
                   B = DataArray(ones(N)))
    for f in elementary_functions
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                  @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end

    # Broadcasting operations between NA's and DataFrames's
    N = 5
    df = DataFrame(A = DataArray(ones(N)),
                   B = DataArray(ones(N)))
    for f in arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert isna(f(df, NA)[i, j])
                @assert isna(f(NA, df)[i, j])
            end
        end
    end

    # Broadcasting operations between scalars and DataFrames's
    N = 5
    df = DataFrame(A = DataArray(ones(N)),
                   B = DataArray(ones(N)))
    for f in arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df, 1)[i, j] == f(df[i, j], 1)
                @assert f(1, df)[i, j] == f(1, df[i, j])
            end
        end
    end

    # Binary operations on pairs of DataFrame's
    # TODO: Test in the presence of in-operable types like Strings
    N = 5
    df = DataFrame(A = DataArray(ones(N)),
                   B = DataArray(ones(N)))
    for f in array_arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert isna(f(df, df)[i, j]) && isna(df[i, j]) ||
                        f(df, df)[i, j] == f(df[i, j], df[i, j])
            end
        end
    end

    # TODO: Columnar operators on DataFrame's

    # Boolean operators on DataFrames's
    N = 5
    df = DataFrame(A = DataArray(falses(N)))
    @assert any(df) == false
    @assert any(!df) == true
    @assert all(df) == false
    @assert all(!df) == true

    df = DataFrame(A = DataArray(falses(N)))
    df[3, 1] = true
    @assert any(df) == true
    @assert all(df) == false

    df = DataFrame(A = DataArray(falses(N)))
    df[2, 1] = NA
    df[3, 1] = true
    @assert any(df) == true
    @assert all(df) == false

    dv = @data [false, NA]
    dv[1] = NA

    df = DataFrame(A = DataArray(falses(N)))
    df[2, 1] = NA
    @assert isna(any(df))
    @assert all(df) == false

    df = DataFrame(A = DataArray(falses(N)))
    df[1, 1] = NA
    @assert isna(any(dv))
    @assert isna(all(dv))

    # Is this a genuine special case?
    @assert isna(NA ^ 2.0)

    #
    # Equality tests
    #
    df = DataFrame({dv})
    alt_dv = @data [false, NA]
    alt_dv[1] = NA
    alt_df = DataFrame({alt_dv})
    # @assert isequal(DataFrame({dv}) .== DataFrame({dv}), DataFrame({DataVector[true, NA]}))

    @assert all(isna(NA .== df))
    @assert all(isna(df .== NA))
end


text:
"""
    AbstractLink

Abstract type defining maps from R^n -> X.
They can be applied by calling `link(x)`.

A series of definitions are given in http://web.pdx.edu/~newsomj/mvclass/ho_link.pdf
"""
abstract type AbstractLink end

"""
	ChainLink(links)

Create a composed chain of different links.
"""
struct ChainLink{Tls} <: AbstractLink
    links::Tls
end

(l::ChainLink)(x) = foldl((x, l) -> l(x), l.ls; init=x)

"""
    Link(f)

General construction for a link with a function `f`.
"""
struct Link{F} <: AbstractLink
    f::F
end

link(f) = Link(f)
link(l::AbstractLink) = l

(l::Link)(x) = l.f(x)

Base.inv(l::Link) = Link(InverseFunctions.inverse(l.f))

"""
    BijectiveSimplexLink(link)

Wrapper to preprocess the inputs by adding a `0` at the end before passing it to
the link `link`.
This is a necessary step to work with simplices.
For example with the [`SoftMaxLink`](@ref), to obtain a `n`-simplex leading to
`n+1` categories for the [`CategoricalLikelihood`](@ref),
one needs to pass `n+1` latent GP.
However, by wrapping the link into a `BijectiveSimplexLink`, only `n` latent are needed.
"""
struct BijectiveSimplexLink{L} <: AbstractLink
    link::L
end

(l::BijectiveSimplexLink)(f::AbstractVector{<:Real}) = l.link(vcat(f, 0))

# alias
const LogLink = Link{typeof(log)}
const ExpLink = Link{typeof(exp)}

const InvLink = Link{typeof(inv)}

const SqrtLink = Link{typeof(sqrt)}
const SquareLink = Link{typeof(InverseFunctions.square)}

const LogitLink = Link{typeof(logit)}
const LogisticLink = Link{typeof(logistic)}

const ProbitLink = Link{typeof(norminvcdf)}
const NormalCDFLink = Link{typeof(normcdf)}

const SoftMaxLink = Link{typeof(softmax)}

"""
    LogLink()

`log` link, f:ℝ⁺->ℝ . Its inverse is the [`ExpLink`](@ref).
"""
LogLink() = Link(log)

"""
    ExpLink()

`exp` link, f:ℝ->ℝ⁺. Its inverse is the [`LogLink`](@ref).
"""
ExpLink() = Link(exp)

"""
    InvLink()

`inv` link, f:ℝ/{0}->ℝ/{0}. It is its own inverse.
"""
InvLink() = Link(inv)

"""
    SqrtLink()

`sqrt` link, f:ℝ⁺∪{0}->ℝ⁺∪{0}. Its inverse is the [`SquareLink`](@ref).
"""
SqrtLink() = Link(sqrt)

"""
    SquareLink()

`^2` link, f:ℝ->ℝ⁺∪{0}. Its inverse is the [`SqrtLink`](@ref).
"""
SquareLink() = Link(InverseFunctions.square)

"""
    LogitLink()

`log(x/(1-x))` link, f:[0,1]->ℝ. Its inverse is the [`LogisticLink`](@ref).
"""
LogitLink() = Link(logit)

"""
    LogisticLink()

`1/(1+exp(-x))` link. f:ℝ->[0,1]. Its inverse is the [`LogitLink`](@ref).
"""
LogisticLink() = Link(logistic)

"""
    ProbitLink()

`ϕ⁻¹(y)` link, where `ϕ⁻¹` is the `invcdf` of a `Normal` distribution, f:[0,1]->ℝ.
Its inverse is the [`NormalCDFLink`](@ref).
"""
ProbitLink() = Link(norminvcdf)

"""
    NormalCDFLink()

`ϕ(y)` link, where `ϕ` is the `cdf` of a `Normal` distribution, f:ℝ->[0,1].
Its inverse is the [`ProbitLink`](@ref).
"""
NormalCDFLink() = Link(normcdf)

(::ChainLink{<:Tuple{LogLink,NormalCDFLink}})(x) = normlogcdf(x) # Specialisation for log + normal cdf

"""
    SoftMaxLink()

`softmax` link, i.e `f(x)ᵢ = exp(xᵢ)/∑ⱼexp(xⱼ)`.
f:ℝⁿ->Sⁿ⁻¹, where Sⁿ⁻¹ is an [(n-1)-simplex](https://en.wikipedia.org/wiki/Simplex)
It has no defined inverse
"""
SoftMaxLink() = Link(softmax)


text:
##### Packaged loglikelihood
using LinearAlgebra
#using Turing
using Random
using AdvancedHMC
using ForwardDiff
include("../../src/faberschauder.jl")
include("../../src/asvar.jl")
include("../../src/pdmp.jl")
# up to constant of proportionality
# π(x) = exp(-ϕ)*C
function ϕ(ξ, L, T; p = T/(1<<L)) # formula (17)
        r = 0.0
        s = p*0.5
        for i in 1:(1<<L)
            x = dotψ(ξ, s, L,  T)
            r += 0.5*p*(b(x)^2 + b′(x))
            s += p
        end
        @assert  T+p*0.5 - 0.0001 < s < T+p*0.5 + 0.0001
    return r + 0.5*norm(ξ)^2
end

function ∇ϕ(ξ, i, L, T; p = T/(1<<L)) # formula (17)
    if i == (2 << L) + 1    # final point
        error("fixed final point")
    elseif i == 1   # initial point
        error("fixed initial point")
    else
        l = lvl(i, L)
        k = (i - 1) ÷ (2 << l)
        # println("index ($l, $k) out of $L level")
        δ = T/(1 << (L-l)) # T/(2^(L-l))
        r = 0.0
        s = δ*k + p*0.5
        # println("T = $T")
        for i in 1:l+1
            # println("integrating at the point $s")
            x = dotψ(ξ, s, L,  T)
            r += 0.5*p*Λ(s, L-l, T)*(2b(x)*b′(x) + b″(x))
            s += p
        end
        @assert δ*(k + 1) + p*0.5 - 0.001 < s < δ*(k + 1) + p*0.5 + 0.001
        return r + ξ[i]
    end
end

# Drift
b(x) = α * sin(x)
# First derivative
b′(x) = α * cos(x)
# Second derivative
b″(x) = -α * sin(x)

# Drift
#α = 1.5
L = 6
T = 50.0
n = (2 << L) + 1
u, v = 0.0, 0.0  # initial and fianl point
D = n - 2

# loglikelihood
function logdensity_f(q)
    ξ = [u/sqrt(T); q;  v/sqrt(T)]
    return -ϕ(ξ, L, T)
end
# gradient of loglikelihood
function grad_f(q)
    ξ = [u / sqrt(T); q; v / sqrt(T)]
    return (-ϕ(ξ, L, T),[-∇ϕ(ξ, i, L, T) for i in 2:(2 << L)])
end

# reshape output from vector of vector to matrix
function vecvec_to_matrix(vecvec)
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(Float64, dim2, dim1)
    for i in 1:dim1
        for j in 1:dim2
            my_array[j,i] = vecvec[i][j]
        end
    end
    return my_array
end

function run_hmc(df)
    q0 = randn(D)
    n_samples, n_adapts, target = 3_000, 2_000, 0.8
    metric = DiagEuclideanMetric(D)
    # h = Hamiltonian(metric, logdensity_f, grad_f)
    h = Hamiltonian(metric, logdensity_f, ForwardDiff)
    eps_init = find_good_stepsize(h, q0)
    int = Leapfrog(eps_init)
    traj = AdvancedHMC.NUTS{MultinomialTS,GeneralisedNoUTurn}(int)
    adaptor = AdvancedHMC.StanHMCAdaptor(
        n_adapts, MassMatrixAdaptor(metric), NesterovDualAveraging(target, eps_init)
        )
    #u, v = 0.0/sqrt(T), 0.0/ sqrt(T)  # initial and fianl point
    hmc_time = @elapsed ((samples, stats) = sample(h, traj, q0, n_samples, adaptor, n_adapts))
    burn1 = 100
     # xT_2 = sqrt(T)/2*(samples[Int((length(q0)+1)/2), burn1+1:niter])
    xT_2 = sqrt(T)/2 .*(getindex.(samples, Int((length(q0)+1)/2)))
    nn = length(xT_2)
    ave = cumsum(xT_2)./(Float64.(1:nn))
    for batches in [50]
        #reshape a vector of vectors in a matrix
        ess = ESS(vecvec_to_matrix(samples)[:,burn1:end], n_batches = batches)
        push!(df, Dict(:sampler => "HMC", :alpha => α, :T => n_samples, :nbatches => batches,
            :stat => "ess_mean", :y => sum(ess[2:end-1])/(length(ess)-2), :runtime =>  hmc_time))
        push!(df, Dict(:sampler => "HMC", :alpha => α, :T => n_samples, :nbatches => batches,
            :stat => "ess_median", :y => median(ess[2:end-1]), :runtime =>  hmc_time))
        push!(df, Dict(:sampler => "HMC", :alpha => α, :T => n_samples, :nbatches => batches,
            :stat => "ess_min", :y => minimum(ess[2:end-1]), :runtime => hmc_time))
        push!(df, Dict(:sampler => "HMC", :alpha => α, :T => n_samples, :nbatches => batches,
            :stat => "ess_x_T2", :y => ess[Int((length(q0)+1)/2)], :runtime => hmc_time))
    end
    return df, ave
end

using DataFrames, Random
function data_collection_hmc(df)
    Random.seed!(0)
    for α′ in [0.1, 0.3, 0.5]
        global α = α′
        println("iteration ")
        df, x = run_hmc(df)
    end
    return df
end

df = DataFrame(sampler = String[], alpha = Float64[], T = Float64[],
    nbatches = Int64[], stat = String[], y = Float64[], runtime = Float64[])
df = data_collection_hmc(df)
using CSV, JLD2
CSV.write("./scripts/compare/benchamrk_hmc.csv", df)

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using Distributions
using Test
using DataFrames
using CSVFiles
using DataStructures

################################################################################
## JULIA
################################################################################
include("../src/utils.jl")
include("../src/sample_sobol.jl")
include("../src/analyze_sobol.jl")
include("../src/test_functions/ishigami.jl")

# define the (uncertain) parameters of the problem and their distributions
data = SobolData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = 1000
)

N = data.N
D = length(data.params)

# sampling
julia_samples = sample(data) |> DataFrame
julia_ishigami = ishigami(convert(Matrix, julia_samples)) |> DataFrame

# analysis
julia_A, julia_B, julia_AB, julia_BA = split_output(convert(Matrix, julia_ishigami), N, D, data.calc_second_order)
julia_results = analyze(data, convert(Matrix, julia_ishigami); num_resamples = 1_000)

################################################################################
## Python
################################################################################

# sampling
py_samples = load("data/py_uniform/py_samples.csv", header_exists=false, colnames = ["x1", "x2", "x3"]) |> DataFrame
py_ishigami = load("data/py_uniform/py_ishigami.csv", header_exists=false) |> DataFrame

# analysis
py_A = load("data/py_uniform/py_A.csv", header_exists=false) |> DataFrame
py_B = load("data/py_uniform/py_B.csv", header_exists=false) |> DataFrame
py_AB = load("data/py_uniform/py_AB.csv", header_exists=false) |> DataFrame
py_BA = load("data/py_uniform/py_BA.csv", header_exists=false) |> DataFrame

py_firstorder = load("data/py_uniform/py_firstorder.csv", header_exists=false) |> DataFrame
py_secondorder = load("data/py_uniform/py_secondorder.csv", header_exists=false) |> DataFrame
py_totalorder = load("data/py_uniform/py_totalorder.csv", header_exists=false) |> DataFrame

py_firstorder_conf = load("data/py_uniform/py_firstorder_conf.csv", header_exists=false) |> DataFrame
py_secondorder_conf = load("data/py_uniform/py_secondorder_conf.csv", header_exists=false) |> DataFrame
py_totalorder_conf = load("data/py_uniform/py_totalorder_conf.csv", header_exists=false) |> DataFrame

################################################################################
## Testing
################################################################################

@testset "Uniform Sampling" begin
    @test convert(Matrix, julia_samples) ≈ convert(Matrix, py_samples) atol = 1e-9
end

@testset "Uniform Analysis" begin
    @test convert(Matrix, julia_ishigami) ≈ convert(Matrix, py_ishigami) atol = 1e-9
    @test julia_A ≈ convert(Matrix, py_A) atol = 1e-9
    @test julia_B ≈ convert(Matrix, py_B) atol = 1e-9
    @test julia_AB ≈ convert(Matrix, py_AB) atol = 1e-9
    @test julia_BA ≈ convert(Matrix, py_BA) atol = 1e-9
    @test julia_results[:firstorder] ≈ convert(Matrix, py_firstorder) atol = 1e-9
    @test julia_results[:totalorder] ≈ convert(Matrix, py_totalorder) atol = 1e-9

    for i = 1:D
        for j = i+1:D
            @test julia_results[:secondorder][i,j] ≈ convert(Matrix, py_secondorder)[i,j] atol = 1e-9
        end
    end

    # TODO Choose proper tolerance for CI comparison
    @test julia_results[:firstorder_conf] ≈ convert(Matrix, py_firstorder_conf) atol = 1e-1
    @test julia_results[:totalorder_conf] ≈ convert(Matrix, py_totalorder_conf) atol = 1e-1

    for i = 1:D
        for j = i+1:D
            @test julia_results[:secondorder_conf][i,j] ≈ convert(Matrix, py_secondorder_conf)[i,j] atol = 1e-1
        end
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export tab

"""
    tab(;kwargs...)
    tab(children::Any;kwargs...)
    tab(children_maker::Function;kwargs...)

A Tab component.
Utility component to pass to Tabs. For more information, see: https://mantine.dev/core/tabs/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Tab content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `disabled` (Bool; optional): A tab can show it is currently unable to be interacted with
- `label` (String; optional): Tab control label
"""
function tab(; kwargs...)
        available_props = Symbol[:children, :id, :className, :disabled, :label]
        wild_props = Symbol[]
        return Component("tab", "Tab", "dash_mantine_components", available_props, wild_props; kwargs...)
end

tab(children::Any; kwargs...) = tab(;kwargs..., children = children)
tab(children_maker::Function; kwargs...) = tab(children_maker(); kwargs...)


text:
# import LinearAlgebra

using LinearAlgebra: BlasInt, Cholesky

## Custom stuff

"""
    λ, C = posdefλ(Q0, Q1)

Provide an initial `λ` such that `Q0 + λ*Q1` is positive-semidefinite (or preferably,
positive-definite).  `λ` may be larger than strictly necessary. `C` is the Cholesky
decomposition of `Q0 + λ*Q1`.
"""
function posdefλ(Q0::SymTridiagonal{T}, Q1::Union{Diagonal{T},UniformScaling};
                 growauto=T(1.0001), growfail=T(1.1)) where T
    λ = zero(T)
    C = cholesky(Q0; check=false)
    if issuccess(C)
        return λ, C
    end
    dv, ev = Q0.dv, Q0.ev
    n = length(dv)
    α = 4*(cos(π/(n+1)))^2
    for i = 1:n-1
        # A symmetric tridiagonal matrix M is positive-definite if
        #     α * M[i,i+1]^2 < M[i,i]*M[i+1,i+1]   ∀ i
        # See Johnson, Neumann, and Tsatsomeros, Conditions for the positivity of
        # determinants, Linear and Multilinear Algebra 40: 241-248 (1996).
        # This criterion is sufficient but not necessary.

        # In this case M[i,i+1] = Q0[i,i+1] and M[i,i] = Q0[i,i] + λ*Q1[i,i].
        # Consequently we can turn this into a quadratic inequality in λ
        #       a*λ^2 + b*λ + c > 0
        c = dv[i]*dv[i+1] - α*ev[i]^2
        if c < 0 || dv[i] < 0
            # This term needs some positive λ
            a = T(Q1[i,i]*Q1[i+1,i+1])
            @assert a >= 0   # Q1 must be possemidef
            b = Q1[i,i]*dv[i+1] + Q1[i+1,i+1]*dv[i]
            s = sqrt(b^2 - 4*a*c)
            # Since the quadratic is convex and below 0 at λ=0, only one root has λ > 0.
            # So we always choose the (-b + s)/(2*a) solution. But write it in a
            # numerically-stable way.
            λ0 = b > 0  ?  -2*c/(s + b)  :  (-b + s)/(2*a)
            λ = max(λ, λ0)
        end
    end
    # Due to roundoff we need to check whether this actually produces a posdef matrix
    λ *= T(growauto)
    C = cholesky(Q0 + λ*Q1; check=false)
    while !issuccess(C)
        λ *= T(growfail)
        C = cholesky(Q0 + λ*Q1; check=false)
    end
    return λ, C
end

function posdefλ(Q0::AbstractMatrix{T}, Q1::Union{AbstractMatrix{T},UniformScaling};
                 growauto=T(1.0001), growfail=T(1.1)) where T
    λ = zero(T)
    C = cholesky(Q0; check=false)
    if issuccess(C)
        return λ, C
    end
    if isa(Q1, UniformScaling)
        D, _ = eigs(Q0; nev=1, which=:SR, ritzvec=false)
        λ = -T(D[1]/T(Q1.λ))
    else
        D, _ = eigs(Q0, Q1; nev=1, which=:SR, ritzvec=false)
        λ = -T(D[1])
    end
    @assert λ > 0
    λ *= T(growauto)
    C = cholesky(Q0 + λ*Q1; check=false)
    while !issuccess(C)
        λ *= T(growfail)
        C = cholesky(Q0 + λ*Q1; check=false)
    end
    return λ, C
end

asmatrix(M::UniformScaling{T}, sz) where T = M.λ * eye(T, sz...)
asmatrix(M::AbstractMatrix, sz) = Matrix(M)

## Stuff that should be in stdlib/LinearAlgebra

function cholesky(A::SymTridiagonal{T}, ::Val{false}=Val(false); check::Bool=true) where T
    n = size(A, 1)
    dv = zeros(T, n)
    ev = zeros(T, n-1)
    L = Bidiagonal{T}(dv, ev, :L)
    if A.dv[1] < 0
        check && error("A is not positive-semidefinte")
        dv[1] = A.dv[1]
        info = -1
        return Cholesky(L, 'L', convert(BlasInt, info))
    end
    dv[1] = sqrt(A.dv[1])
    info = 0
    for i = 1:n-1
        ev[i] = A.ev[i]/dv[i]
        Δ = A.dv[i+1] - ev[i]^2
        if Δ < 0
            check && error("A is not positive-semidefinte")
            info = -1
            dv[i+1] = Δ  # return the magnitude of the violation for later inspection
            break
        end
        dv[i+1] = sqrt(Δ)
    end
    return Cholesky(L, 'L', convert(BlasInt, info))
end

# cholesky(A::AbstractMatrix{T}, ::Val{false}=Val(false); check::Bool=true) where T =
#     cholfact(A)
function cholesky(A::StridedMatrix{T}, ::Val{false}=Val(false); check::Bool=true) where T
    L, info = LinearAlgebra.LAPACK.potrf!('L', copy(A))
    check && @assert info == 0
    return Cholesky(L, 'L', info)
end

# issuccess(C::Cholesky) = C.info == 0
function issuccess(C::Cholesky)
    n = size(C, 1)
    F = C.factors
    for i = 1:n
        F[i,i] >= 0 || return false
    end
    return true
end

function Base.:+(S::SymTridiagonal{T}, D::UniformScaling{T}) where T
    dv, ev = copy(S.dv), copy(S.ev)
    n = length(dv)
    for i = 1:n
        dv[i] += D[i,i]
    end
    return SymTridiagonal(dv, ev)
end

function Base.:\(C::Cholesky{T,Bidiagonal{T}}, b::AbstractVector) where T
    if C.factors.isupper #C.factors.uplo == 'U'
        R = C.factors
        return R \ (R' \ b)
    end
    L = C.factors
    return L' \ (L \ b)
end

LinearAlgebra.eigen(A::AbstractMatrix, B::UniformScaling) = eigen(A, asmatrix(B, size(A)))
LinearAlgebra.eigen(A::AbstractMatrix, B::Diagonal) = eigen(A, asmatrix(B, size(A)))


text:
# This file was generated, do not modify it. # hide
ŷ = predict_mode(mtm, rows=test)
round(misclassification_rate(ŷ, y[test]), sigdigits=3)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
#=

    deterministic.jl

    Contains code to solve the deterministic ODE model

           β     λ    pβ
        T --→ S --→ E --→ S + F

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#

"""
    solve_ode_model(t,[λ,β,p])

Solve the ODE model using an analytical solution obtained in Mathematica.
Returns [A(t),I(t)].

"""
function solve_ode_model(t,θ::Vector)
    λ,β,p = θ
    A = (exp(-t*β)*(exp(-t*((-1+p)*β+λ))*(-1+p)*p*β^3*λ-(-1+p)*(β-λ)*λ*(p*β+λ)^2+exp(t*β)*((-1+p)*β+λ)*(λ^2*(β+λ)+p^2*β*(-λ^2+β^2*(1+t*λ)+β*λ*(1+t*λ))+p*λ*(-λ^2+β^2*(1+t*λ)+β*λ*(1+t*λ)))))/((β+λ)*((-1+p)*β+λ)*(p*β+λ)^2)
    S = (exp(-t*(p*β+λ))*β*((-1+p)*β*λ+exp(t*(p*β+λ))*p*(β+λ)*((-1+p)*β+λ)-exp(t*((-1+p)*β+λ))*(-1+p)*λ*(p*β+λ)))/((β+λ)*((-1+p)*β+λ)*(p*β+λ))
    I = A - S
    return [A,I]
end
solve_ode_model(t,λ,β,p) = solve_ode_model(t,[λ,β,p])

"""
    solve_ode_model_all_vars(t,[λ,β,p])

Solve the ODE model using matrix exponentiation and return all variables.
Returns [T(t),S(t),E(t),F(t)].

"""
function solve_ode_model_all_vars(t,θ::Vector)
    λ,β,p = θ
    M = [-β 0 0 0; β -λ p*β 0; 0 λ -p*β 0; 0 0 p*β 0]
    x₀ = [λ/(λ+β),β/(λ+β),0,0]
    return exp(M*t) * x₀
end
solve_ode_model_all_vars(t,λ,β,p) = solve_ode_model_all_vars(t,[λ,β,p])

"""
    solve_alt_ode_model(t,[γ,λ,β,p])

Solve the alternate ODE model where antibody binds to receptors at constant rate γ.
"""
function solve_alt_ode_model(t,θ::Vector)
    γ,λ,β,p = θ
    M = [-β 0 0 0 0; β -γ 0 p*β 0; 0 γ -λ 0 0; 0 0 λ -p*β 0; 0 0 0 p*β 0]
    x₀ = [λ/(λ+β),β/(λ+β),0,0,0]
    x = exp(M*t) * x₀
    A = sum(x[3:5])
    I = sum(x[4:5])
    return [A,I]
end
solve_alt_ode_model(t,γ,λ,β,p) = solve_alt_ode_model(t,[γ,λ,β,p])

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
#julia uses "using" instead of "import
using ArgParse
#FASTX contains programs to unpack and read fasta/fastq files
using FASTX

using BioSequences
using CodecZlib

#Julia cli setup
cli=ArgParseSettings()
@add_arg_table! cli begin
    "--file"
        help="path to sequence file"
        required=true
        #arg_type=int or float or whatever the type is
        #default=
    #next flag would go here
end

arg=parse_args(ARGS, cli)

function unpacker(file)
    if splitext(file)[2]==".gz"
        myfasta=FASTA.Reader(GzipDecompressorStream(open(file)))
   else myfasta=FASTA.Reader(open(file))
   end
   return myfasta
end

function dnastring(myfasta)

    mydict=Dict{DNA, String}(
        DNA_A => "A",
        DNA_T => "T",
        DNA_C => "C",
        DNA_G => "G",
    )
    all_seqs=String[]
    for pair in myfasta
        dna_seq=FASTA.sequence(pair)
        chars=String[]

        for i in 1:length(dna_seq)
            if haskey(mydict, dna_seq[i])
                nt = mydict[dna_seq[i]]
            else
                nt = "N"
            end
            push!(chars, nt)
        end
        string_seq=join(chars)
        push!(all_seqs, string_seq)
    end
    return all_seqs
end

#get fasta tupes
myfasta=unpacker(arg["file"])

#get string sequences
mydata=dnastring(myfasta)

#I do not know why this is necessary but it is :(
function tester(myfasta)
    contigs = 0
    contig_lengths = []
    nt_length = 0
    nt_counts = Dict()

    for seq in mydata
        len = length(seq)
        contigs +=1
        #appending to contig_lengths
        push!(contig_lengths, len)
        nt_length += len
        for i in 1:len
            n = seq[i]
            holder=get!(nt_counts, n, 0)
            holder+=1
            nt_counts[n] = holder
        end
    end
    println(contigs)
    println("--------")
    println(contig_lengths)
    println("--------")
    println(nt_length)
    println("--------")
    println(keys(nt_counts))
    println("--------")
    println(values(nt_counts))
end

tester(myfasta)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
"""
    UNetDynamic(backbone, inputsize, k_out[; kwargs...])
Create a U-Net model from convolutional `backbone` architecture. After every
downsampling layer (i.e. pooling or strided convolution), a skip connection and
an upsampling block are inserted, resulting in a convolutional network with
the same spatial output dimensions as its input. Outputs an array with `k_out`
channels.
## Keyword arguments
- `upsample`: A *constructor* for an upsampling block callable with `upsample(insize, k_out)`.
    If `insize` is `(h, w, k, b)`, then the output should have size `(2h, 2w, k_out, b)`.
    Defaults to [`FastAI.Models.upsample_block_small`](#).
- `agg`: Aggregation function for skip connection. Default concatenates in the
  channel dimension. Use `+` for summing and see [`Flux.SkipConnection`](#) for more
  details.
- `fdownscale = 0`: Number of upsampling steps to leave out. By default there will be one
    upsampling step for every downsampling step in `backbone`. Hence if the input spatial
    size is `(h, w)`, the output size will be `(h/2^fdownscale, w/2^fdownscale)`, i.e.
    to get outputs at half the resolution, set `fdownscale = 1`.
- `kwargs...`: Other keyword arguments are passed through to `upsample`.
## Examples
```julia
using Metalhead
backbone = Metalhead.ResNet50(pretrain=true).layers[1][1:end-1]
unet = Metalhead.UNetDynamic(backbone, (256, 256, 3, 1); k_out = 10)
Flux.outputsize(unet, (256, 256, 3, 1)) == (256, 256, 10, 1)
unet = Metalhead.UNetDynamic(backbone, (256, 256, 3, 1); fdownscalk_out = 10)
Flux.outputsize(unet, (256, 256, 3, 1)) == (256, 256, 10, 1)
```
"""
function UNetDynamic(
    backbone,
    inputsize,
    k_out::Int;
    final = UNetFinalBlock,
    fdownscale = 0,
    kwargs...,
)
    backbonelayers = collect(iterlayers(backbone))
    unet = unetlayers(
        backbonelayers,
        inputsize;
        m_middle = UNetMiddleBlock,
        skip_upscale = fdownscale,
        kwargs...,
    )
    outsz = Flux.outputsize(unet, inputsize)
    return Chain(unet, final(outsz[end-1], k_out))
end

function catchannels(x1, x2)
    ndims(x1) == ndims(x2) || error("Expected inputs with same number of dimensions!")
    cat(x1, x2; dims = ndims(x1) - 1)
end

function unetlayers(
    layers,
    sz;
    k_out = nothing,
    skip_upscale = 0,
    m_middle = _ -> (identity,),
)
    isempty(layers) && return m_middle(sz[end-1])

    layer, layers = layers[1], layers[2:end]
    outsz = Flux.outputsize(layer, sz)
    does_downscale = sz[1] ÷ 2 == outsz[1]

    if !does_downscale
        # If `layer` does not scale down the spatial dimensions, append
        # it to a Chain
        return Chain(layer, unetlayers(layers, outsz; k_out, skip_upscale)...)
    elseif does_downscale && skip_upscale > 0
        # If `layer` does scale down the spatial dimensions, but we don't
        # to upsample this one, recurse with modified arguments
        return Chain(
            layer,
            unetlayers(layers, outsz; skip_upscale = skip_upscale - 1, k_out)...,
        )
    else
        # `layer` scales down the spatial dimensions and we add an upsampling block
        # and a skip connection that scales the dimensions back up
        childunet = Chain(unetlayers(layers, outsz; skip_upscale)...)
        outsz = Flux.outputsize(childunet, outsz)

        k_in = sz[end-1]
        k_mid = outsz[end-1]
        k_out = isnothing(k_out) ? k_in : k_out
        return FastAI.Vision.Models.UNetBlock(
            Chain(layer, childunet),
            k_in,  # Input channels to upsampling layer
            k_mid,
            k_out,
        )
    end
end

iterlayers(m::Chain) = Iterators.flatten(iterlayers(l) for l in m.layers)
iterlayers(m) = (m,)

"""
    UNetBlock(m, k_in)
Given convolutional module `m` that halves the spatial dimensions
and outputs `k_in` filters, create a module that upsamples the
spatial dimensions and then aggregates features via  a skip connection.
"""
function UNetBlock(m_child, k_in, k_mid, k_out = 2k_in)
    return Chain(
        upsample = SkipConnection(
            Chain(
                child = m_child,                              # Downsampling and processing
                upsample = PixelShuffleICNR(k_mid, k_mid),  # Upsampling
            ),
            Parallel(catchannels, identity, BatchNorm(k_in)),
        ),
        act = xs -> relu.(xs),
        combine = UNetCombineLayer(k_in + k_mid, k_out),  # Data from both branches is combined
    )
end

function PixelShuffleICNR(k_in, k_out; r = 2)
    return Chain(convxlayer(k_in, k_out * (r^2), ks = 1), Flux.PixelShuffle(r))
end

function UNetCombineLayer(k_in, k_out)
    return Chain(convxlayer(k_in, k_out), convxlayer(k_out, k_out))
end

function UNetMiddleBlock(k)
    return Chain(convxlayer(k, 2k), convxlayer(2k, k))
end

function UNetFinalBlock(k_in, k_out)
    return Chain(ResBlock(1, k_in, k_in), convxlayer(k_in, k_out, ks = 1))
end

"""
    upsample_block_small(insize, k_out)
An upsampling block that increases the spatial dimensions of the input by 2
using pixel-shuffle upsampling.
"""
function upsample_block_small(insize, k_out; ks = 3, kwargs...)
    return Chain(Flux.PixelShuffle(2), convxlayer(insize[end-1] ÷ 4, k_out; kwargs...))
end

function conv_final(insize, k_out; ks = 1, kwargs...)
    return convxlayer(insize[end-1], k_out; ks = ks, kwargs...)
end


text:
for i=1:next_x_points
    #println("i:", i, ":")
    lyap_exp = lyaponuv(ts[end-sliding_window:end], J, m, r)
    #println("lyap_exp:", lyap_exp, "\n")
    tasks = Array(Future, sample_size)
    #println("tasks:", tasks, "\n")
    mu = mean(ts[end-sliding_window:end])
    #println("mu:", mu, "\n")

    diff = ts[end-sliding_window+1:end] - ts[end-sliding_window:end-1]
    mu = mean(diff)
    sigma = std(diff)

    sample_values = randn(sample_size) .* sigma .+ ts[end]

    for j=1:sample_size
        tempts=deepcopy(ts[end-sliding_window:end])
        append!(tempts, [sample_values[j]])
        tasks[j] = @spawn lyaponuv(tempts, J, m, r)
    end

    exponents = Array(Float64, sample_size)
    for j=1:sample_size
        exponents[j] = fetch(tasks[j])
    end
    #println("exponent:", mean(lyap_exp))

    exp_diff = abs(exponents .- lyap_exp)
    min_index = findmin(exp_diff)
    best_val = sample_values[min_index[2]]
    append!(ts, [best_val])
    println(i, "   best value:", best_val, "\t")
end

function lyaponuv_next(time_series, J, m, ref, sample_size)
    ts_diff = time_series[2:end] - time_series[1:end-1]
    sigma = std(ts_diff)
    samples = randn(sample_size) * sigma + time_series[end]
    @time norms, lyap_k = lyaponuv_k(time_series, J, m, ref)
    true_exponent = lyaponuv_exp(lyap_k)
    exponents = Array(Float64, sample_size)
    M = size(norms)[1]
    tasks = Array(Future, sample_size)

    for i=1:sample_size
        s = samples[i]
        tasks[i] = @spawn get_next(vcat(time_series, s), m, M, norms, ref, J)
        #@time exponents[i] = get_next(vcat(time_series, s), m, M, norms, ref, J)
        #@printf("process: %d\n", i)
    end

    for i=1:sample_size
        exponents[i]=fetch(tasks[i])
    end

    diff = abs(exponents .-  true_exponent)
    val, idx = findmin(diff)
    println("Next Value:", samples[idx])
    return(samples[idx])
end

function get_next(ts, m, M, norms, ref, J)

    attractor_array = attractor(ts, m, J)
    temp_norms = Array(Float64, M+1, M+1)
    temp_norms[1:M, 1:M] = norms
    col = column_norms(attractor_array, M+1)
    temp_norms[M+1, :] = col
    temp_norms[:, M+1] = col

    pairs=match_pairs(temp_norms)
    lyap_k_temp = follow_points(pairs, temp_norms, ref)
    return(lyaponuv_exp(lyap_k_temp))
end

max_stars_repo_path: saved & not used functions.jl

max_stars_repo_name: guypayeur/Generic-Pred

max_issues_repo_path: saved & not used functions.jl

max_issues_repo_name: guypayeur/Generic-Pred

max_forks_repo_path: saved & not used functions.jl

max_forks_repo_name: guypayeur/Generic-Pred


text:
using Sockets

import PeaceCypher: sign
using PeaceCypher

### Testing Notary

notary = Notary()
signer = newsigner(notary)

value = 19
signature = sign(value, signer)

@show verify(value, signature)
@show id(signer)==id(signature)

### Testing CypherSuite

crypto = CypherSuite(notary)

master = newsigner(crypto.notary)
masterid = id(master)

port = 2019

server = listen(port)

@sync begin
    @async begin
        masters = accept(server)
        masterss = secure(masters, crypto, master)
    end
    sleep(1.)
    slave = connect(port)
    @show slavess = secure(slave, crypto, masterid)
end


text:
function build_app(port=8050)
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])
    app.layout = html_div() do
        #Header
        dcc_markdown("
        #### Inventory Management Simulator
        *A discrete-time multi-product simulator for supply networks.*
        "),
        html_hr(),

        #Bill of Materials
        dcc_upload(
            id="upload_bill_of_materials",
            children=html_div([
                "1. Load bill of materials table: "
                html_a("Select File")
            ]),
            multiple=false
        ),
        dcc_store(id="store_bill_of_materials",data=JSON.json(DataFrame())),
        html_div(id="show_bill_of_materials"),
        html_hr(),

        #Lead Times
        dcc_upload(
            id="upload_lead_times",
            children=html_div([
                "2. Load network lead time table: "
                html_a("Select File")
            ]),
            multiple=false
        ),
        dcc_store(id="store_lead_times",data=JSON.json(DataFrame())),
        html_div(id="show_lead_times"),
        html_hr(),

        #Demand
        dcc_upload(
            id="upload_demand",
            children=html_div([
                "3. Load network demand table: "
                html_a("Select File")
            ]),
            multiple=false
        ),
        dcc_store(id="store_demand",data=JSON.json(DataFrame())),
        html_div(id="show_demand"),
        html_hr(),

        #Specify if tracking inventory position or echelon position
        html_div([
            "4. Select replenishment policy type:",
            dcc_radioitems(
                id="policy_type",
                options = [
                    (label="(r, Q)",value="rQ"),
                    (label="(s, S)",value="sS")
                ],
                value="rQ",
            ),
        ]),
        html_hr(),

        #Specify if tracking inventory position or echelon position
        html_div([
            "5. Select variable for inventory policy:",
            dcc_radioitems(
                id="policy_variable",
                options = [
                    (label="Inventory Position",value="inv_position"),
                    (label="Echelon Position",value="ech_position")
                ],
                value="inv_position",
            ),
        ]),
        html_hr(),

        #Inventory Policy
        dcc_upload(
            id="upload_policy",
            children=html_div([
                "6. Load inventory policy table: "
                html_a("Select File")
            ]),
            multiple=false
        ),
        dcc_store(id="store_policy",data=JSON.json(DataFrame())),
        html_div(id="show_policy"),
        html_hr(),

        #Specify if backlogging
        html_div([
            "7. Select operating mode:",
            dcc_radioitems(
                id="mode",
                options = [
                    (label="Backlog",value="backlog"),
                    (label="Lost Sales",value="lost_sales")
                ],
                value="backlog",
            ),
        ]),
        html_hr(),

        #Run Simulation
        html_div([
            "8. Number of periods to simulate: ",
            dcc_input(
                id="num_periods",
                type="number",
                min=0,
                value=0
            ),
            html_button("Build and Run Simulation Model", id="build_and_run", n_clicks=0),
        ]),
        dcc_confirmdialog(id = "sim_complete", message = ""),
        dcc_loading(html_div(id = "build_and_run_msg")),
        dcc_store(id = "store_inv_on_hand", data=JSON.json(DataFrame())),
        dcc_store(id = "store_inv_level", data=JSON.json(DataFrame())),
        dcc_store(id = "store_inv_position", data=JSON.json(DataFrame())),
        dcc_store(id = "store_ech_position", data=JSON.json(DataFrame())),
        dcc_store(id = "store_inv_pipeline", data=JSON.json(DataFrame())),
        dcc_store(id = "store_market_demand", data=JSON.json(DataFrame())),
        dcc_store(id = "store_replenishments", data=JSON.json(DataFrame())),

        #Download results
        html_div([
            html_button("Download Results", id="download", n_clicks=0),
            dcc_download(id="inv_onhand_result"),
            dcc_download(id="inv_level_result"),
            dcc_download(id="inv_position_result"),
            dcc_download(id="ech_position_result"),
            dcc_download(id="inv_pipeline_result"),
            dcc_download(id="market_demand_result"),
            dcc_download(id="replenishments_result")
        ])
    end

    build_callbacks(app)

    run_server(app, "0.0.0.0", port, debug=true)

    return nothing
end


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 = "FastJet"
version = v"3.4.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://fastjet.fr/repo/fastjet-3.4.0.tar.gz", "ee07c8747c8ead86d88de4a9e4e8d1e9e7d7614973f5631ba8297f7a02478b91"),
    ArchiveSource("http://fastjet.hepforge.org/contrib/downloads/fjcontrib-1.048.tar.gz", "f9989d3b6aeb22848bcf91095c30607f027d3ef277a4f0f704a8f0fc2e766981")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/fastjet-*/
export CXXFLAGS="-O3 -Wall"
export CFLAGS="-O3 -Wall"
if [[ "${target}" == *-freebsd* ]]; then
    # Needed to fix the following errors
    #   undefined reference to `backtrace_symbols'
    #   undefined reference to `backtrace'
    export LDFLAGS="-lexecinfo"
fi
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    # This is needed in order to build the shared library on Windows when we get
    #   libtool: warning: undefined symbols not allowed in x86_64-w64-mingw32 shared libraries; building static only
    FLAGS+=(LDFLAGS="-no-undefined")
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-auto-ptr --disable-static
make -j ${nprocs} "${FLAGS[@]}"
make install
tail -n 340 COPYING > LICENSE
install_license LICENSE
cp include/fastjet/config_win.h ${includedir}/fastjet/

cd $WORKSPACE/srcdir/fjcontrib-*/
for name in ClusteringVetoPlugin ConstituentSubtractor EnergyCorrelator FlavorCone GenericSubtractor JetCleanser JetFFMoments JetsWithoutJets LundPlane Nsubjettiness QCDAwarePlugin RecursiveTools ScJet SoftKiller SubjetCounting ValenciaPlugin VariableR
do
    cd $name
    rm example*
    c++ -fPIC -shared -I${includedir} -O3 -Wall ${LDFLAGS} [A-Z]*.cc -o "${libdir}/lib${name}.${dlext}" -lfastjet -lfastjettools
    cp [A-Z]*.hh ${includedir}/fastjet
    cd ..
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libConstituentSubtractor", :libConstituentSubtractor),
    LibraryProduct("libQCDAwarePlugin", :libQCDAwarePlugin),
    LibraryProduct("libJetFFMoments", :libJetFFMoments),
    LibraryProduct("libLundPlane", :libLundPlane),
    LibraryProduct("libJetsWithoutJets", :libJetsWithoutJets),
    LibraryProduct("libEnergyCorrelator", :libEnergyCorrelator),
    LibraryProduct("libfastjettools", :libfastjettools),
    LibraryProduct("libSoftKiller", :libSoftKiller),
    LibraryProduct("libRecursiveTools", :libRecursiveTools),
    LibraryProduct("libfastjet", :libfastjet),
    LibraryProduct("libValenciaPlugin", :libValenciaPlugin),
    LibraryProduct("libsiscone", :libsiscone),
    LibraryProduct("libNsubjettiness", :libNsubjettiness),
    LibraryProduct("libfastjetplugins", :libfastjetplugins),
    LibraryProduct("libVariableR", :libVariableR),
    LibraryProduct("libsiscone_spherical", :libsiscone_spherical),
    LibraryProduct("libJetCleanser", :libJetCleanser),
    LibraryProduct("libFlavorCone", :libFlavorCone),
    LibraryProduct("libGenericSubtractor", :libGenericSubtractor),
    LibraryProduct("libScJet", :libScJet),
    LibraryProduct("libSubjetCounting", :libSubjetCounting),
    LibraryProduct("libClusteringVetoPlugin", :libClusteringVetoPlugin),
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct HistLLHPrecalulations{N, WT, MPT, BVT, LGT}
    weights::Array{WT, N}
    midpoints::MPT
    volumes::BVT
    logabsgamma::Array{LGT, N}
end

function log_pdf_poisson(λ::T, k::U, logabsgamma::LAT) where {T, U, LAT}
    R = float(promote_type(T,U))
    if λ >= 0 && k >= 0 && isinteger(k)
        result = (iszero(k) ? R(k) : R(log(λ) * k)) - λ - logabsgamma
        R(result)
    else
        R(-Inf)
    end
end

function HistLLHPrecalulations(h::Histogram{<:Any, 1})
    T = promote_type(eltype(h.weights), eltype(h.edges[1]))
    return HistLLHPrecalulations(
        h.weights,
        StatsBase.midpoints(h.edges[1]),
        T[StatsBase.binvolume(h, i) for i in eachindex(h.weights)],
        T[logabsgamma(w + 1)[1] for w in h.weights]
    )
end

struct SpectrumDensity{HP <: HistLLHPrecalulations, D <: UvSpectrumDensity}
    hp::HP
end

@inline DensityInterface.DensityKind(::SpectrumDensity) = IsDensity()
function DensityInterface.logdensityof(object::SpectrumDensity{HP, D}, x) where {HP, D}
    llh = zero(eltype(x))
    for i in eachindex(object.hp.weights)
        llh += log_pdf_poisson(object.hp.volumes[i] * evaluate(D(x), object.hp.midpoints[i]), object.hp.weights[i], object.hp.logabsgamma[i])
    end
    llh
end

"""
    # opt_fit(DT::Type{<:UvSpectrumDensity}, h::Histogram{<:Any, 1},
             p0::AbstractVector, lower_bounds::AbstractVector, upper_bounds::AbstractVector)

Maximum Likelihood Estimation Fit of model density `d` on the histogram `h`.
"""
function opt_fit(DT::Type{<:UvSpectrumDensity}, h::Histogram{<:Any, 1},
             p0::AbstractVector, lower_bounds::AbstractVector, upper_bounds::AbstractVector,
             parshape::ValueShapes.AbstractValueShape = valshape(p0))
    hp = HistLLHPrecalulations(h)
    function f(p::AbstractVector{T}, hp=hp, ps=parshape) where {T}
        x = _par_in_input_form(ps(p))
        d = SpectrumDensity{typeof(hp), DT}(hp)
        return -DensityInterface.logdensityof(d, x)
    end
    opt_result = Optim.optimize( f, promote(lower_bounds, upper_bounds, p0)...,
        Fminbox(BFGS()); autodiff=:forward )
    # return opt_result.minimizer, opt_result
    return DT(_par_in_input_form(parshape(opt_result.minimizer))), opt_result
end

"""
    # opt_fit(DT::Type{<:UvSpectrumDensity}, h::Histogram{<:Any, 1},
             p0::NamedTuple, lower_bounds::NamedTuple, upper_bounds::NamedTuple)

Maximum Likelihood Estimation Fit of model density `d` on the histogram `h`.
"""
function opt_fit(DT::Type{<:UvSpectrumDensity}, h::Histogram{<:Any, 1},
             p0::NamedTuple,
             lower_bounds::NamedTuple,
             upper_bounds::NamedTuple,
             parshape::NamedTupleShape = valshape(p0))
    # flatten parameter and bounds
    flat_p0 = ValueShapes.unshaped(p0, parshape)
    T = eltype(flat_p0)
    rnt_lower = _remove_const_shapes_from_nt(lower_bounds, parshape)
    rnt_upper = _remove_const_shapes_from_nt(upper_bounds, parshape)
    flat_lower::Vector{T} = ValueShapes.unshaped(rnt_lower, valshape(rnt_lower))
    flat_upper::Vector{T} = ValueShapes.unshaped(rnt_upper, valshape(rnt_upper))

    # return parshape(minimizer)[], opt_result
    return opt_fit(DT, Histogram(h.edges[1], T.(h.weights)), flat_p0, flat_lower, flat_upper, parshape)
end


text:
# Softmax Regression/Gradient Descent applied to MNIST Dataset

push!(LOAD_PATH, ".")
push!(LOAD_PATH, "../common")

using MNIST
importall softmax_regression
importall optimization

# parameters
alpha = 0.5
max_iter = 20
min_err = 0.0

# load MNIST data
println("Loading MNIST dataset...")
X_train, y_train = traindata()
X_test, y_test = testdata()
y_train, y_test = y_train', y_test'

# bias feature, mean normalization
X_train = [255*ones(1, size(X_train, 2)); X_train] / 255.0
X_test = [255*ones(1, size(X_test, 2)); X_test] / 255.0

# make data more convenient for further operations
y_train = y_train .+ 1
y_test = y_test .+ 1
y_train = convert(Array{Int64}, y_train)
y_test = convert(Array{Int64}, y_test)

# theta initialization
theta = init_theta(size(X_train, 1), 10)

println("Training...")
tic()
theta, history = gradient_descent(h, J, g(10, y_train), X_train, y_train,
                                  theta, alpha, max_iter, min_err)
toc()

println("train accuracy: ", accuracy(predict(theta, X_train), y_train))
println("test accuracy: ", accuracy(predict(theta, X_test), y_test))

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
using GLPK
using Graphs
using GraphsOptim
using LinearAlgebra
using SparseArrays
using Test

@testset verbose = true "GraphsOptim.jl" begin
    include("minimum_cost_flow.jl")
    include("maximum_weight_matching.jl")
end


text:
using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq #, Plots
using Test

@testset "RC demo" begin
    R = 1.0
    C = 1.0
    V = 1.0
    @parameters t
    @named resistor = Resistor(R=R)
    @named capacitor = Capacitor(C=C)
    @named source = ConstantVoltage(V=V)
    @named ground = Ground()

    rc_eqs = [
            connect(source.p, resistor.p)
            connect(resistor.n, capacitor.p)
            connect(capacitor.n, source.n, ground.g)
            ]

    @named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground])
    sys = structural_simplify(rc_model)
    prob = ODAEProblem(sys, Pair[], (0, 10.0))
    sol = solve(prob, Tsit5())
    #plot(sol)

    @test isapprox(sol[capacitor.v][end], V, atol=1e-2)
end


text:
"""
# Wall

A `Wall` is a pair of `Loc`s with sanity-check style restrictions

"""
immutable Wall
    bottom::Loc
    top::Loc
    center::Loc
    direction::Direction
    adjacent::Vector4{Loc}

    function Wall(a::Loc, b::Loc)

        # cannot be on edges of board
        # TODO I imagine this might eventually have to go away
        # for perf
        @assert ! ((a.x == b.x == BOARD_MAX) ||
                   (a.x == b.x == BOARD_MIN) ||
                   (a.y == b.y == BOARD_MAX) ||
                   (a.y == b.y == BOARD_MIN))

        direction = ifelse(a.x == b.x, Vertical, Horizontal)

        # TODO probably also will have to be removed for perf
        if direction == Vertical
            @assert abs(a.y - b.y) == 2
        else
            @assert abs(a.x - b.x) == 2
        end

        if a < b
            bottom, top = a, b
        else
            bottom, top = b, a
        end

        center = mean(bottom, top)
        across_from_bottom = ifelse(direction == Vertical,
                                    left(bottom), below(bottom))
        across_from_center = ifelse(direction == Vertical,
                                    left(center), below(center))
        adjacent = Vector4(bottom, center,
                           across_from_bottom,
                           across_from_center)
        new(bottom, top, center, direction, adjacent)
    end
end

# have to define this explicitly because inner constructors still dont
# call convert for you? weird
Wall(a::Tuple{Int, Int}, b::Tuple{Int, Int}) = Wall(Loc(a), Loc(b))

center(w::Wall) = w.center
bottom(w::Wall) = w.bottom
top(w::Wall) = w.top
adjacent(w::Wall) = w.adjacent
direction(w::Wall) = w.direction

show(io::IO, w::Wall) = write(io, "($(w.bottom)---$(w.top)")

function separates(wall::Wall, l1::Loc, l2::Loc)
    adjacent_locs = adjacent(wall)
    both_adjacent = l1 in adjacent_locs && l2 in adjacent_locs
    return both_adjacent && (direction(wall) == Vertical && l1.y == l2.y ||
                             direction(wall) == Horizontal && l1.x == l2.x)
end

function any_separates(walls::Vector{Wall}, l1::Loc, l2::Loc)
    # TODO its a shame this has to exist and I can't just use
    # Base.any, which would involve an extra allocation for the
    # intermediate boolean list
    for wall in walls
        if separates(wall, l1, l2)
            return true
        end
    end
    return false
end

function allwalls()
    ret = Wall[]
    for i in 2:9
        for j in 1:8
            bottom_x = Loc(i,j)
            top_x = above(above(bottom_x))
            push!(ret, Wall(bottom_x, top_x))
            bottom_y = Loc(j,i)
            top_y = right(right(bottom_y))
            push!(ret, Wall(bottom_y, top_y))
        end
    end
    return ret
end

ALL_WALLS = allwalls()


text:
using PSFModels: fit

function generate_model(rng, model, params, inds)
    cartinds = CartesianIndices(inds)
    psf = model.(cartinds; params...)
    return psf
end

function test_fitting(rng, model, params, inds; kwargs...)
    psf = generate_model(rng, model, params, inds)
    _keys = keys(params)
    _vals = PSFModels.vector_from_params(Float64, params)
    # perturb starting guess by a little
    _vals .*= 1 .+ 1e-2 .* randn(rng, length(_vals))
    P0 = PSFModels.generate_params(_keys, _vals)
    P, bestfit = fit(model, P0, psf; kwargs...)
    for k in _keys
        if P[k] isa Tuple
            @test P[k][1] ≈ params[k][1] rtol=1e-2
            @test P[k][2] ≈ params[k][2] rtol=1e-2
        else
            @test P[k] ≈ params[k] rtol=1e-2
        end
    end
    @test bestfit.(CartesianIndices(psf)) ≈ psf rtol=1e-2
end

@testset "test fitting synthetic PSFs" begin
    inds = 1:30, 1:30
    @testset "gaussian" begin
        test_fitting(rng, gaussian, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds)
        test_fitting(rng, gaussian, (x=13.5, y=12.3, fwhm=(2.6, 2.4), amp=5), inds)
        test_fitting(rng, gaussian, (x=13.5, y=12.3, fwhm=(2.6, 2.4), theta=12, amp=5), inds)

        # test using L1 loss
        test_fitting(rng, gaussian, (x=13.5, y=12.3, fwhm=(2.6, 2.4), theta=12, amp=5), inds; loss=abs)
        # test using frozen variables
        test_fitting(rng, gaussian, (x=13.5, y=12.3, fwhm=(2.6, 2.4), amp=5), inds; func_kwargs=(;theta=0))

        @test_throws ArgumentError fit(gaussian, (x=12, y=13, fwhm=1, theta=13), randn(10, 10))
    end
    @testset "airydisk" begin
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds)
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=2.6, amp=5, ratio=0.12), inds)
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=(2.6, 2.4), amp=5), inds)
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=(2.6, 2.4), theta=12, amp=5), inds)

        # test using L1 loss
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds; loss=abs)
        # test using frozen variables
        test_fitting(rng, airydisk, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds; func_kwargs=(;ratio=0))

        @test_throws ArgumentError fit(airydisk, (x=12, y=13, fwhm=1, theta=13), randn(10, 10))
    end

    @testset "moffat" begin
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds)
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=2.6, amp=5, alpha=1.2), inds)
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=(2.6, 2.4), amp=5), inds)
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=(2.6, 2.4), amp=5, alpha=1.2), inds)
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=(2.6, 2.4), theta=12, amp=5), inds)
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=(2.6, 2.4), theta=12, amp=5, alpha=1.2), inds)

        # test using L1 loss
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds; loss=abs)
        # test using frozen variables
        test_fitting(rng, moffat, (x=13.5, y=12.3, fwhm=2.6, amp=5), inds; func_kwargs=(;alpha=1))

        @test_throws ArgumentError fit(moffat, (x=12, y=13, fwhm=1, theta=13), randn(10, 10))
    end
end


text:
module Funs

using Base.Meta

import Base: call, eltype

export Fun, Callable
export fcall, @fcall

# TODO: Use Base.return_types
eltype(f::Function) = Any
fcall(f::Function, args...; R::Type=Any) = f(args...)::R

# eltype(f::DataType) = f
fcall(f::DataType, args...; R::Type=f) = f(args...)::R

immutable Fun{T}
    fun::Function
end
eltype{R}(f::Fun{R}) = R
fcall{T}(f::Fun{T}, args...; R::Type=T) = f.fun(args...)::R
# call{T}(f::Fun{T}, args...; R::Type=T) = fcall(f, args..., R=R)
(f::Fun{T}){T}(args...; R::Type=T) = fcall(f, args..., R=R)

typealias Callable Union{Function, DataType, Fun}
# fcall, eltype

iscall(expr) = isexpr(expr, :call)
# Note: The called function is args[1], the arguments are args[2:end]
getcallargs(expr) = (@assert iscall(expr); expr.args)

macro fcall(args...)
    if length(args) == 1
        (expr,) = args
        callargs = getcallargs(expr)
        esc(:(fcall($(callargs...))))
    elseif length(args) == 2
        (R, expr) = args
        callargs = getcallargs(expr)
        esc(:(fcall(R=$R, $(callargs...))))
    else
        error("Expected: @fcall [<type>] <f>(<args>...)")
    end
end

end


text:
module JuliaPackageDemo

using MittagLeffler

export coolfunction1, coolfunction2

"""
    coolfunction1(x)

Computes `2*x + x^2` and returns the value.
"""
function coolfunction1(x)
    2*x + x^2
end

"""
    coolfunction2(x)

Computes `50*mittleff(0.5, -x) - 10` and returns the value.
"""
function coolfunction2(x)
    50*mittleff(0.5, -x) - 10
end

end # module

max_stars_repo_name: moustachio-belvedere/JuliaPackageDemo

max_issues_repo_name: moustachio-belvedere/JuliaPackageDemo

max_forks_repo_name: moustachio-belvedere/JuliaPackageDemo


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -9.75097970893739 <= q <= 8.436779502208893)

                     add_NL_constraint(m, :(softplus(0.5456914178120993 + -0.994137029231358*softplus(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.38451384631825425*softplus(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.9673084611668079*softplus(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + softplus(0.6143810270772536 + 0.32603207961941205*softplus(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.7660966505536861*softplus(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.30918160231413605*softplus(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + softplus(-0.4079478457750749 + 0.06443443960646311*softplus(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.09351478982242645*softplus(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + -0.7258007123919756*softplus(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofRepFAM
include("runtests.jl")
=#
println("Testing...")

using Test
using CytofRepFAM
import CytofRepFAM.MCMC
using BSON
using Random
using Distributions
using Flux

# tests for MCMC
include("MCMC_tests.jl")

# tests for feature allocation model (FAM)
include("Metropolis_tests.jl")

# tests for feature allocation model (FAM)
include("fam_tests.jl")

# tests for repulsive FAM with feature selection
include("repfamFS_tests.jl")

# tests for repulsive FAM with feature selection (parallel tempering version)
include("repfamFSPT_tests.jl")

# tests for repulsive FAM via trained priors
include("repfam_train_prior_tests.jl")

# tests for timing threading for R updates
include("update_r_tests.jl")

# tests for variational inference with FAM
include("vb_fam_tests.jl")


text:
using IterTools

export c_mn, e_i, m_list, n_list

"""
	c_mn(m, n, s_i, r_i, k3, y)

	Calculates the c_mn for general networks

	Uses as input
	m		Vector of integer
	n		Vector of integer
	s_i		stoichiometric coefficient for products
	r_i		stoichiometric coefficient for reactants
	k3		reaction rate for interactions
	y		molecule numbers

	Returns the time dependent c_mn

"""
function c_mn(m::Any, n::Any, s_i::Matrix{Int64}, r_i::Matrix{Int64}, k3::Vector{Float64}, y::Matrix{Float64})

	num_species = size(y,1)		# number of species
	num_int = length(k3)		# number of interaction reactions
	len_ts = size(y,2)			# length of time grid
	res = zeros(len_ts)			# return vector

	for β in 1:num_int

		prod_mu = ones(len_ts)		# product with mu_i
		prod1, prod2 = 1, 1			# initialize products

		for i in 1:num_species

			prod1 *= binomial(s_i[β,i], m[i])
			prod2 *= binomial(r_i[β,i], m[i])

			prod_mu *= binomial(r_i[β,i], n[i])

			if r_i[β,i] == n[i]+1
				prod_mu .*= y[i,:]
			elseif r_i[β,i] == n[i]+2
				prod_mu .*=  y[i,:].^2
			elseif r_i[β,i] > n[i]+2
				prod_mu .*=  y[i,:].^(r_i[β,i]-n[i])
			end

		end

		res += k3[β] * (prod1 - prod2) * prod_mu

	end

	return res

end

"""
	c_mn(m, n, s_i, r_i, k3, y, t)

	Calculates the c_mn for general networks, evaluated at time step t

	Uses as input
	m		Vector of integer
	n		Vector of integer
	s_i		stoichiometric coefficient for products
	r_i		stoichiometric coefficient for reactants
	k3		reaction rate for interactions
	y		molecule numbers
	t		time step

	Returns the c_mn at time step t

"""
function c_mn(m::Any, n::Any, s_i::Matrix{Int64}, r_i::Matrix{Int64}, k3::Vector{Float64}, y::Matrix{Float64}, t::Int64)

	num_species = size(y,1)		# number of species
	num_int = length(k3)		# number of interaction reactions
	len_ts = size(y,2)			# length of time grid
	res = 0.0					# return value

	for β in 1:num_int

		prod_mu = 1.0				# product with mu_i
		prod1, prod2 = 1, 1			# initialize products

		for i in 1:num_species

			prod1 *= binomial(s_i[β,i], m[i])
			prod2 *= binomial(r_i[β,i], m[i])

			prod_mu *= binomial(r_i[β,i], n[i])

			if r_i[β,i] == n[i]+1
				prod_mu *= y[i,t]
			elseif r_i[β,i] == n[i]+2
				prod_mu *=  y[i,t]^2
			elseif r_i[β,i] > n[i]+2
				prod_mu *=  y[i,t]^(r_i[β,i]-n[i])
			end

		end

		res += k3[β] * (prod1 - prod2) * prod_mu

	end

	return res

end

"""
	e_i(n, i)

	Returns a unit vector of length n with
	a 1 at position i

"""
function e_i(n::Int64, i::Int64)

	if n<=0 || i<=0 || i>n
		error("Invalid unit vector!")
	end

	e_i = zeros(Int64, n)
	e_i[i] = 1

	return e_i
end

"""
	calc_prod_resp(n::Tuple, resp_ti::Matrix{Float64})

	Calculates the product of responses with the vector n
	The matrix resp_ti is 'species × time'

	output = prod_i (n_i)! * R_i(tau,:)

"""
function calc_prod_resp(n::Tuple, resp_ti::Matrix{Float64})

	prod_resp = ones(size(resp_ti,2))
	for j in 1:length(n)
		if n[j] == 1
			prod_resp .*= resp_ti[j,:]
		elseif n[j] > 1
			prod_resp .*= factorial(n[j])*(resp_ti[j,:].^n[j])
		end
	end
	return prod_resp
end

"""
	calc_prod_resp(n::Tuple, resp_ti::Vector{Float64})

	Calculates the product of responses with the vector n
	The vector resp_ti is vector of response for each
	species but with specified time tau and tau'

	output = prod_i (n_i)! * R_i(tau,tau')

"""
function calc_prod_resp(n::Tuple, resp_ti::Vector{Float64})

	prod_resp = 1.0
	for j in 1:length(n)
		if n[j] == 1
			prod_resp *= resp_ti[j]
		elseif n[j] > 1
			prod_resp *= factorial(n[j])*(resp_ti[j]^n[j])
		end
	end
	return prod_resp
end

"""
	mn_list(max_mn)

	Generic function to calculate either m_list() or n_list()

"""
function mn_list(max_mn::Vector{Int64})

	num_species = length(max_mn)

	x = [collect(0:max_mn[j]) for j in 1:num_species]
	m_list = collect(Iterators.product(x...))

	return reshape(m_list, length(m_list), 1)

end

"""
	n_list(r_i)

	Returns a list of possible vectors n
	The input r_i is the reactant coefficient matrix

"""
function n_list(r_i)
	max_n = maximum(r_i, dims=1)
	return mn_list(vec(max_n))
end

"""
	m_list(s_i, r_i)

	Returns a list of possible vectors m
	The input s_i and r_i are the stoichiometric
	matrices

"""
function m_list(s_i, r_i)
	max_m = maximum([maximum(r_i, dims=1); maximum(s_i, dims=1)], dims=1)
	return mn_list(vec(max_m))
end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
#Multi-threaded:
#julia -t <number-of-threads> --project=@. main.jl
#Single thread:
#julia --project=@. main.jl
using NEOs
using Dates
using TaylorIntegration
using JLD
using PlanetaryEphemeris
@show Threads.nthreads()

#script parameters (TODO: use ArgParse.jl instead)
varorder = 5 # varorder is the order corresponding to the jet transport perturbation
nv = 7 #number of TaylorN variables
objname = "Bennu"
maxsteps = 10000
nyears_bwd = -13.0 # years in backward integration
nyears_fwd = 10.0 # years in forward integration
dense = false #true
quadmath = false # use quadruple precision
lyap = false #true # compute Lyapunov exponents
neosjlpath = pkgdir(NEOs)
dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!

# integration parameters
order = 25
abstol = 1.0E-20

### observation data files (ra/dec, del/dop)
optfile_bwd = joinpath(neosjlpath, "data", "101955_OPTICAL_1999_2006.dat")
optfile_fwd = joinpath(neosjlpath, "data", "101955_OPTICAL_2011_2018.dat")
radarfile_bwd = joinpath(neosjlpath, "data", "101955_RADAR_1999_2005.dat")
radarfile_fwd = joinpath(neosjlpath, "data", "101955_RADAR_2011.dat")

### path to local Solar System ephemeris file
ss_eph_file = "./sseph343ast016_p56y_et.jld"
ss16asteph_et = JLD.load(ss_eph_file, "ss16ast_eph")

### TaylorN variables setup
if lyap
    ### setup TaylorN variables with order=1, numvars=7
    TNvars = set_variables("δx", order=1, numvars=7)
    # dq corresponding to solution OR7 for Bennu
    dq_OR7 = [7.11289905874202e-7, 3.492170211757043e-7, 1.49735641150018e-8, -3.610309907354978e-10, 7.835479079267429e-11, -2.100961626780984e-10, -4.550232586244751e-14]
    dq = dq_OR7
else
    #### dq: perturbation to nominal initial condition (Taylor1 jet transport)
    #dq = Taylor1.(zeros(7), varorder)
    #dq[end][1] = 1e-14
    # dq: perturbation to nominal initial condition (TaylorN jet transport)
    dq = set_variables("δx", order=varorder, numvars=nv)
    for i in 1:6
        dq[i][1][i] = 1e-8
    end
    if get_numvars() == 7
        dq[7][1][7] = 1e-14
    end
end

### initial conditions from Bennu JPL solution \#97 ### Note: JPL sol \#97 uses nongrav function g(r)=(1 au)/r^2.25, while we use g(r)=(1 au)/r^2
jd0 = datetime2julian(DateTime(2011,1,1)) # JDTDB = 2455562.5
q00 = [-1.1951358208617802, -0.20726185835689961, -0.11201678544935807, 8.881637772597003e-5, -0.013056288090844732, -0.007377624521045638]
q0 = vcat(q00, 0.0) .+ dq ####vcat(q00, -4.614425051841E-14) #.+ dq

####integrator warmup
propagate(objname, dynamics, 1, jd0, nyears_fwd, ss16asteph_et, output=false, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol)
println("*** Finished warmup")

######Full jet transport integration
propagate(objname*"_bwd", dynamics, maxsteps, jd0, nyears_bwd, ss16asteph_et, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol, radarobsfile=radarfile_bwd, opticalobsfile=optfile_bwd, tord=10, niter=5)
propagate(objname*"_fwd", dynamics, maxsteps, jd0, nyears_fwd, ss16asteph_et, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol, radarobsfile=radarfile_fwd, opticalobsfile=optfile_fwd, tord=10, niter=5)
println("*** Finished asteroid ephemeris integration")


text:
mutable struct UniformBlockBuffer
	block::UniformBlock
	buffer::Vector{UInt8}

	# todo: add the option to have own GL buffer object
	# todo: add info for dirty region

	UniformBlockBuffer(block::UniformBlock) = new(block, zeros(UInt8, block.size))
end

mutable struct Material
	shader::Shader
	uniforms::Dict{Symbol, Any}
	blockBuffers::Vector{UniformBlockBuffer}
	states::RenderStateHolder

	Material(shader::Shader) = new(shader, Dict{Symbol, Any}(), Array{UniformBlockBuffer}(undef, length(shader.blocks)), RenderStateHolder())
end

isvalid(mat::Material) = isvalid(mat.shader)

function apply(mat::Material, renderer::Renderer)
	apply(mat.shader)

	for (u, v) in mat.uniforms
		setuniform(mat.shader, u, v)
	end

	for i in 1:length(mat.blockBuffers)
		if isassigned(mat.blockBuffers, i)
			buffer = mat.blockBuffers[i]
			setbufferdata(buffer.block, buffer.buffer)
		end
	end

	apply(mat.states, renderer)
end

resetstate(mat::Material, state::RenderState) = resetstate(mat.states, state)
setstate(mat::Material, state::RenderState) = setstate(mat.states, state)
getstate(mat::Material, ::Type{T}, default) where T <: RenderState = getstate(mat.states, T, default)

function setuniform(mat::Material, uniform::Symbol, value; allowAdd = true)
	if haskey(mat.shader.uniforms, uniform)
		var = mat.shader.uniforms[uniform]
		if inblock(var)
			if !isassigned(mat.blockBuffers, Int(var.blockId))
				if !allowAdd
					return false
				end
				mat.blockBuffers[var.blockId] = UniformBlockBuffer(mat.shader.blocks[var.blockId])
			end
			setvalue(var, mat.blockBuffers[var.blockId].buffer, value)
		else
			if !allowAdd && !haskey(mat.uniforms, uniform)
				return false
			end
			mat.uniforms[uniform] = value
		end
		return true
	end

	@info("setuniform(Material): attempt to set inexistent uniform $uniform")
	return false
end

function getuniform(mat::Material, uniform::Symbol)
	if haskey(mat.shader.uniforms, uniform)
		var = mat.shader.uniforms[uniform]
		if inblock(var)
			if isassigned(mat.blockBuffers, Int(var.blockId))
				return getvalue(var, mat.blockBuffers[var.blockId].buffer)
			end
		else
			return get(mat.uniforms, uniform, nothing)
		end
	end
	return nothing
end

function hasuniform(mat::Material, uniform::Symbol)
	if !haskey(mat.shader.uniforms, uniform)
		return false
	end
	var = mat.shader.uniforms[uniform]
	if inblock(var) && isassigned(mat.blockBuffers, Int(var.blockId))
		return true
	end
	return haskey(mat.uniforms, uniform)
end