text:
import Arpack

"""
    LpNormLMO{T, p}(right_hand_side)

LMO with feasible set being an L-p norm ball:
```
C = {x ∈ R^n, norm(x, p) ≤ right_hand_side}
```
"""
struct LpNormLMO{T,p} <: LinearMinimizationOracle
    right_hand_side::T
end

LpNormLMO{p}(right_hand_side::T) where {T,p} = LpNormLMO{T,p}(right_hand_side)

function compute_extreme_point(lmo::LpNormLMO{T,2}, direction; v = similar(direction), kwargs...) where {T}
    dir_norm = norm(direction, 2)
    n = length(direction)
    # if direction numerically 0
    if dir_norm <= 10eps(eltype(direction))
        @. v = lmo.right_hand_side / sqrt(n)
    else
        @. v = -lmo.right_hand_side * direction / dir_norm
    end
    return v
end

function compute_extreme_point(lmo::LpNormLMO{T,Inf}, direction; v = similar(direction), kwargs...) where {T}
    for idx in eachindex(direction)
        v[idx] = -lmo.right_hand_side * (1 - 2signbit(direction[idx]))
    end
    return v
end

function compute_extreme_point(lmo::LpNormLMO{T,1}, direction; v = nothing, kwargs...) where {T}
    idx = 0
    v = -one(eltype(direction))
    for i in eachindex(direction)
        if abs(direction[i]) > v
            v = abs(direction[i])
            idx = i
        end
    end
    sign_coeff = sign(direction[idx])
    if sign_coeff == 0.0
        sign_coeff -= 1
    end
    return ScaledHotVector(-lmo.right_hand_side * sign_coeff, idx, length(direction))
end

function compute_extreme_point(lmo::LpNormLMO{T,p}, direction; v = similar(direction), kwargs...) where {T,p}
    # covers the case where the Inf or 1 is of another type
    if p == Inf
        v = compute_extreme_point(LpNormLMO{T,Inf}(lmo.right_hand_side), direction, v = v)
        return v
    elseif p == 1
        v = compute_extreme_point(LpNormLMO{T,1}(lmo.right_hand_side), direction)
        return v
    end
    q = p / (p - 1)
    pow_ratio = q / p
    q_norm = norm(direction, q)^(pow_ratio)
    # handle zero_direction first
    # assuming the direction is a vector of 1
    if q_norm < eps()
        one_vec = trues(length(direction))
        @. v = -lmo.right_hand_side * one_vec^(pow_ratio) / oftype(q_norm, 1)
        return v
    end
    @. v = -lmo.right_hand_side * sign(direction) * abs(direction)^(pow_ratio) / q_norm
    return v
end

# temporary oracle for l_1 ball to

struct L1ballDense{T} <: LinearMinimizationOracle
    right_hand_side::T
end

function compute_extreme_point(lmo::L1ballDense{T}, direction; v = zeros(T, length(direction)), kwargs...) where {T}
    idx = 0
    val = -1.0
    for i in eachindex(direction)
        if abs(direction[i]) > val
            val = abs(direction[i])
            idx = i
        end
    end

    v .= 0
    v[idx] = T(-lmo.right_hand_side * sign(direction[idx]))
    return v
end

"""
    KNormBallLMO{T}(K::Int, right_hand_side::T)

LMO with feasible set being the K-norm ball in the sense of
[2010.07243](https://arxiv.org/abs/2010.07243),
i.e., the convex hull over the union of an
L_1-ball with radius τ and an L_∞-ball with radius τ/K:
```
C_{K,τ} = conv { B_1(τ) ∪ B_∞(τ / K) }
```
with `τ` the `right_hand_side` parameter. The K-norm is defined as
the sum of the largest `K` absolute entries in a vector.
"""
struct KNormBallLMO{T} <: LinearMinimizationOracle
    K::Int
    right_hand_side::T
end

function compute_extreme_point(lmo::KNormBallLMO{T}, direction; v = similar(direction), kwargs...) where {T}
    K = max(min(lmo.K, length(direction)), 1)

    oinf = zero(eltype(direction))
    idx_l1 = 0
    val_l1 = -one(eltype(direction))

    @inbounds for (i, dir_val) in enumerate(direction)
        temp = -lmo.right_hand_side / K * sign(dir_val)
        v[i] = temp
        oinf += dir_val * temp
        abs_v = abs(dir_val)
        if abs_v > val_l1
            idx_l1 = i
            val_l1 = abs_v
        end
    end

    v1 = ScaledHotVector(-lmo.right_hand_side * sign(direction[idx_l1]), idx_l1, length(direction))
    o1 = dot(v1, direction)
    if o1 < oinf
        @. v = v1
    end
    return v
end

"""
    NuclearNormLMO{T}(radius)

LMO over matrices that have a nuclear norm less than `radius`.
The LMO returns the rank-one matrix with singular value `radius`.
"""
struct NuclearNormLMO{T} <: LinearMinimizationOracle
    radius::T
end

NuclearNormLMO{T}() where {T} = NuclearNormLMO{T}(one(T))
NuclearNormLMO() = NuclearNormLMO(1.0)

"""
Best rank-one approximation using the greatest singular value computed with Arpack.

Warning: this does not work (yet) with all number types, BigFloat and Float16 fail.
"""
function compute_extreme_point(lmo::NuclearNormLMO{TL}, direction::AbstractMatrix{TD}; tol=1e-8, kwargs...) where {TL, TD}
    T = promote_type(TD, TL)
    Z = Arpack.svds(direction, nsv=1, tol=tol)[1]
    u = -lmo.radius * view(Z.U, :)
    return RankOneMatrix(u::Vector{T}, Z.V[:]::Vector{T})
end

function convert_mathopt(
    lmo::NuclearNormLMO,
    optimizer::OT;
    row_dimension::Integer,
    col_dimension::Integer,
    use_modify=true::Bool,
    kwargs...,
) where {OT}
    MOI.empty!(optimizer)
    x = MOI.add_variables(optimizer, row_dimension * col_dimension)
    (t, _) = MOI.add_constrained_variable(optimizer, MOI.LessThan(lmo.radius))
    MOI.add_constraint(optimizer, [t; x], MOI.NormNuclearCone(row_dimension, col_dimension))
    return MathOptLMO(optimizer, use_modify)
end

"""
    SpectraplexLMO{T,M}(radius::T,gradient_container::M,ensure_symmetry::Bool=true)

Feasible set
```
{X ∈ 𝕊_n^+, trace(X) == radius}
```
`gradient_container` is used to store the symmetrized negative direction.
`ensure_symmetry` indicates whether the linear function is made symmetric before computing the eigenvector.
"""
struct SpectraplexLMO{T,M} <: LinearMinimizationOracle
    radius::T
    gradient_container::M
    ensure_symmetry::Bool
end

function SpectraplexLMO(radius::T, side_dimension::Int, ensure_symmetry::Bool=true) where {T}
    return SpectraplexLMO(
        radius,
        Matrix{T}(undef, side_dimension, side_dimension),
        ensure_symmetry,
    )
end

function SpectraplexLMO(radius::Integer, side_dimension::Int, ensure_symmetry::Bool=true)
    return SpectraplexLMO(float(radius), side_dimension, ensure_symmetry)
end

function compute_extreme_point(lmo::SpectraplexLMO{T}, direction::M; v = nothing, maxiters=500, kwargs...) where {T,M <: AbstractMatrix}
    lmo.gradient_container .= direction
    if !(M <: Union{LinearAlgebra.Symmetric, LinearAlgebra.Diagonal, LinearAlgebra.UniformScaling}) && lmo.ensure_symmetry
        # make gradient symmetric
        @. lmo.gradient_container += direction'
    end
    lmo.gradient_container .*= -1

    _, evec = Arpack.eigs(lmo.gradient_container; nev=1, which=:LR, maxiter=maxiters)
    # type annotation because of Arpack instability
    unit_vec::Vector{T} = vec(evec)
    # scaling by sqrt(radius) so that x x^T has spectral norm radius while using a single vector
    unit_vec .*= sqrt(lmo.radius)
    return FrankWolfe.RankOneMatrix(unit_vec, unit_vec)
end

"""
    UnitSpectrahedronLMO{T,M}(radius::T, gradient_container::M)

Feasible set of PSD matrices with bounded trace:
```
{X ∈ 𝕊_n^+, trace(X) ≤ radius}
```
`gradient_container` is used to store the symmetrized negative direction.
`ensure_symmetry` indicates whether the linear function is made symmetric before computing the eigenvector.
"""
struct UnitSpectrahedronLMO{T,M} <: LinearMinimizationOracle
    radius::T
    gradient_container::M
    ensure_symmetry::Bool
end

function UnitSpectrahedronLMO(radius::T, side_dimension::Int, ensure_symmetry::Bool=true) where {T}
    return UnitSpectrahedronLMO(
        radius,
        Matrix{T}(undef, side_dimension, side_dimension),
        ensure_symmetry,
    )
end

UnitSpectrahedronLMO(radius::Integer, side_dimension::Int) = UnitSpectrahedronLMO(float(radius), side_dimension)

function compute_extreme_point(lmo::UnitSpectrahedronLMO{T}, direction::M; v = nothing, maxiters=500, kwargs...) where {T, M <: AbstractMatrix}
    lmo.gradient_container .= direction
    if !(M <: Union{LinearAlgebra.Symmetric, LinearAlgebra.Diagonal, LinearAlgebra.UniformScaling}) && lmo.ensure_symmetry
        # make gradient symmetric
        @. lmo.gradient_container += direction'
    end
    lmo.gradient_container .*= -1

    e_val::Vector{T}, evec::Matrix{T} = Arpack.eigs(lmo.gradient_container; nev=1, which=:LR, maxiter=maxiters)
    # type annotation because of Arpack instability
    unit_vec::Vector{T} = vec(evec)
    if e_val[1] < 0
        # return a zero rank-one matrix
        unit_vec .*= 0
    else
        # scaling by sqrt(radius) so that x x^T has spectral norm radius while using a single vector
        unit_vec .*= sqrt(lmo.radius)
    end
    return FrankWolfe.RankOneMatrix(unit_vec, unit_vec)
end

function convert_mathopt(
    lmo::Union{SpectraplexLMO{T}, UnitSpectrahedronLMO{T}},
    optimizer::OT;
    side_dimension::Integer,
    use_modify::Bool=true,
    kwargs...,
) where {T, OT}
    MOI.empty!(optimizer)
    X = MOI.add_variables(optimizer, side_dimension * side_dimension)
    MOI.add_constraint(optimizer, X, MOI.PositiveSemidefiniteConeSquare(side_dimension))
    sum_diag_terms = MOI.ScalarAffineFunction{T}([],zero(T))
    # collect diagonal terms of the matrix
    for i in 1:side_dimension
        push!(sum_diag_terms.terms, MOI.ScalarAffineTerm(one(T), X[i + side_dimension * (i-1)]))
    end
    constraint_set = if lmo isa SpectraplexLMO
        MOI.EqualTo(lmo.radius)
    else
        MOI.LessThan(lmo.radius)
    end
    MOI.add_constraint(optimizer, sum_diag_terms, constraint_set)
    return MathOptLMO(optimizer, use_modify)
end


text:
# Utility functions for simulating data
function sim_poisson(trans::AbstractVector{<:Integer})
    # Simulates Poisson noise for a given input vector

    # First, determine the unique distributions from which we need to draw
    λs = unique(trans)
    N = length(λs)
    for λ ∈ λs
        # Generate different Poisson distribution for each
        d = Poisson(λ)
        # Find indices to fill
        ind = findall(==(λ),trans)
        nval = length(ind)
        # Generate nval random number
        vals = rand(d,nval)
        # Put values back into the array
        trans[ind] = vals
    end
    return trans
end

ext: jl


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
using YAActL, Test

t = Ref{Task}()                # this is for debugging
lp = LinkParams(taskref=t)

arg = YAActL.Args(1, 2, c=3, d=4)
@test arg.args == (1, 2)
@test arg.kwargs == pairs((c=3, d=4))

a = b = c = d = 1
e = f = 2

incx(x, by; y=0, z=0) = x+y+z + by
subx(x, y, sub; z=0) = x+y+z - sub

A = Actor(lp, incx, a, y=b, z=c)
sleep(0.1)
@test info(A) == :runnable

# test diag and actor startup, become! (implicitly)
act = YAActL.diag!(A)
sleep(0.1)
@test act.dsp == full
@test act.sta == nothing
@test act.bhv.f == incx
@test act.bhv.args == (1,)
@test act.bhv.kwargs == pairs((y=1,z=1))

# test explicitly become!
become!(A, subx, a, b, z=c)
sleep(0.1)
@test act.bhv.f == subx
@test act.bhv.args == (1,1)
@test act.bhv.kwargs == pairs((z=1,))

# test set!
set!(A, state)
sleep(0.1)
@test act.dsp == state

# test update!
update!(A, (1, 2, 3))
sleep(0.1)
@test act.sta == (1,2,3)
update!(A, Args(2,3, x=1, y=2), s=:arg)
sleep(0.1)
@test act.bhv.args == (2,3)
@test act.bhv.kwargs == pairs((x=1,y=2,z=1))

# test query!
@test query!(A) == (1,2,3)
@test query!(A, :res) == nothing
@test query!(A, :bhv) == subx
@test query!(A, :dsp) == state

# test call!
become!(A, incx, a, y=b, z=c)
set!(A, full)
@test call!(A, 1) == 4
@test query!(A, :res) == 4
@test query!(A) == (1,2,3)
set!(A, state)
update!(A, 1)
update!(A, Args(y=2,z=2), s=:arg)
@test call!(A, 2) == 7
@test query!(A) == 7
@test query!(A, :res) == 7

