text:
## TYPES AND CONSTRUCTOR

mutable struct DeterministicTunedModel{T,M<:Deterministic,A,AR} <: MLJBase.Deterministic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,Vector{<:Real}}
    operation
    range
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::A
    acceleration_resampling::AR
    check_measure::Bool
end

mutable struct ProbabilisticTunedModel{T,M<:Probabilistic,A,AR} <: MLJBase.Probabilistic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    operation
    range
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::A
    acceleration_resampling::AR
    check_measure::Bool
end

const EitherTunedModel{T,M} =
    Union{DeterministicTunedModel{T,M},ProbabilisticTunedModel{T,M}}

MLJBase.is_wrapper(::Type{<:EitherTunedModel}) = true

#todo update:
"""
    tuned_model = TunedModel(; model=nothing,
                             tuning=Grid(),
                             resampling=Holdout(),
                             measure=nothing,
                             weights=nothing,
                             repeats=1,
                             operation=predict,
                             range=nothing,
                             n=default_n(tuning, range),
                             train_best=true,
                             acceleration=default_resource(),
                             acceleration_resampling=CPU1(),
                             check_measure=true)

Construct a model wrapper for hyperparameter optimization of a
supervised learner.

Calling `fit!(mach)` on a machine `mach=machine(tuned_model, X, y)` or
`mach=machine(tuned_model, X, y, w)` will:

- Instigate a search, over clones of `model`, with the hyperparameter
  mutations specified by `range`, for a model optimizing the specified
  `measure`, using performance evaluations carried out using the
  specified `tuning` strategy and `resampling` strategy.

- Fit an internal machine, based on the optimal model
  `fitted_params(mach).best_model`, wrapping the optimal `model`
  object in *all* the provided data `X, y` (or in `task`). Calling
  `predict(mach, Xnew)` then returns predictions on `Xnew` of this
  internal machine. The final train can be supressed by setting
  `train_best=false`.

The `range` objects supported depend on the `tuning` strategy
specified. Query the `strategy` docstring for details. To optimize
over an explicit list `v` of models of the same type, use
`strategy=Explicit()` and specify `model=v[1]` and `range=v`.

The number of models searched is specified by `n`. If unspecified,
then `MLJTuning.default_n(tuning, range)` is used. When `n` is
increased and `fit!(mach)` called again, the old search history is
re-instated and the search continues where it left off.

If `measure` supports weights (`supports_weights(measure) == true`)
then any `weights` specified will be passed to the measure. If more
than one `measure` is specified, then only the first is optimized
(unless `strategy` is multi-objective) but the performance against
every measure specified will be computed and reported in
`report(mach).best_performance` and other relevant attributes of the
generated report.

Specify `repeats > 1` for repeated resampling per model evaluation. See
[`evaluate!](@ref) options for details.

*Important.* If a custom measure `measure` is used, and the measure is
a score, rather than a loss, be sure to check that
`MLJ.orientation(measure) == :score` to ensure maximization of the
measure, rather than minimization. Override an incorrect value with
`MLJ.orientation(::typeof(measure)) = :score`.

*Important:* If `weights` are left unspecified, and `measure` supports
sample weights, then any weight vector `w` used in constructing a
corresponding tuning machine, as in `tuning_machine =
machine(tuned_model, X, y, w)` (which is then used in *training* each
model in the search) will also be passed to `measure` for evaluation.

In the case of two-parameter tuning, a Plots.jl plot of performance
estimates is returned by `plot(mach)` or `heatmap(mach)`.

Once a tuning machine `mach` has bee trained as above, then
`fitted_params(mach)` has these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_fitted_params`| learned parameters of the optimal model

The named tuple `report(mach)` has these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_result`       | corresponding "result" entry in the history
`best_report`       | report generated by fitting the optimal model

plus others specific to the `tuning` strategy, such as `history=...`.

"""
function TunedModel(;model=nothing,
                    tuning=Grid(),
                    resampling=MLJBase.Holdout(),
                    measures=nothing,
                    measure=measures,
                    weights=nothing,
                    operation=predict,
                    ranges=nothing,
                    range=ranges,
                    train_best=true,
                    repeats=1,
                    n=nothing,
                    acceleration=default_resource(),
                    acceleration_resampling=CPU1(),
                    check_measure=true)

    range === nothing && error("You need to specify `range=...`.")
    model == nothing && error("You need to specify model=... .\n"*
                              "If `tuning=Explicit()`, any model in the "*
                              "range will do. ")

    if model isa Deterministic
        tuned_model = DeterministicTunedModel(model, tuning, resampling,
                                       measure, weights, operation, range,
                                              train_best, repeats, n,
                                              acceleration,
                                              acceleration_resampling,
                                              check_measure)
    elseif model isa Probabilistic
        tuned_model = ProbabilisticTunedModel(model, tuning, resampling,
                                       measure, weights, operation, range,
                                              train_best, repeats, n,
                                              acceleration,
                                              acceleration_resampling,
                                              check_measure)
    else
        error("Only `Deterministic` and `Probabilistic` "*
              "model types supported.")
    end

    message = clean!(tuned_model)
    isempty(message) || @info message

    return tuned_model

end

function MLJBase.clean!(tuned_model::EitherTunedModel)
    message = ""
    if tuned_model.measure === nothing
        tuned_model.measure = default_measure(tuned_model.model)
        if tuned_model.measure === nothing
            error("Unable to deduce a default measure for specified model. "*
                  "You must specify `measure=...`. ")
        else
            message *= "No measure specified. "*
            "Setting measure=$(tuned_model.measure). "
        end
    end
    if !(tuned_model.acceleration isa Union{CPU1, CPUProcesses})
        message *= "Supported `acceleration` types are `CPU1` "*
        "and `CPUProcesses`. Setting `acceleration=CPU1()`. "
        tuned_model.acceleration = CPU1()
    end
    return message
end

## FIT AND UPDATE METHODS

# returns a (model, result) pair for the history:
function event(model, resampling_machine, verbosity, tuning, history)
    resampling_machine.model.model = model
    verb = (verbosity == 2 ? 0 : verbosity - 1)
    fit!(resampling_machine, verbosity=verb)
    e = evaluate(resampling_machine)
    r = result(tuning, history, e)

    if verbosity > 2
        println(params(model))
    end
    if verbosity > 1
        println("$r")
    end

    return deepcopy(model), r
end

function assemble_events(models, resampling_machine,
                         verbosity, tuning, history, acceleration::CPU1)
    map(models) do m
        event(m, resampling_machine, verbosity, tuning, history)
    end
end

function assemble_events(models, resampling_machine,
                         verbosity, tuning, history, acceleration::CPUProcesses)
    pmap(models) do m
        event(m, resampling_machine, verbosity, tuning, history)
    end
end

# history is intialized to `nothing` because it's type is not known.
_vcat(history, Δhistory) = vcat(history, Δhistory)
_vcat(history::Nothing, Δhistory) = Δhistory
_length(history) = length(history)
_length(::Nothing) = 0

# builds on an existing `history` until the length is `n` or the model
# supply is exhausted (method shared by `fit` and `update`). Returns
# the bigger history:
function build(history, n, tuning, model::M,
               state, verbosity, acceleration, resampling_machine) where M
    j = _length(history)
    models_exhausted = false
    while j < n && !models_exhausted
        _models = models!(tuning, model, history, state, verbosity)
        models = _models === nothing ? M[] : collect(_models)
        Δj = length(models)
        Δj == 0 && (models_exhausted = true)
        shortfall = n - Δj
        if models_exhausted && shortfall > 0 && verbosity > -1
            @info "Only $j (of $n) models evaluated.\n"*
            "Model supply exhausted. "
        end
        Δj == 0 && break
        shortfall < 0 && (models = models[1:n - j])
        j += Δj

        Δhistory = assemble_events(models, resampling_machine,
                                 verbosity, tuning, history, acceleration)
        history = _vcat(history, Δhistory)
    end
    return history
end

function MLJBase.fit(tuned_model::EitherTunedModel{T,M},
                     verbosity::Integer, data...) where {T,M}
    tuning = tuned_model.tuning
    model = tuned_model.model
    range = tuned_model.range
    n = tuned_model.n === nothing ?
        default_n(tuning, range) : tuned_model.n

    verbosity < 1 || @info "Attempting to evaluate $n models."

    acceleration = tuned_model.acceleration

    state = setup(tuning, model, range, verbosity)

    # instantiate resampler (`model` to be replaced with mutated
    # clones during iteration below):
    resampler = Resampler(model=model,
                          resampling    = tuned_model.resampling,
                          measure       = tuned_model.measure,
                          weights       = tuned_model.weights,
                          operation     = tuned_model.operation,
                          check_measure = tuned_model.check_measure,
                          repeats       = tuned_model.repeats,
                          acceleration  = tuned_model.acceleration_resampling)
    resampling_machine = machine(resampler, data...)

    history = build(nothing, n, tuning, model, state,
                    verbosity, acceleration, resampling_machine)

    best_model, best_result = best(tuning, history)
    fitresult = machine(best_model, data...)

    if tuned_model.train_best
        fit!(fitresult, verbosity=verbosity - 1)
        prereport = (best_model=best_model, best_result=best_result,
                     best_report=MLJBase.report(fitresult))
    else
        prereport = (best_model=best_model, best_result=best_result,
                     best_report=missing)
    end

    report = merge(prereport, tuning_report(tuning, history, state))
    meta_state = (history, deepcopy(tuned_model), state, resampling_machine)

    return fitresult, meta_state, report
end

function MLJBase.update(tuned_model::EitherTunedModel, verbosity::Integer,
                        old_fitresult, old_meta_state, data...)

    history, old_tuned_model, state, resampling_machine = old_meta_state
    acceleration = tuned_model.acceleration

    tuning = tuned_model.tuning
    range = tuned_model.range
    model = tuned_model.model

    # exclamation points are for values actually used rather than
    # stored:
    n! = tuned_model.n === nothing ?
        default_n(tuning, range) : tuned_model.n

    old_n! = old_tuned_model.n === nothing ?
        default_n(tuning, range) : old_tuned_model.n

    if MLJBase.is_same_except(tuned_model, old_tuned_model, :n) &&
        n! >= old_n!

        verbosity < 1 || @info "Attempting to add $(n! - old_n!) models "*
        "to search, bringing total to $n!. "

        history = build(history, n!, tuning, model, state,
                        verbosity, acceleration, resampling_machine)

        best_model, best_result = best(tuning, history)

        fitresult = machine(best_model, data...)

        if tuned_model.train_best
            fit!(fitresult, verbosity=verbosity - 1)
            prereport = (best_model=best_model, best_result=best_result,
                         best_report=MLJBase.report(fitresult))
        else
            prereport = (best_model=best_model, best_result=best_result,
                         best_report=missing)
        end

        _report = merge(prereport, tuning_report(tuning, history, state))

        meta_state = (history, deepcopy(tuned_model), state,
                      resampling_machine)

        return fitresult, meta_state, _report

    else

        return fit(tuned_model, verbosity, data...)

    end

end

MLJBase.predict(tuned_model::EitherTunedModel, fitresult, Xnew) =
    predict(fitresult, Xnew)

function MLJBase.fitted_params(tuned_model::EitherTunedModel, fitresult)
    if tuned_model.train_best
        return (best_model=fitresult.model,
                best_fitted_params=fitted_params(fitresult))
    else
        return (best_model=fitresult.model,
                best_fitted_params=missing)
    end
end

## METADATA

MLJBase.supports_weights(::Type{<:EitherTunedModel{<:Any,M}}) where M =
    MLJBase.supports_weights(M)

MLJBase.load_path(::Type{<:DeterministicTunedModel}) =
    "MLJTuning.DeterministicTunedModel"
MLJBase.package_name(::Type{<:EitherTunedModel}) = "MLJTuning"
MLJBase.package_uuid(::Type{<:EitherTunedModel}) = "MLJTuning"
MLJBase.package_url(::Type{<:EitherTunedModel}) =
    "https://github.com/alan-turing-institute/MLJTuning.jl"
