text:
#############################################################################
# logsumexp.jl
# log of sum of exponentials of each entry in an expression
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

### LogSumExp

# TODO: make this work for a *list* of inputs, rather than just for vector/matrix inputs

struct LogSumExpAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function LogSumExpAtom(x::AbstractExpr)
        if sign(x) == ComplexSign()
            error("The argument should be real but it's instead complex")
        else
            children = (x,)
            return new(:logsumexp, hash(children), children, (1,1))
        end
    end
end

function sign(x::LogSumExpAtom)
    return NoSign()
end

function monotonicity(x::LogSumExpAtom)
    return (Nondecreasing(),)
end

function curvature(x::LogSumExpAtom)
    return ConvexVexity()
end

function evaluate(x::LogSumExpAtom)
    return log(sum(exp.(evaluate(x.children[1]))))
end

logsumexp(x::AbstractExpr) = LogSumExpAtom(x)

function conic_form!(e::LogSumExpAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, e)
        # log(sum(exp(x))) <= t  <=>  sum(exp(x)) <= exp(t) <=> sum(exp(x - t)) <= 1
        t = Variable(e.size)
        z = sum(exp(e.children[1] - t))
        objective = conic_form!(t, unique_conic_forms)
        conic_form!(z, unique_conic_forms)
        conic_form!(z <= 1, unique_conic_forms)

        cache_conic_form!(unique_conic_forms, e, objective)
    end
    return get_conic_form(unique_conic_forms, e)
end

function logisticloss(e::AbstractExpr)
    s = 0
    length(e) == 1 && return logsumexp([e, 0])
    for i = 1:length(e)
        s += logsumexp([e[i]; 0])
    end
    return s
end


text:
module TopoChains

using Functors

export FuncTopo, @functopo_str, @functopo
export TopoChain

include("code.jl")
include("topology.jl")
include("topochain.jl")

Functors.functor(s::TopoChain) = s.models, m -> TopoChain(s.topo, m...)

end

ext: jl


text:
########################################################################################
# Algorithm
########################################################################################
abstract type GaussianODEFilter <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm end
abstract type AbstractEK <: GaussianODEFilter end

"""
    EK0(; prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true)

**Gaussian ODE filtering with zeroth order extended Kalman filtering.**

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[DynamicDiffusion(), DynamicMVDiffusion(), FixedDiffusion(), FixedMVDiffusion()]`.

See also: [`EK1`](@ref)

# References:
- N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective** (2019)
- M. Schober, S. Särkkä, and P. Hennig: **A Probabilistic Model for the Numerical Solution of Initial Value Problems** (2018)
"""
Base.@kwdef struct EK0{DT,IT} <: AbstractEK
    prior::Symbol = :ibm
    order::Int = 3
    diffusionmodel::DT = DynamicDiffusion()
    smooth::Bool = true
    initialization::IT = TaylorModeInit()
end

"""
    EK1(; prior=:ibm, order=3, diffusionmodel=DynamicDiffusion(), smooth=true)

**Gaussian ODE filtering with first order extended Kalman filtering.**

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[DynamicDiffusion(), DynamicMVDiffusion(), FixedDiffusion(), FixedMVDiffusion()]`.

See also: [`EK0`](@ref)

# References:
- N. Bosch, P. Hennig, F. Tronarp: **Calibrated Adaptive Probabilistic ODE Solvers** (2021)
- F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig: **Probabilistic Solutions To Ordinary Differential Equations As Non-Linear Bayesian Filtering: A New Perspective** (2019)
"""
struct EK1{CS,AD,DiffType,DT,IT} <: AbstractEK
    prior::Symbol
    order::Int
    diffusionmodel::DT
    smooth::Bool
    initialization::IT
end
EK1(;
    prior=:ibm,
    order=3,
    diffusionmodel=DynamicDiffusion(),
    smooth=true,
    initialization=TaylorModeInit(),
    chunk_size=0,
    autodiff=true,
    diff_type=Val{:forward},
) = EK1{chunk_size,autodiff,diff_type,typeof(diffusionmodel),typeof(initialization)}(
    prior,
    order,
    diffusionmodel,
    smooth,
    initialization,
)

Base.@kwdef struct EK1FDB{DT,IT} <: AbstractEK
    prior::Symbol = :ibm
    order::Int = 3
    diffusionmodel::DT = DynamicDiffusion()
    smooth::Bool = true
    initialization::IT = TaylorModeInit()
    jac_quality::Int = 1
end


text:
# -------- Miscellaneous ----------------------------------------------------- #

second((a, b)) = b

unzip(a :: Array{<: Tuple}) = map(first, a), map(second, a)

function choose_index(probs :: Vector{Float32}) :: Int

  @assert all(probs .>= 0) && sum(probs) ≈ 1.0 "probability vector not proper"

  r = rand(Float32)
  index = findfirst(x -> r <= x, cumsum(probs))

  isnothing(index) ? length(probs) : index

end

function derive_gametype(players)

  gt = mapreduce(gametype, typeintersect, players, init = Game)

  @assert gt != Union{} "Players do not play compatible games"
  @assert !isabstracttype(gt) "Cannot infere game from abstract type"

  gt

end

isasync(m) = isa(m, Async) ? true : false

# Use the logit of a normal distribution instead of Dirichlet for diluting the
# root prior policy
function logit_normal(n)
  ey = exp.(randn(Float32, n))
  ey / sum(ey)
end

function one_hot(n, k)
  r = zeros(Float32, n)
  r[k] = 1f0
  r
end

get_id(args...) = Int(div(hash(tuple(args...)), Int(1e14)))

# -------- Symmetry ---------------------------------------------------------- #

cind(i, s) = Tuple(CartesianIndices(s)[i])
lind(c, s) = LinearIndices(s)[c...]

hmirror(matrix) = reshape(matrix[end:-1:1, :, :], size(matrix))
hmirror(a :: Tuple, s) = (s[1]+1-a[1], a[2:end]...)
hmirror(a :: Int, s) = lind(hmirror(cind(a, s), s), s)

vmirror(matrix) = reshape(matrix[:, end:-1:1, :], size(matrix))
vmirror(a :: Tuple, s) = (a[1], s[2]+1-a[2], a[3:end]...)
vmirror(a :: Int, s) = lind(vmirror(cind(a, s), s), s)

dmirror(matrix :: Matrix) = permutedims(matrix, (2, 1))
dmirror(matrix :: Array{T,3}) where {T} = permutedims(matrix, (2, 1, 3))
dmirror(a :: Tuple, s) = (a[2], a[1], a[3:end]...)
dmirror(a :: Int, s) = lind(dmirror(cind(a, s), s), s)

# Full symmetry group of square matrices

function apply_dihedral_group(matrix :: Array)
  [
    matrix |> copy,
    matrix |> hmirror,
    matrix |> dmirror,
    matrix |> vmirror,
    matrix |> dmirror |> hmirror,
    matrix |> hmirror |> dmirror,
    matrix |> vmirror |> hmirror,
    matrix |> vmirror |> hmirror |> dmirror
  ]
end

function apply_dihedral_group(action :: Int, s :: Tuple)
  [
    action,
    hmirror(action, s),
    dmirror(action, s),
    vmirror(action, s),
    hmirror(dmirror(action, s), s),
    dmirror(hmirror(action, s), s),
    hmirror(vmirror(action, s), s),
    dmirror(hmirror(vmirror(action, s), s), s)
  ]
end

# Symmetry group for non-square matrices

function apply_klein_four_group(matrix :: Array)
  [
    matrix |> copy,
    matrix |> hmirror,
    matrix |> vmirror,
    matrix |> vmirror |> hmirror
  ]
end

function apply_klein_four_group(action :: Int, s :: Tuple)
  [
    action,
    hmirror(action, s),
    vmirror(action, s),
    vmirror(hmirror(action, s), s)
  ]
end

# -------- Pretty Printing to Monitor Training ------------------------------- #

const gray_crayon = Crayons.Crayon(foreground = :dark_gray)
const default_crayon = Crayons.Crayon(reset = true)

function print_loss_header(loss, use_features)

  names = use_features ? caption(loss) : caption(loss)[1:3]

  components = map(names) do c
    Printf.@sprintf("%10s", string(c)[1:min(end, 10)])
  end

  println(join(["#"; "   epoch"; components; "     total"; "    length"], " "))

end

function print_loss(l, p, epoch, train, test)

  for (set, cr) in [(train, gray_crayon), (test, default_crayon)]

    # Compute the losses and get them as strings
    ls = loss(l, training_model(p), set)
    losses = map(x -> @sprintf("%10.3f", x), ls)

    # Print everything in grey (for train) and white (for test)
    print(cr)
    @printf( "%10d %s %10.3f %10d\n"
           , epoch
           , join(losses, " ")
           , sum(ls)
           , length(set) )

  end

end

function print_ranking(rk)

  # Log that a contest comes next
  print(gray_crayon)
  println("#\n# Contest with $(length(rk.players)) players:\n#")

  # Get the summary of the contest and print it
  s = "# " * replace(summary(rk, true), "\n" => "\n# ") * "\n#"
  println(s)

end

# -------- Feature Compability ----------------------------------------------- #

function check_features(l, model)

  features(l) == features(model) && !isempty(features(model))

end

function check_features(l, model, dataset)

  fl = features(l)
  fm = features(model)
  fd = features(dataset)

  if isempty(fm)        return false
  elseif fm == fl == fd return true
  elseif fm == fd       @warn "Loss does not support model features"
  elseif fm == fl       @warn "Dataset does not support model features"
  else                  @warn "Loss and dataset do not support model features"
  end

  false

end

# ------- Switch Knet Allocator ---------------------------------------------- #

# TODO: This is an ugly hack that should be discussed with the Knet people.
a = true
Knet.cuallocator() = a
switch_knet_allocator() = (global a; a = !a)

# -------- Parallel-Stable Progress Maps ------------------------------------- #

"""
    stepper(description, length; <keyword arguments>)

Return a pair of functions `(step, finish)` that control a progress bar with
title `description` and `length` number of steps. For the keyword arguments, see
`ProgressMeter.Progress`.

The progress bar is running in its own thread and is communicated with through
a `RemoteChannel`. This allows calling the `step` function safely under any
condition. However, you **must** make sure that `step()` is not used after
`finish()` anymore. This means that all tasks or futures have to be fetched
before calling `finish`.

# Examples

```julia
step, finish = stepper("Progress", 15)

# Do some calculations in a parallel fashion
pmap(1:15) do i
  sleep(0.5)
  step()
  i
end

# Since pmap fetches all tasks/futures it creates, we can safely finish
finish()
```
"""
function stepper(description, n :: Int, kwargs...)

  glyphs = ProgressMeter.BarGlyphs("[=>⋅]")
  progress = ProgressMeter.Progress( n + 1
                                   ; dt = 0.1
                                   , desc = description
                                   , barglyphs = glyphs
                                   , kwargs... )

  # Remote channel that can be used to signal a step on any process
  channel = RemoteChannel(() -> Channel{Bool}(0))

  # Thread that manages progressing the progress bar
  thread = @async begin

    while take!(channel)
      ProgressMeter.next!(progress)
    end

    ProgressMeter.printover(progress.output, "")

  end

  step = () -> (put!(channel, true); yield())
  finish = () -> (put!(channel, false); fetch(thread))

  step, finish