# test cast!
update!(A, 2)
update!(A, Args(y=1,z=1), s=:arg)
cast!(A, 2)
@test query!(A) == 6
@test query!(A, :res) == 6
set!(A, full)
update!(A, Args(a, y=3,z=3), s=:arg)
cast!(A, 3)
@test query!(A, :res) == 10
@test query!(A) == 6

# test exec!
@test exec!(A, Func(cos, 2pi)) == 1

# test exit!
exit!(A)
sleep(0.1)
@test info(A).state == :closed


text:
module CheapThreads

using ThreadingUtilities, VectorizationBase
using ArrayInterface: static_length, static_step, static_first
using VectorizationBase: num_threads, cache_linesize, __vload, __vstore!, register_size, False
using StrideArraysCore: object_and_preserve, dereference
import IfElse
using Static
using Requires

export batch, @batch, num_threads

include("request.jl")
include("batch.jl")
include("closure.jl")
include("unsignediterator.jl")

# reset_workers!() = WORKERS[] = UInt128((1 << (num_threads() - 1)) - 1)
dynamic_thread_count() = min((Sys.CPU_THREADS)::Int, Threads.nthreads())
reset_workers!() = WORKERS[] = UInt128((1 << (dynamic_thread_count() - 1)) - 1)
function __init__()
    reset_workers!()
    resize!(STATES, dynamic_thread_count() * cache_linesize())
    STATES .= 0x00
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
end

end


text:
# ---
# title: 654. Maximum Binary Tree
# id: problem654
# author: Indigo
# date: 2021-06-28
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/maximum-binary-tree/description/>
# hidden: true
# ---
#
# Given an integer array with no duplicates. A maximum tree building on this
# array is defined as follow:
#
#   1. The root is the maximum number in the array.
#   2. The left subtree is the maximum tree constructed from left part subarray divided by the maximum number.
#   3. The right subtree is the maximum tree constructed from right part subarray divided by the maximum number.
#
# Construct the maximum tree by the given array and output the root node of this
# tree.
#
# **Example 1:**
#
#
#
#     Input: [3,2,1,6,0,5]
#     Output: return the tree root node representing the following tree:
#
#           6
#         /   \
#        3     5
#         \    /
#          2  0
#            \
#             1
#
#
# **Note:**
#
#   1. The size of the given array will be in the range [1,1000].
#
#
## @lc code=start
using LeetCode

function construct_maximum_binary_tree(nums::AbstractVector{Int})
    isempty(nums) && return nothing
    max_idx = argmax(nums)
    return TreeNode(nums[max_idx],
                    construct_maximum_binary_tree(@view(nums[1:(max_idx - 1)])),
                    construct_maximum_binary_tree(@view(nums[(max_idx + 1):end])))
end
## @lc code=end


text:
using Aqua
using Recalls
using Test

@testset "Aqua" begin
    Aqua.test_all(Recalls)
end

@testset "recall" begin
    @recall 1 + 2
    @recall 3 + 4
    @test Recalls.recall() == 7
    @test Recalls.CALLS[end-1]() == 3
end

@testset "note" begin
    empty!(Recalls.NOTES)
    @note a = 1
    @note a = 2
    notes = Recalls.sortbytime!()
    @test [notes[1].a, notes[2].a] == [1, 2]
end


text:
function budget_test(retired :: Bool)
    lastWorkAge = 30;
    b = init_test_budget(lastWorkAge);
    @test validate_budget(b);

    if retired
        t = lastWorkAge + 1;
        efficiency = 0.0;
    else
        t = lastWorkAge;
        efficiency = [0.9, 2.1];
    end

    k = [1.2, 2.4];

    y = income(b, t, k, efficiency);
    @test size(y) == size(k)

    kp = kprime(b, t, k, 0.0, efficiency);
    @test isapprox(y, kp)

    c = 0.6;
    kp = kprime(b, t, k, c, efficiency);
    @test size(kp) == size(k)

    c2 = consumption(b, t, k, kp, efficiency);
    @test size(c2) == size(k)
    @test all(isapprox.(c, c2))
end

# function retired_test()
#     b = init_test_budget();
#     k = [1.2, 2.4];

#     y = retired_income(b, k);
#     @test size(y) == size(k)

#     kp = retired_kprime(b, k, 0.0);
#     @test isapprox(y, kp)

#     c = 0.6;
#     kp = retired_kprime(b, k, c);
#     @test size(kp) == size(k)

#     c2 = retired_consumption(b, k, kp);
#     @test size(c2) == size(k)
#     @test all(isapprox.(c, c2))
# end

@testset "Budget" begin
    for retired in [true, false]
        budget_test(retired);
    end
    # retired_test();
end

# -----------


text:
using JSON

json(rootNode::RootNode) = map(to_dict, children(rootNode)) |> JSON.json

function to_dict(factGroup::FactsCollectionNode)
  _children = factGroup |> children
  childFactGroups = filter(child -> isa(child, FactsCollectionNode), _children)
  childFacts = filter(child -> isa(child, FactNode), _children)
  childContexts = filter(child -> isa(child, ContextNode), _children)

  Dict(
    "name"       => factGroup |> name,
    "line"       => factGroup |> line,
    "factGroups" => map(to_dict, childFactGroups),
    "facts"      => map(to_dict, childFacts),
    "contexts"   => map(to_dict, childContexts)
  )
end

to_dict(contextNode::ContextNode) = Dict(
    "name"  => contextNode |> name,
    "line"  => contextNode |> line,
    "facts" => map(to_dict, children(contextNode))
  )

to_dict(factNode::FactNode) = Dict(
    "name"         => factNode |> name,
    "line"         => factNode |> line,
    "result"       => factNode |> result |> string,
    "details"      => factNode |> details,
    "stacktrace"   => factNode |> stacktrace
  )


text:
module BenchAppend

using BangBang: Empty, append!!, empty!!
using BangBang.NoBang: SingletonVector
using BenchmarkTools: @benchmarkable, BenchmarkGroup
using StaticArrays: SVector

function repeat_append!(dest, src, n_repeat = 100, n_values = 1000)
    for _ = 1:n_repeat
        x = empty!!(dest)
        for _ in 1:n_values
            append!!(x, src)
        end
    end
end

repeat_append(dest, src) =
    foldl(append!!, ntuple(_ -> src, 10); init = dest)

suite = BenchmarkGroup()
suite["append!!(::Vector, ::SingletonVector)"] =
    @benchmarkable(repeat_append!($(Float64[]), SingletonVector(1.0)))
suite["append!!(::SVector, ::SingletonVector)"] =
    @benchmarkable(repeat_append(SVector{0,Float64}(), SingletonVector(1.0)))
suite["append!!(Empty(SVector), ::SingletonVector)"] =
    @benchmarkable(repeat_append(Empty(SVector), SingletonVector(1.0)))

end  # module
BenchAppend.suite


text:
abstract type TrustRegion end

struct RectangularRegion <: TrustRegion end
struct EllipsoidalRegion <: TrustRegion end

const Φ = 1 / 4  # MathConstants.φ

function bound(::TrustRegion, z, Δz)
    if 2Δz > z
        return Φ * z
    end
    return Δz
end
#
function bound(::RectangularRegion, z::TR, Δz)
    steeper = Δz.r * z.t - Δz.t * z.r > 0
    if 2Δz.r > z.r && steeper
        zʳ = Φ * z.r
        return TR((Δz.t / Δz.r) * zʳ, zʳ)
    end
    if 2Δz.t > z.t && !steeper
        zᵗ = Φ * z.t
        return TR(zᵗ, (Δz.r / Δz.t) * zᵗ)
    end
    return Δz
end
#
function bound(::EllipsoidalRegion, z::TR, Δz)
    if iszero(z.r)
        return TR(min(Δz.t, Φ * z.t), z.r)
    end
    if iszero(z.t)
        return TR(z.t, min(Δz.r, Φ * z.r))
    end
    zᵗ = Δz.t / z.t
    zʳ = Δz.r / z.r
    r = hypot(zᵗ, zʳ)
    if r > 1
        return TR(z.t * (zᵗ / r), z.r * (zʳ / r))
    end
    return Δz
end

function asymptotics_weights!(w, ws, K, S, d, grid)
    # We write `(1 .- inv.(S)).^2 .* S` instead of `(S .- 1).^2 .* inv.(S)`
    # to get `w[j] == Inf` instead of `NaN` whenever `S[j] == Inf`.
    w .= ws .* K.^(d + 1) .* (1 .- inv.(S)).^2 .* S
    return w
end
#
function asymptotics_weights!(w::Vector{U}, ws, K, S, d, grid) where
    {T, U <: MDProjections{2, T}}

    @inbounds for i in eachindex(K)
        Sᵢ = S[i]
        k² = K[i]^2
        wᵣ = ws[i] * k²
        wₜ = wᵣ * k²

        # Projections components of the weights
        t  = wₜ * (1 - inv(Sᵢ.t))^2 * Sᵢ.t
        r₀ = (Sᵢ.r[1] - 1)^2
        r₁ = 3wₜ * r₀ * inv(Sᵢ.r[1])
        r₂ = 6wᵣ * r₀ * inv(Sᵢ.r[2])

        w[i] = MDProjections(t, SVector(r₁, r₂))
    end

    return w
end

memory_term(β, γ) = β * γ
#
function memory_term(β, γ::TR)
    t = β * gett(γ)
    return TR(t, t + 2 * getr(γ))
end

ergodic_param(Sⱼ, Sₑ) = Sⱼ / (Sⱼ + Sₑ)
#
function ergodic_param(Sⱼ::DProjections{2}, Sₑ)
    t = Sⱼ.t / (Sⱼ.t + Sₑ.t)
    @inbounds begin
        r₁ = Sⱼ.r[1] / (Sⱼ.r[1] + Sₑ.r)
        r₂ = Sⱼ.r[2] / (Sⱼ.r[2] + Sₑ.r)
    end
    return DProjections(t, SVector(r₁, r₂))
end
#
function ergodic_param(Sⱼ::LDProjections{1}, Sₑ)
    t = Sⱼ.t / (Sⱼ.t + Sₑ.t)
    @inbounds begin
        r₁ = Sⱼ.r[1] / (Sⱼ.r[1] + Sₑ.r)
    end
    return LDProjections(t, SVector(r₁))
end

function set!(avars::AsymptoticVars, Z, kvars, D₀, ζ)
    @unpack f, fˢ = avars
    @unpack svars, Λ = kvars
    @unpack S, Sˢ, B, w, υ = svars

    γ = D₀ * inv(ζ)

    @inbounds for j in eachindex(Λ)
        Sₑ = memory_term(B[j], γ)
        f[j] = ergodic_param(S[j], Sₑ)
        fˢ[j] = ergodic_param(Sˢ[j], Sₑ)
    end

    avars.ζ∞[] = ζ
    return avars
end

struct FixedPoint{U, V, W} <: Function
    Z::U
    kvars::V
    D₀::W
end
#
(f::FixedPoint)(ζ) = fixedpoint!(f.Z, f.kvars, f.D₀, ζ)

function fixedpoint!(Z, kvars, D₀, ζ)
    @unpack svars, Λ = kvars
    @unpack S, Sˢ, B, w, υ = svars

    γ = D₀ * inv(ζ)

    @inbounds for j in eachindex(Λ)
        Sₑ = memory_term(B[j], γ)
        f = ergodic_param(S[j], Sₑ)
        fˢ = ergodic_param(Sˢ[j], Sₑ)
        Z[j] = product(w[j], f, fˢ)
    end

    return υ * sum(Z)
end


text:
function render_block_trajectory(traj, body, physics)

    anim = Animation()
    T = physics.T
    for t = 1:T
        plot([-2, 2], [-0.015, -0.015], aspect_ratio=:equal, legend=false, linecolor=:black,
            linewidth=4, xlims=[-2,2], ylims=[-0.1,2])
        local_corners = [[-0.2 -0.2 0.2 0.2 -0.2 ];
                         [-0.2 0.2 0.2 -0.2 -0.2];
                         [1. 1. 1. 1. 1.]]

        H = [[cos(traj[3, t]) -sin(traj[3, t]) traj[1, t]];
             [sin(traj[3, t]) cos(traj[3,t]) traj[2, t]];
             [0. 0. 1]]
        world_corners = H * local_corners

        plot!([world_corners[1,:]], [world_corners[2,:]], color=:blue, linewidth=2)
        sleep(0.01)
        frame(anim)
    end
    gif(anim, "media/block_sim.gif", fps=15)

end

function render_block_trajectory(traj, body, physics, use_makie::Bool)
    T = physics.T
    lc = [  [-0.2 -0.2 0.2 0.2 -0.2 ];
            [-0.2 0.2 0.2 -0.2 -0.2];
            [1. 1. 1. 1. 1.]]
    world_corners = Observable([Point2f0(-0.2, -0.2), Point2f0(-0.2, 0.2),
                                Point2f0(0.2, 0.2), Point2f0(0.2, -0.2),
                                Point2f0(-0.2, -0.2)])
    fig = Figure()
    ax = Axis(fig[1,1], aspect=DataAspect(), limits=(-2,2,-0.1,2))
    lines!(ax, [-2,2], [-0.015, -0.015], color=:black, linewidth=4)
    lines!(ax, world_corners, color=:blue, linewidth=2)

    function animstep!(t)
            H = [[cos(traj[3, t]) -sin(traj[3, t]) traj[1, t]];
                    [sin(traj[3, t]) cos(traj[3,t]) traj[2, t]];
                    [0. 0. 1]]
            wc = H * lc
            world_corners[] = [Point2f0(x[1],x[2]) for x in [wc[:,i] for i in 1:size(wc,2)]]
            sleep(0.01)
    end
    record(fig, "media/block_sim_makie.gif", 1:T; framerate=60) do t
        animstep!(t)
    end