MLJBase.is_pure_julia(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.target_scitype(M)


text:
@coeff function zeta(reactor::Reactor)
  B_max_hat ^ 2 / ( 2 * mu_0 * sigma_max_hat )
end

ext: jl


text:
"""
    AccelTurnrate
An action type with a longitudinal acceleration and a turn rate

# Fields
- `a::Float64` longitudinal acceleration [m/s^2]
- `ω::Float64` desired turn rate [rad/sec]
"""
struct AccelTurnrate
    a::Float64
    ω::Float64
end
Base.show(io::IO, a::AccelTurnrate) = @printf(io, "AccelTurnrate(%6.3f,%6.3f)", a.a, a.ω)
Base.length(::Type{AccelTurnrate}) = 2

"""
Convert an input vector containing desired acceleration and turn rate
into the AccelTurnrate type
"""
Base.convert(::Type{AccelTurnrate}, v::Vector{Float64}) = AccelTurnrate(v[1], v[2])

"""
Return a vector containing the desired acceleration and turn rate
"""
function Base.copyto!(v::Vector{Float64}, a::AccelTurnrate)
    v[1] = a.a
    v[2] = a.ω
    v
end

"""
Propagate vehicle forward in time using a desired acceleration and
turn rate
"""
function propagate(veh::Entity{VehicleState,D,I}, action::AccelTurnrate, roadway::Roadway, Δt::Float64; n_integration_steps::Int=4) where {D,I}

    a = action.a # accel
    ω = action.ω # turnrate

    x, y, θ = posg(veh.state)
    v = vel(veh.state)

    δt = Δt / n_integration_steps

    for i in 1 : n_integration_steps
        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end
function Base.get(::Type{AccelTurnrate}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    accel = get(ACC, rec, roadway, vehicle_index, pastframe)
    turnrate = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
    AccelTurnrate(accel, turnrate)
end


text:
# ##  AccuracyAtTopKernels.jl
# This repository is a complementary material to our paper ["General Framework for Nonlinear Binary Classification on Top Samples"]().
#
# # Installation
#
# This package can be installed using pkg REPL as follows
# ```julia
# (v1.2) pkg> add https://github.com/VaclavMacha/AccuracyAtTopKernels_new.jl
# ```

# # Usage
#
# This package provides following methods:
# * Methods:
#     1. `TopPush(C, surrogate)`
#     2. `TopPushK(K, C, surrogate)`
#     3. `PatMat(τ C, surrogate, surrogate)`
#     4. `PatMatNP(τ C, surrogate, surrogate)`
# * Surrogates:
#     1. `Hinge(ϑ)`
#     2. `Quadratic(ϑ)`
# * Problems:
#     1. primal
#     2. dual
# * Solvers:
#     1. `General`
#     2. `Gradient`
#     3. `Coordinate`
#
# ## Simple example
#
# In this example we use our package and StatsPlots to visualize the solution.
using AccuracyAtTopKernels, StatsPlots, Random

Random.seed!(1234)

function random_circle(n::Int; radius::Tuple = (0,1), origin::Tuple = (0,0))
    rmin, rmax = round.(Int, radius.^2)
    r2 = (rmax - rmin) .* rand(n) .+ rmin
    θ  = 2π .* rand(n)

    x = origin[1] .+ sqrt.(r2) .* cos.(θ)
    y = origin[2] .+ sqrt.(r2) .* sin.(θ)
    return x, y
end

nneg, npos = 1000, 200
xneg, yneg = random_circle(nneg; radius = (3, 5));
xpos, ypos = random_circle(npos; radius = (0, 2));

X = [xneg yneg; xpos ypos];
y = 1:(nneg+npos) .> nneg;

#  As an example, we use the following linearly nonseparable data

plt1 = scatter(X[y .== 0, 1], X[y .== 0, 2], label = "positives")
scatter!(X[y .== 1, 1], X[y .== 1, 2], label = "negatives", dpi = 300)
savefig(plt1, "data.png")

# ![Simple example](./scripts/data.png)

# The following example shows, how to compute the primal and dual problem for ToppushK method. To solve, we use our gradient descent based solver for the primal problem and our coordinate descent based solver for the dual problem.

model  = TopPushK(10, 1, Quadratic(1))
data_p = Primal(X, y);
data_d = Dual(model, X, y; kernel = GaussianKernel());

solution_p = solve(Gradient(maxiter = 1000, optimizer = Descent(0.0001)), model, data_p);
solution_d = solve(Coordinate(maxiter = 1000), model, data_d);

# For simple visual verification that methods works, we use density estimates of classification scores.
scores_p = scores(model, data_p, solution_p.w);
scores_d = scores(model, data_d, solution_d.α, solution_d.β);

plt_p = density(scores_p[y .== 0]; trim = true, fill = (0, 0.15), title = "primal", legend = false)
density!(scores_p[y .== 1]; trim = true, fill = (0, 0.15))
plt_d = density(scores_d[y .== 0]; label = "negatives", trim = true, fill = (0, 0.15), title = "dual with kernel")
density!(scores_d[y .== 1]; label = "positives", trim = true, fill = (0, 0.15))
plt = plot(plt_p, plt_d, layout = (1,2), size = (800, 400), dpi = 300)
savefig(plt, "scores.png")

# ![Comparison of scores.](./scripts/scores.png)


text:
module MomentClosure

import Catalyst: species, params, reactions, speciesmap, paramsmap, numspecies,
				 numreactions, numparams, substoichmat, prodstoichmat, netstoichmat

using ModelingToolkit
using Symbolics: value, var_from_nested_derivative, map_subscripts

using SciMLBase, SciMLBase.EnsembleAnalysis
using DiffEqJump
using Random
using Distributions: Geometric

using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint
using SymbolicUtils: Symbolic, Term, Real, expand, simplify, operation,
					 arguments, @rule, @acrule, isnotflat, flatten_term,
					 istree, FnType

using OrderedCollections: OrderedDict
using Combinatorics
using TupleTools: sort
using Cumulants
using Latexify

using DocStringExtensions

export generate_central_moment_eqs, generate_raw_moment_eqs, bernoulli_moment_eqs,
       propensities, ReactionSystemMod, species, params, speciesmap, paramsmap,
	   numspecies, numreactions, numparams, get_S_mat,
       moment_closure, deterministic_IC, JumpProblem,
	   get_raw_moments, get_central_moments, get_cumulants, get_moments_FSP,
	   linear_mapping_approximation

# reexporting from ModelingToolkit & Symbolics
# needed for ReactionSystemMod definition
export @parameters, @variables

include("reaction_systems.jl")
include("moment_equations.jl")
include("symbolic.jl")
include("moment_convert.jl")
include("stochastic_stoichiometry.jl")
include("central_moment_equations.jl")
include("raw_moment_equations.jl")
include("closure_methods/closure.jl")
include("bernoulli.jl")
include("utils.jl")
include("closure_methods/zero_closure.jl")
include("closure_methods/normal_closure.jl")
include("closure_methods/log_normal_closure.jl")
include("closure_methods/poisson_closure.jl")
include("closure_methods/gamma_closure.jl")
include("closure_methods/derivative_matching.jl")
include("closure_methods/conditional_gaussian.jl")
include("closure_methods/conditional_derivative_matching.jl")
include("closure_methods/linear_mapping_approximation.jl")
include("stochastic_simulation.jl")

end


text:
function plot_d_nu(pair, x_de, r̂)
  # extract pair
  d_nu, d_de = pair

  # sort vectors for plotting
  inds = sortperm(x_de)
  xs = x_de[inds]
  r̂s = r̂[inds]

  # actual and estimated d_nu
  ps = pdf.(d_nu, xs)
  p̂s = r̂s .* pdf.(d_de, xs)

  plot(xs, ps, label="d_nu (actual)")
  plot!(xs, p̂s, label="d_nu (estimate)")
end


text:
"""
    module NewtonsCradle

Model of [Newtons' Cradle](https://en.wikipedia.org/wiki/Newton%27s_cradle).

The balls of the Newton Cradle are initially placed in such a way, that there is a small
distance between the balls, so the balls are **not** touching each other. This distance was
determined (by trial-and-error) in such a way, that collisions between the balls
occur in sequence (so if ball1 hits ball2, then ball3 is hit by ball2 only after the
contact between ball1 and ball2 has already ended; this can be easily determined by
setting `simulate!(..., log=true)`). Only in this case, the simulation matches the
real experiment. If the balls are directly touching each other, multiple contacts
occur at the same time instant and the simulation result is then no longer according
to reality.

See also:

M. Payr (2008): **[An Experimental and Theoretical Study of Perfect Multiple
                Contact Collisions in Linear Chains of Balls.](http://www.zfm.ethz.ch/e/dynamics/payr_collisions.htm)**
                Dissertation, ETH Zürich.
"""
module NewtonsCradle

using ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

vmatVisual = VisualMaterial(color="LightBlue", transparency=0.1)
vmatSolids   = VisualMaterial(color="Red"      , transparency=0.0)

diameter   = 0.02
lengthWire = 0.15
rsmall     = diameter/100

Lx = diameter/2
Ly = 5*diameter
Lz = Lx

Pendulum = Model(
    wire = Object3D(feature=Visual(shape=Cylinder(axis=3, diameter=diameter/5, length=lengthWire),
                                   visualMaterial=vmatVisual)),
    frame1 = Object3D(parent=:wire,
                      translation=[0.0, 0.0, lengthWire/2]),
    frame2 = Object3D(parent=:wire,
                      translation=[0.0, 0.0, -lengthWire/2]),
    cylinder = Object3D(parent=:frame1,
                        feature=Visual(shape=Cylinder(axis=1, diameter=0.6*Lx, length=1.4*Lx),
                                       visualMaterial=vmatSolids)),
    sphere = Object3D(parent=:frame2,
                      feature=Solid(shape=Sphere(diameter=diameter),
                                    solidMaterial="Steel",
                                    visualMaterial=vmatSolids,
                                    contactMaterial="BilliardBall", collision=true)),
)

Cradle = Model(
    gravField = UniformGravityField(g=9.81, n=[0, 0,-1]),
    world = Object3D(feature=Scene(gravityField=:gravField,
                                   nominalLength=Ly,
                                   enableContactDetection=true,
                                   elasticContactReductionFactor=1.0,
                                   animationFile="NewtonsCradle.json",
                                   visualizeBoundingBox=true)),
    box = Object3D(parent=:world,
                   feature=Visual(shape=Box(lengthX=Lx, lengthY=Ly, lengthZ=Lz),
                                  visualMaterial=vmatVisual)),

    frame1 = Object3D(parent=:box, translation=[0.0, -2*(diameter+rsmall), 0.0]),
    frame2 = Object3D(parent=:box, translation=[0.0,   -(diameter+rsmall), 0.0]),
    frame3 = Object3D(parent=:box, translation=[0.0,                  0.0, 0.0]),
    frame4 = Object3D(parent=:box, translation=[0.0,    (diameter+rsmall), 0.0]),
    frame5 = Object3D(parent=:box, translation=[0.0,  2*(diameter+rsmall), 0.0]),

    pendulum1 = Pendulum,
    pendulum2 = Pendulum,
    pendulum3 = Pendulum,
    pendulum4 = Pendulum,
    pendulum5 = Pendulum,

    rev1 = Revolute(obj1=:frame1, obj2=:(pendulum1.frame1), axis=1, canCollide=true, phi=Var(init=-pi/3)),
    rev2 = Revolute(obj1=:frame2, obj2=:(pendulum2.frame1), axis=1, canCollide=true),
    rev3 = Revolute(obj1=:frame3, obj2=:(pendulum3.frame1), axis=1, canCollide=true),
    rev4 = Revolute(obj1=:frame4, obj2=:(pendulum4.frame1), axis=1, canCollide=true),
    rev5 = Revolute(obj1=:frame5, obj2=:(pendulum5.frame1), axis=1, canCollide=true)
)

newtonsCradle = @instantiateModel(buildModia3D(Cradle), unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 5.0
tolerance = 1e-8
requiredFinalStates = [-1.180791527568564, -1.033934458288708, 0.04982982739055273, 0.018466029517559813, 0.049844032881367094, 0.018471344490359712, 0.04986127142907837, 0.018476465719447804, 0.04987561416137295, 0.018481494292598495]
simulate!(newtonsCradle, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(newtonsCradle, ["rev1.phi" "rev5.phi"; "rev1.w" "rev5.w"], figure=1)

end


text:
"""
x = nnls(A, b; ...)

Solves non-negative least-squares problem by the active set method
of Lawson & Hanson (1974).

Optional arguments:
    tol: tolerance for nonnegativity constraints
    max_iter: maximum number of iterations (counts inner loop iterations)

References:
    Lawson, C.L. and R.J. Hanson, Solving Least-Squares Problems,
    Prentice-Hall, Chapter 23, p. 161, 1974.
"""
function nnls(A::Matrix{Float64},
              b::Vector{Float64};
              tol::Float64=1e-8,
              max_iter=30*size(A,2),
              kwargs...)

    # dimensions, initialize solution
    m,n = size(A)
    x = zeros(n)

    # P is a bool array storing positive elements of x
    # i.e., x[P] > 0 and x[~P] == 0
    P = zeros(Bool,n)

    # We have reached an optimum when either:
    #   (a) all elements of x are positive (no nonneg constraints activated)
    #   (b) ∂f/∂x = A' * (b - A*x) > 0 for all nonpositive elements of x
    w = A' * (b - A*x)
    iter = 0
    while sum(P)<n && any(w[~P].>tol) && iter < max_iter

        # find i that maximizes w, restricting i to indices not in P
        # Note: the while loop condition guarantees at least one w[~P]>0
        i = indmax(w .* ~P)

        # Move i to P
        P[i] = true

        # Solve least-squares problem, with zeros for columns/elements not in P
        Ap = zeros(m,n)
        Ap[:,P] = A[:,P]
        s = pinv(Ap)*b
        s[~P] = 0.0 # zero out elements not in P

        # Inner loop: deal with negative elements of s
        while any(s[P].<=tol)
            iter += 1

            # find indices in P where s is negative
            ind = (s.<=tol) & P

            # calculate step size, α, to prevent any xᵢ from going negative
            α = minimum(x[ind] ./ (x[ind] - s[ind]))

            # update solution (pushes some xᵢ to zero)
            x += α*(s-x)

            # Remove all i in P where x[i] == 0
            for i = 1:n
                if P[i] && abs(x[i]) < tol
                    # remove i from P
                    P[i] = false
                    # zero out column i of Ap
                    Ap[:,i] *= 0.0
                end
            end

            # Solve least-squares problem again, zeroing nonpositive columns
            s = pinv(Ap)*b
            s[~P] = 0.0 # zero out elements not in P
        end

        # update solution
        x = deepcopy(s)
        w = A' * (b - A*x)
    end
    return x
end

function nnls(A::Matrix{Float64},
              B::Matrix{Float64};
              use_parallel = true,
              kwargs...)

    m,n = size(A)
    k = size(B,2)

    if use_parallel && nprocs()>1
        X = SharedArray(Float64,n,k)
        @sync @parallel for i = 1:k
            X[:,i] = nnls(A, B[:,i]; kwargs...)
        end
        X = convert(Array,X)
    else
        X = Array(Float64,n,k)
        for i = 1:k
            X[:,i] = nnls(A, B[:,i]; kwargs...)
        end
    end

    return X
end


text:
using PyCall
using PyPlot
using SpecialFunctions
# uncomment for non-jupyter use:
#pygui(true)

### FUNCTIONS, used for several subtasks:
# %%

# the Bernoulli PDF:
function bernoulli_pdf(θ, y)
    if 0 <= θ <= 1 && (y==0 || y==1)    # Only admissible values
        return θ^y * (1-θ)^(1-y)
    else
        return 0
    end
end

# Likelihood fct of a (series of) coin tosses with bias:
function likelihood_fct(θ,y)
    # Make array s.t. multiple theta values can be input:
    likelihood = Array{Float64}(undef,length(θ))
    # Do that for every theta value requested:
    for ii in 1:length(θ)
        # Elementwise application of the beronoulli_pdf over outcome vector y, then MULTIPLY together:
        likelihood[ii] = prod(bernoulli_pdf.(θ[ii], y))
    end
    return likelihood
end

#%%

# Log-Bernoulli PDF:
function log_bernoulli_pdf(θ,y)
    if 0 < θ < 1 && (y==0 || y==1)    # Only admissible values
        return y * log(θ) + (1-y) * log(1-θ)
    else
        return -Inf
    end
end

# Log-Likelihood fct of a (series of) coin tosses with bias:
function log_likelihood_fct(θ, y)
    # Make array s.t. multiple theta values can be input:
    likelihood = Array{Float64}(undef,length(θ))
    # Do that for every theta value requested:
    for ii in 1:length(θ)
        # Elementwise application of the log_beronoulli_pdf over outcome vector y, then SUM together:
        likelihood[ii] = sum(log_bernoulli_pdf.(θ[ii],y))
    end
    return likelihood
end

#%%

### Choose and run the subtask:

###########
### Task 2
# Given the parameter, plot outcome probabilities
#%%
θ = 0.5
outcomes = [0, 1]

figure()
bar(outcomes, bernoulli_pdf.(θ,outcomes))
xticks([0, 1])
xlabel("y")
ylabel("θ")
plt.show()

#%%

θ = 0.25

figure()
bar(outcomes, bernoulli_pdf.(θ,outcomes))
xticks([0, 1])
xlabel("y")
ylabel("θ")
plt.show()

#%%

###########
### Task 3
# Given the outcome, plot parameter likelihoods
#%%
θ = 0:0.01:1;
outcome = 1;

figure()
plot(likelihood_fct.(θ,outcome))
plt.show()

#%%
outcome = 0

figure()
plot(likelihood_fct.(θ,outcome))
plt.show()
#%%

###########
### Task 4
# Evaluate likelihood fct after n flips

#%%
θ = 0.5

# a)
println("a) Regular Likelihood:")

for i in 1:3
    n = 10 * 100^(i-1)

    outcome = (rand(n) .>= 0.5) * 1

    print("L(y_1,y_2,...,y_$n | θ = $θ) = ")
    println(likelihood_fct(θ,outcome))
end

println("\nWith rising n, the likelihood -> 0")

#%%
# b) and c)
println("\nb) log-likelihoods, l, and L = exp(l):")

for i in 1:4
    n = 10 * 100^(i-1)

    outcome = (rand(n) .>= 0.5) * 1
    l = log_likelihood_fct(θ,outcome)
    L = exp.(l)

    print("l(y_1,y_2,...,y_$n | θ = $θ) = ")
    println(l)
    print("L(y_1,y_2,...,y_$n | θ = $θ) = ")
    println(L)
end

println("\nYes, n > 100000 could be evaluated without over-/underflow in logarithmic space.")

#%%
# d)
println("\nd) Plotting the likelihood function for a series")
θ = 0:0.01:1;

outcomeSeries = [[1], [1 1], [1 1 0 1]]

for thisOutcome in outcomeSeries
    figure()
    subplot(1,2,1)
    plot(θ,likelihood_fct(θ,thisOutcome),label="L($thisOutcome|θ)")
    grid()
    legend()
    subplot(1,2,2)
    plot(θ,log_likelihood_fct(θ,thisOutcome),label="l($thisOutcome|θ)")
    grid()
    legend()
    plt.show()
end
#%%

###########
### Task 5
# Implement the beta function and create plots

### i)
# Logarithmic version of the beta-distribution
function Log_Beta(θ,a::Any,b::Any)
    #using SpecialFunctions
    if 0 < θ < 1
        # Numerator
        num = log(θ)*(a-1) + log(1-θ)*(b-1)
        # The denominator (Beta-function) can be integrated or made up of Gamma-functions
        # Since we only treat integer values here, Γ(a) = (a-1)!
        #den = log(factorial(a-1) * factorial(b-1) / factorial(a+b-1))
        den = log(SpecialFunctions.beta(a,b))
        return num - den
    else
        return -Inf
    end
end

#%%
### ii) sub-tasks 1-3
function Log_posterior(θ,a,b,y)
    heads = sum(y)
    tails = length(y) - heads
    posterior = Array{Float64}(undef,length(θ))
    # Log addition, the brackets are for converting into a single float:
    for ii in 1:length(θ)
        posterior[ii] = (Log_Beta(θ[ii],a,b) .+ log_likelihood_fct(θ[ii], y) .-
                    log(SpecialFunctions.beta(a+heads,b+tails)))[1];
                    #log(factorial(a+heads-1) * factorial(b+tails-1) / factorial(a+heads+b+tails-1)))[1];
    end
    return posterior
end

#%%

println("\nd) Plotting the likelihood function for a series")
θ = 0:0.001:1;

#%%
# Flat prior:
a = 1;
b = 1;

figure()
subplot(1,2,1)
plot(θ,exp.(Log_Beta.(θ,a,b)), label="beta(θ,$a,$b)")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_Beta.(θ,a,b), label="log(beta(θ,$a,$b))")
grid()
legend()
plt.show()

# %%
y = [1]
a = 1
b = 1

figure()
subplot(1,2,1)
plot(θ,exp.(Log_posterior(θ,a,b,y)), label="posterior for $y")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_posterior(θ,a,b,y), label="log-posterior for $y")
grid()
legend()
plt.show()

# Alternatively directly by Beta:
heads = sum(y)
tails = length(y) - heads

figure()
subplot(1,2,1)
plot(θ,exp.(Log_Beta.(θ,a+heads,b+tails)), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_Beta.(θ,a+heads,b+tails), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
plt.show()

# %%
y = [1 1]
a = 1
b = 1

figure()
subplot(1,2,1)
plot(θ,exp.(Log_posterior(θ,a,b,y)), label="posterior for $y")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_posterior(θ,a,b,y), label="log-posterior for $y")
grid()
legend()
plt.show()

# Alternatively directly by Beta:
heads = sum(y)
tails = length(y) - heads

figure()
subplot(1,2,1)
plot(θ,exp.(Log_Beta.(θ,a+heads,b+tails)), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_Beta.(θ,a+heads,b+tails), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
plt.show()

# %%
y = [1 1 0 1]
a = 1
b = 1

figure()
subplot(1,2,1)
plot(θ,exp.(Log_posterior(θ,a,b,y)), label="posterior for $y")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_posterior(θ,a,b,y), label="log-posterior for $y")
grid()
legend()
plt.show()

# Alternatively directly by Beta:
heads = sum(y)
tails = length(y) - heads

figure()
subplot(1,2,1)
plot(θ,exp.(Log_Beta.(θ,a+heads,b+tails)), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
subplot(1,2,2)
plot(θ,Log_Beta.(θ,a+heads,b+tails), label="beta(θ,$(a+heads),$(b+tails))")
grid()
legend()
plt.show()

#%%

println("The differences are a scaling in the non-log curves and thus
a shifting in the log-curves.
The differences arise because before we only used the LIKELIHOOD function
before, but now we use a normalised POSTERIOR.")

# %%
### iii) More informative prior and reproduce Fig 6.4

# First Column
θ = 0:.01:1
a = 100
b = 100

y = [ones(1,17) zeros(1,3)]

subplot(3,1,1)
plot(θ,exp.(Log_Beta.(θ,a,b)))
grid()
subplot(3,1,2)
plot(θ,exp.(log_likelihood_fct(θ,y)))
grid()
subplot(3,1,3)
plot(θ,exp.(Log_posterior(θ,a,b,y)))
grid()
plt.show()
#%%

# Second Column
θ = 0:.01:1
a = 18.25
b = 6.75

y = [ones(1,17) zeros(1,3)]

subplot(3,1,1)
plot(θ,exp.(Log_Beta.(θ,a,b)))
grid()
subplot(3,1,2)
plot(θ,exp.(log_likelihood_fct(θ,y)))
grid()
subplot(3,1,3)
plot(θ,exp.(Log_posterior(θ,a,b,y)))
grid()
plt.show()

#%%

# Third Column
θ = 0:.001:1
a = 1
b = 1

y = [ones(1,17) zeros(1,3)]

subplot(3,1,1)
plot(θ,exp.(Log_Beta.(θ,a,b)))
grid()
subplot(3,1,2)
plot(θ,exp.(log_likelihood_fct(θ,y)))
grid()
subplot(3,1,3)
plot(θ,exp.(Log_posterior(θ,a,b,y)))
grid()
plt.show()


text:
#=
Here I show the exact correspondence
my spiking implementation and code by
Litwin-Kumar / Auguste Schulz
Checking I<-I and  I<-E connections
For one I neuron connected to one E neuron
=#

push!(LOAD_PATH, abspath(@__DIR__,".."))

using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs

function onesparsemat(w::Real)
  return sparse(fill(w,(1,1)))
end

##

# this is a bare-bone version of the Schulz code
function runsimulation_static_simplified(Ne::Integer,Ni::Integer,
  weights::Array{Float64,2},spiketimes::Array{<:Integer,2},
  v_start::Vector{Float64},input_constant::Vector{Float64}; dt = 0.1, T = 2000)

		# Naming convention
		# e corresponds to excitatory
		# i corresponds to inhibitory
		# x corresponds to external

		#membrane dynamics
		taue = 30 #e membrane time constant ms
		taui = 20 #i membrane time constant ms
		vreste = -60 #e resting potential mV
		vresti = -70 #i resting potential mV
		vpeak = 20 #cutoff for voltage.  when crossed, record a spike and reset mV
		eifslope = 2 #eif slope parameter mV
		C = 300 #capacitance pF
		erev = 0 #e synapse reversal potential mV
		irev = -75 #i synapse reversal potntial mV
		vth0 = -52 #initial spike voltage threshold mV
		thrchange = false # can be switched off to have vth constant at vth0
		ath = 10 #increase in threshold post spike mV
		tauth = 30 #threshold decay timescale ms
		vresete = vreste #reset potential mV
    vreseti = vresti
		taurefrac = 10 #absolute refractory period ms

		# total number of neurons
		Ncells = Ne+Ni

    # constant input
    in_const = input_constant

		# synaptic kernel
		tauerise = 1 #e synapse rise time
		tauedecay = 6 #e synapse decay time
		tauirise = .5 #i synapse rise time
		tauidecay = 2 #i synapse decay time

		# Initialisation of spike arrays ------------------------------
		# spike count and spike times
		# Scope problems - changing arrays are therefore defined here directly

		totalspikes = zeros(Int,Ncells)
		totsp::Int64 = 0; # total numberof spikes
		spmax::Int64 = size(spiketimes,1); # maximum recorded spikes

		# further arrays for storing inputs and synaptic integration
		forwardInputsE = zeros(Float64,Ncells) #sum of all incoming weights from excitatpry inputs both external and EE and IE
		forwardInputsI = zeros(Float64,Ncells) #sum of all incoming weights from inhibitory inputs both II and EI
		forwardInputsEPrev = zeros(Float64,Ncells) #as above, for previous timestep
		forwardInputsIPrev = zeros(Float64,Ncells)

		xerise = zeros(Float64,Ncells) #auxiliary variables for E/I currents (difference of exponentials)
		xedecay = zeros(Float64,Ncells)
		xirise = zeros(Float64,Ncells)
		xidecay = zeros(Float64,Ncells)

		v = zeros(Float64,Ncells) #membrane voltage

		# initialisation of membrane potentials and poisson inputs
    copy!(v,v_start)
		# for cc = 1:Ncells
		# 	v[cc] = vreset # + (vth0-vreset)*rand()
		# end

		vth = vth0*ones(Float64,Ncells) #adaptive threshold
		lastSpike = -100*ones(Ncells) #last time the neuron spiked

		# ---------------------------------- set up storing avg weights -----------------

		# total number of simulation steps, steps when to normalise, and save
		Nsteps = round(Int,T/dt)

		# true time
		t::Float64 = 0.0

		# bool counter if a neuron has just had a spike
		spiked = zeros(Bool,Ncells)
  v_out = Matrix{Float64}(undef,Ncells,Nsteps)
	@time	for tt = 1:Nsteps

				if mod(tt,Nsteps/100) == 1  #print percent complete
					print("\r",round(Int,100*tt/Nsteps))
				end

				forwardInputsE[:] .= 0.
				forwardInputsI[:] .= 0.
				t = dt*tt

				fill!(spiked,zero(Bool)) # reset spike bool without new memory allocation

				for cc = 1:Ncells
					xerise[cc] += -dt*xerise[cc]/tauerise + forwardInputsEPrev[cc]
					xedecay[cc] += -dt*xedecay[cc]/tauedecay + forwardInputsEPrev[cc]
					xirise[cc] += -dt*xirise[cc]/tauirise + forwardInputsIPrev[cc]
					xidecay[cc] += -dt*xidecay[cc]/tauidecay + forwardInputsIPrev[cc]

					if cc <= Ne # excitatory
						if thrchange
						vth[cc] += dt*(vth0 - vth[cc])/tauth;
						end
					end

					if t > (lastSpike[cc] + taurefrac) #not in refractory period
						# update membrane voltage

						ge = (xedecay[cc] - xerise[cc])/(tauedecay - tauerise);
						gi = (xidecay[cc] - xirise[cc])/(tauidecay - tauirise);

						if cc <= Ne #excitatory neuron (eif), has adaptation
							dv = (vreste - v[cc] + eifslope*exp((v[cc]-vth[cc])/eifslope))/taue +
                  ge*(erev-v[cc])/C + gi*(irev-v[cc])/C +
                  in_const[cc]/taue
							v[cc] += dt*dv
							if v[cc] > vpeak
								spiked[cc] = true
							end
						else
							dv = (vresti - v[cc])/taui +
                 ge*(erev-v[cc])/C + gi*(irev-v[cc])/C +
                  in_const[cc]/taui
							v[cc] += dt*dv
							if v[cc] > vth0
								spiked[cc] = true
							end
						end

						if spiked[cc] #spike occurred
							spiked[cc] = true;
							v[cc] = cc<=Ne ? vresete : vreseti;
							lastSpike[cc] = t;
							totalspikes[cc] += 1;
							totsp += 1;
							if totsp < spmax
								spiketimes[totsp,1] = tt; # time index as a sparse way to save spiketimes
								spiketimes[totsp,2] = cc; # cell id
							elseif totsp == spmax
								spiketimes[totsp,1] = tt; # time index
								spiketimes[totsp,2] = cc; # cell id

								totsp = 0 # reset counter total number of spikes
							end

							if cc <= Ne && thrchange # only change for excitatory cells and when thrchange == true
								vth[cc] = vth0 + ath;
							end

							#loop over synaptic projections
							for dd = 1:Ncells # postsynaptic cells dd  - cc presynaptic cells
								if cc <= Ne #excitatory synapse
									forwardInputsE[dd] += weights[cc,dd];
								else #inhibitory synapse
									forwardInputsI[dd] += weights[cc,dd];
								end
							end

						end #end if(spiked)
					end #end if(not refractory)
        v_out[cc,tt] = v[cc]
			end #end loop over cells

			forwardInputsEPrev = copy(forwardInputsE)
			forwardInputsIPrev = copy(forwardInputsI)

		end # tt loop over time
	print("\r")

	println("simulation finished")

	return totalspikes,v_out

end

## two I neurons
const Ttot=1.0

v_start = let n=2,
  (vmin,vmax) = (-60.0,-50.0)
  rand(Uniform(vmin,vmax),n)
end

vtest_ii = let dt=0.1,
  _w = 231.0,
  input = [(70.0 - 50.0) ,0.0],
  T=round(Integer,Ttot*1E3),
  Ne=0,Ni=2,weights=[0.0 _w ;0.0 0.0];
  spiketimes=Matrix{Int32}(undef,1000,2)
  runsimulation_static_simplified(Ne,Ni,weights,spiketimes,
  v_start, input; dt=dt, T=T)[2]
end

_ = let plt = plot()
  ts = (1:size(vtest_ii,2)).*0.1E-3
  plot!(plt,ts,vtest_ii[1,:];linewidth=2)
  plot!(plt,ts,vtest_ii[2,:];linewidth=2)
end

##
vtest_ie = let dt=0.1,
  _w = 300.0,
  input = [(60.0 - 50.0) , 0.0],
  T=round(Integer,Ttot*1E3),
  Ne=1,Ni=1,weights=[0.0 _w ;0.0 0.0];
  spiketimes=Matrix{Int32}(undef,1000,2)
  runsimulation_static_simplified(Ne,Ni,weights,spiketimes,
  v_start, input; dt=dt, T=T)[2]
end

_ = let plt = plot()
  ts = (1:size(vtest_ie,2)).*0.1E-3
  plot!(plt,ts,vtest_ie[1,:];linewidth=2)
  plot!(plt,ts,vtest_ie[2,:];linewidth=2)
end

##
# time in seconds, voltage in mV
dt = 0.1E-3
myτe = 30E-3 # seconds
myτi = 20E-3 # seconds
vreste = -60.0
vresti = -70.0
τrefr= 10E-3 # refractoriness
vth_e = 20.   # mV
vthexp = -52.0 # actual threshold for spike-generation
vth_i = vthexp
eifslope = 2.0
v_reset_e = vreste
v_reset_i = vresti
v_rev_e = 0.0
v_rev_i = -75.0
v_leak_e = vreste
v_leak_i = vresti
Cap = 300.0 #capacitance mF

in_const_e = (-50.0 - vreste)*Cap/myτe
in_const_i = (-50.0 - vresti)*Cap/myτi

# synaptic kernel
tauerise = 1E-3 #e synapse rise time
tauedecay = 6E-3 #e synapse decay time
#taueplus,taueminus  = tauerise, tauedecay
taueplus,taueminus  = tauedecay, tauerise
tauirise = 0.5E-3 #i synapse rise time
tauidecay = 2E-3 #i synapse decay time
tauiplus,tauiminus = tauidecay,tauirise

##
nt_e = let sker = S.SKExpDiff(taueplus,taueminus)
  sgen = S.SpikeGenEIF(vthexp,eifslope)
  S.NTLIFConductance(sker,sgen,myτe,Cap,
   vth_e,v_reset_e,vreste,τrefr,v_rev_e)
end
nt_i = let sker = S.SKExpDiff(tauiplus,tauiminus)
  sgen = S.SpikeGenNone()
  S.NTLIFConductance(sker,sgen,myτi,Cap,
   vth_i,v_reset_i,vresti,τrefr,v_rev_i)
end
ps_e = S.PSLIFConductance(nt_e,1)
ps_i1 = S.PSLIFConductance(nt_i,1)
ps_i2 = S.PSLIFConductance(nt_i,1)

# one static input
ps_input1e = S.PSSimpleInput(S.InputSimpleOffset(in_const_e))
ps_input1i = S.PSSimpleInput(S.InputSimpleOffset(in_const_i))

# I to E connection
conn_ii = let  _w = 231.0
  wmat = sparse(fill(_w,(1,1)))
  S.ConnGeneralIF2(wmat)
end

pop_i1 = S.Population(ps_i1,(S.FakeConnection(),ps_input1i))
pop_i2 = S.Population(ps_i2,(conn_ii,ps_i1))

##
# that's it, let's make the network
myntw = S.RecurrentNetwork(dt,pop_i1,pop_i2)

# record spiketimes and internal potential
krec = 1
rec_state_i1 = S.RecStateNow(ps_i1,krec,dt,Ttot)
rec_state_i2 = S.RecStateNow(ps_i2,krec,dt,Ttot)
rec_spikes_i1 = S.RecSpikes(ps_i1,100.0,Ttot)
rec_spikes_i2 = S.RecSpikes(ps_i2,100.0,Ttot)

## Run

times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_i1,rec_spikes_i1])
S.reset!.([rec_state_i2,rec_spikes_i2])
S.reset!.([ps_e,ps_i1,ps_i2])
S.reset!(conn_ii)
# initial conditions
ps_i1.state_now .= v_start[1]
ps_i2.state_now .= v_start[2]

@time begin
  for (k,t) in enumerate(times)
    rec_state_i1(t,k,myntw)
    rec_state_i2(t,k,myntw)
    rec_spikes_i1(t,k,myntw)
    rec_spikes_i2(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
end

S.add_fake_spikes!(1.0vth_i,rec_state_i1,rec_spikes_i1)
S.add_fake_spikes!(1.0vth_i,rec_state_i2,rec_spikes_i2)
##

_ = let plt=plot()
  plot!(plt,rec_state_i1.times,rec_state_i1.state_now[1,:];linewidth=2,leg=false)
  ts = (1:size(vtest_ii,2)).*0.1E-3
  plot!(plt,ts,vtest_ii[1,:];linewidth=2,linestyle=:dash)
end

_ = let plt=plot()
  plot!(plt,rec_state_i2.times,rec_state_i2.state_now[1,:];linewidth=2,leg=false)
  ts = (1:size(vtest_ii,2)).*0.1E-3
  plot!(plt,ts,vtest_ii[2,:];linewidth=2,linestyle=:dash)
end

##

# I to E connection
conn_ie = let  _w = 300.0
  wmat = sparse(fill(_w,(1,1)))
  S.ConnGeneralIF2(wmat)
end

ps_e = S.PSLIFConductance(nt_e,1)
ps_i = ps_i1

pop_e = S.Population(ps_e,(S.FakeConnection(),ps_input1e))
pop_i = S.Population(ps_i,(conn_ie,ps_e))

myntw = S.RecurrentNetwork(dt,pop_e,pop_i)

krec = 1
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot)
rec_state_i = S.RecStateNow(ps_i,krec,dt,Ttot)
rec_spikes_e = S.RecSpikes(ps_e,100.0,Ttot)
rec_spikes_i = S.RecSpikes(ps_i,100.0,Ttot)

##
times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_state_e,rec_spikes_e])
S.reset!.([rec_state_i,rec_spikes_i])
S.reset!.([ps_e,ps_i])
S.reset!(conn_ie)
# initial conditions
ps_e.state_now .= v_start[1]
ps_i.state_now .= v_start[2]