end

# -------- Random Branching -------------------------------------------------- #

prepare(; steps = 0) = game -> random_turns!(copy(game), steps)

function branch(; prob = 0., steps = 1)
  game -> rand() < prob ? random_turns!(copy(game), steps) : nothing
end

# -------- Neural Network Head Creation -------------------------------------- #

function prepare_head(head, s, l, gpu)

  if isnothing(head)

    head = Dense(prod(s), l, gpu = gpu)

  else

    @assert valid_insize(head, s) "Head incompatible with trunk."
    @assert prod(outsize(head, s)) == l "Head incompatible with game."
    head = (on_gpu(head) == gpu) ? head : swap(head)

  end

  head

end

max_stars_repo_name: roSievers/jtac


text:
hangman1 = """
                 _____
                |     |
                      |
                      |
                      |
                     _|_"""

hangman2 = """
                 _____
                |     |
                O     |
                      |
                      |
                     _|_"""

hangman3 = """
                 _____
                |     |
                O     |
                |     |
                |     |
                     _|_"""

hangman4 = """
                 _____
                |     |
                O     |
               /|     |
                |     |
                     _|_"""

hangman5 = """
                 _____
                |     |
                O     |
               /|\\    |
                |     |
                     _|_"""

hangman6 = """
                 _____
                |     |
                O     |
               /|\\    |
                |     |
               /     _|_"""

hangman7 = """
                 _____
                |     |
                O     |
               /|\\    |
                |     |
               / \\   _|_"""

hangman = [hangman1, hangman2, hangman3, hangman4, hangman5, hangman6, hangman7]

max_stars_repo_name: fungktt/JuliaGames

max_issues_repo_name: fungktt/JuliaGames

max_forks_repo_name: fungktt/JuliaGames


text:
# include("C:/Users/joaquimgarcia/.julia/dev/MathOptInterface/test/Utilities/DoubleDicts.jl")
# include("C:/Users/joaquimgarcia/.julia/dev/MathOptInterface/src/Utilities/DoubleDicts.jl")

using MathOptInterface, Test
const MOI = MathOptInterface
const CI{F,S} = MOI.ConstraintIndex{F,S}

const DoubleDicts = MathOptInterface.Utilities.DoubleDicts

const CI_I = CI{MOI.SingleVariable, MOI.Integer}
const CI_B = CI{MOI.SingleVariable, MOI.ZeroOne}

function test_iterator(dict)
    kk, vv = [], []
    for (k,v) in dict
        push!(kk, k)
        push!(vv, v)
    end
    @test length(kk) == length(keys(dict))
    @test isempty(setdiff(kk, keys(dict)))
    @test length(vv) == length(values(dict))
    @test isempty(setdiff(vv, values(dict)))
end

function basic_functionality(dict, k_values, v_values)
    @test_throws ErrorException sizehint!(dict, 1)

    @test isempty(dict)
    @test length(dict) == 0

    dict[k_values[1]] = v_values[1]
    delete!(dict, k_values[1])
    test_iterator(dict)

    dict[k_values[3]] = v_values[3]
    test_iterator(dict)

    empty!(dict)
    @test isempty(dict)
    @test !haskey(dict, k_values[1])

    dict[k_values[1]] = v_values[1]
    @test haskey(dict, k_values[1])
    @test values(dict) == [v_values[1]]
    @test keys(dict) == [k_values[1]]
    @test dict[k_values[1]] == v_values[1]

    test_iterator(dict)

    delete!(dict, k_values[1])
    @test !haskey(dict, k_values[1])
    @test isempty(dict)

    dict[k_values[1]] = v_values[1]
    @test length(dict) == 1

    for (k,v) in zip(k_values, v_values)
        dict[k] = v
    end

    test_iterator(dict)

    empty!(dict)
    @test isempty(dict)
    @test length(dict) == 0

    dict[k_values[1]] = v_values[1]
    idict = DoubleDicts.with_type(dict, MOI.SingleVariable, MOI.Integer)
    bdict = DoubleDicts.with_type(dict, MOI.SingleVariable, MOI.ZeroOne)
    idict_ = dict[MOI.SingleVariable, MOI.Integer]
    @test idict.dict === idict_.dict
    sizehint!(idict, 2)
    @test length(idict) == 1
    @test length(bdict) == 0
    @test values(bdict) == typeof(v_values[3])[]
    @test keys(bdict) == typeof(k_values[3])[]

    @test haskey(idict, k_values[1])
    @test !haskey(idict, k_values[2])
    @test !haskey(bdict, k_values[3])

    @test values(idict) == [v_values[1]]
    @test keys(idict) == [k_values[1]]
    @test idict[k_values[1]] == v_values[1]

    idict[k_values[2]] = v_values[2]
    @test haskey(idict, k_values[2])
    test_iterator(idict)

    delete!(idict, k_values[2])
    @test !haskey(idict, k_values[2])
    test_iterator(idict)

    @test !isempty(idict)
    @test isempty(bdict)
    empty!(idict)
    empty!(bdict)
    @test isempty(idict)
    @test isempty(bdict)

    bdict[k_values[3]] = v_values[3]
    length(bdict) == 1

    edict = DoubleDicts.with_type(dict, MOI.SingleVariable, MOI.EqualTo{Bool})
    ek = CI{MOI.SingleVariable, MOI.EqualTo{Bool}}(1)
    delete!(edict, ek)
    @test_throws KeyError edict[ek]
    sizehint!(edict, 0)
    @test length(edict) == 0
    @test_throws KeyError edict[ek]
    delete!(edict, ek)
    test_iterator(edict)
end

@testset "DoubleDict" begin
    dict = DoubleDicts.DoubleDict{Float64}()
    keys = [
        CI_I(1),
        CI_I(2),
        CI_B(1),
    ]
    vals = [
        1.0,
        2.0,
        1.0,
    ]
    basic_functionality(dict, keys, vals)
end

@testset "IndexDoubleDict" begin
    dict = DoubleDicts.IndexDoubleDict()
    keys = [
        CI_I(1),
        CI_I(2),
        CI_B(1),
    ]
    vals = keys
    basic_functionality(dict, keys, vals)

    src = DoubleDicts.IndexDoubleDict()
    for (k,v) in zip(keys, vals)
        dict[k] = v
    end
    dest = DoubleDicts.IndexDoubleDict()
    MOIU._reverse_dict(dest, src)
    for (k, v) in src
        @test dest[v] == k
    end
end

@testset "FunctionSetDoubleDict" begin
    dict = DoubleDicts.FunctionSetDoubleDict()
    keys = [
        CI_I(1),
        CI_I(2),
        CI_B(1),
    ]
    vals = [
        (MOI.SingleVariable(MOI.VariableIndex(1)), MOI.Integer()),
        (MOI.SingleVariable(MOI.VariableIndex(2)), MOI.Integer()),
        (MOI.SingleVariable(MOI.VariableIndex(1)), MOI.ZeroOne()),
    ]
    basic_functionality(dict, keys, vals)
end


text:
using Test

@testset "Element" begin
    @test configuration(n"Fe") == "1s² 2s² 2p⁶ 3s² 3p⁶ 4s² 3d⁶"
    @test element(20) == n"Ca"
    @test a(n"Ca")==40.0784
    @test all(zz->z(elements[zz])==zz,1:92)
    @test symbol(n"Pd")=="Pd"
    @test density(n"Pb")==11.34
    @test atomic_weight[n"Pb"].lo == 206.14
    @test atomic_weight[n"Pb"].hi == 207.94
    @test isapprox(atomic_weight[n"Ca"], UncertainValue(40.078, 0.004), atol=0.00001)
end


text:
using Knet: @diff, Param, value, grad
using Knet: sigm, tanh

relu(x) = max(0,x)

mutable struct Layer
    wf1::Param
    wf2::Param
    wfs::Param
    wi1::Param
    wi2::Param
    wk1::Param
    wk2::Param
    wks::Param
end

Layer(in_size1, in_size2, layer_size) =
begin
    sq = sqrt(2/(in_size1+in_size2+layer_size))
    wf1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wf2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wfs = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    wi1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wi2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wk1 = Param(2*sq .* randn(in_size1,layer_size)   .-sq)
    wk2 = Param(2*sq .* randn(in_size2,layer_size)   .-sq)
    wks = Param(2*sq .* randn(layer_size,layer_size) .-sq)
    layer = Layer(wf1, wf2, wfs, wi1, wi2, wk1, wk2, wks)
layer
end

(layer::Layer)(state, in_1, in_2) =
begin
    attention = sigm.(in_1 * layer.wf1 + in_2 * layer.wf2 + state * layer.wfs)
    remember  = sigm.(in_1 * layer.wk1 + in_2 * layer.wk2 + state * layer.wks)
    interm    = tanh.(in_1 * layer.wi1 + in_2 * layer.wi2 + state .* attention)
    state = remember .* state + (1 .- remember) .* interm
[state, state]
end

struct IS
    layer_1::Layer
    layer_out::Layer
end