end

function render_ball_trajectory(traj, body, physics)

    anim = Animation()
    T = physics.T
    r = body.l/2.
    circle(x,y,r) = (x .+ r*sin.(LinRange(0, 2*π, 500)), y .+ r*cos.(LinRange(0, 2*π, 500)))
    for t = 1:T
        plot([-2, 2], [-0.015, -0.015], aspect_ratio=:equal,
            legend=false, linecolor=:black,
            linewidth=4, xlims=[-2,2], ylims=[-0.1,2])

        local_center = [0.0; 0.0; 1.0]

        H = [[cos(traj[3, t]) -sin(traj[3, t]) traj[1, t]];
             [sin(traj[3, t]) cos(traj[3,t]) traj[2, t]];
             [0. 0. 1]]
        world_center = H * local_center
        x = world_center[1]
        y = world_center[2]

        # plot!(circle(x,y,r), seriestype=[:shape], color=:blue, linewidth=2)
        plot!(circle(x,y,r), color=:blue, linewidth=2)
        sleep(0.01)
        frame(anim)
    end
    gif(anim, "media/ball_sim.gif", fps=15)

end

function render_ball_trajectory(traj, body, physics, use_makie::Bool)
    T = physics.T
    r = body.l/2.
    lc = [0.0; 0.0; 1.0]
    center = Observable(Point2f0(0.,0.))

    fig = Figure()
    ax = Axis(fig[1,1], aspect=DataAspect(), limits=(-2,2,-0.1,2))
    lines!(ax, [-2,2], [-0.015, -0.015], color=:black, linewidth=4)
    Makie.scatter!(ax, center; marker=:circle,  linewidth=4, strokecolor=:blue, markersize= 80)

    function animstep!(t)
            H = [[cos(traj[3, t]) -sin(traj[3, t]) traj[1, t]];
                    [sin(traj[3, t]) cos(traj[3,t]) traj[2, t]];
                    [0. 0. 1]]
            wc = H * lc
            center[] = Point2f0(wc[1], wc[2])
            sleep(0.01)
    end
    record(fig, "media/ball_sim_makie.gif", 1:T; framerate=15) do t
        animstep!(t)
    end
end

function render_walking_trajectory(ts, zs, walker, steps, fps)
    zs_tp = hcat(zs[:,1], zs[:,3], zs[:,5], zs[:,6])
    num_var = 4
    total_frames = Int(round(ts[end]*fps))
    t = 0:ts[end]/total_frames:ts[end]
    z = zeros((total_frames+1,num_var))
    for i = 1:num_var
        interp = LinearInterpolation(ts, zs_tp[:,i])
        z[:,i] = interp.(t)
    end

    l = walker.l
    c = walker.c

    mm = size(z, 1)
    min_xh = min(z[:,3]...); max_xh = max(z[:,3]...)

    window_xmin = -1*l; window_xmax = 1*l
    window_ymin = -0.1; window_ymax = 1.1*l
    anim = Plots.Animation()
    rampref = [min_xh-1 max_xh+1; -0.015 -0.015]

    plot(rampref[1,:], rampref[2,:], aspect_ratio=:equal, legend=false, linecolor=:black, linewidth=4,xlims=[window_xmin, window_xmax], ylims=[window_ymin, window_ymax])

    for i=1:mm
        plot(rampref[1,:], rampref[2,:], aspect_ratio=:equal, legend=false, linecolor=:black, linewidth=4,xlims=[window_xmin, window_xmax], ylims=[window_ymin, window_ymax])

        θ₁ = z[i,1]; θ₂ = z[i,2]
        xh = z[i,3]; yh = z[i,4]

        hinge = [xh; yh]
        stance_foot = [xh+l*sin(θ₁); yh-l*cos(θ₁)]
        stance_leg_com = [xh+c*sin(θ₁); yh-c*cos(θ₁)]
        swing_foot = [xh+l*sin(θ₁ +θ₂); yh-l*cos(θ₁+θ₂)]
        swing_leg_com = [xh+c*sin(θ₁+θ₂); yh-c*cos(θ₁+θ₂)]

        plot!([hinge[1]], [hinge[2]], markershape=:square, markercolor=:black, markersize=15)

        plot!([hinge[1], stance_foot[1]], [hinge[2], stance_foot[2]], color=:red, linewidth=2)
        plot!([hinge[1], swing_foot[1]], [hinge[2], swing_foot[2]], color=:blue, linewidth=2)

        fl = 0.05
        stance_feet = [[stance_foot[1] - (fl/2.)*cos(θ₁), stance_foot[2] - (fl/2.)*sin(θ₁)], [stance_foot[1] + (fl/2.)*cos(θ₁), stance_foot[2] + (fl/2.)*sin(θ₁)]]

        swing_feet = [[swing_foot[1] - (fl/2.)*cos(θ₁+θ₂), swing_foot[2] - (fl/2.)*sin(θ₁+θ₂)], [swing_foot[1] + (fl/2.)*cos(θ₁+θ₂), swing_foot[2] + (fl/2.)*sin(θ₁+θ₂)]]

        plot!([stance_feet[1][1], stance_feet[2][1]], [stance_feet[1][2], stance_feet[2][2]], color=:red, linewidth=4)

        plot!([swing_feet[1][1], swing_feet[2][1]], [swing_feet[1][2], swing_feet[2][2]], color=:blue, linewidth=4)

        sleep(0.01)
        frame(anim)
    end
    if walker.use_controller
        gif(anim, "media/walker_sim_ctrl.gif", fps=15)
    else
        gif(anim, "media/walker_sim_no_ctrl.gif", fps=15)
    end
end

function render_walking_trajectory(ts, zs, walker, steps, fps, use_makie::Bool)
    zs_tp = hcat(zs[:,1], zs[:,3], zs[:,5], zs[:,6])
    num_var = 4
    total_frames = Int(round(ts[end]*fps))
    t = 0:ts[end]/total_frames:ts[end]
    z = zeros((total_frames+1,num_var))
    for i = 1:num_var
        interp = LinearInterpolation(ts, zs_tp[:,i])
        z[:,i] = interp.(t)
    end

    l = walker.l
    c = walker.c
    fl = 0.05

    T = size(z, 1)
    min_xh = min(z[:,3]...); max_xh = max(z[:,3]...)

    window_xmin = -1*l; window_xmax = 1*l
    window_ymin = -0.1; window_ymax = 1.1*l
    rampref = [min_xh-1 max_xh+1; -0.015 -0.015]

    fig = Figure()
    ax = Axis(fig[1,1], aspect=DataAspect(), limits=(window_xmin,window_xmax,window_ymin,window_ymax))
    lines!(ax, rampref[1,:], rampref[2,:], color=:black, linewidth=10)

    hinge = Observable(Point2f0(0.0, 0.0))
    stance_foot = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    swing_foot = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    stance_feet = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    swing_feet = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])

    Makie.scatter!(ax, hinge; marker=:rect,  linewidth=4, strokecolor=:black, color=:black, markersize= 80)
    lines!(ax, stance_foot, color=:red, linewidth=6)
    lines!(ax, swing_foot, color=:blue, linewidth=6)
    lines!(ax, stance_feet, color=:red, linewidth=6)
    lines!(ax, swing_feet, color=:blue, linewidth=6)

    function animstep!(i)
        θ₁ = z[i,1]; θ₂ = z[i,2]
        xh = z[i,3]; yh = z[i,4]
        hinge[] = Point2f0(xh, yh)
        stance_foot[] = [Point2f0(xh, yh), Point2f0(xh+l*sin(θ₁), yh-l*cos(θ₁))]
        swing_foot[] = [Point2f0(xh, yh), Point2f0(xh+l*sin(θ₁ +θ₂), yh-l*cos(θ₁+θ₂))]

        stance_feet[] = [Point2f0(xh+l*sin(θ₁) - (fl/2.)*cos(θ₁), yh-l*cos(θ₁) - (fl/2.)*sin(θ₁)), Point2f0(xh+l*sin(θ₁) + (fl/2.)*cos(θ₁), yh-l*cos(θ₁) + (fl/2.)*sin(θ₁))]
        swing_feet[] = [Point2f0(xh+l*sin(θ₁ +θ₂) - (fl/2.)*cos(θ₁+θ₂), yh-l*cos(θ₁+θ₂) - (fl/2.)*sin(θ₁+θ₂)), Point2f0(xh+l*sin(θ₁ +θ₂) + (fl/2.)*cos(θ₁+θ₂), yh-l*cos(θ₁+θ₂) + (fl/2.)*sin(θ₁+θ₂))]

        sleep(0.001)
    end

    record(fig, "media/walker_sim_makie.gif", 1:T; framerate=60) do t
        for j=1:2
            animstep!(t)
        end
    end

end

function render_kick_trajectory(ts, zs, qs, walker, steps, fps, body, physics,
                                use_makie::Bool; trail=true)
    zs_tp = hcat(zs[:,1], zs[:,3], zs[:,5], zs[:,6])
    num_var = 4
    total_frames = Int(round(ts[end]*fps))
    t = 0:ts[end]/total_frames:ts[end]
    z = zeros((total_frames+1,num_var))
    for i = 1:num_var
        interp = LinearInterpolation(ts, zs_tp[:,i])
        z[:,i] = interp.(t)
    end

    l = walker.l
    c = walker.c
    fl = 0.05
    ind = 1

    T = size(z, 1)
    N = physics.T
    min_xh = min(z[:,3]...); max_xh = max(z[:,3]...)

    window_xmin = -1*l; window_xmax = 1*l
    window_ymin = -0.1; window_ymax = 1.1*l
    rampref = [min_xh-1 max_xh+1; -0.015 -0.015]

    fig = Figure()
    ax = Axis(fig[1,1], aspect=DataAspect(), limits=(-2,4,-0.1,2))
    lines!(ax, [-2,4], [-0.015, -0.015], color=:black, linewidth=10)

    hinge = Observable(Point2f0(0.0, 0.0))
    stance_foot = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    swing_foot = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    stance_feet = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    swing_feet = Observable([Point2f0(0.0, 0.0),Point2f0(0.0, 0.0)])
    tail = 300
    traj = CircularBuffer{Point2f0}(tail)
    fill!(traj, Point2f0(0.0, 0.2))
    traj = Observable(traj)

    c = to_color(:purple)
    tailcol = [RGBAf0(c.r, c.g, c.b, (i/tail)^2) for i in 1:tail]

    lc = [0.0; 0.0; 1.0]
    center = Observable(Point2f0(0.1,0.19))

    Makie.scatter!(ax, hinge; marker=:rect,  linewidth=4, strokecolor=:black, color=:black, markersize= 50)
    Makie.scatter!(ax, center; marker=:circle,  linewidth=4, strokecolor=:blue, markersize= 50)
    lines!(ax, stance_foot, color=:red, linewidth=6)
    lines!(ax, swing_foot, color=:blue, linewidth=6)
    lines!(ax, stance_feet, color=:red, linewidth=6)
    lines!(ax, swing_feet, color=:blue, linewidth=6)
    lines!(ax, traj; linewidth = 3, color = tailcol)

    function ballanimstep!(t)
        if t <= N
            t = Int(round(t))
            H = [[cos(qs[3, t]) -sin(qs[3, t]) qs[1, t]];
                    [sin(qs[3, t]) cos(qs[3,t]) qs[2, t]];
                    [0. 0. 1]]
            wc = H * lc
            center[] = Point2f0(wc[1], wc[2])
            push!(traj[], Point2f0(wc[1], wc[2]))
            traj[] = traj[]
            sleep(0.001)
        end
    end

    function animstep!(i)
        if i<T
            θ₁ = z[i,1]; θ₂ = z[i,2]
            xh = z[i,3]; yh = z[i,4]
            hinge[] = Point2f0(xh, yh)
            stance_foot[] = [Point2f0(xh, yh), Point2f0(xh+l*sin(θ₁), yh-l*cos(θ₁))]
            swing_foot[] = [Point2f0(xh, yh), Point2f0(xh+l*sin(θ₁ +θ₂), yh-l*cos(θ₁+θ₂))]

            stance_feet[] = [Point2f0(xh+l*sin(θ₁) - (fl/2.)*cos(θ₁), yh-l*cos(θ₁) - (fl/2.)*sin(θ₁)), Point2f0(xh+l*sin(θ₁) + (fl/2.)*cos(θ₁), yh-l*cos(θ₁) + (fl/2.)*sin(θ₁))]
            swing_feet[] = [Point2f0(xh+l*sin(θ₁ +θ₂) - (fl/2.)*cos(θ₁+θ₂), yh-l*cos(θ₁+θ₂) - (fl/2.)*sin(θ₁+θ₂)), Point2f0(xh+l*sin(θ₁ +θ₂) + (fl/2.)*cos(θ₁+θ₂), yh-l*cos(θ₁+θ₂) + (fl/2.)*sin(θ₁+θ₂))]
            if xh >= -0.252
                # println("attack angle ",θ₂)
                ballanimstep!(ind)
                ind+=1
            end
        else
            # println("from T ")
            ballanimstep!(i+ind-T-1)
        end

        sleep(0.001)
    end

    record(fig, "media/kick100.mp4", 1:T+N; framerate=60) do t
        animstep!(t)
    end

end


text:
import DelayEmbeddings.Dataset

abstract type DatasetTriangulationPartition{D, T} <: AbstractTriangulationPartition{D, T} end
abstract type DatasetTriangulationPartitionFull{D, T} <: AbstractTriangulationPartitionFull{D, T} end