@time begin
  for (k,t) in enumerate(times)
    rec_state_e(t,k,myntw)
    rec_state_i(t,k,myntw)
    rec_spikes_e(t,k,myntw)
    rec_spikes_i(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
end

S.add_fake_spikes!(1.0vth_i,rec_state_i1,rec_spikes_i1)
S.add_fake_spikes!(1.0vth_i,rec_state_i2,rec_spikes_i2)

##
_ = let plt=plot()
  plot!(plt,rec_state_e.times,rec_state_e.state_now[1,:];linewidth=2,leg=false)
  ts = (1:size(vtest_ie,2)).*0.1E-3
  plot!(plt,ts,vtest_ie[1,:];linewidth=2,linestyle=:dash)
end

_ = let plt=plot()
  plot!(plt,rec_state_i.times,rec_state_i.state_now[1,:];linewidth=2,leg=false)
  ts = (1:size(vtest_ie,2)).*0.1E-3
  plot!(plt,ts,vtest_ie[2,:];linewidth=2,linestyle=:dash)
end


text:
using CpuId

register_size = simdbytes()
register_count = cpufeature(CpuId.AVX512F) ? 32 : 16
register_capacity = register_size * register_count
cache_size = CpuId.cachesize()
num_cores = CpuId.cpucores()

cpu_info_string = """
const REGISTER_SIZE = $register_size
const REGISTER_COUNT = $register_count
const FP256 = $(cpufeature(CpuId.FP256)) # Is AVX2 fast?
const CACHELINE_SIZE = $(cachelinesize())
const CACHE_SIZE = $cache_size
const NUM_CORES = $num_cores
const FMA3 = $(cpufeature(CpuId.FMA3))
const AVX2 = $(cpufeature(CpuId.AVX2))
const AVX512F = $(cpufeature(CpuId.AVX512F))
const AVX512ER = $(cpufeature(CpuId.AVX512ER))
const AVX512PF = $(cpufeature(CpuId.AVX512PF))
const AVX512VL = $(cpufeature(CpuId.AVX512VL))
const AVX512BW = $(cpufeature(CpuId.AVX512BW))
const AVX512DQ = $(cpufeature(CpuId.AVX512DQ))
const AVX512CD = $(cpufeature(CpuId.AVX512CD))
const SIMD_NATIVE_INTEGERS = $(cpufeature(CpuId.AVX2))
"""


text:
mutable struct VAEvamp{T<:AbstractFloat,V<:Val}
	q  # encoder (inference modul)
    g  # decoder (generator)
    centroids::AbstractMatrix{T}
    labels::Vector{Int}
	variant::V
end

VAEvamp(q, g, c, l, s::Symbol = :unit) = VAEvamp(q, g, c, l, Val(s))

Flux.@treelike(VAEvamp)

function add_labeled_centroid!(m::VAEvamp, c, l)
    m.centroids = hcat(m.centroids, c)
    push!(m.labels, l)
end

function sample_vamp(m::VAEvamp{T,V}, n) where {T,V<:Val{:unit}}
    components = rand(1:size(m.centroids, 2), n)
	μz = m.q(m.centroids[:, components])
	return gaussiansample(μz, 1f0)
end

function sample_vamp(m::VAEvamp{T,V}, n0, n1) where {T,V<:Val{:unit}}
    normal_ids = findall(m.labels .== 0)
    anomalous_ids = findall(m.labels .== 1)
    components_0 = rand(1:length(normal_ids), n0)
    components_1 = rand(1:length(anomalous_ids), n1)
    μz0 = m.q(m.centroids[:, components_0])
    μz1 = m.q(m.centroids[:, components_1])
	return gaussiansample(μz0, 1f0), gaussiansample(μz1, 1f0)
end

function wloss(m::VAEvamp{T,V}, x, d, β) where {T,V<:Val{:unit}}
    μz = m.q(x)
    μx = m.g(μz)
    prior = sample_vamp(m, size(x, 2))
    Ω = d(μz, prior)
	-mean(log_normal(x, μx)) + β * Ω
end

function wloss(m::VAEvamp{T,V}, x, l, d, β, α) where {T,V<:Val{:unit}}
    normal_ids = findall(l .== 0)
    anomalous_ids = findall(l .== 1)
    if length(anomalous_ids) == 0
        return wloss(m, x, d, β)
    end

    μz = m.q(x)
    μx = m.g(μz)
    prior_0, prior_1 = sample_vamp(m, length(normal_ids), length(anomalous_ids))
    Ω0 = d(μz[:, normal_ids], prior_0)
    Ω1 = d(μz[:, anomalous_ids], prior_1)
    -mean(log_normal(x, μx)) + β * ((1f0 - α) * Ω0 + α * Ω1)
end

function log_pxexpectedz(m::VAEvamp{T, V}, x) where {T, V<:Val{:unit}}
	μz = m.q(x)
	log_normal(x, m.g(μz))
end

function pz(m::VAEvamp{T,V}, x) where {T,V<:Val{:unit}}
	if size(x, 2) > 1
		xs = [x[:, i] for i in 1:size(x, 2)]
		return hcat(map(x -> collect(pz_singleinstance(m, x)), xs)...)
	else
		return pz_singleinstance(m, x)
	end
end

function pz_singleinstance(m::VAEvamp{T,V}, x) where {T,V<:Val{:unit}}
	@assert size(x, 2) == 1
	μz = m.q(x).data
    centroids_0 = m.q(m.centroids[:, m.labels .== 0]).data
    centroids_1 = m.q(m.centroids[:, m.labels .== 1]).data
	p0 = 0
    for i in 1:size(centroids_0, 2)
		p0 += exp(log_normal(μz, centroids_0[:, i])[1])
	end
    p0 /= size(centroids_0, 2)
    p1 = 0
	for i in 1:size(centroids_1, 2)
		p1 += exp(log_normal(μz, centroids_1[:, i])[1])
	end
    p1 /= size(centroids_1, 2)
    return p0, p1
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
module ExtendedVariantModeColorGradeTrigger

using ..Ahorn, Maple

@mapdef Trigger "ExtendedVariantMode/ColorGradeTrigger" ColorGradeTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	colorGrade::String="none", revertOnDeath::Bool=true, onlyOnce::Bool=false)