IS(storage_size, hm_vectors, vector_size, layer_sizes, out_size) =
begin
    is = IS(
        Layer(storage_size, hm_vectors*vector_size, layer_sizes[1]),
        Layer(storage_size, layer_sizes[end], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[is, state]
end

(is::IS)(states, in_1, in_2) =
begin
    out_1, state_1 = is.layer_1(states[1], in_1, in_2)
    out, state_out = is.layer_out(states[end], in_1, out_1)
[out, [state_1, state_out]]
end

struct GS
    layer_1::Layer
    layer_out::Layer
end

GS(in_size1, in_size2, layer_sizes, out_size) =
begin
    gs = GS(
        Layer(in_size1, in_size2, layer_sizes[1]),
        Layer(in_size1, layer_sizes[1], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[gs, state]
end

(gs::GS)(states, in_1, in_2) =
begin
    out_1, state_1 = gs.layer_1(states[1], in_1, in_2)
    out, state_out = gs.layer_out(states[2], in_1, out_1)
[out, [state_1, state_out]]
end

struct GO
    layer_1::Layer
    layer_out::Layer
end

GO(in_size1, in_size2, layer_sizes, out_size) =
begin
    go = GO(
        Layer(in_size1, in_size2, layer_sizes[1]),
        Layer(in_size1, layer_sizes[1], out_size),
    )
    state = [zeros(1, layer_sizes[1]),
             zeros(1, out_size)]
[go, state]
end

(go::GO)(states, in_1, in_2) =
begin
    out_1, state_1 = go.layer_1(states[1], in_1, in_2)
    out, state_out = go.layer_out(states[2], in_1, out_1)
[out, [state_1, state_out]]
end

struct ENCODER
    is::IS
    gs::GS
    go::GO
end

ENCODER(hm_vectors, vector_size, storage_size, is_layers, gs_layers, go_layers) =
begin
    is, is_state = IS(storage_size, hm_vectors, vector_size, is_layers, storage_size)
    gs, gs_state = GS(storage_size, storage_size, gs_layers, storage_size)
    go, go_state = GO(storage_size, vector_size, go_layers, vector_size)
    encoder = ENCODER(is, gs, go)
    state   = [is_state,gs_state,[go_state for i in 1:hm_vectors]]
[encoder, state]
end

(encoder::ENCODER)(state, storage, vectors) =
begin
    is_out, is_state = encoder.is(state[1], storage, hcat(vectors...))
    gs_out, gs_state = encoder.gs(state[2], storage, is_out)
    go_out, go_state = [], []
    for (vector, state_v) in zip(vectors, state[3])
        go_out_v, go_state_v = encoder.go(state_v, is_out, vector)
        push!(go_out, go_out_v)
        push!(go_state, go_state_v)
    end
    output = [gs_out, go_out]
    state  = [is_state, gs_state, go_state]
[output, state]
end

struct DECODER
    is::IS
    gs::GS
    go::GO
end

DECODER(hm_vectors, vector_size, storage_size, is_layers, gs_layers, go_layers) =
begin
    is, is_state = IS(storage_size, hm_vectors, vector_size, is_layers, storage_size)
    gs, gs_state = GS(storage_size, storage_size, gs_layers, storage_size)
    go, go_state = GO(storage_size, vector_size, go_layers, vector_size)
    decoder = DECODER(is, gs, go)
    state   = [is_state,gs_state,[go_state for i in 1:hm_vectors]]
[decoder, state]
end

(decoder::DECODER)(state, storage, vectors, enc_storages, enc_vectors) =
begin
    attended = attend(storage, enc_storages, enc_vectors)
    is_out, is_state = decoder.is(state[1], storage, hcat(attended...))
    gs_out, gs_state = decoder.gs(state[2], storage, is_out)
    go_out, go_state = [], []
    for (vector, state_v) in zip(vectors, state[3])
        go_out_v, go_state_v = decoder.go(state_v, is_out, vector)
        push!(go_out, go_out_v)
        push!(go_state, go_state_v)
    end
    output = [gs_out, go_out]
    state  = [is_state, gs_state, go_state]
[output, state]
end

attend(dec_storage, enc_storages, enc_vectors) =
begin
    similarities = soft([sum(dec_storage .* enc_storage) for enc_storage in enc_storages])
    # similarities = [sum(dec_storage .* enc_storage) for enc_storage in enc_storages]
# enc_vectors[argmax(reshape(similarities, length(similarities)))]
sum([vec .* sim for (vec, sim) in zip(enc_vectors,similarities)])
end

transfer_states!(enc_state, dec_state) =
begin
    for (dec_module, enc_module) in zip(dec_state, enc_state)
        for (dec_layer, enc_layer) in zip(dec_module, enc_module)
            i = 0
            for enc_neuron in enc_layer
                i +=1
                dec_layer[i] = enc_neuron
            end
        end
    end
end

propogate(enc, dec, enc_state, x, length_y; dec_state=nothing) =
begin
    storage_size = length(enc_state[2][end])
    enc_outs = [zeros(1, storage_size), nothing]
    enc_gs_time, enc_go_time = [], []
    for timestep in x
        enc_outs, enc_state = enc(enc_state, enc_outs[1], timestep)
        push!(enc_gs_time, enc_outs[1])
        push!(enc_go_time, enc_outs[2])
    end

    # transfer_states!(enc_state, dec_state)

    dec_outs = enc_outs
    # dec_outs = [zeros(1, storage_size), [zeros(1, vector_size) for i in 1:hm_vectors]]
    if dec_state == nothing dec_state = enc_state end # [zeros(1, storage_size), [zeros(1, vector_size) for i in 1:hm_vectors]]
    dec_go_time = []
    for timestep in 1:length_y
        dec_outs, dec_state = dec(dec_state, dec_outs[1], dec_outs[2], enc_gs_time, enc_go_time)
        push!(dec_go_time, dec_outs[2])
    end
dec_go_time
end

# using Statistics: norm

grads(result, encoder, decoder) =
begin
    grads = []
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    gradient = grad(result, getfield(layer, lfield))
                    # if gradient > 0
                    #     gradient = 1.0
                    # end
                    push!(grads, gradient)
                    # @show norm(gradient)
                end
            end
        end
    end
grads
end

upd!(encoder, decoder, grads, lr) =
begin
    i = 0
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    i +=1
                    setfield!(layer, lfield, Param(getfield(layer, lfield) - grads[i].*lr))
                end
            end
        end
    end
end

upd_rms!(encoder, decoder, grads, lr, accu_grads; alpha=0.9) =
begin
    i = 0
    for model in [encoder, decoder]
        for mfield in fieldnames(typeof(model))
            net = getfield(model, mfield)
            for nfield in fieldnames(typeof(net))
                layer = getfield(net, nfield)
                for lfield in fieldnames(typeof(layer))
                    i +=1
                    accu_grads[i] = alpha .* accu_grads[i] .+ (1 - alpha) .* grads[i].^2
                    setfield!(layer, lfield, Param(getfield(layer, lfield) - grads[i].*lr./sqrt.(accu_grads[i].+1e-8)))
                end
            end
        end
    end
end

soft(arr) = ((arr) -> arr/sum(arr))(exp.(arr))

sqe(out, y) = sum(out .- y) ^2

sequence_loss(out, y) =
begin
    sum([sum(sum([sqe(out_e, y_e) for (out_e,y_e) in zip(out_t,y_t)])) for (out_t,y_t) in zip(out, y)])
end

drop(x; rate=0.1) = rand() < rate ? x : 0.0
begin
    # for i in 1:length(array):
    #     if rand() < rate
    #         array[i]
#     new_arr = [rand() < rate ? 0.0 : e for e in array]
# reshape(new_arr, size(array))
# array
end

max_stars_repo_name: developersfox/GSTM

max_issues_repo_name: developersfox/GSTM

max_forks_repo_name: developersfox/GSTM


text:
module SingleFieldFEFunctionsTests

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -48.987102771054815 <= q <= 45.974621434893855)

                     add_NL_constraint(m, :(gelu(-0.368688824745 + 0.553104331767496*gelu(-0.34679271547764046 + -0.03387429483407178*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.17925952214788987*gelu(0.6682158984675199 + 0.6233864552377666*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4698400936211189*gelu(0.8534714115378583 + 0.6703763846356168*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.7276938776175426*gelu(-0.09327714623589056 + -0.24359337663872616*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4202042464647322*gelu(-0.09609361556421225 + -0.8712983668552394*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + gelu(-0.8761821794177092 + 0.4562914248397121*gelu(-0.34679271547764046 + -0.03387429483407178*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3827259206071312*gelu(0.6682158984675199 + 0.6233864552377666*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.20062384412841983*gelu(0.8534714115378583 + 0.6703763846356168*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.28064116456943733*gelu(-0.09327714623589056 + -0.24359337663872616*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4949171063694977*gelu(-0.09609361556421225 + -0.8712983668552394*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + gelu(0.7716665307334409 + 0.20199447224764766*gelu(-0.34679271547764046 + -0.03387429483407178*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3172295658101083*gelu(0.6682158984675199 + 0.6233864552377666*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.10629001518557368*gelu(0.8534714115378583 + 0.6703763846356168*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4930771030063288*gelu(-0.09327714623589056 + -0.24359337663872616*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.5476037221052898*gelu(-0.09609361556421225 + -0.8712983668552394*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + gelu(0.6820199557891846 + 0.4693680154187745*gelu(-0.34679271547764046 + -0.03387429483407178*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.6552396777765495*gelu(0.6682158984675199 + 0.6233864552377666*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.535415419107649*gelu(0.8534714115378583 + 0.6703763846356168*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.022291893207070146*gelu(-0.09327714623589056 + -0.24359337663872616*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.04141935716867051*gelu(-0.09609361556421225 + -0.8712983668552394*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + gelu(-0.779255880286021 + 0.8816016628889711*gelu(-0.34679271547764046 + -0.03387429483407178*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3409151924818685*gelu(0.6682158984675199 + 0.6233864552377666*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.5165976261581497*gelu(0.8534714115378583 + 0.6703763846356168*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.06878582640645536*gelu(-0.09327714623589056 + -0.24359337663872616*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.8469129987678201*gelu(-0.09609361556421225 + -0.8712983668552394*gelu(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*gelu(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*gelu(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*gelu(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*gelu(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# "Derivative of Hankel function of zeroth kind"
# function diffhankelh1(n, z)
#     if n != 0
#         (hankelh1(-1 + n, z) - hankelh1(1 + n, z)) / 2
#     else
#         -hankelh1(one(n), z)
#     end
# end
#
# "Derivative of Bessel function of first kind"
# function diffbesselj(n, z)
#     if n != 0
#         (besselj(-1 + n, z) - besselj(1 + n, z)) / 2
#     else
#         -besselj(one(n), z)
#     end
# end


text:
# Chapter 8. Strings
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap08
module ch8
'x'
🍌 = 5 # \:banana:
🐷 = "Yi-Kun Yang" # \:pig:
@show(🐷)

# A String is a sequence
fruit = "banana"
fruit[1]
fruit[end]

# nextind, firstind, sizeof
fruits = "🍌🍎🍐🍑"
sizeof(fruits)
nextind(fruits, 1)
fruits[5]
nextind(fruits, 5)
fruits[9]

i = 13

z = 1
i = 1
while i < 3
    global i
    i += 1
end

#@show(z)

# Traversal
# Use the one below which is way easier
index = firstindex(fruits)
while index <= sizeof(fruits)
    letter = fruits[index]
    println(letter)
    global index = nextind(fruits, index)
end

# Way easier
for letter in fruits
    println(letter)
end

# Exercise 8-1
function printbackward(word)
    wordLength = length(word)
    for i in wordLength:-1:1
        println(word[i])
    end
end

printbackward("ABC")

# Exercise 8-2

# String Slices
str = "Julius Caesar"

str[1:6]
str[:]
#fruits[1:5]

# Strings are immutable

# String Interpolation
greet = "Hello"
whom = "World"
"$greet, $(whom)"
"1 + 2 = $(1 +2)"

# Searching
# Exercise 8-4
function find(word, letter, index = 1)::Int64
    if !(letter isa Char)
        println("letter must be a Char")
        return -1
    end

    #index = firstindex(word)

    while index <= sizeof(word)
        if word[index] == letter
            return index
        end
        index = nextind(word, index)
    end
    -1
end

find("abc", 'a', 2)

function count(word, char)::Int64
    if !(char isa Char)
        println("char must be a Char")
        return -1
    end
    count = 0
    for letter in word
        if letter == char
            count += 1
        end
    end

    return count
end
count("banana", 'a')

# String Library
uppercase("Hello, World!")
ind = findfirst("a", "banana")
ind.start
ind.stop
ind = findfirst("na", "banana")
ind.start
ind.stop

# The ∈ Operator (\in)
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_the_operator
# The IN Operator
@show(∈)
#@show(ϵ)

'a' ∈ "banana"

function inboth(word1, word2)
    for char in word1
        if char ∈ word2
            print(char, " ")
        end
    end
end

inboth("apple", "orange")

# String Comparison
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_string_comparison

# ==

# Debugging
# Exercise 8-6

end # module ch8


text:
module TestBioSequences

using Test
using Random
using LinearAlgebra: normalize
import BioSymbols
using BioSequences
using IntervalTrees: IntervalValue
using StatsBase
using YAML
using BioCore

const codons = [
    "AAA", "AAC", "AAG", "AAU",
    "ACA", "ACC", "ACG", "ACU",
    "AGA", "AGC", "AGG", "AGU",
    "AUA", "AUC", "AUG", "AUU",
    "CAA", "CAC", "CAG", "CAU",
    "CCA", "CCC", "CCG", "CCU",
    "CGA", "CGC", "CGG", "CGU",
    "CUA", "CUC", "CUG", "CUU",
    "GAA", "GAC", "GAG", "GAU",
    "GCA", "GCC", "GCG", "GCU",
    "GGA", "GGC", "GGG", "GGU",
    "GUA", "GUC", "GUG", "GUU",
    "UAA", "UAC", "UAG", "UAU",
    "UCA", "UCC", "UCG", "UCU",
    "UGA", "UGC", "UGG", "UGU",
    "UUA", "UUC", "UUG", "UUU",
    # translatable ambiguities in the standard code
    "CUN", "CCN", "CGN", "ACN",
    "GUN", "GCN", "GGN", "UCN"
]

function random_translatable_rna(n)
    probs = fill(1.0 / length(codons), length(codons))
    cumprobs = cumsum(probs)
    r = rand()
    x = Vector{String}(undef, n)
    for i in 1:n
        x[i] = codons[searchsorted(cumprobs, rand()).start]
    end

    return string(x...)
end

function get_bio_fmt_specimens()
    path = joinpath(dirname(@__FILE__), "BioFmtSpecimens")
    if !isdir(path)
        run(`git clone --depth 1 https://github.com/BioJulia/BioFmtSpecimens.git $(path)`)
    end
end

# The generation of random test cases...

function random_array(n::Integer, elements, probs)
    cumprobs = cumsum(probs)
    x = Vector{eltype(elements)}(undef, n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end

# Return a random DNA/RNA sequence of the given length.
function random_seq(n::Integer, nts, probs)
    cumprobs = cumsum(probs)
    x = Vector{Char}(undef, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return String(x)
end

function random_seq(::Type{A}, n::Integer) where A<:Alphabet
    nts = BioSymbols.alphabet(A)
    probs = Vector{Float64}(undef, length(nts))
    fill!(probs, 1 / length(nts))
    return BioSequence{A}(random_seq(n, nts, probs))
end

function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

function intempdir(fn::Function, parent=tempdir())
    dirname = mktempdir(parent)
    try
        cd(fn, dirname)
    finally
        rm(dirname, recursive=true)
    end
end

function random_dna_kmer(len)
    return random_dna(len, [0.25, 0.25, 0.25, 0.25])
end

function random_rna_kmer(len)
    return random_rna(len, [0.25, 0.25, 0.25, 0.25])
end

function random_dna_kmer_nucleotides(len)
    return random_array(len, [DNA_A, DNA_C, DNA_G, DNA_T],
                        [0.25, 0.25, 0.25, 0.25])
end

function random_rna_kmer_nucleotides(len)
    return random_array(len, [RNA_A, RNA_C, RNA_G, RNA_U],
                        [0.25, 0.25, 0.25, 0.25])
end

function dna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    for (i, c) in enumerate(seq)
        if c     ==   'A'
            seqc[i] = 'T'
        elseif c ==   'C'
            seqc[i] = 'G'
        elseif c ==   'G'
            seqc[i] = 'C'
        elseif c ==   'T'
            seqc[i] = 'A'
        else
            seqc[i] = seq[i]
        end
    end
    return String(seqc)
end

function rna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    for (i, c) in enumerate(seq)
        if c == 'A'
            seqc[i] = 'U'
        elseif c == 'C'
            seqc[i] = 'G'
        elseif c == 'G'
            seqc[i] = 'C'
        elseif c == 'U'
            seqc[i] = 'A'
        else
            seqc[i] = seq[i]
        end
    end
    return String(seqc)
end

function random_interval(minstart, maxstop)
    start = rand(minstart:maxstop)
    return start:rand(start:maxstop)
end

include("symbols.jl")

@testset "Sequences" begin
    a = dna"A-CG-G"; b = rna"A-CG-G"; c = aa"AK-MV-";
    @test ungap(a) == dna"ACGG"
    @test ungap(b) == rna"ACGG"
    @test ungap(c) == aa"AKMV"
    @test ungap!(a) === a && a == dna"ACGG"
    @test ungap!(b) === b && b == rna"ACGG"
    @test ungap!(c) === c && c == aa"AKMV"
end

@testset "BioSequences" begin
    include("bioseq/conversion.jl")
    include("bioseq/basics.jl")
    include("bioseq/hashing.jl")
    include("bioseq/iteration.jl")
    include("bioseq/subseq.jl")
    include("bioseq/mutability.jl")
    include("bioseq/print.jl")
    include("bioseq/transformations.jl")
    include("bioseq/mutability.jl")
    include("bioseq/predicates.jl")
    include("bioseq/find.jl")
    include("bioseq/counting.jl")
    include("bioseq/gc_content.jl")
    include("bioseq/ambiguous.jl")
    include("bioseq/shuffle.jl")
end

@testset "ReferenceSequences" begin
    include("refseq/conversion.jl")
    include("refseq/basics.jl")
    include("refseq/long.jl")
    include("refseq/print.jl")
    include("refseq/random.jl")
end

include("composition.jl")

@testset "Kmers" begin
    include("kmers/conversion.jl")
    include("kmers/comparisons.jl")
    include("kmers/length.jl")
    include("kmers/arithmetic.jl")
    include("kmers/access.jl")
    include("kmers/random.jl")
    include("kmers/find.jl")
    include("kmers/print.jl")
    include("kmers/transformations.jl")
    include("kmers/mismatches.jl")
    include("kmers/eachkmer.jl")
    include("kmers/debruijn_neighbors.jl")
    include("kmers/shuffle.jl")
end

@testset "Search" begin
    include("search/exact.jl")
    include("search/approximate.jl")
    include("search/regex.jl")
    include("search/pwm.jl")
end

include("translation.jl")

include("demultiplexer.jl")

@testset "Reading and Writing" begin
    include("io/FASTA.jl")
    include("io/FASTQ.jl")
    include("io/twobit.jl")
    include("io/abif.jl")
end

end


text:
@testset "event" begin
    seed!(42)
    y = rand(1:10, 10)
    l = PoissonLikelihood(3.0)
    @test AGP.treat_labels!(y, l) == y
    @test_throws ErrorException AGP.treat_labels!(rand(10), l)
end


text:
module PotentialCalculation

using Reexport

include("SubModules/identical.jl")
@reexport using .IdenticalTools

include("SubModules/atoms.jl")
@reexport using .Atoms

include("SubModules/molecules.jl")
@reexport using .Molecules

include("SubModules/clusters.jl")
@reexport using .Clusters

include("SubModules/unitconversions.jl")
@reexport using .UnitConversions

include("SubModules/fileaccess.jl")
@reexport using .Fileaccess

include("SubModules/calculators.jl")
@reexport using .Calculators

include("SubModules/sample.jl")
@reexport using .Sample

include("SubModules/restarttools.jl")
@reexport using .Restarttools

include("SubModules/psi4.jl")

end # module


text:
export Distributor
export createFromSends, createFromRecvs
export resolve, resolveReverse, resolveWaits
export resolvePosts, resolveReversePosts, resolveReverseWaits

# methods (and docs) are currently based straight off Epetra_Distributor to match Comm

"""
The base type for gather/scatter setup.
All subtypes must have the following methods, with DistributorImpl standing
in for the subtype:

    createFromSends(dist::DistributorImpl,exportPIDs::AbstractArray{PID})::Integer where PID <:Integer
Sets up the Distributor object using a list of process IDs to which we
export.  Returns the number of IDs this processor will be receiving.

    createFromRecvs(dist::DistributorImpl, remoteGIDs::AbstractArray{GID}, remotePIDs::AbstractArray{PID})::Tuple{AbstractArray{GID}, AbstractArray{PID}} where GID <: Integer where PID <: Integer
Sets up the Distributor object using a list of remote global IDs and corresponding
PIDs.  Returns a tuple with the global IDs and their respective processor IDs
being sent by this processor.

    resolvePosts(dist::DistributorImpl, exportObjs::AbstractArray)
Post buffer of export objects.  Other, local work can be done before resolving
the waits.  Otherwise, as [`resolve`](@ref).

    resolveWaits(dist::DistributorImpl)::AbstractArray
Wait on a set of posts

    resolveReversePosts(dist::DistributorImpl, exportObjs::AbstractArray)
Do reverse post of buffer of export objects Other, local work can be done before
resolving the waits.  Otherwise, as
[`resolveReverse`](@ref).

    resolveReverseWaits(dist::DistributorImpl)::AbstractArray
Wait on a set of reverse posts.

"""
abstract type Distributor{GID <: Integer, PID <: Integer, LID <: Integer}
end

"""
    createFromSends(dist::Distributor, exportPIDs::AbstractArray{<:Integer})::Integer

Sets up the Distributor object using a list of process IDs to which we
export and the number of IDs being exported.  Returns the number of
IDs this processor will be receiving
"""
function createFromSends(dist::Distributor{GID, PID, LID}, exportPIDs::AbstractArray{<:Integer}) where GID <: Integer where PID <: Integer where LID <: Integer
    createFromSends(dist, Array{PID}(exportPIDs))
end

"""
    createFromRecvs(dist::Distributor{GID, PID, LID}, remoteGIDs::AbstractArray{<:Integer}, remotePIDs::AbstractArray{<:Integer})::Tuple{AbstractArray{GID}, AbstractArray{PID}}

Sets up the Distributor object using a list of remote global IDs and
corresponding PIDs.  Returns a tuple with the global IDs and their
respective processor IDs being sent by me.
"""
function createFromRecvs(dist::Distributor{GID, PID, LID}, remoteGIDs::AbstractArray{<:Integer},
        remotePIDs::AbstractArray{<:Integer}) where GID <: Integer where PID <: Integer where LID <: Integer
    createFromRecvs(dist, Array{GID}(remoteGIDs), Array{PID}(remotePIDs))
end

"""
    resolve(dist::Distributor, exportObjs::AbstractArray{T})::AbstractArray{T}

Execute the current plan on buffer of export objects and return the
objects set to this processor
"""
function resolve(dist::Distributor, exportObjs::AbstractArray{T})::AbstractArray{T} where T
    resolvePosts(dist, exportObjs)
    resolveWaits(dist)::AbstractArray{T}
end

"""
    resolveReverse(dist::Distributor, exportObjs::AbstractArray{T})::AbstractArray{T}

Execute the reverse of the current plan on buffer of export objects and
return the objects set to this processor
"""
function resolveReverse(dist::Distributor, exportObjs::AbstractArray{T})::AbstractArray{T} where T
    resolveReversePosts(dist, exportObjs)
    resolveReverseWaits(dist)
end

#### required method documentation stubs ####

"""
    resolvePosts(dist::Distributor, exportObjs::AbstractArray)

Post buffer of export objects.  Other, local work can be done before resolving
the waits.  Otherwise, as [`resolve`](@ref).
"""
function resolvePosts end

"""
    resolveWaits(dist::Distributor)::AbstractArray

Wait on a set of posts.
"""
function resolveWaits end

"""
    resolveReversePosts(dist::Distributor, exportObjs::AbstractArray)
Do reverse post of buffer of export objects Other, local work can be done before
resolving the waits.  Otherwise, as
[`resolveReverse`](@ref).
"""
function resolveReversePosts end

"""
    resolveReverseWaits(dist::Distributor)::AbstractArray

Wait on a set of reverse posts.
"""
function resolveReverseWaits end


text:
module Numerics

  using Reexport

  @reexport using DataFrames
  @reexport using DataStructures
  @reexport using LinearAlgebra
  @reexport using SparseArrays
  @reexport using Statistics
  @reexport using StatsBase
  @reexport using Random

end


text:
# Julia wrapper for header: /usr/local/include/libavutil/tree.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_tree_node_alloc,
    av_tree_find,
    av_tree_insert,
    av_tree_destroy,
    av_tree_enumerate

function av_tree_node_alloc()
    ccall((:av_tree_node_alloc,libavutil),Ptr{AVTreeNode},())
end

function av_tree_find(root,key,cmp,next::NTuple{2,Ptr{Void}})
    ccall((:av_tree_find,libavutil),Ptr{Void},(Ptr{AVTreeNode},Ptr{Void},Ptr{Void},NTuple{2,Ptr{Void}}),root,key,cmp,next)
end

function av_tree_insert(rootp,key,cmp,next)
    ccall((:av_tree_insert,libavutil),Ptr{Void},(Ptr{Ptr{AVTreeNode}},Ptr{Void},Ptr{Void},Ptr{Ptr{AVTreeNode}}),rootp,key,cmp,next)
end

function av_tree_destroy(t)
    ccall((:av_tree_destroy,libavutil),Void,(Ptr{AVTreeNode},),t)
end

function av_tree_enumerate(t,opaque,cmp,enu)
    ccall((:av_tree_enumerate,libavutil),Void,(Ptr{AVTreeNode},Ptr{Void},Ptr{Void},Ptr{Void}),t,opaque,cmp,enu)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


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

using JuliaFEM, Test

# Test that if 3d continuum problem have some elements we don't know how to
# deal with, raise error with clear message

elements = [Element(Seg3, (1, 2, 3))]
problem = Problem(Elasticity, "test seg3", 3)
add_elements!(problem, elements)
@test_throws ErrorException assemble!(problem, 0.0)


text:
using Test
using CuboidGrids
Lar = CuboidGrids

@testset "LarGrid-0-1 Tests" begin
   @testset "grid_0" begin
      @test Lar.grid_0(1) == [0  1]
      @test Lar.grid_0(2) == [0  1  2]
      @test Lar.grid_0(3) == [0  1  2  3]
   end

   @testset "Grid_1 Tests" begin
      @test repr(Lar.grid_1(1)) == "[0; 1]"
      @test Lar.grid_1(2) == [0 1; 1 2]
      @test Lar.grid_1(3) == [0 1 2; 1 2 3]
   end

   @testset "Grid Tests" begin
      @test repr(Lar.larGrid(1)(0)) == "[0 1]"
      @test repr(Lar.larGrid(1)(1)) == "[0; 1]"
      @test Lar.larGrid(2)(0) == [0  1  2]
      @test Lar.larGrid(2)(1) == [0 1; 1 2]
      @test Lar.larGrid(3)(0) == [0  1  2  3]
      @test Lar.larGrid(3)(1) == [0  1  2; 1  2  3]
   end
end

@testset "LarVertProd" begin
      vertexDomain(n) = hcat([k for k in 0:n-1]...)

      @testset "LarVertProd 1D" begin
            shape = [3]
	      vertLists = [vertexDomain(k+1) for k in shape]
            @test typeof(Lar.larVertProd(vertLists))==Array{Int64,2}
            @test size(Lar.larVertProd(vertLists))==(1, 4)
            @test Lar.larVertProd(vertLists)[:,1]==[0]
            @test Lar.larVertProd(vertLists)[:,4]==[3]
      end

      @testset "LarVertProd 2D" begin
            shape = [3,2]
	      vertLists = [vertexDomain(k+1) for k in shape]
            @test typeof(Lar.larVertProd(vertLists))==Array{Int64,2}
            @test size(Lar.larVertProd(vertLists))==(2, 12)
            @test Lar.larVertProd(vertLists)[:,1]==[0;0]
            @test Lar.larVertProd(vertLists)[:,12]==[3;2]
      end

      @testset "LarVertProd 3D" begin
            shape = [3,2,1]
	      vertLists = [vertexDomain(k+1) for k in shape]
            @test typeof(Lar.larVertProd(vertLists))==Array{Int64,2}
            @test size(Lar.larVertProd(vertLists))==(3, 24)
            @test Lar.larVertProd(vertLists)[:,1]==[0;0;0]
            @test Lar.larVertProd(vertLists)[:,24]==[3;2;1]
      end
end

@testset "Index2addr Tests" begin
   @testset "Shape 1D Tests" begin
      @test Lar.index2addr([10])([0])==1
      @test Lar.index2addr([10])([9])==10
      @test [Lar.index2addr([10])([index]) for index in collect(0:9)]==collect(1:10)
   end

   @testset "Shape 2d Tests" begin
      aa = Lar.cart([[0;1;2],[0;1]])
      bb = Lar.cart([collect(0:9),collect(0:1)])
      cc = Lar.cart([collect(0:2),collect(0:2)])
      dd = "Tuple{Int64,Int64}[(0,0),(0,1),(1,0),(1,1),(2,0),(2,1)]"
      ee = "Tuple{Int64,Int64}[(0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1)]"
      ff = repr( Lar.cart([collect(0:2),collect(0:1)]) )
      @test [ Lar.index2addr([3,2])(collect(index)) for index in aa ]==collect(1:6)
      @test [ Lar.index2addr([10,2])(collect(index)) for index in bb ] == collect(1:20)
      @test [ Lar.index2addr([3,3])(collect(index)) for index in cc] == collect(1:9)
   end

   @testset "Shape 3d Tests" begin
      aaa = Lar.cart([collect(0:3),collect(0:2),collect(0:1)])
      bbb = Lar.cart([[0;1;2],[0;1],[0;1]])
      ccc = Lar.cart([[0;1],[0;1],[0;1;2]])
      @test Lar.index2addr([3,2,1])([0,0,0]) == 1
      @test [Lar.index2addr([4,3,2])(collect(index))
      	for index in aaa ] == collect(1:24)
      @test [Lar.index2addr([3,2,2])(collect(index))
      	for index in bbb ] == collect(1:12)
      @test [Lar.index2addr([2,2,3])(collect(index))
      	for index in ccc ] == collect(1:12)
   end
end

@testset "LarCellProd Tests" begin
	 @testset "3D shape" begin
		@test typeof(Lar.larGridSkeleton([3,2,1])(3)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3,2,1])(2)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3,2,1])(1)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3,2,1])(0)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([1,1,1])(3)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([1,1,1])(2)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([1,1,1])(1)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([1,1,1])(0)) ==
		   		Array{Array{Int64,1},1}
	 end
	 @testset "2D shape" begin
		@test typeof(Lar.larGridSkeleton([3,2])(2)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3,2])(1)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3,2])(0)) ==
		   		Array{Array{Int64,1},1}
	 end
	 @testset "1D shape" begin
		@test typeof(Lar.larGridSkeleton([3])(1)) ==
		   		Array{Array{Int64,1},1}
		@test typeof(Lar.larGridSkeleton([3])(0)) ==
		   		Array{Array{Int64,1},1}
	 end
end

@testset "FilterByOrder Tests" begin
      term = "000"
      bit = '0'
      theTerm = Vector{Char}(term)
      @test typeof(theTerm) == Array{Char,1}
      @test parse(Int8,bit) == 0
      @test [parse(Int8,bit) for bit in theTerm] == zeros(3)
      out = hcat([[parse(Int8,bit) for bit in term] for term in
      	Lar.binaryRange(3)]...)
      @test typeof(out) == Array{Int8,2}
      @test size(out) == (3,8)
      @test repr(out) == "Int8[0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1]"

      @testset "$n" for n in 1:4
            data = [Lar.filterByOrder(n)[k] for (k,el) in enumerate(Lar.filterByOrder(n))]
            @test sum(map(length,data)) == 2^n
      end

      @testset "$n,$k" for n in 1:4, k in 0:n
            @test length(Lar.filterByOrder(n)[k+1]) == binomial(n,k)
      end
end

@testset "LarImageVerts Tests" "$shape" for shape in [[3,2,1],[3,2],[10,10,10]]
          @test size(Lar.larImageVerts(shape)) == (length(shape),prod(shape .+ 1))
       end

@testset "LarCuboids Tests" begin
      @testset "Simple shape" begin
            shape = [3,2,1]
            cubes = Lar.cuboidGrid(shape,true)
            verts, cells = cubes
            @test typeof(verts) == Array{Float64,2}
            VV,EV,FV,CV = cells
            @testset "$basis" for basis in [VV,EV,FV,CV]
              @test typeof(basis) == Array{Array{Int64,1},1}
            end
      end

      @testset "$shape" for shape in [[3,2,1],[1,1,1],[3,3,10]]
            endShape = shape .+ 1
            cubes = Lar.cuboidGrid(shape,true);
            verts, cells = cubes;
            VV,EV,FV,CV = cells;
            @test size(Lar.cuboidGrid(shape)[1],2) == prod(collect(endShape))
            @test length(Lar.cuboidGrid(shape)[2]) == prod(shape)
            @test typeof(verts) == Array{Float64,2}
	      @test typeof(VV) == Array{Array{Int64,1},1}
	      @test typeof(EV) == Array{Array{Int64,1},1}
	      @test typeof(FV) == Array{Array{Int64,1},1}
	      @test typeof(FV) == Array{Array{Int64,1},1}
      end

      @testset "$shape" for shape in [[3,2],[1,1],[3,10]]
            endShape = shape .+ 1
            @test size(Lar.cuboidGrid(shape)[1],2) == prod(collect(endShape))
            @test length(Lar.cuboidGrid(shape)[2]) == prod(shape)
            cubes = Lar.cuboidGrid(shape,true)
            verts, cells = cubes
            @test typeof(verts) == Array{Float64,2}
            VV,EV,FV = cells
	      @test typeof(VV) == Array{Array{Int64,1},1}
	      @test typeof(EV) == Array{Array{Int64,1},1}
	      @test typeof(FV) == Array{Array{Int64,1},1}
      end

      @testset "$shape" for shape in [[3],[1],[10]]
            endShape = shape .+ 1
            @test size(Lar.cuboidGrid(shape)[1],2) == prod(collect(endShape))
            @test length(Lar.cuboidGrid(shape)[2]) == prod(shape)
            cubes = Lar.cuboidGrid(shape,true)
            verts, cells = cubes
            @test typeof(verts) == Array{Float64,2}
            VV,EV = cells
	      @test typeof(VV) == Array{Array{Int64,1},1}
	      @test typeof(EV) == Array{Array{Int64,1},1}
      end
end

@testset "LarModelProduct Tests" begin
      geom_0 = hcat([[x] for x=0.:10]...)
      topol_0 = [[i,i+1] for i=1:9]
      geom_1 = hcat([[0.],[1.],[2.]]...)
      topol_1 = [[1,2],[2,3]]
      model_0 = (geom_0, topol_0)
      model_1 = (geom_1, topol_1)
      model_2 = Lar.larModelProduct(model_0, model_1)
      model_3 = Lar.larModelProduct(model_2, model_1)

      @testset "Result size" begin
            @test size(model_2[1],2) == size(model_0[1],2) * size(model_1[1],2)
            @test length(model_2[2]) == length(model_0[2]) * length(model_1[2])
            @test size(model_3[1],2) == size(model_2[1],2) * size(model_1[1],2)
            @test length(model_3[2]) == length(model_2[2]) * length(model_1[2])
      end

      @testset "Model 0 type" begin
            m0_V,m0_EV = model_0
            @test typeof(m0_V) == Array{Float64,2}
            @test typeof(m0_EV) == Array{Array{Int64,1},1}
      end

      @testset "Model 1 type" begin
            m1_V,m1_EV = model_1
            @test typeof(m1_V) == Array{Float64,2}
            @test typeof(m1_EV) == Array{Array{Int64,1},1}
      end

      @testset "Model 2 type" begin
            m2_V,m2_EV = model_2
            @test typeof(m2_V) == Array{Float64,2}
            @test typeof(m2_EV) == Array{Array{Int64,1},1}
      end

      @testset "Model 3 type" begin
            m3_V,m3_EV = model_3
            @test typeof(m3_V) == Array{Float64,2}
            @test typeof(m3_EV) == Array{Array{Int64,1},1}
      end

      @testset "Simple product" begin
            modelOne = ([0 0 1 1; 0 1 0 1], Array{Int64,1}[[1, 2, 3, 4]])
            modelTwo = ([0 1], Array{Int64,1}[[1, 2]])
            @test Lar.larModelProduct(modelOne, modelTwo) == Lar.larModelProduct(modelTwo, modelOne)
            @test Lar.larModelProduct(modelTwo, modelTwo)[1] == [0  0  1  1; 0  1  0  1]
            @test Lar.larModelProduct(modelOne,modelOne)[2][1] == Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
      end
end


text:
abstract type AbstractUpdateSelector end
abstract type AbstractUpdate end
abstract type WholespaceUpdate <: AbstractUpdate end
abstract type ManifoldUpdate <: AbstractUpdate end
abstract type AbstractUpdateState end

@with_kw struct PartlySmoothOptimizer{UpSel, WhUp, ManUp} <: Optimizer
    update_selector::UpSel = AlternatingUpdateSelector()
    wholespace_update::WhUp = WholespaceProximalGradient()
    manifold_update::ManUp = ManifoldIdentity()
end

mutable struct PartlySmoothOptimizerState{Tambiant, Tmanpoint, Tmanvec} <: OptimizerState
    it::Int64
    x::Point{Tambiant, Tmanpoint}
    x_old::Point{Tambiant, Tmanpoint}
    M::Manifold
    M_old::Manifold
    f_x::Float64
    g_x::Float64
    ∇f_x::Tambiant
    temppoint_amb::Tambiant
    tempvec_man::Tmanvec
    selected_update
    previous_update
    update_to_updatestate::Dict
    ncalls_f::Int64
    ncalls_g::Int64
    ncalls_∇f::Int64
    ncalls_proxg::Int64
    ncalls_gradₘF::Int64
    ncalls_HessₘF::Int64
    ncalls_retr::Int64
    niter_CG::Int64
    niter_manls::Int64
end

function PartlySmoothOptimizerState(
    o::PartlySmoothOptimizer,
    x_amb::T,
    g::R,
    M::Manifold
) where {T,R}
    x_man = project(M, x_amb)
    return PartlySmoothOptimizerState(
        -1,
        Point(x_amb, x_man, M, ambiant_repr),
        Point(x_amb, x_man, M, ambiant_repr),
        M,
        M,
        -Inf,
        -Inf,
        zero(x_amb),
        zero(x_amb),
        zero_tangent_vector(M, x_man),
        nothing,
        nothing,
        Dict(),
        0, 0, 0, 0, 0, 0, 0, 0, 0
    )
end

Base.summary(o::PartlySmoothOptimizer) = string(
    "PSOpt - ", summary(o.update_selector), " - ",
    summary(o.wholespace_update), " - ", summary(o.manifold_update)
)

function print_header(pso::PartlySmoothOptimizer)
    println("---------------------")
    println("--- Partly Smooth Optimizer")
    println(" - update selection    ", summary(pso.update_selector))
    println(" - wholespace update   ", summary(pso.wholespace_update))
    println(" - manifold update     ", summary(pso.manifold_update))
end

function update_fg∇f!(state::PartlySmoothOptimizerState, pb)
    x = get_repr(state.x)

    state.f_x = f(pb, x); state.ncalls_f += 1
    state.g_x = g(pb, x); state.ncalls_g += 1
    ∇f!(pb, state.∇f_x, x); state.ncalls_∇f += 1
    state.it += 1
    return
end

function initial_state(o::PartlySmoothOptimizer, x, reg; manifold=wholespace_manifold(reg, x))
    initstate = PartlySmoothOptimizerState(o, x, reg, manifold)
    initstate.update_to_updatestate[o.wholespace_update] = initial_state(o.wholespace_update, x, reg)
    initstate.update_to_updatestate[o.manifold_update] = initial_state(o.manifold_update, x, reg)
    return initstate
end

function convert_point_repr!(x::Point{Tamb, Tman}, ::WholespaceUpdate) where {Tamb, Tman}
    if x.repr == manifold_repr
        x.amb_repr = embed(x.M, x.man_repr)
        x.repr = ambiant_repr
    end
    return
end

function convert_point_repr!(x::Point{Tamb, Tman}, ::ManifoldUpdate) where {Tamb, Tman}
    if x.repr == ambiant_repr
        x.man_repr = project(x.M, x.amb_repr)
        x.repr = manifold_repr
    end
    return
end

function update_iterate!(ostate::PartlySmoothOptimizerState, pb, optimizer::PartlySmoothOptimizer)
    ostate.M = ostate.M_old

    @timeit_debug "updateselect" select_update!(optimizer.update_selector, ostate, optimizer, pb)

    # A composite update expects points in ambiant space, while a manifold one in manifold representation.
    @timeit_debug "convpointrepr" convert_point_repr!(ostate, optimizer, pb)

    @timeit_debug "updateit" update_iterate!(ostate, pb, ostate.selected_update)

    return
end

function build_optimstate(::PartlySmoothOptimizer, state, iteration, time, normstep, minsubgradient_tan, minsubgradient_norm, optimstate_extensions)
    return OptimizationState(
        it = iteration,
        time = time,
        f_x = state.f_x,
        g_x = state.g_x,
        norm_step = normstep,
        norm_minsubgradient_tangent = minsubgradient_tan,
        norm_minsubgradient_normal = minsubgradient_norm,
        ncalls_f = state.ncalls_f,
        ncalls_g = state.ncalls_g,
        ncalls_∇f = state.ncalls_∇f,
        ncalls_proxg = state.ncalls_proxg,
        ncalls_gradₘF = state.ncalls_gradₘF,
        ncalls_HessₘF = state.ncalls_HessₘF,
        ncalls_retr = state.ncalls_retr,
        niter_CG = state.niter_CG,
        niter_manls = state.niter_manls,
        additionalinfo = (;
            zip(
                [osextension.key for osextension in optimstate_extensions],
                [
                    deepcopy(osextension.getvalue(state))
                    for osextension in optimstate_extensions
                ],
            )...
        ),
    )
end

#
### Printing and logging
#
str_updatelog(::AbstractUpdate, ::OptimizerState) = ""
str_updatelog(::Nothing, ::OptimizerState) = ""

display_logs_header_post(o::PartlySmoothOptimizer) = print("Update")

function display_logs_pre(::PartlySmoothOptimizer, state, pb)
    return isa(state.selected_update,WholespaceUpdate) ? "90" : "0"
end

function display_logs_post(::PartlySmoothOptimizer, state, pb)
    @printf "%-13s\t" summary(state.selected_update)

    if !isnothing(state.selected_update)
        print(str_updatelog(state.selected_update, state.update_to_updatestate[state.selected_update]))
    end
    return
end
display_logs_header_post(::PartlySmoothOptimizerState) = print("      Update")


text:
# Test the get/set mean functions
@test all(getmean(m) .== means)
setmean(m, ones(flex_data.numRVs))
@test all(getmean(m) .== ones(flex_data.numRVs))
@test_throws ErrorException setmean(m, ones(flex_data.numRVs - 1))
setmean(m, means)
@test getmean(T[1]) == 620


text:
# Bi-objective unidimensionnal 01 knapsack problem (biukp)
#
# Exercise 10.2 page 290 of
# Multicriteria Optimization (2nd edt), M. Ehrgott, Springer 2005.

# ---- Packages to use
using vOptGeneric, JuMP, GLPK

# ---- Values of the instance to solve
p1 = [10, 3,  6, 8, 2]  # coefficients's vector of the objective 1
p2 = [12, 9, 11, 5, 6]  # coefficients's vector of the objective 2
w  = [ 4, 5,  2, 5, 6]  # coefficients's vector of weights
c  = 17                 # nominal capacity
n  = length( p1 )       # number of items

# ---- setting the model
biukp = vModel( GLPK.Optimizer ) #; JuMP.set_silent( biukp )
@variable( biukp, x[1:n], Bin )
@addobjective( biukp, Max, sum( p1[j]*x[j] for j=1:n ) )
@addobjective( biukp, Max, sum( p2[j]*x[j] for j=1:n ) )
@constraint( biukp, sum( w[j]*x[j] for j=1:n ) <= c )

# ---- Invoking the solver (dichotomic method)
vSolve( biukp, method=:dichotomy, verbose=false )

# ---- Querying the results
Y_N = getY_N( biukp )

# ---- Displaying the results (X_{SE} and Y_{SN})
for i = 1:length(Y_N)
    X = value.(x, i)
    print("X = ", findall(elt -> elt ≈ 1, X))
    println(" | Z = ",Y_N[i])
end

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
import .DataFrames.DataFrame

function tabular(io::IO, A::DataFrame; alignment::String="", rounding::Int = 0)
    (r,c) = size(A)

    if length(alignment)==0
        alignment = ALIGN^c
    end

    println(io, "\\begin{tabular}{", alignment, "}")
    println(io, "\\hline")
    for (i,name) in enumerate(names(A))
        print(io, name)
        if i < length(names(A))
            print(io, " & ")
        end
    end
    println(io, get_eol(false))
    println(io, "\\hline")
    for a=1:r
        for b=1:c
            print(io, "\$",latex_form((rounding > 0) && typeof(A[a,b])<:Number ? round(A[a,b],rounding) : A[a,b]),"\$")
            if b<c
                print(io, " & ")
            else
                if a<r
                    println(io, "\\\\")
                end
            end
        end
    end
    println(io, "\n\\end{tabular}")
end

tabular(A::DataFrame; kwargs...) = tabular(stdout, A; kwargs...)


text:
## Loading packages.
using MicroHH
using Tullio
using Statistics

## Loading settings.
include("drycbl_settings.jl")

## Initialize the model.
n_domains = 1
m = Model("drycbl", n_domains, settings, float_type)

## Create the initials fields.
f = m.fields[1]; g = m.grid[1]
s = @view f.s[g.is:g.ie, g.js:g.je, g.ks:g.ke]
z = @view g.z[g.ks:g.ke]
rand2d = rand(g.itot, g.jtot)
rand2d .-= mean(rand2d)
s[:, :, 1] .+= rand2d[:, :]
f.s_gradbot[:, :] .= - 0.1 / settings[1]["fields"]["visc"]
f.s_gradtop[:, :] .= 0.003
@tullio s[i, j, k] += 0.003 * z[k]

## Save the restart files.
save_model(m)


text:
mutable struct Particle <: AbstractParticle
    pose_::Vector{Float64}
    weight_::Float64
    function Particle(pose::Vector{Float64}, weight::Float64)
        @assert length(pose) == 3
        new(copy(pose), weight)
    end
end

function Base.copy(p::Particle)
    p_ = Particle(copy(p.pose_), p.weight_)
    return p_
end

function motion_update(
    p::AbstractParticle,
    v::Float64,
    ω::Float64,
    dt::Float64,
    mv::MvNormal{Float64},
)
    noises = rand(mv)
    noised_v = v + noises[1] * sqrt(abs(v) / dt) + noises[2] * sqrt(abs(ω) / dt)
    noised_ω = ω + noises[3] * sqrt(abs(v) / dt) + noises[4] * sqrt(abs(ω) / dt)
    p.pose_ = state_transition(p.pose_, noised_v, noised_ω, dt)
end

function observation_update(
    p::Particle,
    observation::Vector{Vector{Float64}},
    envmap::Map,
    distance_dev_rate::Float64,
    direction_dev::Float64,
)
    N = size(observation)[1]
    for i = 1:N
        obsv = observation[i]
        obs_pos = obsv[1:2]
        obs_id = convert(Int64, obsv[3])

        pos_on_map = envmap.landmarks_[obs_id+1].pos_
        particle_suggest_pos = observation_function(p.pose_, pos_on_map)

        distance_dev = distance_dev_rate * particle_suggest_pos[1]
        cov = [distance_dev^2 0.0; 0.0 direction_dev^2]
        p.weight_ *= pdf(MvNormal(particle_suggest_pos, cov), obs_pos)
    end
end

mutable struct Mcl <: AbstractMcl
    particles_::Vector{Particle}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::Particle
    pose_::Vector{Float64}
    function Mcl(
        initial_pose::Vector{Float64},
        num::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05;
        glob = false, # GlobalMcl
        xlim = [-5.0, 5.0],
        ylim = [-5.0, 5.0],
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])
        if glob
            pose_distrib = PoseUniform(xlim, ylim)
            initial_pose = uniform(pose_distrib)
            new(
                [Particle(uniform(pose_distrib), 1.0 / num) for i = 1:num],
                MvNormal([0.0, 0.0, 0.0, 0.0], cov),
                distance_dev_rate,
                direction_dev,
                Particle(initial_pose, 0.0),
                copy(initial_pose),
            )
        else
            new(
                [Particle(initial_pose, 1.0 / num) for i = 1:num],
                MvNormal([0.0, 0.0, 0.0, 0.0], cov),
                distance_dev_rate,
                direction_dev,
                Particle(initial_pose, 0.0),
                copy(initial_pose),
            )
        end
    end
end

mutable struct ResetMcl <: AbstractResetMcl
    particles_::Vector{Particle}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::Particle
    pose_::Vector{Float64}
    αs_::Dict{Int64,Vector{Float64}}
    α_threshold::Float64
    reset_distrib::PoseUniform
    function ResetMcl(
        initial_pose::Vector{Float64},
        num::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05;
        α_threshold = 0.001,
        xlim = [-5.0, 5.0],
        ylim = [-5.0, 5.0],
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])
        new(
            [Particle(initial_pose, 1.0 / num) for i = 1:num],
            MvNormal([0.0, 0.0, 0.0, 0.0], cov),
            distance_dev_rate,
            direction_dev,
            Particle(initial_pose, 0.0),
            copy(initial_pose),
            Dict{Int64,Vector{Float64}}(),
            α_threshold,
            PoseUniform(xlim, ylim),
        )
    end
end

mutable struct AMcl <: AbstractResetMcl
    particles_::Vector{Particle}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::Particle
    pose_::Vector{Float64}
    reset_distrib::PoseUniform
    amcl_params::Dict{String,Float64}
    fast_α_term::Float64
    slow_α_term::Float64
    function AMcl(
        initial_pose::Vector{Float64},
        num::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05;
        xlim = [-5.0, 5.0],
        ylim = [-5.0, 5.0],
        amcl_params = Dict("slow" => 0.001, "fast" => 0.1, "v" => 3.0),
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])
        new(
            [Particle(initial_pose, 1.0 / num) for i = 1:num],
            MvNormal([0.0, 0.0, 0.0, 0.0], cov),
            distance_dev_rate,
            direction_dev,
            Particle(initial_pose, 0.0),
            copy(initial_pose),
            PoseUniform(xlim, ylim),
            amcl_params,
            1.0,
            1.0,
        )
    end
end

function set_ml(mcl::AbstractMcl)
    ind = findmax([p.weight_ for p in mcl.particles_])[2]
    mcl.ml_ = copy(mcl.particles_[ind])
    mcl.pose_ = copy(mcl.ml_.pose_)
end

function motion_update(
    mcl::Union{Mcl,ResetMcl,AMcl},
    v::Float64,
    ω::Float64,
    dt::Float64;
    kwargs...,
)
    N = length(mcl.particles_)
    for i = 1:N
        motion_update(mcl.particles_[i], v, ω, dt, mcl.motion_noise_rate_pdf)
    end
end

function observation_update(
    mcl::Mcl,
    observation::Vector{Vector{Float64}},
    envmap::Map;
    kwargs...,
)
    resample = true
    if haskey(kwargs, :resample)
        resample = convert(Bool, kwargs[:resample])
    end
    N = length(mcl.particles_)
    for i = 1:N
        observation_update(
            mcl.particles_[i],
            observation,
            envmap,
            mcl.distance_dev_rate,
            mcl.direction_dev,
        )
    end
    set_ml(mcl)
    if resample
        resampling(mcl)
    end
end

function observation_update(
    mcl::ResetMcl,
    observation::Vector{Vector{Float64}},
    envmap::Map;
    kwargs...,
)
    resample = true
    if haskey(kwargs, :resample)
        resample = convert(Bool, kwargs[:resample])
    end
    sensor_reset = true
    if haskey(kwargs, :sensor_reset)
        sensor_reset = convert(Bool, kwargs[:sensor_reset])
    end
    N = length(mcl.particles_)
    for i = 1:N
        observation_update(
            mcl.particles_[i],
            observation,
            envmap,
            mcl.distance_dev_rate,
            mcl.direction_dev,
        )
    end

    set_ml(mcl)
    α = sum([p.weight_ for p in mcl.particles_])
    if α < mcl.α_threshold
        if sensor_reset
            sensor_resetting(mcl, observation, envmap)
        end
    else
        if resample
            resampling(mcl)
        end
    end
end

function observation_update(
    mcl::AMcl,
    observation::Vector{Vector{Float64}},
    envmap::Map;
    kwargs...,
)
    N = length(mcl.particles_)
    for i = 1:N
        observation_update(
            mcl.particles_[i],
            observation,
            envmap,
            mcl.distance_dev_rate,
            mcl.direction_dev,
        )
    end
    set_ml(mcl)
    adaptive_resetting(mcl, observation, envmap)
end

function resampling(mcl::AbstractMcl)
    ws = cumsum([p.weight_ for p in mcl.particles_])

    if ws[end] < 1e-100
        ws = [e + 1e-100 for e in ws]
    end

    step = ws[end] / length(ws)
    r = rand(Uniform(0.0, step))
    cur_pos = 1
    ps = []

    while length(ps) < length(mcl.particles_)
        if r < ws[cur_pos]
            push!(ps, copy(mcl.particles_[cur_pos]))
            r += step
        else
            cur_pos += 1
        end
    end

    mcl.particles_ = copy(ps)
    for i = 1:length(ps)
        mcl.particles_[i].weight_ = 1.0 / length(ps)
    end
end

function draw(mcl::AbstractMcl, p::Plot{T}) where {T}
    xs = [p.pose_[1] for p in mcl.particles_]
    ys = [p.pose_[2] for p in mcl.particles_]
    vxs =
        [cos(p.pose_[3]) * 0.5 * p.weight_ * length(mcl.particles_) for p in mcl.particles_]
    vys =
        [sin(p.pose_[3]) * 0.5 * p.weight_ * length(mcl.particles_) for p in mcl.particles_]
    p = quiver!(xs, ys, quiver = (vxs, vys), color = "blue", alpha = 0.5)
    pose = mcl.pose_
    annota = "($(round(pose[1], sigdigits=3)), $(round(pose[2], sigdigits=3)), $(round(pose[3], sigdigits=3)))"
    p = annotate!(pose[1] + 1.0, pose[2] + 1.0, text(annota, 10))
end

function sensor_resetting(
    mcl::AbstractResetMcl,
    observation::Vector{Vector{Float64}},
    envmap::Map,
)
    d_obs = findmin([obsv[1] for obsv in observation])
    nearest_ind = d_obs[2]
    d_obs = d_obs[1]
    values = observation[nearest_ind][1:2]
    landmark_id = convert(Int64, observation[nearest_ind][3] + 1)
    @assert values[1] == d_obs
    ψ_obs = values[2]

    N = length(mcl.particles_)
    for i = 1:N
        sensor_resetting_draw(
            mcl,
            mcl.particles_[i],
            envmap[landmark_id].pos_,
            d_obs,
            ψ_obs,
        )
    end
end

function sensor_resetting_draw(
    mcl::AbstractResetMcl,
    p::Particle,
    landmark_pos::Vector{Float64},
    d_obs::Float64,
    ϕ_obs::Float64,
)
    ψ = (rand() - 0.5) * (2 * pi) # ∈ [-π, π]
    d = rand(Normal(d_obs, (mcl.distance_dev_rate * d_obs)^2))
    p.pose_[1] = landmark_pos[1] + d * cos(ψ)
    p.pose_[2] = landmark_pos[2] + d * sin(ψ)

    ϕ = rand(Normal(ϕ_obs, mcl.direction_dev^2))
    p.pose_[3] = atan(landmark_pos[2] - p.pose_[2], landmark_pos[1] - p.pose_[1]) - ϕ
    p.weight_ = 1.0 / length(mcl.particles_)
end

function adaptive_resetting(mcl::AMcl, observation::Vector{Vector{Float64}}, envmap::Map)
    if length(observation) == 0
        return
    end

    amcl_params = mcl.amcl_params
    α = sum([p.weight_ for p in mcl.particles_])
    mcl.slow_α_term += amcl_params["slow"] * (α - mcl.slow_α_term)
    mcl.fast_α_term += amcl_params["fast"] * (α - mcl.fast_α_term)
    sl_num = 1.0 - amcl_params["v"] * mcl.fast_α_term / mcl.slow_α_term
    sl_num = length(mcl.particles_) * maximum([0.0, sl_num])

    resampling(mcl)
    nearest_obs = findmin([obsv[1] for obsv in observation])
    ind = nearest_obs[2]
    obs_d, obs_ψ = observation[ind][1], observation[ind][2]
    landmark_id = convert(Int64, observation[ind][3] + 1)

    N = length(mcl.particles_)
    for n = 1:convert(Int64, floor(sl_num))
        ind = rand(1:N)
        p = mcl.particles_[ind]
        sensor_resetting_draw(mcl, p, envmap[landmark_id].pos_, obs_d, obs_ψ)
    end
end

mutable struct KldMcl <: AbstractMcl
    particles_::Vector{Particle}
    motion_noise_rate_pdf::MvNormal{Float64}
    distance_dev_rate::Float64
    direction_dev::Float64
    ml_::Particle
    pose_::Vector{Float64}
    reso::Vector{Float64}
    ϵ::Float64
    δ::Float64
    maxnum::Int64
    binnum_::Int64
    function KldMcl(
        initial_pose::Vector{Float64},
        maxnum::Int64,
        motion_noise_stds = Dict("vv" => 0.19, "vω" => 0.001, "ωv" => 0.13, "ωω" => 0.2),
        distance_dev_rate = 0.14,
        direction_dev = 0.05,
        reso = [0.2, 0.2, pi / 18],
        ϵ = 0.1,
        δ = 0.01,
    )
        v = motion_noise_stds
        cov = Diagonal([v["vv"]^2, v["vω"]^2, v["ωv"]^2, v["ωω"]^2])

        new(
            [Particle(initial_pose, 1.0)],
            MvNormal([0.0, 0.0, 0.0, 0.0], cov),
            distance_dev_rate,
            direction_dev,
            Particle(initial_pose, 1.0),
            copy(initial_pose),
            reso,
            ϵ,
            δ,
            maxnum,
            0,
        )
    end
end

function motion_update(mcl::KldMcl, v::Float64, ω::Float64, dt::Float64; kwargs...)
    ws = [p.weight_ for p in mcl.particles_]

    if sum(ws) < 1e-100
        ws = [e + 1e-100 for e in ws]
    end

    new_particles = Vector{Particle}(undef, 0)
    bins = Set{Vector{Int}}()
    for cnt = 1:mcl.maxnum
        p = sample(mcl.particles_, Weights(ws))
        motion_update(p, v, ω, dt, mcl.motion_noise_rate_pdf)
        bin = [convert(Int64, ceil(p.pose_[i] / mcl.reso[i])) for i = 1:length(p.pose_)]
        # local variables must be copied in Julia??
        push!(bins, copy(bin))
        push!(new_particles, copy(p))

        mcl.binnum_ = (length(bins) > 1) ? length(bins) : 2
        if length(new_particles) >
           ceil(chisqinvcdf(mcl.binnum_ - 1, 1.0 - mcl.δ) / (2 * mcl.ϵ))
            break
        end
    end

    mcl.particles_ = new_particles
    N = length(new_particles)
    for i = 1:N
        mcl.particles_[i].weight_ = 1.0 / N
    end
end

function observation_update(
    mcl::KldMcl,
    observation::Vector{Vector{Float64}},
    envmap::Map;
    kwargs...,
)
    resample = true
    if haskey(kwargs, :resample)
        resample = convert(Bool, kwargs[:resample])
    end
    N = length(mcl.particles_)
    for i = 1:N
        observation_update(
            mcl.particles_[i],
            observation,
            envmap,
            mcl.distance_dev_rate,
            mcl.direction_dev,
        )
    end
    set_ml(mcl)
end

function draw(mcl::KldMcl, p::Plot{T}) where {T}
    xs = [p.pose_[1] for p in mcl.particles_]
    ys = [p.pose_[2] for p in mcl.particles_]
    vxs =
        [cos(p.pose_[3]) * 0.5 * p.weight_ * length(mcl.particles_) for p in mcl.particles_]
    vys =
        [sin(p.pose_[3]) * 0.5 * p.weight_ * length(mcl.particles_) for p in mcl.particles_]
    p = quiver!(xs, ys, quiver = (vxs, vys), color = "blue", alpha = 0.5)
    pose = mcl.pose_
    annota = "($(round(pose[1], sigdigits=3)), $(round(pose[2], sigdigits=3)), $(round(pose[3], sigdigits=3)), particle: $(length(mcl.particles_)), bins: $(mcl.binnum_))"
    p = annotate!(pose[1] + 1.0, pose[2] + 1.0, text(annota, 10))
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
@testset "align minimum support sampler" begin

	cc = dummyChoiceContext()

	@testset "constructor, getparams, numparams, paramranges, setparams" begin

		subsA = DataGenerators.GeometricSampler()
		s = DataGenerators.AlignMinimumSupportSampler(subsA)

	    @testset "numparams and paramranges" begin
			@test DataGenerators.numparams(s) == DataGenerators.numparams(subsA)
			@test DataGenerators.paramranges(s) == DataGenerators.paramranges(subsA)
	    end

	    @testset "default params" begin
			@test DataGenerators.getparams(s) == DataGenerators.getparams(subsA)
	    end

		@testset "setparams"  begin
			DataGenerators.setparams!(s, [0.4])
			@test DataGenerators.getparams(subsA) == [0.4]
			@test DataGenerators.getparams(s) == [0.4]
		end

	end

	@testset "sampling" begin

		subsA = DataGenerators.GeometricSampler([0.3])
		s = DataGenerators.AlignMinimumSupportSampler(subsA)

		@mtestset "support is applied for support $cpsupport" reps=Main.REPS alpha=Main.ALPHA for cpsupport in [(17, typemax(Int)), (-4, -4), (0, 10),]
			x, trace = DataGenerators.sample(s, cpsupport, cc)
			@test cpsupport[1] <= x # this sampler does not guarantee x <= cpsupport[2]
			@mtest_distributed_as Geometric(0.3) x-cpsupport[1]
		end

	end

	@testset "estimateparams" begin

 			subs1A = DataGenerators.GeometricSampler()
			s1 = DataGenerators.AlignMinimumSupportSampler(subs1A)
			params = [0.6]
			DataGenerators.setparams!(s1, params)

			subs2A = DataGenerators.GeometricSampler()
			s2 = DataGenerators.AlignMinimumSupportSampler(subs2A)
			otherparams = [0.3]
			DataGenerators.setparams!(s2, otherparams)

			traces = map(1:100) do i
				cpsupport =[rand(-100:100) for j in 1:2]
				x, trace = DataGenerators.sample(s1, (minimum(cpsupport),maximum(cpsupport)), cc)
				trace
			end
			# note: varying support in traces - to check if re-estimation is based on underlying sample from distribution
			# and not the adjustment

			estimateparams!(s2, traces)

			@mtestset "consistent with geometric" reps=Main.REPS alpha=Main.ALPHA begin
				x, trace = DataGenerators.sample(s2, (0,1), cc)
				@mtest_distributed_as Geometric(params[1]) x
			end

	end

end


text: activate_github("QuantEcon/QuantEconLectureAllPackages", tag = "v0.9.8", add_default_environment = true)

max_stars_repo_name: sbacelar/quantecon-notebooks-jl

max_issues_repo_name: sbacelar/quantecon-notebooks-jl

max_forks_repo_name: sbacelar/quantecon-notebooks-jl