abstract type MutableDatasetTriangulationPartition{D, T} <: AbstractMutableTriangulationPartition{D, T} end
abstract type MutableDatasetTriangulationPartitionFull{D, T} <: AbstractMutableTriangulationPartitionFull{D, T} end

struct DatasetTriangulation{D, T} <: DatasetTriangulationPartition{D, T}
	data::Dataset
	simplexindices::DelaunayTriangulation
end

mutable struct MutableDatasetTriangulation{D, T} <: MutableDatasetTriangulationPartition{D, T}
	data::Dataset
	simplexindices::DelaunayTriangulation
end

# Full info
struct DatasetTriangulationFull{D, T} <: DatasetTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

# Full info
struct MutableDatasetTriangulationFull{D, T} <: MutableDatasetTriangulationPartitionFull{D, T}
	data::AbstractArray{T, 2}
	simplexindices::DelaunayTriangulation
    simplices::Vector{Simplex}
    radii::Vector{Float64}
    orientations::Vector{Float64}
    volumes::Vector{Float64}
    centroids::Vector{Vector{Float64}}
end

PT = Union{DatasetTriangulationPartition, DatasetTriangulationPartitionFull,
    MutableDatasetTriangulationPartition, MutableDatasetTriangulationPartitionFull}
getpoints(pt::PT) = [pt.data[i] for i = 1:length(pt.data)]

getsimplices(pt::Union{DatasetTriangulationPartitionFull, MutableDatasetTriangulationPartitionFull}) =
    pt.simplices
getsimplices(pt::Union{DatasetTriangulationPartition, MutableDatasetTriangulationPartition}) =
    [Simplex(pt.data[pt.simplexindices[i]]) for i = 1:length(pt.simplexindices)]

####################
# Pretty printing
####################
function summarise(tp::Union{DatasetTriangulation, MutableDatasetTriangulation})
    _type = typeof(tp)
    n_simplices = nsimplices(tp)
    n_pts = length(tp.data)
    D = dimension(tp)
    summary = "$_type from $n_simplices simplices and $n_pts Datasets"
end

Base.show(io::IO, tp::Union{DatasetTriangulation, MutableDatasetTriangulation}) =
    println(io, summarise(tp))

export
DatasetTriangulation,
MutableDatasetTriangulation,
DatasetTriangulationFull,
MutableDatasetTriangulationFull


text:
"""
# ZeroConf.jl

[DNS Service Discovery](https://en.wikipedia.org/wiki/Zero-configuration_networking#DNS-based_service_discovery)
(DNS-SD, ZeroConf, Bonjour) interface for Julia.

On macOS (or BSD) this package uses the `dns-sd` tool (see `man dns-sd`).
On Linux the [Avahi](https://www.avahi.org) tools are used instead
(`apt-get install avahi-utils`).

"""
module ZeroConf

export dns_service_browse, register_dns_service, unregister_dns_service

"""
name => (host, port)
"""
const DNSService = Pair{String, Tuple{String, Int16}}

"""
    dns_service_browse([type = "_http._tcp"]) -> Channel{DNSService}

Open a channel to browse for DNS Services of a certian `type`.
Reading from channel yeilds pairs: `name => (host, port)`.

e.g.

    julia> c = dns_service_browse()
    julia> take!(c)
    "AxisCamera" => ("axis-ptz.local.", 80)
"""
dns_service_browse() = dns_service_browse("_http._tcp")

"""
    register_dns_service(name, service_type, port)

Register a DNS Service.
"""
register_dns_service(name, service_type, port)

dns_process = Dict{Tuple{String,String},Base.Process}()

"""
    unregister_dns_service(name, service_type)

Cancel a service registered by `register_dns_service`.
"""
function unregister_dns_service(name, service_type)
    global dns_process
    kill(dns_process[(name, service_type)])
end

const dnssd_is_installed = Sys.which("dns-sd") != nothing
const avahi_is_installed = Sys.which("avahi-publish") != nothing

@static if dnssd_is_installed
    include("dnssd.jl")
elseif avahi_is_installed
    include("avahi.jl")
else
    @error """
        Can't find `dns-sd` or `avahi-publish`.
        Try `apt-get install avahi-utils`?
        """
end

readme() = join([
    Docs.doc(@__MODULE__),
    Docs.doc(dns_service_browse),
    Docs.doc(register_dns_service),
    Docs.doc(unregister_dns_service)
   ], "\n\n")

end # module


text:
using ArgParse, ApplicationBuilder

const julia_v07 = VERSION > v"0.7-"

s = ArgParseSettings(autofix_names = true)  # turn "-" into "_" for arg names.

@add_arg_table s begin
    "juliaprog_main"
        arg_type = String
        required = true
        help = "Julia program to compile -- must define julia_main()"
    "appname"
        arg_type = String
        default = nothing
        help = "name to call the generated .app bundle; defaults to basename(juliaprog_main)"
    "builddir"
        arg_type = String
        default = "builddir"
        help = "build directory for output."
    "--verbose", "-v"
        action = :store_true
        help = "increase verbosity"
    "--resource", "-R"
        arg_type = String
        action = :append_arg  # Can specify multiple
        dest_name = "resources"  # each -R appends to resources.
        default = String[]
        metavar = "<resource>"
        help = """specify files or directories to be copied to
                  MyApp.app/Contents/Resources/. This should be done for all
                  resources that your app will need to have available at
                  runtime. Can be repeated.

                  NOTE: following the system conventions, -R /path/dir will copy
                  "dir" to Resources/, but -R /path/dir/ will copy all *contents*
                  of `dir/*` to Resources/."""
    "--lib", "-L"
        arg_type = String
        action = :append_arg  # Can specify multiple
        dest_name = "libraries"  # each -L appends to libraries.
        default = String[]
        metavar = "<file>"
        help = """specify user library files to be copied to
                  MyApp.app/Contents/Libraries/. This should be done for all
                  libraries that your app will need to reference at
                  runtime. Can be repeated."""
    "--bundle-identifier"
        arg_type = String
        default = nothing
        metavar = "com.user.appname"
        help = "the bundle identifier for this app. Default: 'com.<USER>.<appname>'."
    "--icns"
        arg_type = String
        default = nothing
        dest_name = "icns_file"
        metavar = "<file>"
        help = ".icns file to be used as the app's icon"
    "--app-version"
        arg_type = String
        default = "0.1"
        #range_tester = (x -> r"^[0-9]+(\.[0-9]+)*$"(x))  # can the version have other characters in it? idk..
        metavar = "0.0.1"
        help = ".icns file to be used as the app's icon"
    "--certificate"
        arg_type = String
        default = nothing
        metavar = "<cert_name>"
        help = "name of the certificate to use to sign the app"
    "--entitlements"
        arg_type = String
        default = nothing
        metavar = "<file>"
        help = ".entitlements file. Must set --certificate to add entitlements."
end
s.epilog = """
    examples:\n
    \ua0\ua0 # Build HelloApp.app from hello.jl\n
    \ua0\ua0 build.jl hello.jl HelloApp\n
    \ua0\ua0 # Build MyGame, and copy in imgs/, mus.wav and all files in libs/\n
    \ua0\ua0 build.jl -R imgs -R mus.wav -L lib/* main.jl MyGame
    """

parsed_args = parse_args(ARGS, s; as_symbols=true)

juliaprog_main = pop!(parsed_args, :juliaprog_main)
if julia_v07
    filter!(kv -> kv.second ∉ (nothing, false), parsed_args)
else
    filter!((k, v) -> v ∉ (nothing, false), parsed_args)
end

ApplicationBuilder.build_app_bundle(juliaprog_main; parsed_args...)


text:
#
#   Posituve Definite Karcher Mean (Matlab Manopt)
#
using Manopt, Manifolds, ManifoldsBase, LinearAlgebra, Random
Random.seed!(42)
n = 5
m = 5
M = SymmetricPositiveDefinite(n)
x = random_point(M)
A = [random_point(M) for _ in 1:m]
A = [Symmetric(a) for a in A]
F(X::Array{Float64,2}) = sum([distance(M, X, B)^2 for B in A]) / (2 * m)
∇F(X::Array{Float64,2}) = -sum([log(M, X, B) for B in A]) / m

@time quasi_Newton(
    M,
    F,
    ∇F,
    x;
    memory_size=100,
    stopping_criterion=StopWhenGradientNormLess(norm(M, x, ∇F(x)) * 10^(-6)),
    debug=[:Iteration, " ", :Cost, "\n", 1, :Stop],
)

# B1 = quasi_Newton(M,F,∇F,x; memory_size = 100, debug = [:Iteration, " ", :Cost, "\n", 1, :Stop])
# B2 = mean(M,A)


text:
function descendants(
    node, key;
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true, # like continue in the R package, but actually the opposite
    self = false,
    filter_fun = nothing,
    recursivity_level = -1,
    ignore_nothing::Bool = false,
    type::Union{Union,DataType} = Any)

    # Check the filters once, and then compute the descendants recursively using `descendants_`
    check_filters(node, scale = scale, symbol = symbol, link = link)

    # Change the filtering function if we also want to remove nodes with nothing values:
    filter_fun_ = filter_fun_nothing(filter_fun, ignore_nothing, key)

    val = Array{type,1}()

    if self
        keep = is_filtered(node, scale, symbol, link, filter_fun_)

        if keep
            val_ = unsafe_getindex(node, key)
            push!(val, val_)
        elseif !all
            # We don't keep the value and we have to stop at the first filtered-out value
            return val
        end
    end

    descendants_(node, key, scale, symbol, link, all, filter_fun_, val, recursivity_level)
    return val
end

function descendants_(node, key, scale, symbol, link, all, filter_fun, val, recursivity_level)

    if !isleaf(node) && recursivity_level != 0
        for chnode in ordered_children(node)
            # Is there any filter happening for the current node? (FALSE if filtered out):
            keep = is_filtered(chnode, scale, symbol, link, filter_fun)

            if keep
                val_ = unsafe_getindex(chnode, key)
                push!(val, val_)
                # Only decrement the recursivity level when the current node is not filtered-out
                recursivity_level -= 1
            end

            # If we want to continue even if the current node is filtered-out
            if all || keep
                descendants_(chnode, key, scale, symbol, link, all, filter_fun, val, recursivity_level)
            end
        end
    end
end

function descendants!(
    node, key;
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true, # like continue in the R package, but actually the opposite
    self = false,
    filter_fun = nothing,
    recursivity_level = -1,
    ignore_nothing = false,
    type::Union{Union,DataType} = Any)

    # Check the filters once, and then compute the descendants recursively using `descendants_`
    check_filters(node, scale = scale, symbol = symbol, link = link)

    key_cache = cache_name(key, scale, symbol, link, all, self, filter_fun, type)
    val = Array{type,1}()

    # Change the filtering function if we also want to remove nodes with nothing values:
    filter_fun_ = filter_fun_nothing(filter_fun, ignore_nothing, key)

    if self
        keep = is_filtered(node, scale, symbol, link, filter_fun_)

        if keep
            val_ = unsafe_getindex(node, key)
            push!(val, val_)
        elseif !all
            # We don't keep the value and we have to stop at the first filtered-out value
            return val
        end
    end

    if node[key_cache] === nothing

        descendants_!(node, key, scale, symbol, link, all, filter_fun_, val, recursivity_level, key_cache)

        # Caching the result into a cache attribute named after the SHA of the function arguments:
        node[key_cache] = val
    else
        append!(val, node[key_cache])
    end

    return val
end

"""
Fast version of descendants_ that mutates the mtg nodes to cache the information.
"""
function descendants_!(node, key, scale, symbol, link, all, filter_fun, val, recursivity_level, key_cache)

    val_i = Array{eltype(val),1}()
    if !isleaf(node) && recursivity_level != 0
        if node[key_cache] === nothing # Is there any cached value? If so, do not recompute
            for chnode in ordered_children(node)
                # Is there any filter happening for the current node? (FALSE if filtered out):
                keep = is_filtered(chnode, scale, symbol, link, filter_fun)

                if keep
                    val_key = unsafe_getindex(chnode, key)
                    push!(val_i, val_key)
                    # Only decrement the recursivity level when the current node is not filtered-out
                    recursivity_level -= 1
                    # chnode[key_cache] = val_key
                end
                # If we want to continue even if the current node is filtered-out
                if all || keep
                    descendants_!(chnode, key, scale, symbol, link, all, filter_fun, val_i, recursivity_level, key_cache)
                end
            end
            node[key_cache] = val_i
            append!(val, val_i)
        else
            append!(val, copy(node[key_cache]))
            # node[key_cache]
        end
    end
end

"""
    clean_cache!(mtg)

Clean the cached variables in the mtg, usually added from [`descendants!`](@ref).
"""
function clean_cache!(mtg)
    cached_vars = find_cached_vars(mtg)
    for i in cached_vars
        traverse!(mtg, x -> pop!(x, i))
    end
end

function find_cached_vars(node)
    collect(keys(node.attributes))[findall(x -> occursin("_cache_", x), String.(keys(node.attributes)))]
end