const colorGrades = String["none", "oldsite", "panicattack", "templevoid", "reflection", "credits", "cold", "hot", "feelingdown", "golden",
	"max480/extendedvariants/celsius/tetris",
	"max480/extendedvariants/greyscale", "max480/extendedvariants/sepia", "max480/extendedvariants/inverted",
	"max480/extendedvariants/rgbshift1", "max480/extendedvariants/rgbshift2", "max480/extendedvariants/hollys_randomnoise"]

const placements = Ahorn.PlacementDict(
	"Color Grade Trigger (Extended Variant Mode)" => Ahorn.EntityPlacement(
		ColorGradeTrigger,
		"rectangle"
	)
)

Ahorn.editingOptions(trigger::ColorGradeTrigger) = Dict{String, Any}(
	"colorGrade" => colorGrades
)

end


text:
@testset "Robots: hopper" begin
    q0 = nominal_configuration(hopper)

    @test norm(RoboDojo.kinematics_body(hopper, q0) - q0[1:2]) < 1.0e-8
    @test norm(RoboDojo.kinematics_body_jacobian(hopper, q0) - [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0]) < 1.0e-8
    @test norm(RoboDojo.kinematics_foot(hopper, q0) - [0.0; 0.5]) < 1.0e-8
    @test norm(RoboDojo.kinematics_foot_jacobian(hopper, q0) - [1.0 0.0 (q0[4] * cos(q0[3])) sin(q0[3]);
                                                                 0.0 1.0 (q0[4] * sin(q0[3])) (-1.0 * cos(q0[3]))]) < 1.0e-8
    @test lagrangian(hopper, zeros(hopper.nq), zeros(hopper.nq)) ≈ 0.0
    @test norm(RoboDojo.input_jacobian(hopper, q0) - transpose([0.0 -sin(q0[3]); 0.0 cos(q0[3]); 1.0 0.0; 0.0 1.0])) < 1.0e-8
    @test norm(RoboDojo.contact_jacobian(hopper, q0) - [1.0 0.0 0.0 0.0;
                                                         0.0 1.0 0.0 0.0;
                                                         1.0 0.0 (q0[4] * cos(q0[3])) sin(q0[3]);
                                                         0.0 1.0 (q0[4] * sin(q0[3])) (-1.0 * cos(q0[3]));
                                                         0.0 0.0 0.0 1.0;
                                                         0.0 0.0 0.0 -1.0]) < 1.0e-8

    # visualizer
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, hopper, [q0], Δt=0.1);
end