"""
    descendants(node::Node,key,<keyword arguments>)
    descendants!(node::Node,key,<keyword arguments>)

Get attribute values from the descendants (acropetal). The mutating version (`descendants!`)
cache the results in a cached variable named after the hash of the function call. This version
is way faster for large trees, but require to clean the chache sometimes (see [`clean_cache!`](@ref)).
It also only works for trees with attributes of subtype of `AbstractDict`.

# Arguments

## Mandatory arguments

- `node::Node`: The node to start at.
- `key`: The key, or attribute name. Make it a `Symbol` for faster computation time.

## Keyword Arguments

- `scale = nothing`: The scale to filter-in (i.e. to keep). Usually a Tuple-alike of integers.
- `symbol = nothing`: The symbol to filter-in. Usually a Tuple-alike of Strings.
- `link = nothing`: The link with the previous node to filter-in. Usually a Tuple-alike of Char.
- `all::Bool = true`: Return all filtered-in nodes (`true`), or stop at the first node that
is filtered out (`false`).
- `self = false`: is the value for the current node needed ?
- `filter_fun = nothing`: Any filtering function taking a node as input, e.g. [`isleaf`](@ref).
- `recursivity_level = -1`: The maximum number of recursions allowed (considering filters).
*E.g.* to get the first level children only: `recursivity_level = 1`, for children +
grand-children: `recursivity_level = 2`. If a negative value is provided (the default), the
function returns all valid values from the node to the leaves.
- `ignore_nothing = false`: filter-out the nodes with `nothing` values for the given `key`
- `type::Union{Union,DataType}`: The type of the attribute. Makes the function run much
faster if provided (≈4x faster).

# Tips

To get the values of the leaves use [`isleaf`](@ref) as the filtering function, e.g.:
`descendants(mtg, :Width; filter_fun = isleaf)`.

# Note

In most cases, the `type` argument should be given as a union of `Nothing` and the data type
of the attribute to manage missing or inexistant data, e.g. measurements made at one scale
only. See examples for more details.

# Examples

```julia
# Importing the mtg from the github repo:
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

descendants(mtg, :Length) # Short to write, but slower to execute

# Fast version, note that we pass a union of Nothing and Float64 because there are some nodes
# without a `Length` attribute:
descendants(mtg, :Length, type = Union{Nothing,Float64})

# Filter by scale:
descendants(mtg, :XX, scale = 1, type = Float64)
descendants(mtg, :Length, scale = 3, type = Float64)

# Filter by symbol:
descendants(mtg, :Length, symbol = "Leaf")
descendants(mtg, :Length, symbol = ("Leaf","Internode"))

# Filter by function, e.g. get the values for the leaves only:
descendants(mtg, :Width; filter_fun = isleaf)

# You can also ask for different attributes by passing them as a vector:
descendants(mtg, [:Width, :Length]; filter_fun = isleaf)
# The output is an array of arrays of length of the attributes you asked for.

# It is possible to cache the results in the mtg. This is wqy faster when using
# `@mutate_mtg` (note the `!` at the end of the function name):
@mutate_mtg!(mtg, subtree_length = sum(descendants!(node, :Length, symbol = "Internode")), symbol = "Internode")

function compute_subtree_length(x)
    length_descendants = filter(x -> x !== nothing, descendants(x, :Length, symbol = "Internode", self = true))
    length(length_descendants) > 0 ? sum(length_descendants) : nothing
end

function compute_subtree_length!(x)
    length_descendants = filter(x -> x !== nothing, descendants!(x, :Length, symbol = "Internode", self = true))
    length(length_descendants) > 0 ? length_descendants : nothing
end

@mutate_mtg!(mtg, subtree_length = compute_subtree_length(node), symbol = "Internode")

file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
@mutate_mtg!(mtg, subtree_length = compute_subtree_length!(node))

file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
compute_subtree_length!(mtg)
DataFrame(mtg, [:Length, :_cache_9011cfa452383c48086b78014718eeebab7b12b9])

# You can then clean the cach to avoid using too much memory:
clean_cache!(mtg)
mtg[1][1][1].attributes
mtg[1][1][1][2].attributes
DataFrame(mtg, [:Length, :subtree_length])

file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
node = mtg
key = :Length
scale = nothing
symbol = nothing
link = nothing
all= true
self = false
filter_fun = nothing
recursivity_level = -1
type = Any

# compute_subtree_length!(mtg)

descendants!(mtg, :Length, self = true)
DataFrame(mtg, [:Length, :subtree_length, :_cache_9011cfa452383c48086b78014718eeebab7b12b9])
DataFrame(mtg, [:Length, :subtree_length, :_cache_403abd80258f45cfa2a64226edcf3c39c44a3302])

```
"""
descendants!, descendants


text:
function pressure(model::EoSModel, V, T, z=SA[1.])
    return -∂f∂V(model,V,T,z)
end

function VT_entropy(model::EoSModel, V, T, z=SA[1.])
    return -∂f∂T(model,V,T,z)
end

function VT_entropy_res(model::EoSModel, V, T, z=SA[1.])
    fun(x) = eos_res(model,V,x,z)
    return -ForwardDiff.derivative(fun,T)
end

function VT_chemical_potential(model::EoSModel, V, T, z=SA[1.])
    fun(x) = eos(model,V,T,x)
    return ForwardDiff.gradient(fun,z)
end

function VT_chemical_potential_res(model::EoSModel, V, T, z=SA[1.])
    fun(x) = eos_res(model,V,T,x)
    return ForwardDiff.gradient(fun,z)
end

function VT_internal_energy(model::EoSModel, V, T, z=SA[1.])
    dA, A = ∂f(model,V,T,z)
    ∂A∂V, ∂A∂T = dA
    return A - T*∂A∂T
end

function VT_enthalpy(model::EoSModel, V, T, z=SA[1.])
    dA, A = ∂f(model,V,T,z)
    ∂A∂V, ∂A∂T = dA
    return A - V*∂A∂V - T*∂A∂T
end

function VT_gibbs_free_energy(model::EoSModel, V, T, z=SA[1.])
    dA, A = ∂f(model,V,T,z)
    ∂A∂V, ∂A∂T = dA
    return A - V*∂A∂V
end

function VT_helmholtz_free_energy(model::EoSModel, V, T, z=SA[1.])
    return eos(model,V,T,z)
end

function VT_isochoric_heat_capacity(model::EoSModel, V, T, z=SA[1.])
    A(x) = eos(model,V,x,z)
    ∂A∂T(x) = ForwardDiff.derivative(A,x)
    ∂²A∂T²(x) = ForwardDiff.derivative(∂A∂T,x)
    return -T*∂²A∂T²(T)
end

function VT_isobaric_heat_capacity(model::EoSModel, V, T, z=SA[1.])
    d²A = f_hess(model,V,T,z)
    ∂²A∂V∂T = d²A[1,2]
    ∂²A∂V² = d²A[1,1]
    ∂²A∂T² = d²A[2,2]
    return -T*(∂²A∂T² - ∂²A∂V∂T^2/∂²A∂V²)
end

function VT_isothermal_compressibility(model::EoSModel, V, T, z=SA[1.])
    p0,dpdV = p∂p∂V(model,V,T,z)
    return -1/V*dpdV^-1
end

function VT_isentropic_compressibility(model::EoSModel, V, T, z=SA[1.])
    d²A = f_hess(model,V,T,z)
    ∂²A∂V∂T = d²A[1,2]
    ∂²A∂V² = d²A[1,1]
    ∂²A∂T² = d²A[2,2]
    return 1/V*(∂²A∂V²-∂²A∂V∂T^2/∂²A∂T²)^-1
end

function VT_speed_of_sound(model::EoSModel, V, T, z=SA[1.])
    Mr = molecular_weight(model,z)
    d²A = f_hess(model,V,T,z)
    ∂²A∂V∂T = d²A[1,2]
    ∂²A∂V² = d²A[1,1]
    ∂²A∂T² = d²A[2,2]
    return V*sqrt((∂²A∂V²-∂²A∂V∂T^2/∂²A∂T²)/Mr)
end

function VT_isobaric_expansivity(model::EoSModel, V, T, z=SA[1.])
    d²A = f_hess(model,V,T,z)
    ∂²A∂V∂T = d²A[1,2]
    ∂²A∂V² = d²A[1,1]
    ∂²A∂T² = d²A[2,2]
    return ∂²A∂V∂T/(V*∂²A∂V²)
end

function VT_joule_thomson_coefficient(model::EoSModel, V, T, z=SA[1.])
    d²A = f_hess(model,V,T,z)
    ∂²A∂V∂T = d²A[1,2]
    ∂²A∂V² = d²A[1,1]
    ∂²A∂T² = d²A[2,2]
    return -(∂²A∂V∂T - ∂²A∂V²*((T*∂²A∂T² + V*∂²A∂V∂T) / (T*∂²A∂V∂T + V*∂²A∂V²)))^-1
end

"""
    second_virial_coefficient(model::EoSModel, T, z=SA[1.])

Calculates the second virial coefficient `B`, defined as:

```julia
B = lim(V->∞)[ V^2/RT *  (∂Aᵣ∂V + V*∂²Aᵣ∂V²) ]
```
where `Aᵣ` is the residual helmholtz energy.
"""
function second_virial_coefficient(model::EoSModel, T, z=SA[1.])
    TT = promote_type(eltype(z),typeof(T))
    V = 1/sqrt(eps(TT))
    fAᵣ(x) = eos_res(model,x,T,z)
    Aᵣ,∂Aᵣ∂V,∂²Aᵣ∂V² = Solvers.f∂f∂2f(fAᵣ,V)
    return V^2/(R̄*T)*(∂Aᵣ∂V+V*∂²Aᵣ∂V²)
end

function VT_compressibility_factor(model::EoSModel, V, T, z=SA[1.])
    p = pressure(model,V,T,z)
    return p*V/(R̄*T)
end

"""
    pip(model::EoSModel,V,T,z=[1.0])

Phase identification parameter `Π`. as described in _1_. If `Π > 1`, then the phase is clasified as a liquid or a liquid-like vapor, being a vapor or vapor-like liquid otherwise.

This identification parameter fails at temperatures and pressures well aboVe the critical point.

Calculated as:
```
Π = V*((∂²p/∂V∂T)/(∂p/∂T) - (∂²p/∂V²)/(∂p/∂V))
```

1.  G. Venkatarathnama, L.R. Oellrich, Identification of the phase of a fluid using partial derivatives of pressure, volume,and temperature without reference to saturation properties: Applications in phase equilibria calculations, Fluid Phase Equilibria 301 (2011) 225–233

"""
function pip(model::EoSModel, V, T, z=SA[1.0])
    _∂2p = ∂2p(model,V,T,z)
    hess_p, grad_p, _ = _∂2p
    Π = V*(hess_p[1,2]/grad_p[2]  - hess_p[1,1]/grad_p[1])
end

export second_virial_coefficient,pressure


text:
struct HydraSummary
    F_report
    F
    G
    mod_mat
end

function dblcen(D)
    A = D .- mean(D,dims = 1)
    return A .- mean(A,dims = 2)
end

function QRfit(x,y)
    Q, R = LinearAlgebra.qr(x)
    β = (R^(-1)*Q' *y)'
    fitted =   x * β'#sum(x .* β, dims = 2)
    residuals = y .-fitted
    return sum(fitted .^2), sum(Diagonal(residuals) .^2)
end

function QRfit(Q,R,x,y)
    β = (R^(-1)*Q' *y)'
    fitted =   x * β'#sum(x .* β, dims = 2)
    residuals = y .-fitted
    return sum(fitted .^2)/ sum(Diagonal(residuals) .^2)
end

function hydra2(data,D, formula = @formula(1~1); pairs = false)

    mod_mat = StatsModels.modelmatrix(formula,data)
    n,m = size(mod_mat)
    G = -0.5 * dblcen(D .^2)
    Gfit, Gres = QRfit(mod_mat,G)
    F_report = (Gfit/(m-1)) / (Gres/(n-m))
    F = Gfit / Gres

    if pairs == true
        pair_array = Array{HydraSummary}(undef,m+1,m+1)
        bool_mat = Bool.(hcat(sum(mod_mat, dims = 2) .==0, mod_mat))
        for i in 2:m+1
            for j in 1:i-1
                bool_mask = bool_mat[:,i] .| bool_mat[:,j]
                mm = float.(bool_mat[bool_mask,[i,j]])
                g = G[bool_mask,bool_mask]
                pair_array[i,j] = hydra2(g,mm)
            end
        end
        pair_array[1,m+1] = HydraSummary(F_report,F, G,mod_mat)
        return pair_array
    else
    return HydraSummary(F_report,F, G,mod_mat)
    end
end

function hydra2(data,M,metric ::DataType, formula = @formula(1~1); pairs = false)
    D = pairwise(metric(),M',M')
return  hydra2(data,D, formula ; pairs = pairs)
end

function hydra2(G ::Matrix,mod_mat ::Matrix)
    n,m = size(mod_mat)
    Gfit, Gres = QRfit(mod_mat,G)
    F_report = (Gfit/(m-1)) / (Gres/(n-m))
    F = Gfit / Gres
    return HydraSummary(F_report,F, G,mod_mat)
end

function permute(H ::HydraSummary,n_perm = 1000)
    mod_mat = H.mod_mat
    G = H.G
    Q,R = qr(mod_mat)
    inds = collect(1:size(mod_mat)[1])
    F = Vector{Float64}(undef,n_perm)
     for i in 1:n_perm
       F[i] = QRfit(Q,R,mod_mat,view(G,shuffle!(inds),shuffle!(inds)))
    end
    return sum(H.F .< F) /n_perm
end

function permute(H ::Matrix{HydraSummary},n_perm = 1000)
    n = size(H)[1]
    P = zeros(n,n)

    Threads.@threads for i in 2:n
                        for j in 1:i-1
                            P[i,j] = permute(H[i,j],n_perm)
                        end
                    end
                    P[1,n] = permute(H[1,n])
    return P
end

function permute(H ::Matrix{HydraSummary},levels ::Vector,n_perm = 1000)

    n = size(H)[1]
    P = zeros(n,n)
    if length(levels) >n
        levels = unique(sort(levels))
    end

    Threads.@threads for i in 2:n
                        for j in 1:i-1
                            P[i,j] = permute(H[i,j],n_perm)
                        end
                    end
                    P[1,n] = permute(H[1,n])
    return NamedArray(P,(levels,levels))
end

max_stars_repo_name: EvoArt/Hydra

max_issues_repo_name: EvoArt/Hydra

max_forks_repo_name: EvoArt/Hydra


text:
module Oracle

using QuantExQASM.GateOps
using QuantExQASM.NCU
using QuantExQASM.Circuit

function apply_x!(c::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)
    Circuit.add_gatecall!(c, GateOps.pauli_x(tgt, reg))
end

"""
    bitstring_ncu(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx, U::GateOps.GateLabel)

Takes bitstring as the binary pattern and indices as the qubits to operate upon. Applies the appropriate PauliX gates to the control lines to call the NCU with the given matrix
"""
function bitstring_ncu(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx, U::GateOps.GateLabel)
    bitmask =  0x1
    aux_idx = typeof(ctrl_indices)()

    # Filter qubit values to mark specified pattern
    for idx in collect(0:length(ctrl_indices))
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                apply_x!(cct, ctrl_indices[idx+1])
            else
                apply_x!(cct, tgt_idx)
            end
        end
    end
    NCU.apply_ncu!(cct, ctrl_indices, aux_idx, tgt_idx, U)
    for idx in collect(0:length(ctrl_indices))
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                apply_x!(cct, ctrl_indices[idx+1])
            else
                apply_x!(cct, tgt_idx)
            end
        end
    end
end

"""
    bitstring_ncu(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, aux_indices::Vector, tgt_idx, U::GateOps.GateLabel)

Takes bitstring as the binary pattern and indices as the qubits to operate upon. Applies the appropriate PauliX gates to the control lines to call the NCU with the given matrix. Uses aux_indices to reduce circuit depth by expanding width.
"""
function bitstring_ncu(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, aux_indices::Vector, tgt_idx, U::GateOps.GateLabel)
    bitmask =  0x1

    # Filter qubit values to mark specified pattern
    for idx in collect(0:length(ctrl_indices))
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                apply_x!(cct, ctrl_indices[idx+1])
            else
                apply_x!(cct, tgt_idx)
            end
        end
    end
    NCU.apply_ncu!(cct, ctrl_indices, aux_indices, tgt_idx, U)
    for idx in collect(0:length(ctrl_indices))
        if ~( (bitstring & (bitmask << idx) ) != 0)
            if idx < length(ctrl_indices)
                apply_x!(cct, ctrl_indices[idx+1])
            else
                apply_x!(cct, tgt_idx)
            end
        end
    end
end

"""
    bitstring_phase_oracle(bitstring::Unsigned, ctrl_indices::Vector, tgt_idx::Unsigned)

Applies PauliX gates to the appropriate lines in the circuit, then applies a n-controlled PauliZ to mark the state.
"""
function bitstring_phase_oracle(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx)
    return bitstring_ncu(cct, bitstring, ctrl_indices, tgt_idx, GateOps.GateLabel(:z) )
end

"""
    bitstring_phase_oracle(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, tgt_idx)

Applies PauliX gates to the appropriate lines in the circuit, then applies a n-controlled PauliZ to mark the state. Uses aux qubits to reduce circuit depth.
"""
function bitstring_phase_oracle(cct::Circuit.Circ, bitstring::Integer, ctrl_indices::Vector, aux_indices::Vector, tgt_idx)
    return bitstring_ncu(cct, bitstring, ctrl_indices, aux_indices, tgt_idx, GateOps.GateLabel(:z) )
end

end


text:
using JuMP
using CPLEX

include("cluster.jl")

"""
Structure représentant un arbre de décision
"""
mutable struct Tree

    D::Int64 # Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
    a::Matrix{Float64} # Vecteurs caractérisant les règles de branchement (a.x < b)
    b::Vector{Float64} # Idem
    c::Vector{Int64} # c[t] : classe du noeud d'indice t (-1 si aucune classe n'est attribuée)

    function Tree()
        return new()
    end
end

"""
Constructeur par défaut d'un arbre
(les séparations sont celles fournies par le solveur, sans recentrage ni post-traitement)

Entrée :
- D : Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
- a : Vecteurs caractérisant les règles de branchement (a.x < b)
- b : Idem
- c : c[t] est la classe du noeud d'indice t (-1 si aucune classe n'est attribuée)

Sortie :
- this::Tree : l'arbre correspondant
"""
function Tree(D::Int64, a::Matrix{Float64}, b::Vector{Float64}, c::Vector{Int64})
    this = Tree()
    this.D = D
    this.a = a
    this.b = b
    this.c = c

    return this
end

"""
Constructeur d'arbre univarié où la valeur de b est ajustée pour recentrer la séparation par rapport aux données

Entrées :
- D : Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
- a : Vecteurs caractérisant les règles de branchement (a.x < b)
- c : c[t] est la classe du noeud d'indice t (-1 si aucune classe n'est attribuée)
- u : flot des données dans l'arbre
- x : caractéristiques des données d'entraînement

Sortie :
- this::Tree : un arbre
"""
function Tree(D::Int64, a::Matrix{Float64}, c::Vector{Int64}, u::Matrix{Int64}, x::Matrix{Float64})

    this = Tree()
    this.D = D
    this.a = a
    this.c = c

    sepCount = 2^D - 1
    leavesCount = 2^D
    dataCount = length(x[:, 1])
    featuresCount = length(x[1, :])

    # Borne supérieure de b pour chaque séparation
    ub_b = ones(Float64, sepCount)

    # Borne inférieure de b pour chaque séparation
    lb_b = zeros(Float64, sepCount)
    ## Fixe lb_b et ub_b

    # Pour chaque donnée
    for i in 1:dataCount

        # Trouve la feuille où elle arrive
        # (en partant des sommets les plus bas dans l'arbre, trouver le premier sommet atteint par la donnée)
        t = sepCount+leavesCount
        while u[i, t] == 0
            t -= 1
            if t == 0
                break
            end
        end

        # Si cette donnée a un flot dans l'arbre (i.e., si elle est correctement classifiée)
        if t != 0

            # Pour chaque sommet parent a_t et tant que la racine n'est pas atteinte
            while t>1
                a_t = t ÷ 2
                if t % 2 == 0
                    lb_b[a_t] = max(lb_b[a_t], sum(a[j, a_t]*x[i, j] for j in 1:featuresCount))
                else
                    ub_b[a_t] = min(ub_b[a_t], sum(a[j, a_t]*x[i, j] for j in 1:featuresCount))
                end
                t = a_t
            end
        end
    end

    this.b = zeros(Float64, sepCount)
    for t in 1:sepCount
        if c[t] == -1
            this.b[t] = (ub_b[t] + lb_b[t])/2
        end
    end

    return this
end

"""
FONCTION SIMILAIRE A LA PRECEDENTE UTILISEE UNIQUEMENT SI VOUS FAITES DES REGROUPEMENTS

Constructeur d'arbre univarié où la valeur de b est ajustée pour recentrer la séparation par rapport aux données

Entrées :
- D : Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
- a : Vecteurs caractérisant les règles de branchement (a.x < b)
- c : c[t] est la classe du noeud d'indice t (-1 si aucune classe n'est attribuée)
- u : flot des données dans l'arbre
- clusters : partition des données d'entraînement

Sortie :
- this::Tree : un arbre
"""
function Tree(D::Int64, a::Matrix{Float64}, c::Vector{Int64}, u::Matrix{Int64}, clusters::Vector{Cluster})

    this = Tree()
    this.D = D
    this.a = a
    this.c = c

    sepCount = 2^D - 1
    leavesCount = 2^D
    dataCount = length(clusters)
    featuresCount = length(clusters[1].lBounds)

    # Upper bound of b for each separation
    ub_b = ones(Float64, sepCount)

    # Lower bound of b for each separation
    lb_b = zeros(Float64, sepCount)
    ## Fixe lb_b et ub_b

    # Pour chaque cluster
    for i in 1:dataCount

        # Trouver la feuille où il arrive
        # (en partant des sommets les plus bas dans l'arbre, trouver le premier sommet atteint par le cluster)
        t = sepCount+leavesCount
        while u[i, t] == 0
            t -= 1
            if t == 0
                break
            end
        end

        # Si cette donnée a un flot dans l'arbre (i.e., si elle est correctement classifiée)
        if t != 0

            # Pour chaque sommet parent a_t et tant que la racine n'est pas atteinte
            while t>1
                a_t = t ÷ 2

                ## Warning: improve this later. Here the cluster is considered as a box. It would be better to get the distance to the separation for each data point in the cluster

                # If the cluster goes to the left of at
                if t % 2 == 0
                    lb_b[a_t] = max(lb_b[a_t], sum(a[j, a_t]*clusters[i].uBounds[j] for j in 1:featuresCount))

                # If the cluster goes to the right of at
                else
                    ub_b[a_t] = min(ub_b[a_t], sum(a[j, a_t]*clusters[i].lBounds[j] for j in 1:featuresCount))
                end
                t = a_t
            end
        end
    end

    this.b = zeros(Float64, sepCount)
    for t in 1:sepCount
        if c[t] == -1
            this.b[t] = (ub_b[t] + lb_b[t])/2
        end
    end

    return this
end

"""
Constructeur d'arbre multivarié où les valeurs de a et b sont ajustées pour recentrer la séparation par rapport aux données

Entrées :
- D : Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
- c : c[t] est la classe de la donnée d'indice t (-1 si aucune classe n'est attribuée)
- u : flot des données dans l'arbre
- s_model : s[j, t] == 1 ssi le coefficient de la caractéristique j est non nul dans la séparation du sommet t
- x : caractéristiques des données d'entraînement

Sortie :
- this::Tree : un arbre
"""
function Tree(D::Int64, c::Vector{Int64}, u::Matrix{Int64}, s_model::Matrix{Int64}, x::Matrix{Float64})
    this = Tree()
    this.D = D
    this.c = c

    sepCount = 2^D - 1
    dataCount = length(x[:, 1])
    featuresCount = length(x[1, :])

    this.a = zeros(Float64, featuresCount, sepCount)
    this.b = zeros(Float64, sepCount)

    # Pour chaque sommet interne de l'arbre
    for t in 1:sepCount

        # Si une séparation y est effectuée
        if c[t] == -1

            # Indice des données allant à droite ou à gauche
            I_R = Int64[]
            I_L = Int64[]
            for i in 1:dataCount
                if u[i, t*2] == 1
                    push!(I_L, i)
                elseif u[i, t*2+1] == 1
                    push!(I_R, i)
                end
            end

            len_l = length(I_L)
            len_r = length(I_R)

            m = Model(CPLEX.Optimizer)
            set_silent(m)
            @variable(m, a[1:featuresCount], base_name = "a_j")
            @variable(m, s[1:featuresCount], Bin, base_name = "s_j")
            @variable(m, b, base_name = "b")
            @variable(m, e[1:(len_l + len_r)], base_name = "ecart") # Ecart entre chaque variable et la séparation
            @variable(m, e_min, base_name = "ecart_min")

            @constraint(m, -1 <= b)
            @constraint(m, b <= 1)
            @constraint(m, [j in 1:featuresCount], -s[j] <= a[j])
            @constraint(m, [j in 1:featuresCount], a[j] <= s[j])
            @constraint(m, sum(s[j]  for j in 1:featuresCount) <= sum(s_model[j, t] for j in 1:featuresCount)) # on ne veut pas augmenter la "complexité" de la séparation
            @constraint(m, [i in 1:(len_l + len_r)], e[i] >= e_min)
            @constraint(m, [i in 1:len_l], e[i] == b - sum(a[j]*x[I_L[i], j] for j in 1:featuresCount))
            @constraint(m, [i in 1:len_r], e[i+len_l] == - b + sum(a[j]*x[I_R[i], j] for j in 1:featuresCount))

            @objective(m, Max, e_min)
            optimize!(m)

            this.b[t] = value.(b)
            for j in 1:featuresCount
                this.a[j, t] = value.(a[j])
            end
        end
    end

    return this
end

"""
FONCTION SIMILAIRE A LA PRECEDENTE UTILISEE UNIQUEMENT SI VOUS FAITES DES REGROUPEMENTS

Constructeur d'arbres multivariés où les valeurs de a et b sont ajustées pour recentrer la séparation par rapport aux données

Entrées :
- D : Nombre maximal de séparations d'une branche (profondeur de l'arbre - 1)
- c : c[t] est la classe du cluster d'indice t (-1 si aucune classe n'est attribuée)
- u : flot des données dans l'arbre
- s_model : s[j, t] == 1 ssi le coefficient de la caractéristique j est non nul dans la séparation du sommet t
- clusters : liste des clusters regroupant les données d'entraînement

Sortie :
- this::Tree : un arbre
"""
function Tree(D::Int64, c::Vector{Int64}, u::Matrix{Int64}, s_model::Matrix{Int64}, clusters::Vector{Cluster})

    this = Tree()
    this.D = D
    this.c = c

    sepCount = 2^D - 1
    dataCount = length(clusters)
    featuresCount = length(clusters[1].lBounds)

    this.a = zeros(Float64, featuresCount, sepCount)
    this.b = zeros(Float64, sepCount)

    # Pour chaque sommet interne de l'arbre
    for t in 1:sepCount

        # Si une sépration y est effectuée
        if c[t] == -1

            # Déterminer les données qui vont à gauche ou droite
            rightData = Matrix{Float64}(undef, 0, featuresCount)
            leftData = Matrix{Float64}(undef, 0, featuresCount)

            for i in 1:dataCount
                if u[i, t*2] == 1
                    leftData = vcat(leftData, clusters[i].x)
                elseif u[i, t*2+1] == 1
                    rightData = vcat(rightData, clusters[i].x)
                end
            end

            len_l = size(leftData, 1)
            len_r = size(rightData, 1)

            m = Model(CPLEX.Optimizer)
            set_silent(m)
            @variable(m, a[1:featuresCount], base_name = "a_j")
            @variable(m, s[1:featuresCount], Bin, base_name = "s_j")
            @variable(m, b, base_name = "b")
            @variable(m, e[1:len_l+len_r], base_name = "ecart") # Ecart entre chaque variable et la séparation
            @variable(m, e_min, base_name = "ecart_min")

            @constraint(m, -1 <= b)
            @constraint(m, b <= 1)
            @constraint(m, [j in 1:featuresCount], -s[j] <= a[j])
            @constraint(m, [j in 1:featuresCount], a[j] <= s[j])
            @constraint(m, sum(s[j]  for j in 1:featuresCount) <= sum(s_model[j, t] for j in 1:featuresCount)) # on ne veut pas augmenter la "complexité" de la séparation
            @constraint(m, [i in 1:len_l + len_r], e[i] >= e_min)
            @constraint(m, [i in 1:len_l], e[i] == b - sum(a[j]*leftData[i, j] for j in 1:featuresCount))
            @constraint(m, [i in 1:len_r], e[i+len_l] == - b + sum(a[j]*rightData[i, j] for j in 1:featuresCount))

            @objective(m, Max, e_min)
            optimize!(m)

            this.b[t] = value.(b)
            for j in 1:featuresCount
                 this.a[j, t] = value.(a[j])
            end
        end
    end

    return this
end

max_stars_repo_name: boumpteryx/RODM

max_issues_repo_name: boumpteryx/RODM

max_forks_repo_name: boumpteryx/RODM


text:
#=
 The BUGS Book
 Example 2.4.1 Transformation

 Cf ~/jags/bugs_book_2_4_1.R
 model {
  Z ~ dnorm(0, 1)
  Y <- pow(2*Z + 1, 3)
  P10 <- step(Y - 10)
 }
       Mean      SD Naive SE Time-series SE
P10  0.2874  0.4526 0.001848       0.001839
Y   13.3067 39.9025 0.162901       0.162901
Z    0.0112  0.9993 0.004080       0.004025

  ~/webppl/bugs_book_2_4_2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_4_1()

    z ~ Normal(0,1)
    y = (2.0*z + 1.0)^3.0
    p10 ~ y > 10.0 ? flip(1.0) : flip(0.0)

    return p10
end

model = bugs_book_2_4_1()

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)
chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

display(chains)
# display(plot(chains))

gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module MyPkg

export testHelloWorld, addOne

"""
    testHelloWorld()

Prints 'Hello World!'
"""
function testHelloWorld()
    return "Hello World!"
end

"""
    addOne(int)

Adds one to a given integer
"""
function addOne(n)
    return n+1
end

end # module

max_stars_repo_name: jannittayao/MyPkg

max_issues_repo_name: jannittayao/MyPkg

max_forks_repo_name: jannittayao/MyPkg


text:
using LinearAlgebra

export euler_step

"""
    euler_step(x0, params, tspan, plf::Plefka, s_i, r_i, invFunc)

    Initializes the fields and runs the dynamics.
    Uses an Euler step integrator.

    Uses as input
    x0		Initial condition
    params	Reaction rates
    tspan	time grid vector
    plf		Plefka struct
    s_i,r_i Stoichiometric matrix, reactions in rows
    invFunc	Function, on which the expansion is based.

    If no function `invFunc` is specified, the default f(x)=1+x function
    will be used. Note that the defined function has to have the same
    Taylor expansion as f(x)=1+x.

    Returns a 3-element list with trajectories, responses and fields

"""
function euler_step(x0, params, tspan, plf::Plefka, s_i, r_i, invFunc)

	global num_species = length(params[1])

	len_ts = length(tspan)
	dt = tspan[2] - tspan[1]

	# Initialize trajectory
	y = zeros(num_species, len_ts)
	y[:,1] = x0

	# Initialize response functions
	resp = zeros(num_species, len_ts, len_ts)
	for i in 1:num_species
		resp[i,:,:] = Matrix{Int}(I, len_ts, len_ts)
	end

	# Initialize fields (depending on Plefka theory)
	fields = initialize_fields(plf, num_species, len_ts)

	# run dynamics
	for i in 1:len_ts-1

		# update fields
		if plf.alphaOrder == 1
			update_fields!(fields, s_i, r_i, params[3], i, dt, y)
		else
			update_fields!(fields, s_i, r_i, params[3], i, dt, y, resp)
		end

		# update response functions
		if plf.orderParam == "linear"
			update_responses!(resp, params, dt, i)
		else
			update_responses!(resp, params, dt, i, plf.α, fields, invFunc)
		end

		# Euler step
		y[:,i+1] = y[:,i] + dt * y_derivative(y, params, i, dt, plf.α, fields, invFunc)

	end

	trajectories = Result(tspan, y)
	responses = Responses(tspan, resp)

	return [trajectories, responses, fields]

end

function euler_step(x0, params, tspan, plf::Plefka, s_i, r_i)
	@. f_theta(x) = 1 + x
	return euler_step(x0, params, tspan, plf, s_i, r_i, f_theta)
end

function euler_step(p::Parameters, plf::Plefka)
	return euler_step(p.x0, p.k, tspan(p), plf, p.s_i, p.r_i)
end

function euler_step(p::Parameters, plf::Plefka, f_theta)
	return euler_step(p.x0, p.k, tspan(p), plf, p.s_i, p.r_i, f_theta)
end

"""
	initialize_fields(plf::Plefka, num_species, len_ts)

	Initializes the fields depending on the Plefka theory

	Returns a struct with the needed fields

"""
function initialize_fields(plf::Plefka, num_species, len_ts)

	if plf.orderParam == "linear" && plf.alphaOrder == 1
		hatTheta1 = zeros(num_species, len_ts)
		fields = Fields_lin1(hatTheta1)

	elseif plf.orderParam == "linear" && plf.alphaOrder == 2
		hatTheta1 = zeros(num_species, len_ts)
		hatTheta2 = zeros(num_species, len_ts)
		fields = Fields_lin2(hatTheta1, hatTheta2)

	elseif plf.orderParam == "quad" && plf.alphaOrder == 1
		hatTheta1 = zeros(num_species, len_ts)
		hatR1 = zeros(num_species, len_ts)
		fields = Fields_quad1(hatTheta1, hatR1)

	elseif plf.orderParam == "quad" && plf.alphaOrder == 2
		hatTheta1 = zeros(num_species, len_ts)
		hatTheta2 = zeros(num_species, len_ts)
		hatR1 = zeros(num_species, len_ts)
		hatR2 = zeros(num_species, len_ts, len_ts)
		fields = Fields_quad2(hatTheta1, hatTheta2, hatR1, hatR2)
	end

	return fields

end

"""
	update_fields!(fields::Fields_lin1, s_i, r_i, k3, t_i, dt, y)

	Updates the fields according to the theoretic equations
	for general reaction networks in the linear-1 case

	Returns the updated fields

"""
function update_fields!(fields::Fields_lin1, s_i, r_i, k3, t_i, dt, y)

	N = size(y,1)		# number of species
	len_ts = size(y,2)	# length of time grid
	z = zeros(Int, N)

	for i in 1:N
		fields.hatTheta1[i,t_i] = -c_mn(e_i(N,i), z, s_i, r_i, k3, y, t_i)
	end

end

"""
	update_fields!(fields::Fields_lin2, s_i, r_i, k3, t_i, dt, y, resp)

	Updates the fields according to the theoretic equations
	for general reaction networks in the linear-2 case

	Returns the updated fields

"""
function update_fields!(fields::Fields_lin2, s_i, r_i, k3, t_i, dt, y, resp)

	N = size(y,1)		# number of species
	len_ts = size(y,2)	# length of time grid
	z = zeros(Int, N)

	for i in 1:N

		# update Theta1
		fields.hatTheta1[i,t_i] = -c_mn(e_i(N,i), z, s_i, r_i, k3, y, t_i)

		# update Theta2
		sum_n = zeros(len_ts)
		for n in n_list(r_i)
			if sum(n) > 1
				prod_resp = calc_prod_resp(n, resp[:,t_i,:])
				sum_n += c_mn(e_i(N,i), n, s_i, r_i, k3, y, t_i)*c_mn(n, z, s_i, r_i, k3, y).*prod_resp
			end
		end
		fields.hatTheta2[i,t_i] = -2*dt*sum(sum_n)

	end

end

"""
	update_fields!(fields::Fields_quad1, s_i, r_i, k3, t_i, dt, y)

	Updates the fields according to the theoretic equations
	for general reaction networks in the quad-1 case

	Returns the updated fields

"""
function update_fields!(fields::Fields_quad1, s_i, r_i, k3, t_i, dt, y)

	N = size(y,1)		# number of species
	len_ts = size(y,2)	# length of time grid
	z = zeros(Int, N)

	for i in 1:N
		c_eiei = c_mn(e_i(N,i), e_i(N,i), s_i, r_i, k3, y, t_i)
		fields.hatTheta1[i,t_i] = -c_mn(e_i(N,i), z, s_i, r_i, k3, y, t_i) + c_eiei*y[i,t_i]
		fields.hatR1[i,t_i] = -c_eiei/dt
	end

end

"""
	update_fields!(fields::Fields_quad2, s_i, r_i, k3, t_i, dt, y, resp)

	Updates the fields according to the theoretic equations
	for general reaction networks in the quad-2 case

	Returns the updated fields

"""
function update_fields!(fields::Fields_quad2, s_i, r_i, k3, t_i, dt, y, resp)

	N = size(y,1)		# number of species
	len_ts = size(y,2)	# length of time grid
	z = zeros(Int, N)

	for i in 1:N

		# update Theta1 and hatR1
		c_eiei = c_mn(e_i(N,i), e_i(N,i), s_i, r_i, k3, y, t_i)
		fields.hatTheta1[i,t_i] = -c_mn(e_i(N,i), z, s_i, r_i, k3, y, t_i) + c_eiei*y[i,t_i]
		fields.hatR1[i,t_i] = -c_eiei/dt

		# update Theta2 and hatR2
		sum_T2 = 0.0
		local_R2 = 0.0
		sum_R2 = zeros(t_i)

		for n in n_list(r_i)
			c_ein = c_mn(e_i(N,i), n, s_i, r_i, k3, y, t_i)
			c_einei = c_mn(e_i(N,i), n.+e_i(N,i), s_i, r_i, k3, y, t_i)
			for t in 1:t_i
				prod_resp = calc_prod_resp(n, resp[:,t_i,t+1])
				if sum(n) > 1
					c_n0 = c_mn(n, z, s_i, r_i, k3, y, t)
					# hatTheta2
					sum_T2 += (c_ein - c_einei*y[i,t_i])*c_n0*prod_resp
					# local hatR2
					local_R2 += c_einei*c_n0*prod_resp
				end
				if sum(n) > 1 || (sum(n) == 1 && n[i] != 1)
					c_nei = c_mn(n, e_i(N,i), s_i, r_i, k3, y, t)
					# hatTheta2 with cein-cnei
					sum_T2 += -c_ein*c_nei*y[i,t]*prod_resp
					# Non-local hatR2
					sum_R2[t] += c_ein*c_nei*prod_resp
				end
			end
		end
		fields.hatTheta2[i,t_i] = -2*dt*sum_T2
		fields.hatR2[i,t_i+1,1:t_i] = -2*sum_R2
		fields.hatR2[i,t_i+1,t_i] = -2*local_R2
	end

end

"""
	update_responses!(resp, params, delta_t, t_i)

	Updates the responses for linear order parameters

	Returns the updated response functions

"""
function update_responses!(resp, params, delta_t, t_i)

	resp[:,t_i+1,:] = (1 .- params[2]*delta_t).*resp[:,t_i,:]
	resp[:,t_i+1,t_i+1] .= 1.0

end

"""
	update_responses!(resp, params, delta_t, t_i, α, fields::Fields_quad1, invFunc)

	Updates the responses for quadratic order parameters (alpha lin)

	Returns the updated response functions

"""
function update_responses!(resp, params, delta_t, t_i, α, fields::Fields_quad1, invFunc)

	for i in 1:size(resp,1)
		resp[i,t_i+1,:] = (1 - params[2][i]*delta_t - α*delta_t^2*fields.hatR1[i,t_i]) .* resp[i,t_i,:]
	end
	resp[:,t_i+1,t_i+1] .= 1.0

end

"""
	update_responses!(resp, params, delta_t, t_i, α, fields::Fields_quad2, invFunc)

	Updates the responses for quadratic order parameters (alpha quad)

	Returns the updated response functions

"""
function update_responses!(resp, params, delta_t, t_i, α, fields::Fields_quad2, invFunc)

	for i in 1:size(resp,1)
		for k in t_i:-1:1
			x = α*delta_t^2*fields.hatR1[i,t_i]*resp[i,t_i,k]
			xx = 0.5*α^2*delta_t^2*sum(fields.hatR2[i,t_i+1,:].*resp[i,:,k])
			x == 0 ? corr = 0 : corr = x * invFunc(xx/x)
			resp[i,t_i+1,k] = (1 - params[2][i]*delta_t)*resp[i,t_i,k] - corr
		end
		resp[i,t_i+1,t_i+1] = 1.0
	end

end

"""
	y_derivative(y, params, t_i, delta_t, α, fields::Fields_lin1, invFunc)

	calculates the dy/dt derivative for the euler integration

	returns dydt

"""
function y_derivative(y, params, t_i, delta_t, α, fields::Fields_lin1, invFunc)

	dydt = params[1] .- params[2] .* y[:,t_i] .- α*fields.hatTheta1[:,t_i]
	return dydt

end

"""
	y_derivative(y, params, t_i, delta_t, α, fields::Fields_lin2, invFunc)

	calculates the dy/dt derivative for the euler integration

	returns dydt

"""
function y_derivative(y, params, t_i, delta_t, α, fields::Fields_lin2, invFunc)

	dydt = zeros(size(y,1))
	for i in 1:size(y,1)
		t1 = α*fields.hatTheta1[i,t_i]
		t2 = 0.5*α^2*fields.hatTheta2[i,t_i]
		t1 == 0.0 ? x = 0 : x = t1 * invFunc(t2/t1)
		dydt[i] = params[1][i] - params[2][i] * y[i,t_i] - x
	end
	return dydt

end

"""
	y_derivative(y, params, t_i, delta_t, α, fields::Fields_quad1, invFunc)

	calculates the dy/dt derivative for the euler integration

	returns dydt

"""
function y_derivative(y, params, t_i, delta_t, α, fields::Fields_quad1, invFunc)

	t1 = α*fields.hatTheta1[:,t_i]
	r1 = α*delta_t*y[:,t_i].*fields.hatR1[:,t_i]
	dydt = params[1] .- params[2] .* y[:,t_i] .- t1 .- r1
	return dydt

end

"""
	y_derivative(y, params, t_i, delta_t, α, fields::Fields_quad2, invFunc)

	calculates the dy/dt derivative for the euler integration

	returns dydt

"""
function y_derivative(y, params, t_i, delta_t, α, fields::Fields_quad2, invFunc)

	dydt = zeros(size(y,1))
	for i in 1:size(y,1)
		t1 = α*fields.hatTheta1[i,t_i]
		t2 = 0.5*α^2*fields.hatTheta2[i,t_i]
		r1 = α*delta_t*y[i,t_i].*fields.hatR1[i,t_i]
		r2 = 0.5*α^2*delta_t*sum(y[i,:].*fields.hatR2[i,t_i+1,:])
		t1+r1 == 0.0 ? x = 0 : x = (t1+r1) * invFunc((t2+r2)/(t1+r1))
		dydt[i] = params[1][i] - params[2][i] * y[i,t_i] - x
	end
	return dydt

end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
using ProductizeBayes, Test

TestDir = @__DIR__

test_versioned = [
    "test_versioned_model/test_versioned_model.jl"
  ]

  @testset "VersionedModel basic tests" begin
    for test in test_versioned
      println("\nTesting: $test.")
      include(joinpath(TestDir, test))
    end
    println()
  end


text:
@static if VERSION >= v"1.5-"
    function _foreach_thread(f)
        ids = zeros(Threads.nthreads())
        Threads.@threads :static for i in 1:Threads.nthreads()
            ids[i] = Threads.threadid()
            f()
        end
        return check_threadids(ids)
    end
else
    function _foreach_thread(f)
        ids = zeros(Threads.nthreads())
        Threads.@threads for i in 1:Threads.nthreads()
            ids[i] = Threads.threadid()
            f()
        end
        return check_threadids(ids)
    end
end

check_threadids(ids) = length(Set(ids)) == length(ids)

function nonsticky!(task)
    task.sticky = false
    return task
end

const SIMDFlag = Union{Bool, Symbol, Val{true}, Val{false}, Val{:ivdep}}

struct Err{T}
    value::T
end

struct Ok{T}
    value::T
end

Base.getindex(err::Err) = throw(err.value)
Base.getindex(ok::Ok) = ok.value

struct Promise
    value::Base.RefValue{Any}
    isset::Threads.Atomic{Bool}
    notify::Threads.Condition
end

Promise() = Promise(Ref{Any}(), Threads.Atomic{Bool}(false), Threads.Condition())

tryfetch(::Nothing) = nothing
tryfetch(p::Promise) = p.isset[] ? Some(p.value[]) : nothing

function tryput!(p::Promise, value)
    p.isset[] && return Some(p.value[])
    lock(p.notify) do
        p.isset[] && return Some(p.value[])
        p.value[] = value
        p.isset[] = true
        notify(p.notify)
        return nothing
    end
end

function Base.fetch(p::Promise)
    p.isset[] && return p.value[]
    lock(p.notify) do
        while !p.isset[]
            wait(p.notify)
        end
        return p.value[]
    end
end

struct Future
    f::Any
    p::Promise
end

Future(@nospecialize f) = Future(f, Promise())

function run!(future::Future)
    ans = try
        Ok(future.f())
    catch err
        @debug("FUTURE FAILED", objectid(future), future.f, exception = (err, catch_backtrace()))
        Err(err)
    end
    @assert tryput!(future.p, ans) === nothing
end

Base.fetch(future::Future) = fetch(future.p)[]

runner(f::Future) = function run_future()
    run!(f)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

# Let's drop here all integration schemes and some defaults for different element types maybe parse from txt file ..?

### Gauss quadrature rules for one dimension

function get_integration_points(::Type{Val{1}})
    return [2.0], [0.0]
end

function get_integration_points(::Type{Val{2}})
    return [1.0, 1.0], sqrt(1.0/3.0)*[-1.0, 1.0]
end

function get_integration_points(::Type{Val{3}})
    return 1.0/9.0*[5.0, 8.0, 5.0], sqrt(3.0/5.0)*[-1.0, 0.0, 1.0]
end

function get_integration_points(::Type{Val{4}})
    weights = 1.0/36.0*[
        18.0-sqrt(30.0),
        18.0+sqrt(30.0),
        18.0+sqrt(30.0),
        18.0-sqrt(30.0)]
    points = [
        -sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)),
        -sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
         sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
         sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0))]
    return weights, points