text:
# consoleutil.jl (v0.1)

export cin

function putmenu()

	f = open("consolemenu.txt");

	println(read(f, String));

	close(f);

	choice::Int = 0;

end

# cin reads from STDIN. It reads the DataType t.
# lines specifies how many lines one wishes to read. One line is standard.
# -1 lines reads until "__stop" is given.
# prompt is shown when cin is called.
function cin(t::DataType, lines::Int = 1, prompt::String = "")

	line::String = "";
	ret::Array{t, 1} = [];

	print(prompt);

	if lines == -1

		println("Send \"__stop\" to end reading.");
		stop::Bool = false;

		while !stop
			line = readline();

			if line == "__stop"
				stop = true;
				continue;
			else
				try
					push!(ret, parse(t, line));
				catch err
					if err isa MethodError && t == String
						push!(ret, line);
					else
						println(err);
					end
				end
			end
		end

	elseif lines == 1

		if t == String
			return readline();
		else
			while true
				try
					return parse(t, readline());
				catch err
					println(err);
				end
			end
		end

	else
		if t == String
			for i = 1:lines
				push!(ret, readline());
			end
		else
			for i = 1:lines
				push!(ret, cin(t));
			end
		end

	end # end if lines == -1

	return ret;
end


text:
import OptimizingIR
const OIR = OptimizingIR
using Test

@testset "Graph" begin
    include("test_graph.jl")
end

foreign_fun_binary(a, b) = a^2 + b
foreign_fun_3args(a, b, c) = a^3 + b^2 + c
foreign_fun_4args(a, b, c, d) = a^4 + b^3 + c^2 + d
foreign_fun_5args(a, b, c, d, e) = a^5 + b^4 + c^3 + d^2 + e
foreign_fun_6args(a, b, c, d, e, f) = a^6 + b^5 + c^4 + d^3 + e^2 + f
foreign_fun_7args(a, b, c, d, e, f, g) = a + b + c + d + e + f + g
julia_basic_block_test_function(x::Number) = (((-((10.0 * 2.0 + x) / 1.0) + (x + 10.0 * 2.0) + 1.0) * 1.0 / 2.0) + (0.0 * x) + 1.0) * 1.0
julia_native_test_function(x::Number, y::Number, z::Number) = x^3 + y^2 + z

function push_into_two_args(num::Integer, a::Vector, b::Vector)
    push!(a, num)
    push!(b, num)
end

const OP_SUM = OIR.Op(+, pure=true, commutative=true, hasleftidentity=true, hasrightidentity=true, identity_element=0)
const OP_SUB = OIR.Op(-, pure=true, hasrightidentity=true, identity_element=0)
const OP_MUL = OIR.Op(*, pure=true, commutative=true, hasleftidentity=true, hasrightidentity=true, identity_element=1)
const OP_DIV = OIR.Op(/, pure=true, hasrightidentity=true, identity_element=1)
const OP_POW = OIR.Op(^, pure=true, hasrightidentity=true, identity_element=1)
const OP_ABS = OIR.Op(abs, pure=true)
const OP_FOREIGN_FUN_BINARY = OIR.Op(foreign_fun_binary, pure=true)
const OP_FOREIGN_FUN_3ARGS = OIR.Op(foreign_fun_3args, pure=true)
const OP_FOREIGN_FUN_4ARGS = OIR.Op(foreign_fun_4args, pure=true)
const OP_FOREIGN_FUN_5ARGS = OIR.Op(foreign_fun_5args, pure=true)
const OP_FOREIGN_FUN_6ARGS = OIR.Op(foreign_fun_6args, pure=true)
const OP_FOREIGN_FUN_7ARGS = OIR.Op(foreign_fun_7args, pure=true)
const OP_ZEROS = OIR.Op(zeros)

const OP_GETINDEX = OIR.Op(Base.getindex, pure=true)
const OP_SETINDEX = OIR.Op(Base.setindex!, pure=false, mutable_arg=1)

const OP_PUSH_INTO_TWO_ARGS = OIR.Op(push_into_two_args, mutable_arg=(2, 3))

@testset "LookupTable" begin
    table = OIR.LookupTable{Int}()
    @test isempty(table)
    @test OIR.addentry!(table, 10) == 1
    @test !isempty(table)
    @test OIR.addentry!(table, 20) == 2
    @test length(table) == 2
    @test OIR.addentry!(table, 10) == 1
    @test OIR.addentry!(table, 20) == 2
    @test length(table) == 2
    @test 10 ∈ table
    @test 20 ∈ table
    @test 1 ∉ table
    @test lastindex(table) == 2
    @test OIR.indexof(table, 10) == 1
    @test OIR.indexof(table, 20) == 2
    @test table[1] == 10
    @test table[2] == 20

    for (i, item) in enumerate(table)
        if i == 1
            @test item == 10
        elseif i == 2
            @test item == 20
        else
            @test false
        end
    end

    for item in table
        @test item == 10 || item == 20
    end

    @test collect(table) == [10, 20]
    @test filter(i -> i > 10, table) == [20]

    @testset "Vector constructor" begin
        v = collect(1:3)
        tb = OIR.LookupTable(v)

        for (i, v) in enumerate(v)
            @test OIR.indexof(tb, v) == i
        end
    end
end

@testset "OptimizationRule" begin
    @test OIR.is_pure(OP_SUM)
    @test !OIR.is_pure(OP_ZEROS)
    @test !OIR.is_impure(OP_SUM)
    @test OIR.is_impure(OP_ZEROS)
    @test !OIR.has_left_identity_property(OP_ZEROS)
    @test OIR.has_left_identity_property(OP_SUM)
    @test !OIR.has_right_identity_property(OP_ZEROS)
    @test OIR.has_right_identity_property(OP_SUM)
    @test OIR.has_identity_property(OP_SUB)
    @test !OIR.has_identity_property(OP_ZEROS)
    @test OIR.has_identity_element(OP_SUM)
end

@testset "AbstractValue" begin
    @test OIR.is_immutable(OIR.NullPointer())
    @test OIR.is_immutable(OIR.SSAValue(1))
    @test OIR.is_immutable(OIR.Const(1.0))
    @test OIR.is_immutable(OIR.ImmutableVariable(:sym))
    @test !OIR.is_immutable(OIR.MutableVariable(:sym))
    @test OIR.is_immutable(OIR.Variable{OIR.Immutable})

    @test !OIR.is_mutable(OIR.NullPointer())
    @test !OIR.is_mutable(OIR.SSAValue(1))
    @test !OIR.is_mutable(OIR.Const(1.0))
    @test !OIR.is_mutable(OIR.ImmutableVariable(:sym))
    @test OIR.is_mutable(OIR.MutableVariable(:sym))
end

@testset "call" begin

    @testset "Call Foreign Fun Binary" begin
        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(10.0)
        arg4 = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_BINARY, var_x, arg2))
        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg4)

        input = 20.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ 20.0^2 + 10.0
    end

    @testset "Call3Args" begin

        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(10.0)
        arg3 = OIR.constant(2.0)
        arg4 = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_3ARGS, var_x, arg2, arg3))

        var_x_2nd_instance = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x_2nd_instance)
        arg6 = OIR.constant(10.0)
        arg7 = OIR.constant(2.0)
        arg8 = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_3ARGS, var_x_2nd_instance, arg6, arg7))
        arg9 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg4, arg8))
        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg9)

        println(bb)
        println(OIR.CompiledBasicBlock(bb))
        @test length(bb.instructions) == 3

        input = 20.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ 2 * (input^3 + 10^2 + 2)

        fnative = OIR.compile(OIR.Native, bb)
        @test fnative(input) ≈ 2 * (input^3 + 10^2 + 2)
    end

    @testset "Call4Args" begin

        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(2.0)
        arg3 = OIR.constant(3.0)
        arg4 = OIR.constant(4.0)
        arg_result = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_4ARGS, var_x, arg2, arg3, arg4))

        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg_result)

        println(bb)
        println(OIR.CompiledBasicBlock(bb))
        @test length(bb.instructions) == 2

        input = 5.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ foreign_fun_4args(input, 2.0, 3.0, 4.0)

        fnative = OIR.compile(OIR.Native, bb)
        @test fnative(input) ≈ foreign_fun_4args(input, 2.0, 3.0, 4.0)
    end

    @testset "Call5Args" begin

        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(2.0)
        arg3 = OIR.constant(3.0)
        arg4 = OIR.constant(4.0)
        arg5 = OIR.constant(5.0)
        arg_result = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_5ARGS, var_x, arg2, arg3, arg4, arg5))

        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg_result)

        println(bb)
        println(OIR.CompiledBasicBlock(bb))
        @test length(bb.instructions) == 2

        input = 5.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ foreign_fun_5args(input, 2.0, 3.0, 4.0, 5.0)

        fnative = OIR.compile(OIR.Native, bb)
        @test fnative(input) ≈ foreign_fun_5args(input, 2.0, 3.0, 4.0, 5.0)
    end

    @testset "Call6Args" begin

        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(2.0)
        arg3 = OIR.constant(3.0)
        arg4 = OIR.constant(4.0)
        arg5 = OIR.constant(5.0)
        arg6 = OIR.constant(6.0)
        arg_result = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_6ARGS, var_x, arg2, arg3, arg4, arg5, arg6))

        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg_result)

        println(bb)
        println(OIR.CompiledBasicBlock(bb))
        @test length(bb.instructions) == 2

        input = 5.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ foreign_fun_6args(input, 2.0, 3.0, 4.0, 5.0, 6.0)

        fnative = OIR.compile(OIR.Native, bb)
        @test fnative(input) ≈ foreign_fun_6args(input, 2.0, 3.0, 4.0, 5.0, 6.0)
    end

    @testset "Call7Args" begin

        bb = OIR.BasicBlock()
        var_x = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, var_x)
        arg2 = OIR.constant(2.0)
        arg3 = OIR.constant(3.0)
        arg4 = OIR.constant(4.0)
        arg5 = OIR.constant(5.0)
        arg6 = OIR.constant(6.0)
        arg7 = OIR.constant(7.0)

        # currently supports up to 6 args
        @test_throws MethodError arg_result = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_7ARGS, var_x, arg2, arg3, arg4, arg5, arg6, arg7))
#=
        var_output = OIR.ImmutableVariable(:output)
        OIR.addoutput!(bb, var_output)
        OIR.assign!(bb, var_output, arg_result)

        println(bb)
        println(OIR.CompiledBasicBlock(bb))
        @test length(bb.instructions) == 2

        input = 5.0
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(input) ≈ foreign_fun_7args(input, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)

        fnative = OIR.compile(OIR.Native, bb)
        @test fnative(input) ≈ foreign_fun_7args(input, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
=#
    end

    @testset "ImpureCall" begin
        bb = OIR.BasicBlock()
        arg1 = OIR.constant(Float64)
        arg2 = OIR.constant(3)
        arg3 = OIR.addinstruction!(bb, OIR.call(OP_ZEROS, arg1, arg2))
        arg4 = OIR.addinstruction!(bb, OIR.call(OP_ZEROS, arg1, arg2))
        @test length(bb.instructions) == 2
    end
end

@testset "OP_GETINDEX" begin
    bb = OIR.BasicBlock()
    var_x = OIR.ImmutableVariable(:x)
    OIR.addinput!(bb, var_x)
    arg2 = OIR.constant(2)
    arg3 = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_x, arg2)) # reads x[2]
    arg4 = OIR.constant(5.0)
    arg5 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg4, arg3))
    arg6 = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_x, arg2)) # reads x[2]
    arg7 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg6, arg5))
    var_output = OIR.ImmutableVariable(:output)
    OIR.addoutput!(bb, var_output)
    OIR.assign!(bb, var_output, arg7)
    @test length(bb.instructions) == 4

    # println(bb)

    finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
    @test finterpreter([5.0, 6.0, 7.0]) ≈ 6.0 * 5.0 * 6.0

    fnative = OIR.compile(OIR.Native, bb)
    @test fnative([5.0, 6.0, 7.0]) ≈ 6.0 * 5.0 * 6.0
end

@testset "OP_SETINDEX" begin
    bb = OIR.BasicBlock()
    var_x = OIR.ImmutableVariable(:x)
    OIR.addinput!(bb, var_x)
    arg1 = OIR.constant(Float64)
    arg3 = OIR.addinstruction!(bb, OIR.call(OP_ZEROS, arg1, var_x))
    var_vec = OIR.MutableVariable(:vec)
    OIR.assign!(bb, var_vec, arg3)
    arg4 = OIR.constant(1)
    arg_inspect = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_vec, arg4))
    var_inspect1 = OIR.ImmutableVariable(:inspect1)
    OIR.assign!(bb, var_inspect1, arg_inspect)
    arg_input_value = OIR.constant(10)
    arg5 = OIR.addinstruction!(bb, OIR.call(OP_SETINDEX, var_vec, arg_input_value, arg4))
    arg_inspect = OIR.addinstruction!(bb, OIR.call(OP_GETINDEX, var_vec, arg4))
    var_inspect2 = OIR.ImmutableVariable(:inspect2)
    OIR.assign!(bb, var_inspect2, arg_inspect)
    arg6 = OIR.addinstruction!(bb, OIR.call(OP_MUL, OIR.constant(2.0), arg_inspect))
    var_inspect3 = OIR.ImmutableVariable(:inspect3)
    OIR.assign!(bb, var_inspect3, arg6)
    arg7 = OIR.addinstruction!(bb, OIR.call(OP_MUL, OIR.constant(1.0), arg6))
    arg8 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg7, OIR.constant(1.0)))
    arg9 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg8, OIR.constant(0.0)))
    arg10 = OIR.addinstruction!(bb, OIR.call(OP_SUM, OIR.constant(0.0), arg9))
    arg11 = OIR.addinstruction!(bb, OIR.call(OP_SUB, arg10, OIR.constant(0.0)))
    arg12 = OIR.addinstruction!(bb, OIR.call(OP_DIV, arg11, OIR.constant(1.0)))
    var_inspect4 = OIR.ImmutableVariable(:inspect4)
    OIR.assign!(bb, var_inspect4, arg12)
    @test length(bb.instructions) == 10

    OIR.addoutput!(bb, var_vec)
    OIR.addoutput!(bb, var_inspect1)
    OIR.addoutput!(bb, var_inspect2)
    OIR.addoutput!(bb, var_inspect3)
    OIR.addoutput!(bb, var_inspect4)

    println(bb)

    let
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        out_vec, out_inspect1, out_inspect2, out_inspect3, out_inspect4 = finterpreter(3)
        @test out_inspect1 ≈ 0.0
        @test out_inspect2 ≈ 10.0
        @test isa(out_inspect2, Float64)
        @test out_inspect3 ≈ 20.0
        @test out_inspect4 ≈ 20.0
        @test out_vec == [ 10.0, 0.0, 0.0 ]
    end

    let
        fnative = OIR.compile(OIR.Native, bb)
        out_vec, out_inspect1, out_inspect2, out_inspect3, out_inspect4 = fnative(3)
        @test out_inspect1 ≈ 0.0
        @test out_inspect2 ≈ 10.0
        @test isa(out_inspect2, Float64)
        @test out_inspect3 ≈ 20.0
        @test out_inspect4 ≈ 20.0
        @test out_vec == [ 10.0, 0.0, 0.0 ]
    end
end

@testset "Basic Block" begin
    bb = OIR.BasicBlock()
    x = OIR.ImmutableVariable(:x)
    @test !OIR.has_symbol(bb, :x)
    OIR.addinput!(bb, x)
    @test OIR.has_symbol(bb, :x)
    arg1 = OIR.constant(10.0)
    arg2 = OIR.constant(2.0)
    arg3 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg1, arg2))
    arg4 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg3, x))
    arg5 = OIR.constant(1.0)
    arg6 = OIR.addinstruction!(bb, OIR.call(OP_DIV, arg4, arg5))
    arg7 = OIR.addinstruction!(bb, OIR.call(OP_SUB, arg6))
    arg8 = OIR.addinstruction!(bb, OIR.call(OP_SUM, x, arg3))
    arg9 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg8, arg7))
    arg10 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg9, arg5))
    arg11 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg10, arg5))
    arg12 = OIR.addinstruction!(bb, OIR.call(OP_DIV, arg11, arg2))
    arg13 = OIR.constant(0.0)
    arg14 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg13, x))
    arg15 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg14, arg12))
    arg16 = OIR.constant(1.0)
    arg17 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg16, arg15))
    arg18 = OIR.addinstruction!(bb, OIR.call(OP_MUL, arg16, arg17))
    @test !OIR.has_symbol(bb, :output)
    var_output = OIR.MutableVariable(:output)
    OIR.addoutput!(bb, var_output)
    @test OIR.has_symbol(bb, :output)
    OIR.assign!(bb, var_output, arg18)
    @test length(bb.instructions) == 9

    # println(bb)

    let
        input = 10.0
        f = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test f(input) ≈ julia_basic_block_test_function(input)
    end

    let
        input = 10.0
        f = OIR.compile(OIR.Native, bb)
        @test f(input) ≈ julia_basic_block_test_function(input)
    end

    @testset "Multiply by zero" begin
        last_instruction_ssavalue = OIR.SSAValue(lastindex(bb.instructions))
        arg_zero = OIR.constant(0.0)
        arg_result = OIR.addinstruction!(bb, OIR.call(OP_MUL, last_instruction_ssavalue, arg_zero))
        OIR.assign!(bb, var_output, arg_result)

        let
            f = OIR.compile(OIR.BasicBlockInterpreter, bb)
            @test f(10.0) == 0.0
        end

        let
            f = OIR.compile(OIR.Native, bb)
            @test f(10.0) == 0.0
        end
    end
end

@testset "Inputs" begin
    bb = OIR.BasicBlock()
    x = OIR.ImmutableVariable(:x)
    y = OIR.ImmutableVariable(:y)
    z = OIR.ImmutableVariable(:z)
    @test OIR.addinput!(bb, x) == 1
    @test OIR.addinput!(bb, y) == 2
    @test OIR.addinput!(bb, z) == 3
    @test OIR.addinput!(bb, OIR.ImmutableVariable(:y)) == 2
    out = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_3ARGS, x, y, z))
    var_result = OIR.ImmutableVariable(:result)
    @test OIR.addoutput!(bb, var_result) == 1
    @test OIR.addoutput!(bb, OIR.ImmutableVariable(:result)) == 1
    OIR.assign!(bb, var_result, out)

    # cannot be optimized to a constant since it depends on the inputs
    @test length(bb.instructions) == 2

    # println(bb)

    let
        f = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test f(30.0, 20.0, 10.0) ≈ foreign_fun_3args(30.0, 20.0, 10.0)
    end

    let
        f = OIR.compile(OIR.Native, bb)
        @test f(30.0, 20.0, 10.0) ≈ foreign_fun_3args(30.0, 20.0, 10.0)
    end
end

@testset "Passes" begin
    @testset "Constant Propagation" begin
        @testset "1ARG" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(-30.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_ABS, x))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ 30.0
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ 30.0
            end
        end

        @testset "2ARGS" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(30.0)
            y = OIR.constant(20.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_SUM, x, y))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ 50.0
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ 50.0
            end
        end

        @testset "3ARGS" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(30.0)
            y = OIR.constant(20.0)
            z = OIR.constant(10.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_3ARGS, x, y, z))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ foreign_fun_3args(30.0, 20.0, 10.0)
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ foreign_fun_3args(30.0, 20.0, 10.0)
            end
        end

        @testset "4ARGS" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(30.0)
            y = OIR.constant(20.0)
            z = OIR.constant(10.0)
            w = OIR.constant(100.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_4ARGS, x, y, z, w))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ foreign_fun_4args(30.0, 20.0, 10.0, 100.0)
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ foreign_fun_4args(30.0, 20.0, 10.0, 100.0)
            end
        end

        @testset "5ARGS" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(30.0)
            y = OIR.constant(20.0)
            z = OIR.constant(10.0)
            w = OIR.constant(100.0)
            t = OIR.constant(200.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_5ARGS, x, y, z, w, t))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ foreign_fun_5args(30.0, 20.0, 10.0, 100.0, 200.0)
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ foreign_fun_5args(30.0, 20.0, 10.0, 100.0, 200.0)
            end
        end

        @testset "6ARGS" begin
            bb = OIR.BasicBlock()
            x = OIR.constant(30.0)
            y = OIR.constant(20.0)
            z = OIR.constant(10.0)
            w = OIR.constant(100.0)
            t = OIR.constant(200.0)
            p = OIR.constant(300.0)
            out = OIR.addinstruction!(bb, OIR.call(OP_FOREIGN_FUN_6ARGS, x, y, z, w, t, p))
            var_result = OIR.ImmutableVariable(:result)
            OIR.addoutput!(bb, var_result)
            OIR.assign!(bb, var_result, out)
            @test isa(out, OIR.Const)
            @test length(bb.instructions) == 1

            let
                f = OIR.compile(OIR.BasicBlockInterpreter, bb)
                @test f() ≈ foreign_fun_6args(30.0, 20.0, 10.0, 100.0, 200.0, 300.0)
            end

            let
                f = OIR.compile(OIR.Native, bb)
                @test f() ≈ foreign_fun_6args(30.0, 20.0, 10.0, 100.0, 200.0, 300.0)
            end
        end
    end