end

function get_integration_points(::Type{Val{5}})
    weights = [
        1.0/900.0*(322.0-13.0*sqrt(70.0)),
        1.0/900.0*(322.0+13.0*sqrt(70.0)),
        128.0/225.0,
        1.0/900.0*(322.0+13.0*sqrt(70.0)),
        1.0/900.0*(322.0-13.0*sqrt(70.0))]
    points = [
        -1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0)),
        -1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)),
         0.0,
         1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)),
         1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0))]
    return weights, points
end

function get_integration_points(order::Int64)
    if order <= 5
        return get_integration_points(Val{order})
    else
        points, weights = Base.QuadGK.gauss(Float64, order)
        return weights, points
    end
end

### "cartesian" elements, integration rules comes from tensor product

function get_integration_points(element::Poi1)
    [ (1.0, [] ) ]
end

typealias CartesianLineElement Union{Seg2, Seg3, NSeg}
typealias CartesianSurfaceElement Union{Quad4, Quad8, Quad9, NSurf}
typealias CartesianVolumeElement Union{Hex8, Hex20, Hex27, NSolid}

function get_integration_points(element::CartesianLineElement, order::Int64)
    w, xi = get_integration_points(order)
    [ (w[i], [xi[i]]) for i=1:order ]
end

function get_integration_points(element::CartesianSurfaceElement, order::Int64)
    w, xi = get_integration_points(order)
    [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:order, j=1:order ]
end

function get_integration_points(element::CartesianVolumeElement, order::Int64)
    w, xi = get_integration_points(order)
    [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:order, j=1:order, k=1:order ]
end

### triangular and tetrahedral elements

# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html

typealias TriangularElement Union{Tri3, Tri6, Tri7}

function get_integration_points(element::TriangularElement, ::Type{Val{1}})
    weights = [0.5]
    points = Vector{Float64}[1.0/3.0*[1.0, 1.0]]
    return zip(weights, points)
end

function get_integration_points(element::TriangularElement, ::Type{Val{2}})
    weights = 1.0/6.0*[1.0, 1.0, 1.0]
    points = Vector{Float64}[
        [2.0/3.0, 1.0/6.0],
        [1.0/6.0, 2.0/3.0],
        [1.0/6.0, 1.0/6.0]]
    return zip(weights, points)
end

function get_integration_points(element::TriangularElement, ::Type{Val{3}})
    weights = 0.5*[
        -0.5625,
         0.5208333333333333,
         0.5208333333333333,
         0.5208333333333333]
    points = Vector{Float64}[
        [1.0/3.0, 1.0/3.0],
        [0.2, 0.2],
        [0.2, 0.6],
        [0.6, 0.2]]
    return zip(weights, points)
end

function get_integration_points(element::TriangularElement, ::Type{Val{4}})
    weights = 0.5*[
        0.22338158967801,
        0.22338158967801,
        0.22338158967801,
        0.10995174365532,
        0.10995174365532,
        0.10995174365532]
    points = Vector{Float64}[
        [0.44594849091597, 0.44594849091597],
        [0.44594849091597, 0.10810301816807],
        [0.10810301816807, 0.44594849091597],
        [0.09157621350977, 0.09157621350977],
        [0.09157621350977, 0.81684757298046],
        [0.81684757298046, 0.09157621350977]]
    return zip(weights, points)
end

function get_integration_points(element::TriangularElement, ::Type{Val{5}})
    weights = 0.5*[
        0.22500000000000,
        0.13239415278851,
        0.13239415278851,
        0.13239415278851,
        0.12593918054483,
        0.12593918054483,
        0.12593918054483]
    points = Vector{Float64}[
        [0.33333333333333, 0.33333333333333],
        [0.47014206410511, 0.47014206410511],
        [0.47014206410511, 0.05971587178977],
        [0.05971587178977, 0.47014206410511],
        [0.10128650732346, 0.10128650732346],
        [0.10128650732346, 0.79742698535309],
        [0.79742698535309, 0.10128650732346]]
    return zip(weights, points)
end

### 3d elements

typealias TetrahedralElement Union{Tet4, Tet10}

function get_integration_points(element::TetrahedralElement, ::Type{Val{1}})
    weights = 1.0/6.0*[1.0]
    points = Vector{Float64}[
        1.0/4.0*[1.0, 1.0, 1.0]]
    return zip(weights, points)
end

function get_integration_points(element::TetrahedralElement, ::Type{Val{2}})
    a = (5.0+3.0*sqrt(5.0))/20.0
    b = (5.0-sqrt(5.0))/20.0
    w = 1.0/24.0
    weights = [w, w, w, w]
    points = Vector{Float64}[
        [a, b, b],
        [b, a, b],
        [b, b, a],
        [b, b, b]]
    return zip(weights, points)
end

function get_integration_points(element::TetrahedralElement, ::Type{Val{3}})
    a = 1.0/4.0
    b = 1.0/6.0
    c = 1.0/2.0
    weights = [-2.0/15.0, 3.0/40.0, 3.0/40.0, 3.0/40.0, 3.0/40.0]
    points = Vector{Float64}[
        [a, a, a],
        [b, b, b],
        [b, b, c],
        [b, c, b],
        [c, b, b]]
    return zip(weights, points)
end

function get_integration_points(element::TetrahedralElement, ::Type{Val{4}})
    a = 0.25
    b1 = 1.0/34.0*(7.0 + sqrt(15.0))
    b2 = 1.0/34.0*(7.0 - sqrt(15.0))
    c1 = 1.0/34.0*(13.0 - 3.0*sqrt(15.0))
    c2 = 1.0/34.0*(13.0 + 3.0*sqrt(15.0))
    d = 1.0/20.0*(5.0 - sqrt(15.0))
    e = 1.0/20.0*(5.0 + sqrt(15.0))
    w1 = 8.0/405.0
    w2 = (2665.0 - 14.0*sqrt(15.0))/226800.0
    w3 = (2665.0 + 14.0*sqrt(15.0))/226800.0
    w4 = 5.0/567.0
    weights = [w1, w2, w2, w2, w2, w3, w3, w3, w3, w4, w4, w4, w4, w4]
    points = Vector{Float64}[
        [a, a, a],
        [b1, b1, b1],
        [b1, b1, c1],
        [b1, c1, b1],
        [c1, b1, b1],
        [b2, b2, b2],
        [b2, b2, c2],
        [b2, c2, b2],
        [c2, b2, b2],
        [d, d, e],
        [d, e, d],
        [e, d, d],
        [d, e, e],
        [e, d, e],
        [e, e, d]]
    return zip(weights, points)
end

typealias PrismaticElement Union{Wedge6}

function get_integration_points(element::PrismaticElement, ::Type{Val{2}})
    weights = 1/6*[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    points = Vector{Float64}[
        [0.5, 0.0, -1.0/sqrt(3)],
        [0.0, 0.5, -1.0/sqrt(3)],
        [0.5, 0.5, -1.0/sqrt(3)],
        [0.5, 0.0,  1.0/sqrt(3)],
        [0.0, 0.5,  1.0/sqrt(3)],
        [0.5, 0.5,  1.0/sqrt(3)]]
    return zip(weights, points)
end

function get_integration_points(element::Union{TriangularElement,
            TetrahedralElement, PrismaticElement}, order::Int64)
    return get_integration_points(element, Val{order})
end

### default number of integration points for each element
### 2 for linear elements, 3 for quadratic

typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Wedge6, Hex8}

typealias QuadraticElement Union{Seg3, Tri6, Tri7, Tet10, Quad8, Quad9, Hex20, Hex27}

function get_integration_order(element::LinearElement)
    return 2
end

function get_integration_order(element::QuadraticElement)
    return 3
end

function get_integration_points(element::LinearElement)
    order = get_integration_order(element)
    get_integration_points(element, order)
end

function get_integration_points(element::QuadraticElement)
    order = get_integration_order(element)
    get_integration_points(element, order)
end