end

@testset "Variables" begin
    bb = OIR.BasicBlock()
    x = OIR.ImmutableVariable(:x)
    OIR.addinput!(bb, x)
    z = OIR.constant(1.0)
    var_cnst = OIR.ImmutableVariable(:cnst)
    OIR.assign!(bb, var_cnst, z)
    var_slot = OIR.MutableVariable(:slot)
    OIR.assign!(bb, var_slot, z)
    out = OIR.addinstruction!(bb, OIR.call(OP_SUM, var_slot, x))
    OIR.assign!(bb, var_slot, out)
    out = OIR.addinstruction!(bb, OIR.call(OP_SUM, var_slot, x))
    var_output = OIR.ImmutableVariable(:output)
    OIR.assign!(bb, var_output, out)
    @test length(bb.instructions) == 6

    OIR.addoutput!(bb, var_cnst)
    OIR.addoutput!(bb, var_slot)
    OIR.addoutput!(bb, var_output)

    # println(bb)

    input = 10.0

    let
        f = OIR.compile(OIR.BasicBlockInterpreter, bb)
        cnst, slot, output = f(input)
        @test slot == 11.0
        @test output == 1.0 + 10.0 + 10.0
        @test cnst == 1.0
    end

    let
        f = OIR.compile(OIR.BasicBlockInterpreter, OIR.CompiledBasicBlock(bb))
        cnst, slot, output = f(input)
        @test slot == 11.0
        @test output == 1.0 + 10.0 + 10.0
        @test cnst == 1.0
    end

    let
        mem_buff = Vector{Any}(undef, 10_000)
        input_buff = Vector{Any}(undef, 1_000)

        f = OIR.BasicBlockInterpreter(bb, mem_buff, input_buff)
        cnst, slot, output = f(input)
        @test slot == 11.0
        @test output == 1.0 + 10.0 + 10.0
        @test cnst == 1.0

        @test length(mem_buff) == 10_000
        @test length(input_buff) == 1_000
    end

    let
        f = OIR.compile(OIR.Native, bb)
        cnst, slot, output = f(input)
        @test slot == 11.0
        @test output == 1.0 + 10.0 + 10.0
        @test cnst == 1.0
    end

    let
        f = OIR.compile(OIR.Native, OIR.CompiledBasicBlock(bb))
        cnst, slot, output = f(input)
        @test slot == 11.0
        @test output == 1.0 + 10.0 + 10.0
        @test cnst == 1.0
    end
end

@testset "has_symbol" begin
    @testset "input" begin
        bb = OIR.BasicBlock()
        @test !OIR.has_symbol(bb, :x)
        OIR.addinput!(bb, OIR.ImmutableVariable(:x))
        @test OIR.has_symbol(bb, :x)
    end

    @testset "locals" begin
        bb = OIR.BasicBlock()
        @test !OIR.has_symbol(bb, :x)
        OIR.assign!(bb, OIR.ImmutableVariable(:x), OIR.constant(1.0))
        @test OIR.has_symbol(bb, :x)
        @test !OIR.has_symbol(bb, :y)
        OIR.assign!(bb, OIR.MutableVariable(:y), OIR.constant(1.0))
        @test OIR.has_symbol(bb, :y)
    end

    @testset "output mut" begin
        bb = OIR.BasicBlock()
        @test !OIR.has_symbol(bb, :z)
        OIR.addoutput!(bb, OIR.MutableVariable(:z))
        @test OIR.has_symbol(bb, :z)
    end

    @testset "output immut" begin
        bb = OIR.BasicBlock()
        @test !OIR.has_symbol(bb, :z)
        OIR.addoutput!(bb, OIR.ImmutableVariable(:z))
        @test OIR.has_symbol(bb, :z)
    end
end

@testset "gensym" begin
    bb = OIR.BasicBlock()
    sym = OIR.generate_unique_variable_symbol(bb)
    OIR.addinput!(bb, OIR.ImmutableVariable(sym))
    new_sym = OIR.generate_unique_variable_symbol(bb)
    @test new_sym != sym
end

@testset "Native" begin
    bb = OIR.BasicBlock()
    in1 = OIR.ImmutableVariable(:x)
    in2 = OIR.ImmutableVariable(:y)
    in3 = OIR.ImmutableVariable(:z)
    OIR.addinput!(bb, in1)
    OIR.addinput!(bb, in2)
    OIR.addinput!(bb, in3)
    c3 = OIR.constant(3)
    c2 = OIR.constant(2)
    arg1 = OIR.addinstruction!(bb, OIR.call(OP_POW, in1, c3))
    arg2 = OIR.addinstruction!(bb, OIR.call(OP_POW, in2, c2))
    arg3 = in3
    s1 = OIR.addinstruction!(bb, OIR.call(OP_SUM, arg1, arg2))
    s2 = OIR.addinstruction!(bb, OIR.call(OP_SUM, s1, arg3))
    var_result = OIR.ImmutableVariable(:result)
    OIR.assign!(bb, var_result, s2)
    OIR.addoutput!(bb, var_result)

    let
        finterpreter = OIR.compile(OIR.BasicBlockInterpreter, bb)
        @test finterpreter(30.0, 20.0, 10.0) == julia_native_test_function(30.0, 20.0, 10.0)
    end

    let
        f = OIR.compile(OIR.Native, bb)
        @test f(30.0, 20.0, 10.0) == julia_native_test_function(30.0, 20.0, 10.0)
    end
end

@testset "World Age Problem" begin

    function gen_and_run()
        bb = OIR.BasicBlock()
        input_var = OIR.ImmutableVariable(:x)
        OIR.addinput!(bb, input_var)
        c2 = OIR.constant(2)
        arg1 = OIR.addinstruction!(bb, OIR.call(OP_POW, input_var, c2))
        var_result = OIR.ImmutableVariable(:result)
        OIR.assign!(bb, var_result, arg1)
        OIR.addoutput!(bb, var_result)

        f = OIR.compile(OIR.Native, bb)
        @test Base.invokelatest(f, 2) == 4
    end

    gen_and_run()
end

@testset "can't mutate input" begin
    bb = OIR.BasicBlock()
    var_input = OIR.ImmutableVariable(:x)
    OIR.addinput!(bb, var_input)
    cnst_index = OIR.constant(1)
    cnst_val = OIR.constant(2)
    @test_throws AssertionError OIR.addinstruction!(bb, OIR.call(OP_SETINDEX, var_input, cnst_val, cnst_index))
    @test_throws AssertionError OIR.addinstruction!(bb, OIR.call(OP_PUSH_INTO_TWO_ARGS, cnst_val, cnst_val, cnst_val))

    var_mutable_a = OIR.MutableVariable(:a)
    var_mutable_b = OIR.MutableVariable(:b)
    OIR.addinstruction!(bb, OIR.call(OP_PUSH_INTO_TWO_ARGS, cnst_val, var_mutable_a, var_mutable_b))
end

#=
function fif(x)
   y = 0
   if x > 1
       y = y + 1
   end
   return y * 2
end

#=
1 ─      y = 0
│   %2 = x > 1
└──      goto #3 if not %2
2 ─      y = y + 1
3 ┄ %5 = y * 2
└──      return %5
=#

@testset "CFG" begin
    cfg = OIR.CFG()
    bb = cfg.start
    OIR.assign!(bb, OIR.Variable(:y), OIR.constant(0))
end
=#

@testset "Benchmarks" begin
    include("benchmarks.jl")
end


text:
## NOT USED AT PRESENT

push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Documenter, Jchemo

makedocs(
    modules = [Jchemo],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Matthieu Lesnoff",
    sitename = "Jchemo.jl",
    pages = Any["index.md",
        "regression.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/mlesnoff/Jchemo.jl.git",
)


text:
export normcdf, norminv, normcdfCody

const DBL_EPSILON = eps()
const DBL_MIN = floatmin() #different from math.SmallestNon0Float64

using SpecialFunctions

const sqrt2 = sqrt(2)
const OneOverSqrtTwo = 1 / sqrt(2)

# normcdf(z::Real) =  normcdfCody(z) #
normcdf(z::Float64) = erfc(-z * OneOverSqrtTwo)/2
norminv(z::Float64) = -sqrt2*erfcinv(2*z)

normcdf(z::T) where {T} = erfc(-z/sqrt(T(2)))/2
norminv(z::T) where {T} = -sqrt(T(2))*erfcinv(2*z)

#const 1.0 = 1.0
#const 0.5 = 0.5
#const 0 = 0.0
const _normal_sixten = 1.6
const thrsh = 0.66291
const root32 = 5.656854248
const OneOverSqrt2Pi = 1 / sqrt(2 * pi)

 #Reference: Cody, W.D. (1993). ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of Special Function Routines and Test Drivers" ACM Transactions on Mathematical Software. 19, 22-32.
 # This function evaluates the normal distribution function: The main computation evaluates near-minimax approximations derived from those in "Rational Chebyshev approximations for the error function" by W. J. Cody, Math. Comp., 1969,
 # 631-637. This transportable program uses rational functions that theoretically approximate the normal distribution function to at least 18 significant decimal digits. The accuracy achieved depends on the arithmetic system, the compiler,
 # the intrinsic functions, and proper selection of the machine-dependent constants.

#
#normcdf(z::Float64) = normcdfCody(z)

const c0 = 0.39894151208813466764
const c1 = 8.8831497943883759412
const c2 = 93.506656132177855979
const c3 = 597.27027639480026226
const c4 = 2494.5375852903726711
const c5 = 6848.1904505362823326
const c6 = 11602.651437647350124
const c7 = 9842.7148383839780218
const c8 = 1.0765576773720192317e-8

const d0 = 22.266688044328115691
const d1 = 235.38790178262499861
const d2 = 1519.377599407554805
const d3 = 6485.558298266760755
const d4 = 18615.571640885098091
const d5 = 34900.952721145977266
const d6 = 38912.003286093271411
const d7 = 19685.429676859990727

const p0 = 0.21589853405795699
const p1 = 0.1274011611602473639
const p2 = 0.022235277870649807
const p3 = 0.001421619193227893466
const p4 = 2.9112874951168792e-5
const p5 = 0.02307344176494017303

const q0 = 1.28426009614491121
const q1 = 0.468238212480865118
const q2 = 0.0659881378689285515
const q3 = 0.00378239633202758244
const q4 = 7.29751555083966205e-5

const a0 = 2.2352520354606839287
const a1 = 161.02823106855587881
const a2 = 1067.6894854603709582
const a3 = 18154.981253343561249
const a4 = 0.065682337918207449113

const b0 = 47.20258190468824187
const b1 = 976.09855173777669322
const b2 = 10260.932208618978205
const b3 = 45507.789335026729956

function normcdfCody(x::Float64)::Float64

    xden = 0.0
    temp = 0.0
    xnum = 0.0
    result = 0.0
    ccum = 0.0
    min = 0.0
    eps = 0.0
    xsq = 0.0
    y = 0.0

    eps = DBL_EPSILON * 0.5
    min = DBL_MIN
    y = abs(x)
    if y <= thrsh
        # Evaluate pnorm for |z| <= 0.66291 */
        xsq = 0
        if y > eps
            xsq = x^2
        end
        xnum = a4 * xsq
        xden = xsq
        xnum = (xnum + a0) * xsq
        xden = (xden + b0) * xsq
        xnum = (xnum + a1) * xsq
        xden = (xden + b1) * xsq
        xnum = (xnum + a2) * xsq
        xden = (xden + b2) * xsq
        result = x * (xnum + a3) / (xden + b3)
        temp = result
        result = 0.5 + temp
        ccum = 0.5 - temp
    elseif y <= root32
        #/* Evaluate pnorm for 0.66291 <= |z| <= sqrt(32) */
        xnum = c8 * y
        xden = y
        xnum = (xnum + c0) * y
        xden = (xden + d0) * y
        xnum = (xnum + c1) * y
        xden = (xden + d1) * y
        xnum = (xnum + c2) * y
        xden = (xden + d2) * y
        xnum = (xnum + c3) * y
        xden = (xden + d3) * y
        xnum = (xnum + c4) * y
        xden = (xden + d4) * y
        xnum = (xnum + c5) * y
        xden = (xden + d5) * y
        xnum = (xnum + c6) * y
        xden = (xden + d6) * y
        result = (xnum + c7) / (xden + d7)
        result *= pdfDenormalized(y)
        #            xsq = fint(y * _normal_sixten) / _normal_sixten;
        #            del = (y - xsq) * (y + xsq);
        #            result = math.Exp(-xsq * xsq * 0.5) * math.Exp(-del * 0.5) * result;
        ccum = 1.0 - result
        if x > 0
            temp = result
            result = ccum
            ccum = temp
        end
    else
        # Evaluate pnorm for |z| > sqrt(32) */
        result = 0
        xsq = 1.0 / (x * x)
        xnum = p5 * xsq
        xden = xsq
        xnum = (xnum + p0) * xsq
        xden = (xden + q0) * xsq
        xnum = (xnum + p1) * xsq
        xden = (xden + q1) * xsq
        xnum = (xnum + p2) * xsq
        xden = (xden + q2) * xsq
        xnum = (xnum + p3) * xsq
        xden = (xden + q3) * xsq
        result = xsq * (xnum + p4) / (xden + q4)
        result = (OneOverSqrt2Pi - result) / y
        result *= pdfDenormalized(y)
        #            xsq = fint(x * _normal_sixten) / _normal_sixten;
        #            del = (x - xsq) * (x + xsq);
        #           result = math.Exp(-xsq * xsq * 0.5) * math.Exp(-del * 0.5) * result;
        ccum = 1.0 - result
        if x > 0
            temp = result
            result = ccum
            ccum = temp
        end
    end
    if result < min
        result = 0.0
    end
    if ccum < min
        ccum = 0.0
    end
    return result
end

function pdfDenormalized(x::Float64)::Float64
	return exp(-0.5 * x * x)
end


text:
module Clusters

export AbstractCluster,
       AbstractClusterWithSymbols,
       center_cluster!,
       center_coordinates,
       Cluster,
       cluster_angle,
       ClusterNoSymbols,
       dihedral_angle,
       distances,
       move!,
       print_xyz,
       rotate_x!,
       rotate_y!,
       rotate_z!,
       rotate_randomly!

using ..Atoms
using Distances: Euclidean, pairwise, euclidean
using LinearAlgebra
using Rotations

import Base.==, Base.+

abstract type AbstractCluster end
abstract type AbstractClusterWithSymbols{T} <: AbstractCluster where T<:AbstractAtom end

"""
    ClusterNoSymbols <: AbstractCluster

Basic cluster has only location of atoms.

# Fields
- `xyz::Array{Float64,2}` : location of atoms in 3d space, first index is x, y, z coordinate
"""
mutable struct ClusterNoSymbols <: AbstractCluster
    xyz::Array{Float64,2}
    function ClusterNoSymbols(xyz::AbstractArray{<:Real,2})
        if size(xyz,1) != 3
            throw(DimensionMismatch("ClusterNoSymbols - xyz has wrong dimensions"))
        end
        new(xyz)
    end
    function ClusterNoSymbols(xyz::AbstractArray{<:Real,1})
        if length(xyz) != 3
            throw(DimensionMismatch("ClusterNoSymbols - atoms has wrong dimensions length=$(size(xyz))"))
        end
        new(reshape(xyz,3,1))
    end
end

"""
    Cluster{T} <: AbstractClusterWithSymbols{T} where T<:AbstractAtom

Structure to hold location data of clusters/molecules

# Fields
- `xyz::Array{Float64,2}` : location of atoms in 3d space, first index is x, y, z coordinate
- `atoms::Vector{T}` : atom type information
"""
mutable struct Cluster{T} <: AbstractClusterWithSymbols{T}
    "Location of atoms"
    xyz::Array{Float64,2}
    "Symbols for atoms"
    atoms::Vector{T}
    function Cluster(xyz::AbstractArray{<:Real,2}, atoms::AbstractVector{T}) where T<:AbstractAtom
        if size(xyz,2) != length(atoms)
            throw(DimensionMismatch("Cluster has different sizes for atoms $(size(atoms)) and xyz $(size(xyz))"))
        elseif size(xyz,1) != 3
            throw(DimensionMismatch("Cluster - xyz has wrong dimensions"))
        end
        new{T}(xyz,atoms)
    end
    function Cluster(xyz::AbstractArray{<:Real,1}, atoms::AbstractVector{T}) where T<:AbstractAtom
        if length(atoms) != 1
            throw(DimensionMismatch("Cluster has different number for atoms $(length(atoms)) and xyz 1"))
        elseif length(xyz) != 3
            throw(DimensionMismatch("Cluster - xyz has wrong dimensions"))
        end
        new{T}(reshape(xyz,3,1),atoms)
    end
    function Cluster(xyz::AbstractVector{<:Real}, atom::T) where T<:AbstractAtom
        if length(xyz) != 3
            throw(DimensionMismatch("Cluster - xyz has wrong dimensions size=$(size(xyz))"))
        end
        new{T}(reshape(xyz,3,1), [atom])
    end
end

Base.show(io::IO, C::AbstractCluster) = print(io,"$(typeof(C)) of ", size(C), " atoms")

Base.size(a::AbstractCluster) = size(a.xyz,2)

Base.length(a::AbstractCluster) = size(a)

Base.lastindex(a::AbstractCluster) = length(a)

function (+)(c1::Cluster{T}, c2::Cluster{T}) where T
    return Cluster(hcat(c1.xyz,c2.xyz),vcat(c1.atoms,c2.atoms))
end

function (+)(c1::ClusterNoSymbols, c2::ClusterNoSymbols)
    return ClusterNoSymbols(hcat(c1.xyz,c2.xyz))
end

function Base.getindex(C::Cluster, i::Int)
    return Cluster(C.xyz[:,i], [C.atoms[i]])
end

function Base.getindex(C::ClusterNoSymbols, i::Int)
    return ClusterNoSymbols(C.xyz[:,i])
end

function Base.getindex(C::Cluster, i::AbstractUnitRange)
    return Cluster(C.xyz[:,i], C.atoms[i])
end

function Base.getindex(C::ClusterNoSymbols, i::AbstractUnitRange)
    return ClusterNoSymbols(C.xyz[:,i])
end

function Base.print(io::IO, C::AbstractClusterWithSymbols)
    for i in 1:size(C)
        println(io, C.atoms[i].id, "    ",C.xyz[1,i], "  ", C.xyz[2,i], "  ",C.xyz[3,i])
    end
end

"""
    distances(c1::AbstractCluster, c2::AbstractCluster)

Return distances between atoms of given clusters
"""
function distances(c1::AbstractCluster, c2::AbstractCluster)
    return pairwise(Euclidean(),c1.xyz,c2.xyz,dims=2)
end

"""
    distances(c::AbstractCluster, ar1::AbstractRange, ar2::AbstractRange)

Returns distances between atoms in given unit ranges
"""
function distances(c::AbstractCluster, ur1::AbstractUnitRange, ur2::AbstractUnitRange)
    return pairwise(Euclidean(),c.xyz[:,ur1],c.xyz[:,ur2],dims=2)
end

"""
    distances(c::AbstractCluster, i, j)

Returns distance of atoms `i` and `j`
"""
function distances(c::AbstractCluster, i, j)
    return euclidean(c.xyz[:,i], c.xyz[:,j])
end

"""
    distances(c1::AbstractCluster, i, c2::AbstractCluster, j)

Return distance between atom `i` in `c1` and atom `j` in `c2`
"""
function distances(c1::AbstractCluster, i, c2::AbstractCluster, j)
    return euclidean(c1.xyz[:,i], c2.xyz[:,j])
end

"""
    distances(c1::AbstractCluster, ur1::AbstractUnitRange,
              c2::AbstractCluster, ur2::AbstractUnitRange)

Return distance between atoms `ur1` in `c1` and atoms `ur2` in `c2`
"""
function distances(c1::AbstractCluster, ur1::AbstractUnitRange,
                   c2::AbstractCluster, ur2::AbstractUnitRange)
    return pairwise(Euclidean(),c1.xyz[:,ur1],c2.xyz[:,ur2],dims=2)
end

"""
    center_coordinates(c::AbstractCluster)

Gives coordinates to aricmetric mean of clusters atoms
"""
function center_coordinates(c::AbstractCluster)
    return sum(c.xyz, dims=2) / length(c)
end

"""
    move!(c::AbstractCluster,r)

Moves cluster by `r`
"""
function move!(c::AbstractCluster,r)
    c.xyz .+= r
end

"""
    center_cluster!(c::AbstractCluster)

Centers cluster to origin of coordinates
"""
function center_cluster!(c::AbstractCluster)
    move!(c, -center_coordinates(c))
end

"""
    rotate_x!(c::AbstractCluster, θ)

Rotates cluster around x-axis by angle `θ`
"""
function rotate_x!(c::AbstractCluster, θ)
    c.xyz = RotX(θ) * c.xyz
end

"""
    rotate_y!(c::AbstractCluster, θ)

Rotates cluster around y-axis by angle `θ`
"""
function rotate_y!(c::AbstractCluster, θ)
    c.xyz = RotY(θ) * c.xyz
end

"""
    rotate_z!(c::AbstractCluster, θ)

Rotates cluster around z-axis by angle `θ`
"""
function rotate_z!(c::AbstractCluster, θ)
    c.xyz = RotZ(θ) * c.xyz
end

"""
    rotate_randomly!(c::AbstractCluster)

Rotate cluster by random angle and axis
"""
function rotate_randomly!(c::AbstractCluster)
    c.xyz = rand(RotMatrix{3}) * c.xyz
end

"""
    print_xyz(io::IO, c::AbstractClusterWithSymbols, note=""; printheader=true)

Prints cluster in xyz file format

# Arguments
- `io::IO` : stream where writing is done
- `c::AbstractClusterWithSymbols` : cluster that is writen
- `note=""` : message writen on note line
- `printheader=true` : wheather or not header is writen (number of atoms and note)
"""
function print_xyz(io::IO, c::AbstractClusterWithSymbols, note=""; printheader=true)
    if printheader
        println(io, "    ",length(c))
        println(io, note)
    end
    for i in 1:length(c)
        println(io, c.atoms[i].id, "   ", c.xyz[1,i], "  ", c.xyz[2,i], "  ", c.xyz[3,i])
    end
end

"""
cluster_angle(c::AbstractCluster, i, j, k)

Calculates angle (radians) between atoms i,j,k in cluster
"""
function cluster_angle(c::AbstractCluster, i, j, k)
    r1 = c.xyz[:,i] - c.xyz[:,j]
    r2 = c.xyz[:,k] - c.xyz[:,j]
    return acos(dot(r1,r2)/sqrt(dot(r1,r1)*dot(r2,r2)))
end

"""
cluster_angle(c1::AbstractCluster, i, j, c2::AbstractCluster, k)

Calculates angle (radians) between atons in different clusters

# Arguments
- `c1::AbstractCluster` : first cluster
- `i` : index in `c1`
- `j` : index in `c1`
- `c2::AbstractCluster` : second cluster
- `k` : index in `c2`
"""
function cluster_angle(c1::AbstractCluster, i, j, c2::AbstractCluster, k)
    r1 = c1.xyz[:,i] - c1.xyz[:,j]
    r2 = c2.xyz[:,k] - c1.xyz[:,j]
    return acos(dot(r1,r2)/sqrt(dot(r1,r1)*dot(r2,r2)))
end

"""
cluster_dihedral_angle(c::AbstractCluster, i, j, k, m)
"""
function dihedral_angle(c::AbstractCluster, i, j, k, m)
    r1 = c.xyz[:,j] - c.xyz[:,i]
    r2 = c.xyz[:,k] - c.xyz[:,j]
    r3 = c.xyz[:,m] - c.xyz[:,k]
    t1 = cross(cross(r1,r2), cross(r2,r3))
    t2 = dot(cross(r1,r2), cross(r2,r3))
    return atan( dot(t1, r2./norm(r2)), t2 )
end

end #module Clusters

# This is for compability for older versions
module clusters
    using ..Clusters
    cluster = Clusters.Cluster
end #module clusters


text:
struct PoissonSolver{DT <: Real}
    nx::Int
    Δx::DT
    xgrid::Vector{DT}
    ρ::Vector{DT}
    ϕ::Vector{DT}

    function PoissonSolver{DT}(nx::Int) where {DT}
        Δx = 1/nx
        xgrid = collect(0:Δx:1)
        new(nx, Δx, xgrid, zeros(DT, nx), zeros(DT, nx))
    end
end

function solve!(p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    h = fit(Histogram, mod.(x, 1), p.xgrid)
    p.ρ .= h.weights ./ length(x)
    ρ̂ = rfft(p.ρ)
    k² = [(i-1)^2 for i in eachindex(ρ̂)]
    ϕ̂ = - ρ̂ ./ k²
    ϕ̂[1] = 0
    p.ϕ .= irfft(ϕ̂, length(p.ρ))
    return p
end

function eval_field(p::PoissonSolver{DT}, x::DT) where {DT}
    y = mod(x, one(x))
    i1 = floor(Int, y / p.Δx) + 1
    i2 = mod( ceil(Int, y / p.Δx), p.nx) + 1
    i1 == i2 && (i1 = i1-1)
    i1 == 0 && (i1 = lastindex(p.ϕ))
    return - (p.ϕ[i2] - p.ϕ[i1]) / p.Δx
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments