text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(model::InfiniteModel,
                               index::ReducedVariableIndex
                               )::ReducedVariableRef
    return ReducedVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::VariableData{<:ReducedVariable}
                          )::ReducedVariableIndex
    return MOIUC.add_item(model.reduced_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel, ::Type{ReducedVariable}
    )::MOIUC.CleverDict{ReducedVariableIndex, VariableData{ReducedVariable{GeneralVariableRef}}}
    return model.reduced_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(vref::ReducedVariableRef
    )::MOIUC.CleverDict{ReducedVariableIndex, VariableData{ReducedVariable{GeneralVariableRef}}}
    return JuMP.owner_model(vref).reduced_vars
end

# Extend _data_object
function _data_object(vref::ReducedVariableRef
    )::VariableData{ReducedVariable{GeneralVariableRef}}
  object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
  object === nothing && error("Invalid point variable reference, cannot find " *
                        "corresponding variable in the model. This is likely " *
                        "caused by using the reference of a deleted variable.")
  return object
end

"""
    internal_reduced_variable(vref::ReducedVariableRef,
                              key::Val{:my_ext_key})::ReducedVariable

Return the reduced variable object of `vref` assuming it is an internal variable
made during measure expansion within an optimizer model. This will apply to
optimizer model extensions that utilize `add_measure_variable` in combination
with `expand_measure`.
"""
function internal_reduced_variable end

# Extend _core_variable_object
function _core_variable_object(vref::ReducedVariableRef
    )::ReducedVariable{GeneralVariableRef}
    if !haskey(_data_dictionary(vref), JuMP.index(vref))
        model = JuMP.owner_model(vref)
        return internal_reduced_variable(vref, Val(optimizer_model_key(model)))
    else
        return _data_object(vref).variable
    end
end

# Extend _object_numbers
function _object_numbers(vref::ReducedVariableRef)::Vector{Int}
    return _core_variable_object(vref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(vref::ReducedVariableRef)::Vector{Int}
    return _core_variable_object(vref).parameter_nums
end

################################################################################
#                             DEFINITION METHODS
################################################################################
"""
    JuMP.build_variable(_error::Function, ivref::GeneralVariableRef,
                        eval_supports::Dict{Int, Float64}; [check::Bool = true]
                        )::ReducedVariable{GeneralVariableRef}

Extend the `JuMP.build_variable` function to build a reduced infinite variable
based on the infinite variable/derivative/parameter function `ivref` with reduction support `eval_supports`.
Will check that input is appropriate if `check = true`. Errors if `ivref` is
not an infinite variable, `eval_supports` violate infinite parameter domains, or
if the support dimensions don't match the infinite parameter dimensions of `ivref`.
This is intended an internal method for use in evaluating measures.
"""
function JuMP.build_variable(_error::Function, ivref::GeneralVariableRef,
                             eval_supports::Dict{Int, Float64};
                             check::Bool = true
                             )::ReducedVariable{GeneralVariableRef}
    # check the inputs
    dvref = dispatch_variable_ref(ivref)
    if check && !(dvref isa Union{InfiniteVariableRef, DerivativeRef, ParameterFunctionRef})
         _error("Must specify an infinite variable/derivative/parameter function dependency.")
    elseif check && maximum(keys(eval_supports)) > length(parameter_list(dvref))
        _error("Support evaluation dictionary indices do not match the infinite " *
               "parameter dependencies of $(ivref).")
    end
    prefs = parameter_list(dvref)
    if check
        for (index, value) in eval_supports
            pref = prefs[index]
            if JuMP.has_lower_bound(pref) && !supports_in_set(value, infinite_set(pref))
                _error("Evaluation support violates infinite parameter domain(s).")
            end
        end
    end
    # get the parameter object numbers of the dependencies
    object_nums = Int[]
    for i in eachindex(prefs)
        if !haskey(eval_supports, i)
            union!(object_nums, _object_number(prefs[i]))
        end
    end
    # get the parameter numbers
    orig_nums = _parameter_numbers(ivref)
    param_nums = [orig_nums[i] for i in eachindex(orig_nums)
                  if !haskey(eval_supports, i)]
    # round the support values in accordance with the significant digits
    for (k, v) in eval_supports
        eval_supports[k] = round(v, sigdigits = significant_digits(prefs[k]))
    end
    # build the variable
    return ReducedVariable(ivref, eval_supports, param_nums, object_nums)
end

"""
    JuMP.add_variable(model::InfiniteModel, var::ReducedVariable,
                      [name::String = ""])::GeneralVariableRef

Extend the [`JuMP.add_variable`](@ref JuMP.add_variable(::JuMP.Model, ::JuMP.ScalarVariable, ::String))
function to accomodate `InfiniteOpt` reduced variable types. Adds `var` to the
infinite model `model` and returns a [`GeneralVariableRef`](@ref).
Primarily intended to be an internal function used in evaluating measures.
"""
function JuMP.add_variable(model::InfiniteModel, var::ReducedVariable,
                           name::String = "")::GeneralVariableRef
    ivref = dispatch_variable_ref(var.infinite_variable_ref)
    JuMP.check_belongs_to_model(ivref, model)
    data_object = VariableData(var, name)
    vindex = _add_data_object(model, data_object)
    push!(_reduced_variable_dependencies(ivref), vindex)
    gvref = GeneralVariableRef(model, vindex.value, typeof(vindex))
    return gvref
end

################################################################################
#                          PARAMETER REFERENCE METHODS
################################################################################
"""
    infinite_variable_ref(vref::ReducedVariableRef)::GeneralVariableRef

Return the infinite variable/derivative/parameter function reference associated
with the reduced infinite variable `vref`.

**Example**
```julia-repl
julia> infinite_variable_ref(vref)
g(t, x)
```
"""
function infinite_variable_ref(vref::ReducedVariableRef)::GeneralVariableRef
    return _core_variable_object(vref).infinite_variable_ref
end

"""
    eval_supports(vref::ReducedVariableRef)::Dict{Int, Float64}

Return the evaluation supports associated with the reduced infinite variable
`vref`.

**Example**
```julia-repl
julia> eval_supports(vref)
Dict{Int64,Float64} with 1 entry:
  1 => 0.5
```
"""
function eval_supports(vref::ReducedVariableRef)::Dict{Int, Float64}
    return _core_variable_object(vref).eval_supports
end

"""
    raw_parameter_refs(vref::ReducedVariableRef)::VectorTuple{GeneralVariableRef}

Return the raw [`VectorTuple`](@ref) of the parameter references that `vref`
depends on. This is primarily an internal method where
[`parameter_refs`](@ref parameter_refs(vref::ReducedVariableRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(vref::ReducedVariableRef
                            )::VectorTuple{GeneralVariableRef}
    orig_prefs = raw_parameter_refs(infinite_variable_ref(vref))
    eval_supps = eval_supports(vref)
    delete_indices = [haskey(eval_supps, i) for i in eachindex(orig_prefs)]
    return deleteat!(copy(orig_prefs), delete_indices)
end

"""
    parameter_refs(vref::ReducedVariableRef)::Tuple

Return the infinite parameter references associated with the reduced infinite variable
`vref`. This is formatted as a `Tuple` of containing the parameter references as
they were inputted to define the untranscripted infinite variable except, the
evaluated parameters are excluded.

**Example**
```julia-repl
julia> parameter_refs(vref)
(t, [x[1], x[2]])
```
"""
function parameter_refs(vref::ReducedVariableRef)::Tuple
    return Tuple(raw_parameter_refs(vref))
end

"""
    parameter_list(vref::ReducedVariableRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `vref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(vref::ReducedVariableRef))
is intended as the preferred user function.
"""
function parameter_list(vref::ReducedVariableRef)::Vector{GeneralVariableRef}
    orig_prefs = parameter_list(infinite_variable_ref(vref))
    eval_supps = eval_supports(vref)
    return [orig_prefs[i] for i in eachindex(orig_prefs) if !haskey(eval_supps, i)]
end

################################################################################
#                            VARIABLE DEPENDENCIES
################################################################################
# Extend _derivative_dependencies
function _derivative_dependencies(vref::ReducedVariableRef
                                  )::Vector{DerivativeIndex}
    return _data_object(vref).derivative_indices
end

"""
    used_by_derivative(vref::ReducedVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(vref::ReducedVariableRef)::Bool
    return !isempty(_derivative_dependencies(vref))
end

"""
    is_used(vref::ReducedVariableRef)::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```julia-repl
julia> is_used(vref)
true
```
"""
function is_used(vref::ReducedVariableRef)::Bool
    if used_by_measure(vref) || used_by_constraint(vref)
        return true
    end
    for dindex in _derivative_dependencies(vref)
        if is_used(DerivativeRef(JuMP.owner_model(vref), dindex))
            return true
        end
    end
    return false
end

################################################################################
#                            VARIABLE INFO METHODS
################################################################################
"""
    JuMP.has_lower_bound(vref::ReducedVariableRef)::Bool

Extend [`JuMP.has_lower_bound`](@ref) to return a `Bool` whether the original
infinite variable of `vref` has a lower bound.

**Example**
```julia-repl
julia> has_lower_bound(vref)
true
```
"""
function JuMP.has_lower_bound(vref::ReducedVariableRef)::Bool
    return JuMP.has_lower_bound(infinite_variable_ref(vref))
end

"""
    JuMP.lower_bound(vref::ReducedVariableRef)::Float64

Extend [`JuMP.lower_bound`](@ref) to return the lower bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a lower bound.

**Example**
```julia-repl
julia> lower_bound(vref)
0.0
```
"""
function JuMP.lower_bound(vref::ReducedVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_lower_bound(ivref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return JuMP.lower_bound(ivref)
end

# Extend to return the index of the lower bound constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._lower_bound_index(vref::ReducedVariableRef)::ConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_lower_bound(ivref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return InfiniteOpt._lower_bound_index(ivref)
end

"""
    JuMP.LowerBoundRef(vref::ReducedVariableRef)::InfOptConstraintRef

Extend [`JuMP.LowerBoundRef`](@ref) to extract a constraint reference for the
lower bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = LowerBoundRef(vref)
var >= 0.0
```
"""
function JuMP.LowerBoundRef(vref::ReducedVariableRef)::InfOptConstraintRef
    return JuMP.LowerBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.has_upper_bound(vref::ReducedVariableRef)::Bool

Extend [`JuMP.has_upper_bound`](@ref) to return a `Bool` whether the original
infinite variable of `vref` has an upper bound.

**Example**
```julia-repl
julia> has_upper_bound(vref)
true
```
"""
function JuMP.has_upper_bound(vref::ReducedVariableRef)::Bool
    return JuMP.has_upper_bound(infinite_variable_ref(vref))
end

"""
    JuMP.upper_bound(vref::ReducedVariableRef)::Float64

Extend [`JuMP.upper_bound`](@ref) to return the upper bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a upper bound.

**Example**
```julia-repl
julia> upper_bound(vref)
0.0
```
"""
function JuMP.upper_bound(vref::ReducedVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_upper_bound(ivref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return JuMP.upper_bound(ivref)
end

# Extend to return the index of the upper bound constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._upper_bound_index(vref::ReducedVariableRef)::ConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_upper_bound(ivref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return InfiniteOpt._upper_bound_index(ivref)
end

"""
    JuMP.UpperBoundRef(vref::ReducedVariableRef)::InfOptConstraintRef

Extend [`JuMP.UpperBoundRef`](@ref) to extract a constraint reference for the
upper bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = UpperBoundRef(vref)
var <= 1.0
```
"""
function JuMP.UpperBoundRef(vref::ReducedVariableRef)::InfOptConstraintRef
    return JuMP.UpperBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_fixed(vref::ReducedVariableRef)::Bool

Extend [`JuMP.is_fixed`](@ref) to return `Bool` whether the original infinite
variable of `vref` is fixed.

**Example**
```julia-repl
julia> is_fixed(vref)
true
```
"""
function JuMP.is_fixed(vref::ReducedVariableRef)::Bool
    return JuMP.is_fixed(infinite_variable_ref(vref))
end

"""
    JuMP.fix_value(vref::ReducedVariableRef)::Float64

Extend [`JuMP.fix_value`](@ref) to return the fix value of the original infinite
variable of `vref`. Errors if variable is not fixed.

**Example**
```julia-repl
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::ReducedVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_fixed(ivref)
        error("Variable $(vref) is not fixed.")
    end
    return JuMP.fix_value(ivref)
end

# Extend to return the index of the fix constraint associated with the original
# infinite variable of `vref`.
function InfiniteOpt._fix_index(vref::ReducedVariableRef)::ConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_fixed(ivref)
        error("Variable $(vref) is not fixed.")
    end
    return InfiniteOpt._fix_index(ivref)
end

"""
    JuMP.FixRef(vref::ReducedVariableRef)::InfOptConstraintRef

Extend [`JuMP.FixRef`](@ref) to return the constraint reference of the fix
constraint associated with the original infinite variable of `vref`. Errors
`vref` is not fixed.

**Examples**
```julia-repl
julia> cref = FixRef(vref)
var == 1.0
```
"""
function JuMP.FixRef(vref::ReducedVariableRef)::InfOptConstraintRef
    return JuMP.FixRef(infinite_variable_ref(vref))
end

# Fallback dispatch
function JuMP.start_value(vref::ReducedVariableRef)
    return JuMP.start_value(infinite_variable_ref(vref))
end

"""
    start_value_function(vref::ReducedVariableRef)::Union{Nothing, Function}

Return the function that is used to generate the start values of `vref` for
particular support values. Returns `nothing` if no start behavior has been
specified.

**Example**
```julia-repl
julia> start_value_func(vref)
my_func
```
"""
function start_value_function(vref::ReducedVariableRef)::Union{Nothing, Function}
    return start_value_function(infinite_variable_ref(vref))
end

"""
    JuMP.is_binary(vref::ReducedVariableRef)::Bool

Extend [`JuMP.is_binary`](@ref) to return `Bool` whether the original infinite
variable of `vref` is binary.

**Example**
```julia-repl
julia> is_binary(vref)
true
```
"""
function JuMP.is_binary(vref::ReducedVariableRef)::Bool
    return JuMP.is_binary(infinite_variable_ref(vref))
end

# Extend to return the index of the binary constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._binary_index(vref::ReducedVariableRef)::ConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_binary(ivref)
        error("Variable $(vref) is not binary.")
    end
    return InfiniteOpt._binary_index(ivref)
end

"""
    JuMP.BinaryRef(vref::ReducedVariableRef)::InfOptConstraintRef

Extend [`JuMP.BinaryRef`](@ref) to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be binary.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(vref::ReducedVariableRef)::InfOptConstraintRef
    return JuMP.BinaryRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_integer(vref::ReducedVariableRef)::Bool

Extend [`JuMP.is_integer`](@ref) to return `Bool` whether the original infinite
variable of `vref` is integer.

**Example**
```julia-repl
julia> is_integer(vref)
true
```
"""
function JuMP.is_integer(vref::ReducedVariableRef)::Bool
    return JuMP.is_integer(infinite_variable_ref(vref))
end

# Extend to return the index of the integer constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._integer_index(vref::ReducedVariableRef)::ConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_integer(ivref)
        error("Variable $(vref) is not an integer.")
    end
    return InfiniteOpt._integer_index(ivref)
end

"""
    JuMP.IntegerRef(vref::ReducedVariableRef)::InfOptConstraintRef

Extend [`JuMP.IntegerRef`](@ref) to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be integer.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(vref::ReducedVariableRef)::InfOptConstraintRef
    return JuMP.IntegerRef(infinite_variable_ref(vref))
end

################################################################################
#                                  DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::ReducedVariableRef)::Nothing
    # remove mapping to infinite variable
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    filter!(e -> e != JuMP.index(vref), _reduced_variable_dependencies(ivref))
    # delete associated derivative variables and mapping
    model = JuMP.owner_model(vref)
    for index in _derivative_dependencies(vref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    return
end


text:
push!(LOAD_PATH, "../src")

using Quante
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Quante.jl" begin

#include("C60.jl")
include("SzaboOstlundReference.jl") # further references from Szabo and Ostlund
#include("vrrspeedtest.jl") # speed test of two-centre integrals

function microtest() # collection of microscopic tests, within scope of functions
    Quante.test_utils()
    Quante.test_pgbf()
    Quante.test_cgbf()
    Quante.test_overlap()
    Quante.test_kinetic()
    Quante.test_a_terms()
    Quante.test_gamma()
    Quante.test_na()
    Quante.test_fgamma()
    Quante.test_one()
    Quante.test_na2()
    Quante.test_two_terms()
    Quante.test_coul1()
    Quante.test_vrr()
    Quante.test_hrr()
    Quante.test_geo_basis()
end

microtest()

# OK! Let's try the full blown restricted Hartree Fock...
function test_h2() # via PyQuante
    @time Energy, E, U = rhf(Quante.h2,verbose=true)
    @test Energy ≈ -1.1170996
end

function test_lih() # via PyQuante
    @time Energy, E, U = rhf(Quante.lih,verbose=true)
    @test Energy ≈ -7.860745582058085 #-7.86073270525799) <-- old PyQuante version
end

function test_h2o() # via PyQuante
    @time Energy,E,U = rhf(Quante.h2o)
    @test Energy ≈ -74.9598566070913 #-74.9597609118851) <-- old PyQuante version
end

println("Testing simple molecules, sto-3g energy...")
test_h2()
test_h2o()
test_lih()

end # testset
println("All tests passed! 8^)")

ext: jl


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            ДАНО: Робот - рядом с горизонтальной перегородкой (под ней), бесконечно продолжающейся в обе стороны,
                    в которой имеется проход шириной в одну клетку.

            РЕЗУЛЬТАТ: Робот - в клетке под проходом
"""

 """Решение задачи 8"""
 function solve_8!(r)
    escape_find(r, Nord)
    show(r)
end


text:
# This file is a part of TextSearch.jl

export Vocabulary, occs, ndocs, weight, token, vocsize, trainsize

struct Vocabulary
    token::Vector{String}
    occs::Vector{Int32}
    ndocs::Vector{Int32}
    weight::Vector{Float32}
    token2id::Dict{String,UInt32}
    corpuslen::Int
end

"""
    Vocabulary(n::Integer)

Creates a `Vocabulary` struct
"""
function Vocabulary(n::Integer)
    # n == 0 means unknown
    voc = Vocabulary(String[], Int32[], Int32[], Float32[], Dict{String,UInt32}(), n)
    vocsize = ceil(Int, n^0.6)  # approx based on Heaps law
    sizehint!(voc.token, vocsize)
    sizehint!(voc.occs, vocsize)
    sizehint!(voc.ndocs, vocsize)
    sizehint!(voc.weight, vocsize)
    sizehint!(voc.token2id, vocsize)
    voc
end

function locked_tokenize_and_push(voc, textconfig, doc, buff, l)
    empty!(buff)

    for token in tokenize(textconfig, doc, buff)
        lock(l)
        try
            push!(voc, token, 1, 1, 0f0)
        finally
            unlock(l)
        end
    end
end

"""
    Vocabulary(textconfig, corpus; minbatch=0)

Computes a vocabulary from a corpus using the TextConfig `textconfig`
"""
function Vocabulary(textconfig::TextConfig, corpus::AbstractVector; minbatch=0)
    n = length(corpus)
    voc = Vocabulary(n)
    l = Threads.SpinLock()
    minbatch = getminbatch(minbatch, n)

    Threads.@threads for i in 1:n
        doc = corpus[i]

        buff = take!(CACHES)
        try
            if doc isa AbstractString
                locked_tokenize_and_push(voc, textconfig, doc, buff, l)
            else
                for text in doc
                    locked_tokenize_and_push(voc, textconfig, text, buff, l)
                end
            end
        finally
            put!(CACHES, buff)
        end
    end

    voc
end

"""
    Vocabulary(corpus)

Computes a vocabulary from an already tokenized corpus
"""
function Vocabulary(corpus)
    voc = Vocabulary(length(corpus))

    for tokens in corpus
        for token in tokens
            push!(voc, token, 1, 1, 0f0)
        end
    end

    voc
end

Base.length(voc::Vocabulary) = length(voc.occs)
Base.eachindex(voc::Vocabulary) = eachindex(voc.occs)
vocsize(voc::Vocabulary) = length(voc)
trainsize(voc::Vocabulary) = voc.corpuslen
ndocs(voc::Vocabulary, tokenID::Integer) = tokenID == 0 ? zero(eltype(voc.ndocs)) : voc.ndocs[tokenID]
occs(voc::Vocabulary, tokenID::Integer) = tokenID == 0 ? zero(eltype(voc.occs)) : voc.occs[tokenID]
weight(voc::Vocabulary, tokenID::Integer) = tokenID == 0 ? zero(eltype(voc.weight)) : voc.weight[tokenID]
token(voc::Vocabulary, tokenID::Integer) = tokenID == 0 ? "" : voc.token[tokenID]

function Base.push!(voc::Vocabulary, token::String, occs::Integer, ndocs::Integer, weight::Real)
    id = get(voc.token2id, token, zero(UInt32))

    if id == 0
        id = length(voc) + 1
        push!(voc.token, token)
        push!(voc.occs, occs)
        push!(voc.ndocs, ndocs)
        push!(voc.weight, weight)
        voc.token2id[token] = id
    else
        voc.occs[id] += occs
        voc.ndocs[id] += ndocs
        voc.weight[id] = weight
    end

    voc
end

Base.get(voc::Vocabulary, token::String, default)::UInt32 = get(voc.token2id, token, default)

function Base.getindex(voc::Vocabulary, token::String)
    getindex(voc, get(voc, token, 0))
end

function Base.getindex(voc::Vocabulary, tokenID::Integer)
    id = convert(UInt32, tokenID)

    if id == 0
        (; id=id, occs=zero(eltype(voc.occs)), ndocs=zero(eltype(voc.ndocs)), weight=zero(eltype(voc.weight)))
    else
        (; id=id, occs=voc.occs[id], ndocs=voc.ndocs[id], weight=voc.weight[id])
    end
end

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel


text:
export heisenberg, iter_groundstate!, itime_groundstate!, vqe_solve!

"""
    heisenberg(nbit::Int; periodic::Bool=true)

heisenberg hamiltonian, for its ground state, refer `PRB 48, 6141`.
"""
function heisenberg(nbit::Int; periodic::Bool=true)
    sx = i->put(nbit, i=>X)
    sy = i->put(nbit, i=>Y)
    sz = i->put(nbit, i=>Z)
    mapreduce(i->(j=i%nbit+1; sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j)), +, 1:(periodic ? nbit : nbit-1))
end

"""
    iter_groundstate!({reg::AbstractRegister}, h::AbstractBlock; niter::Int=100) -> AbstractRegister

project wave function to ground state by iteratively apply -h.
"""
iter_groundstate!(h::AbstractBlock; niter::Int=100) = reg -> iter_groundstate!(reg, h, niter=niter)
function iter_groundstate!(reg::AbstractRegister, h::AbstractBlock; niter::Int=100)
    for i = 1:niter
        reg |> h
        i%5 == 0 && reg |> normalize!
    end
    reg |> normalize!
end

"""
    itime_groundstate!({reg::AbstractRegister}, h::AbstractBlock; τ::Int=20, tol=1e-4) -> AbstractRegister

Imaginary time evolution method to get ground state, i.e. by projecting wave function to ground state by exp(-hτ). `tol` is for `expmv`.
"""
itime_groundstate!(h::AbstractBlock; τ::Real=20, tol=1e-4) = reg -> itime_groundstate!(reg, h; τ=τ, tol=tol)
function itime_groundstate!(reg::AbstractRegister, h::AbstractBlock; τ::Int=20, tol=1e-4)
    span = 1.0
    te = time_evolve(h, -im*span)
    for i = 1:τ÷span
        reg |> te |> normalize!
    end
    if τ%span != 0
        reg |> time_evolve(h, τ%span) |> normalize!
    end
    reg
end

"""
    vqe_solve!(circuit::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) -> circuit

variational quantum eigensolver, faithful simulation with optimizer Adam(lr=0.01).
"""
function vqe_solve!(circuit::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N
    optimizer = Adam(lr=0.01)
    dbs = collect_blocks(AbstractDiff, circuit)
    params = parameters(circuit)
    for i = 1:niter
        grad = opdiff.(()->zero_state(N) |> circuit, dbs, Ref(hamiltonian))
        dispatch!(circuit, update!(params, grad, optimizer))
        println("Step $i, Energy = $(expect(hamiltonian, zero_state(N) |> circuit))")
    end
    circuit
end


text:
using DandelionWebSockets: Frame
import Base: write, close

export WriterProxy

abstract type WriteRequest end;

struct WriteFrame <: WriteRequest
    frame::Frame
end

struct CloseSocket <: WriteRequest end

struct WriterProxy <: IO
    channel::Channel{WriteRequest}
    writer::IO

    function WriterProxy(writer::IO)
        proxy = new(Channel{WriteRequest}(Inf), writer)
        @async run_writerproxy(proxy)
        proxy
    end
end

function run_writerproxy(w::WriterProxy)
    try
        for writerequest in w.channel
            if writerequest == CloseSocket()
                break
            end
            write(w.writer, writerequest.frame)
        end
    catch ex
    finally
        close(w.writer)
    end
end

write(w::WriterProxy, frame::Frame) = put!(w.channel, WriteFrame(frame))
stopproxy(w::WriterProxy) = close(w.channel)
close(w::WriterProxy) = put!(w.channel, CloseSocket())


text:
struct FastReadIO{I <: IO,T} <: IO
    io::I
    ref::Base.RefValue{T}
end

FastReadIO(io::IO) = FastReadIO(io, Ref(0x00))
FastReadIO(o::FastReadIO) = o

for f in [:seekstart, :seekend, :position, :eof, :close]
    @eval Base.$f(o::FastReadIO) = Base.$f(o.io)
end
Base.seek(o::FastReadIO, pos) = seek(o.io, pos)

function read_(o::FastReadIO, ::Type{UInt8})
    read!(o.io, o.ref)
    o.ref[]
end

function read_(o::FastReadIO, ::Type{UInt32})
    b4 = UInt32(read_(o, UInt8)) << 0
    b3 = UInt32(read_(o, UInt8)) << 8
    b2 = UInt32(read_(o, UInt8)) << 16
    b1 = UInt32(read_(o, UInt8)) << 24
    b1 + b2 + b3 + b4
end

read_(o::FastReadIO, ::Type{Nothing}) = nothing
read_(o::FastReadIO, ::Type{Float32}) = reinterpret(Float32, read_(o, UInt32))
read_(o::FastReadIO, ::Type{Int32}) = reinterpret(Int32, read_(o, UInt32))
read_(o::FastReadIO, ::Type{Int8}) = reinterpret(Int8, read_(o, UInt8))
read_(o::FastReadIO, ::Type{Char}) = Char(read_(o, UInt8))

read_(io::IO, ::Type{T}) where {T} = read(io, T)
read_(io::IO, ::Type{Nothing}) = nothing

@generated function read_(io::IO, ::Type{T}) where {T <: Tuple}
    args = [:(read_(io, $Ti)) for Ti ∈ T.parameters]
    Expr(:call, :tuple, args...)
end

function bytelength(io::IO)
    # this dose the same as filesize
    # except that filesize does not work
    # on IOBuffer
    init_pos = Base.position(io)
    seekstart(io)
    start_pos = Base.position(io)
    seekend(io)
    end_pos = Base.position(io)
    seek(io, init_pos)
    end_pos - start_pos
end

function getbit(x::Integer, i)
    (x & (1 << i)) == (1 << i)
end

function setbit(x::Integer, val::Bool, i)
    T = typeof(x)
    newbit = T(val)
    mask = (T(-newbit) ⊻ x) & (T(1) << i)
    x ⊻ mask
end

function unsafe_getbyoffset(::Type{T}, s::S, offset::Int) where {T, S}
    rt = Ref{T}()
    rs = Ref{S}(s)
    GC.@preserve rt rs begin
        pt = Ptr{UInt8}(Base.unsafe_convert(Ref{T}, rt))
        ps = Ptr{UInt8}(Base.unsafe_convert(Ref{S}, rs)) + offset
        Base._memcpy!(pt, ps, sizeof(T))
    end
    return rt[]
end

function check_getbyoffset(T, s, offset)
    S = typeof(s)
    isbitstype(T) || throw(ArgumentError("Can only cast into bitstype."))
    isbitstype(S) || throw(ArgumentError("Can only cast from bitstype."))
    @assert offset >= 0
    sizeof(T) + offset <= sizeof(S) || throw(ArgumentError("Can only cast between types of equal size."))
end

function getbyoffset(T, s, offset)
    @boundscheck check_getbyoffset(T, s, offset)
    unsafe_getbyoffset(T, s, offset)
end

if VERSION < v"1.4"
    function only(itr)
        @argcheck length(itr) == 1
        first(itr)
    end
end


text:
"""
circuit optimisation
"""
module Optimise
using SimpleTraits
using YaoBlocks, YaoBlocks.ConstGate

include("to_basictypes.jl")

export replace_block
"""
    replace_block(actor, tree::AbstractBlock) -> AbstractBlock
    replace_block(pair::Pair{Type{ST}, TT}, tree::AbstractBlock) -> AbstractBlock

replace blocks in a circuit, where `actor` is a function that given input block,
returns the block to replace, or `nothing` for skip replacing and visit sibling.
If `pair` is provided, then replace original block with type `ST` with new block (`pair.second`).
"""
function replace_block(actor, tree::AbstractBlock)
    res = actor(tree)
    if res === tree || res === nothing # not replaced
        return chsubblocks(tree, replace_block.(Ref(actor), subblocks(tree)))
    else
        return res
    end
end

function replace_block(
    pair::Pair{ST,TT},
    tree::AbstractBlock,
) where {ST<:AbstractBlock,TT<:AbstractBlock}
    replace_block(x -> (x == pair.first ? pair.second : nothing), tree)
end

export is_pauli, IsPauliGroup

"""
    IsPauliGroup{X}

Trait to check if `X` is an element of Pauli group.
"""
@traitdef IsPauliGroup{X}
IsPauliGroup(x) = IsPauliGroup{typeof(x)}()

"""
    is_pauli(x)

Check if `x` is an element of pauli group.

!!! note
    this function is just a binding of `SimpleTraits.istrait`, it will not work
    if the type is not registered as a trait with `@traitimpl`.
"""
is_pauli(x::T) where {T} = SimpleTraits.istrait(IsPauliGroup{T})
is_pauli(xs...) = all(is_pauli, xs)

for G in [:I2, :X, :Y, :Z]
    ImG = Symbol(:Im, G)
    nImG = Symbol(:nIm, G)
    nG = Symbol(:n, G)

    @eval const $ImG = im * $G
    @eval const $nImG = -im * $G
    @eval const $nG = -$G

    @eval @traitimpl IsPauliGroup{typeof($G)}
    @eval @traitimpl IsPauliGroup{typeof($ImG)}
    @eval @traitimpl IsPauliGroup{typeof($nImG)}
    @eval @traitimpl IsPauliGroup{typeof($nG)}
end

export merge_pauli
merge_pauli(x) = x
merge_pauli(x::AbstractBlock, y::AbstractBlock) = x * y

function merge_pauli(ex::ChainBlock{1})
    L = length(ex)
    new_ex = chain(1)

    # find all contiguous pauli and merge them
    # note we need to iterate in inverse order
    iterm = L
    while iterm > 0
        if iterm > 1 && is_pauli(ex[iterm], ex[iterm-1])
            pushfirst!(new_ex, merge_pauli(ex[iterm], ex[iterm-1]))
            iterm = iterm - 2
        else
            # search next
            pushfirst!(new_ex, ex[iterm])
            iterm -= 1
        end
    end
    return new_ex
end

merge_pauli(::XGate, ::YGate) = ImZ
merge_pauli(::XGate, ::ZGate) = -ImY
merge_pauli(::YGate, ::XGate) = -ImZ
merge_pauli(::YGate, ::ZGate) = ImX
merge_pauli(::ZGate, ::XGate) = ImY
merge_pauli(::ZGate, ::YGate) = nImX

for G in [:X, :Y, :Z]
    GT = Symbol(G, :Gate)

    @eval merge_pauli(::I2Gate, x::$GT) = x
    @eval merge_pauli(x::$GT, ::I2Gate) = x
    @eval merge_pauli(::$GT, ::$GT) = I2
end

merge_pauli(::I2Gate, ::I2Gate) = I2

export eliminate_nested
eliminate_nested(ex::AbstractBlock) = ex

# TODO: eliminate nested expr e.g chain(X, chain(X, Y))
function eliminate_nested(ex::T) where {T<:Union{ChainBlock,Add}}
    _flatten(x) = (x,)
    _flatten(x::T) = subblocks(x)

    isone(length(ex)) && return first(subblocks(ex))
    return chsubblocks(ex, Iterators.flatten(map(_flatten, subblocks(ex))))
end

# temporary utils
_unscale(x::AbstractBlock) = x
_unscale(x::Scale) = content(x)
merge_alpha(alpha, x::AbstractBlock) = alpha
merge_alpha(alpha, x::Scale) = alpha * x.alpha
merge_alpha(alpha, x::Scale{Val{S}}) where {S} = alpha * S

# since we don't have T in blocks, this is a workaround
# to get correct identity in type stable term
merge_alpha(::Nothing, x::Scale) = x.alpha
merge_alpha(::Nothing, x::Scale{Val{S}}) where {S} = S
merge_alpha(::Nothing, x::AbstractBlock) = nothing

merge_scale(ex::AbstractBlock) = ex

# a simple function to find one for Val and Number
_one(x) = one(x)
_one(::Type{Val{S}}) where {S} = one(S)
_one(::Val{S}) where {S} = one(S)

export merge_scale

function merge_scale(ex::Union{Scale{S,N},ChainBlock{N}}) where {S,N}
    alpha = nothing
    for each in subblocks(ex)
        alpha = merge_alpha(alpha, each)
    end
    ex = chsubblocks(ex, map(_unscale, subblocks(ex)))
    if alpha === nothing
        return ex
    else
        return alpha * ex
    end
end

export combine_similar

combine_similar(ex::AbstractBlock) = ex

combine_alpha(alpha, x) = alpha
combine_alpha(alpha, x::AbstractBlock) = alpha + 1
combine_alpha(alpha, x::Scale) = alpha + x.alpha
combine_alpha(alpha, x::Scale{Val{S}}) where {S} = alpha + S

function combine_similar(ex::Add{N}) where {N}
    table = zeros(Bool, length(ex))
    list = []
    p = 1
    while p <= length(ex)
        if table[p] == true
            # checked term, skip
            p += 1
        else
            # check similar term
            term = ex[p]
            table[p] = true # mark it in the table
            alpha = 1
            for (k, each) in enumerate(ex)
                if table[k] == true # checked term, skip
                    continue
                else
                    # check if unscaled term is the same
                    # merge them if they are
                    if _unscale(term) == _unscale(each)
                        alpha = combine_alpha(alpha, each)
                        # mark checked term in the table
                        table[k] = true
                    end
                end
            end

            # eliminate zeros
            if alpha != 0
                alpha = imag(alpha) == 0 ? real(alpha) : alpha
                alpha = isinteger(alpha) ? Integer(alpha) : alpha
                push!(list, alpha * term)
            end
        end
    end

    if isempty(list)
        return Add{N}()
    else
        return Add(list...)
    end
end

export simplify

const __default_simplification_rules__ = Function[
    merge_pauli,
    eliminate_nested,
    merge_scale,
    combine_similar,
]

# Inspired by MasonPotter/Symbolics.jl
"""
    simplify(block[; rules=__default_simplification_rules__])

Simplify a block tree accroding to given rules, default to use
[`__default_simplification_rules__`](@ref).
"""
function simplify(ex::AbstractBlock; rules = __default_simplification_rules__)
    out1 = simplify_pass(rules, ex)
    out2 = simplify_pass(rules, out1)
    counter = 1
    while (out1 isa AbstractBlock) && (out2 isa AbstractBlock) && (out2 != out1)
        out1 = simplify_pass(rules, out2)
        out2 = simplify_pass(rules, out1)
        counter += 1
        if counter > 1000
            @warn "possible infinite loop in simplification rules. Breaking"
            return out2
        end
    end
    return out2
end

function simplify_pass(rules, ex)
    ex = chsubblocks(ex, map(x -> simplify_pass(rules, x), subblocks(ex)))

    for rule in rules
        ex = rule(ex)
    end
    return ex
end

end


text:
using MathOptInterface
const MOI = MathOptInterface

# Some tests are excluded because UniversalFallback accepts absolutely
# everything.

MOI.Test.runtests(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    ),
    MOI.Test.Config(),
    exclude = [
        "test_model_ScalarFunctionConstantNotZero",
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
        "test_model_supports_constraint_VariableIndex_EqualTo",
        "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
    ],
    warn_unsupported = true,
)

# Run the previously excluded tests, this time without UniversalFallback.

MOI.Test.runtests(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.Model{Float64}(),
        scalar_function_constant_non_zero = true,
    ),
    MOI.Test.Config(),
    include = [
        "test_model_ScalarFunctionConstantNotZero",
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
        "test_model_supports_constraint_VariableIndex_EqualTo",
        "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
    ],
)

# Test for Issue #1757

MOI.Test.test_model_ScalarFunctionConstantNotZero(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.Model{Float64}(),
        scalar_function_constant_non_zero = false,
    ),
    MOI.Test.Config(exclude = Any[MOI.ConstraintFunction]),
)

# Test exclude_tests_after. This should work despite no methods being added for IncompleteOptimizer
# because every test should get skipped.

struct IncompleteOptimizer <: MOI.AbstractOptimizer end

MOI.Test.runtests(
    IncompleteOptimizer(),
    MOI.Test.Config();
    exclude_tests_after = v"0.0.1",
)

# Non-Float64 tests

MOI.Test.runtests(
    MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{BigFloat}()),
        BigFloat,
    ),
    MOI.Test.Config(BigFloat),
    exclude = [
        # ========================= Expected failures ==========================
        # UniversalFallback supports these tests.
        "test_model_copy_to_UnsupportedAttribute",
        "test_model_copy_to_UnsupportedConstraint",
        "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
        "test_model_supports_constraint_VariableIndex_EqualTo",
        "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
    ],
)


text:
#-------------------------------------------------------------------------------
# Entry point utilities
#
# These make it easy for users to open `DataSet`s and map them into types
# understood by their program.

function extract_dtypes(call)
    dtypes = []
    jtypes = []
    argnames = []
    for ex in call.args[2:end]
        @assert ex.head == :call && ex.args[1] == :(=>)
        @assert ex.args[2].head == :(::) && length(ex.args[2].args) == 2
        push!(argnames, ex.args[2].args[1])
        push!(dtypes, ex.args[2].args[2])
        push!(jtypes, ex.args[3])
    end
    argnames, dtypes, jtypes
end

"""
    @datafunc function f(x::DT=>T, y::DS=>S...)
        ...
    end

Define the function `f(x::T, y::S, ...)` and add data dispatch rules so that
`f(x::DataSet, y::DataSet)` will open datasets matching dataset types `DT,DS`
as Julia types `T,S`.
"""
macro datafunc(func_expr)
    @assert func_expr.head == :function
    call = func_expr.args[1]
    body = func_expr.args[2]
    funcname = call.args[1]
    argnames, dtypes, jtypes = extract_dtypes(call)
    real_args = [:($n::$t) for (n,t) in zip(argnames, jtypes)]
    table_name = Symbol("#_$(funcname)_datasets_dispatch")
    esc_funcname = esc(funcname)
    esc_table_name = esc(table_name)
    func_expr.args[1].args[2:end] = real_args
    quote
        if !$(esc(:(@isdefined($table_name))))
            function $esc_funcname(ds::DataSet...)
                _run($esc_funcname, $esc_table_name, ds...)
            end
            const $esc_table_name = Dict()
        end
        push!($esc_table_name, tuple($(map(string, dtypes)...)) =>
                               tuple($(map(esc, jtypes)...)))
        $(esc(func_expr))
    end
end

function datarun(proj::AbstractDataProject, func::Function, data_names::AbstractString...)
    ds = map(n->dataset(proj, n), data_names)
    func(ds...)
end

"""
    @datarun [proj] func(args...)

Run `func` with the named `DataSet`s from the list `args`.

# Example

Load `DataSet`s named a,b as defined in Data.toml, and pass them to `f()`.
```
proj = DataSets.load_project("Data.toml")
@datarun proj f("a", "b")
```
"""
macro datarun(args...)
    if length(args) == 2
        proj, call = args
        esc_proj = esc(proj)
    elseif length(args) == 1
        esc_proj = :PROJECT
        call = args[1]
    else
        throw(ArgumentError("@datarun macro expects one or two arguments"))
    end
    esc_funcname = esc(call.args[1])
    esc_funcargs = esc.(call.args[2:end])
    quote
        datarun($esc_proj, $esc_funcname, $(esc_funcargs...))
    end
end

"""
    dataset_type(dataset)

Get a string representation of the "DataSet type", which represents the type of
the data *outside* Julia.

A given DataSet type may be mapped into many different Julia types. For example
consider the "Blob" type which is an array of bytes (commonly held in a file).
When loaded into Julia, this may be represented as a
    * IO             — via open())
    * String         — via open() |> read(_,String)
    * Vector{UInt8}  — via mmap)
    * Path
"""
function dataset_type(d::DataSet)
    # TODO: Enhance this once maps can be applied on top of the storage layer
    # Should we use MIME type? What about layering?
    d.storage["type"]
end

function _openall(func, opened, (dataset,T), to_open...)
    open(T, dataset) do newly_opened
        _openall(func, (opened..., newly_opened), to_open...)
    end
end

function _openall(func, opened)
    func(opened...)
end

# Match `dataset_type` of `ds` against `dispatch_table`, using the match to
# determine the appropriate Julia types we will open.
function _run(func, dispatch_table, ds::DataSet...)
    # For now, uses a simplistic exact matching strategy. We don't use Julia's
    # builtin dispatch here because
    # a) It seems wasteful to create a pile of tag types just for the purposes
    #    of matching some strings
    # b) It seems like a good idea to separate the declarative "data
    #    typesystem" (implicitly defined outside Julia) from Julia's type
    #    system and dispatch rules.
    dtypes = dataset_type.(ds)
    if !haskey(dispatch_table, dtypes)
        table = join(string.(keys(dispatch_table)), "\n")
        throw(ArgumentError("""No matching function $func for DataSet types $dtypes.

                            The types must match one of the following:
                            $table
                            """))
    end
    julia_types = dispatch_table[dtypes]
    to_open = Pair.(ds, julia_types)
    _openall(func, (), to_open...)
end


text:
module NeuroJulia
    export NeuroData,NeuroAdmin

    using JSON
    using Requests
    using MbedTLS

    global token = ENV["JUPYTER_TOKEN"]
    global domain = ENV["NV_DOMAIN"]
    domain = if contains(domain,"prd")
        #"https://neuroverse.com.au"
        "https://15ded47f-ef38-4ee3-b989-685820ca3d36.cloudapp.net"
    elseif contains(domain,"tst")
        "https://neuroqa.d3s.com.au"
    elseif contains(domain,"sit")
        "https://neurosit.d3s.com.au"
    elseif contains(domain,"dev")
        "https://neurodev.d3s.com.au"
    else
        "http://localhost"
    end
    global homedir = "/home/jovyan/session/"

    prepo = LibGit2.GitRepo(abspath(Pkg.dir(), "NeuroJulia"))
    phead = LibGit2.head(prepo)
    global branchname = LibGit2.shortname(phead)

    function neurocall(port,service,method,requestbody;timeout=1200)
        url = domain * ":8080/NeuroApi/" * port * "/" * service * "/api/" * replace(lowercase(service),"service","") * "/" * method
        if domain=="http://localhost"
            url = domain * ":8082/NeuroApi/" * port * "/" * service * "/api/" * replace(lowercase(service),"service","") * "/" * method
        end
        msgdata = nothing
        msgdatalength = 0
        if requestbody!=nothing
            msgdata = JSON.json(requestbody)
            msgdatalength = length(msgdata)
        end
        headers = Dict("Content-Length" => string(msgdatalength), "Token" => token, "Accept" => "application/json", "Content-Type" => "application/json")
        response = post(url; headers=headers, data=msgdata, timeout=timeout, tls_conf=MbedTLS.SSLConfig(false))
        if response.status != 200
            if response.status == 401
                error("Session has expired: Log into Neuroverse and connect to your Notebooks session or reload the Notebooks page in Neuroverse")
            elseif response.status == 404
                error(replace(lowercase(service),"service","") * "/" * method * " does not exist")
            else
                error("Neuroverse connection error: Http code " * string(response.status))
            end
        end
        responseobj = JSON.parse(readstring(response))
        if responseobj["ErrorCode"] != 0
            error("Neuroverse Error: " * responseobj["Error"])
        end
        return responseobj
    end

    function neurocall(service,method,requestbody;timeout=1200)
        url = domain * ":8080/NeuroApi/8080/" * service * "/api/" * replace(lowercase(service),"service","") * "/" * method
        if domain=="http://localhost"
            url = domain * ":8082/NeuroApi/8080/" * service * "/api/" * replace(lowercase(service),"service","") * "/" * method
        end
        msgdata = nothing
        msgdatalength = 0
        if requestbody!=nothing
            msgdata = JSON.json(requestbody)
            msgdatalength = length(msgdata)
        end
        headers = Dict("Content-Length" => string(msgdatalength), "Token" => token, "Accept" => "application/json", "Content-Type" => "application/json")
        response = post(url; headers=headers, data=msgdata, timeout=timeout, tls_conf=MbedTLS.SSLConfig(false))
        if response.status != 200
            if response.status == 401
                error("Session has expired: Log into Neuroverse and connect to your Notebooks session or reload the Notebooks page in Neuroverse")
            else
                error("Neuroverse connection error: Http code " * string(response.status))
            end
        end
        responseobj = JSON.parse(readstring(response))
        if responseobj["ErrorCode"] != 0
            error("Neuroverse Error: " * responseobj["Error"])
        end
        return responseobj
    end

    function get_notebook(filename)
        directory=homedir * "00_NeuroTemplates"
        try
            run(`mkdir $directory`)
        end
        filename=split(filename,'.')[1] * ".ipynb"
        output=directory * "/" * split(filename,'.')[1] * "_" * replace(split(string(Dates.now()),'.')[1],':','_') * ".ipynb"
        run(`curl https://raw.githubusercontent.com/SnowdenNeuroverse/NeuroNotebooks/$branchname/Notebooks/$filename --output $output`)
    end

    include(Pkg.dir() * "/NeuroJulia/src/NeuroData.jl")
    include(Pkg.dir() * "/NeuroJulia/src/NeuroAdmin.jl")
end # module

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


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

mutable struct Dependency <: SwaggerModel
    dependsOn::Any # spec type: Union{ Nothing, Vector{BasicDependency} } # spec name: dependsOn
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    resourceType::Any # spec type: Union{ Nothing, String } # spec name: resourceType
    resourceName::Any # spec type: Union{ Nothing, String } # spec name: resourceName

    function Dependency(;dependsOn=nothing, id=nothing, resourceType=nothing, resourceName=nothing)
        o = new()
        validate_property(Dependency, Symbol("dependsOn"), dependsOn)
        setfield!(o, Symbol("dependsOn"), dependsOn)
        validate_property(Dependency, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(Dependency, Symbol("resourceType"), resourceType)
        setfield!(o, Symbol("resourceType"), resourceType)
        validate_property(Dependency, Symbol("resourceName"), resourceName)
        setfield!(o, Symbol("resourceName"), resourceName)
        o
    end
end # type Dependency

const _property_map_Dependency = Dict{Symbol,Symbol}(Symbol("dependsOn")=>Symbol("dependsOn"), Symbol("id")=>Symbol("id"), Symbol("resourceType")=>Symbol("resourceType"), Symbol("resourceName")=>Symbol("resourceName"))
const _property_types_Dependency = Dict{Symbol,String}(Symbol("dependsOn")=>"Vector{BasicDependency}", Symbol("id")=>"String", Symbol("resourceType")=>"String", Symbol("resourceName")=>"String")
Base.propertynames(::Type{ Dependency }) = collect(keys(_property_map_Dependency))
Swagger.property_type(::Type{ Dependency }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Dependency[name]))}
Swagger.field_name(::Type{ Dependency }, property_name::Symbol) =  _property_map_Dependency[property_name]

function check_required(o::Dependency)
    true
end

function validate_property(::Type{ Dependency }, name::Symbol, val)
end


text:
using IterativeSolvers
using Test
using Random
using LinearAlgebra

@testset "IDR(s)" begin

n = 10
m = 6
Random.seed!(1234567)

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n) + n * I
    b = rand(T, n)
    tol = √eps(real(T))

    @testset "Without residual smoothing" begin
        x, history = idrs(A, b, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ tol
    end

    # with smoothing
    @testset "With residual smoothing" begin
        x, history = idrs(A, b; smoothing=true, log=true)
        @test history.isconverged
        @test norm(A*x - b) / norm(b) ≤ tol
    end
end

@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
    A = sprand(T, n, n, 0.5) + n * I
    b = rand(T, n)
    tol = √eps(real(T))

    x, history = idrs(A, b, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ tol
end

@testset "Maximum number of iterations" begin
    x, history = idrs(rand(5, 5), rand(5), log=true, maxiter=2)
    @test history.iters == 2
    @test length(history[:resnorm]) == 2
end
end


text:
function random_greedy_color(g::AbstractGraph{T}, reps::Integer) where T <: Integer
    best = @distributed (LightGraphs.best_color) for i in 1:reps
        seq = shuffle(vertices(g))
        LightGraphs.perm_greedy_color(g, seq)
    end

    return convert(LightGraphs.coloring{T}, best)
end

greedy_color(g::AbstractGraph{U}; sort_degree::Bool=false, reps::Integer=1) where {U <: Integer} =
sort_degree ? LightGraphs.degree_greedy_color(g) : Parallel.random_greedy_color(g, reps)


text:
@info "Caesar.jl is loading tools using ImageMagick.jl."

export toFormat
export fetchDataImage

"""
    $SIGNATURES

Convert image to PNG bytestream.

Notes
- `using ImageMagick, FileIO`

DevNotes
- TODO somehow consolidate with `MIME"image/png"`.

See also: [`makeImage`](@ref)
"""
function toFormat(format::DataType,
                  img::AbstractMatrix{<:Colorant} )
  #
  io = IOBuffer()
  pngSm = Stream(format, io)
  save(pngSm, img)  # think FileIO is required for this
  take!(io)
end

toFormat(img::AbstractMatrix{<:Colorant}) = toFormat(format"PNG", img)

"""
    $SIGNATURES

`Data: Entry => Blob` helper function to load images stored in standard (png, jpg, jpeg) format from supported a DFG data blob store.

Example
```julia
# Skip if AbstractBlobStore is already set up
  ## if a FolderStore already exists (assumed inside `getLogPath(fg)`)
  storeDir = joinLogPath(fg,"data")
  datastore = FolderStore{Vector{UInt8}}(:default_folder_store, storeDir)
  addBlobStore!(fg, datastore)

# Fetch the image
img = fetchDataImage(fg, :x4, :KEYFRAME_IMG)
```

Notes
- https://juliarobotics.org/Caesar.jl/latest/concepts/interacting_fgs/#Retrieving-a-Data-Blob
"""
function fetchDataImage(dfg::AbstractDFG,
                        varLbl::Symbol,
                        dataLbl::Symbol,
                        getDataLambda::Function = (g,vl,dl) -> getData(g,vl,dl),
                        checkMimeType::Bool=true )
#
  imgEntry, imgBytes = getDataLambda(dfg, varLbl, dataLbl)
  allowList = ["image/png"; "image/jpg"; "image/jpeg"]
  checkMimeType && (@assert imgEntry.mimeType in allowList "Should we allow image format DataBlobEntry.mimeType=$(imgEntry.mimeType), current allowList=$allowList")
  ImageMagick.readblob(imgBytes)
end
fetchDataImage(dfg::AbstractDFG,datastore::AbstractBlobStore,varLbl::Symbol,dataLbl::Symbol,checkMimeType::Bool=true) = fetchDataImage(dfg, varLbl, dataLbl, (g,vl,dl) -> getData(g,datastore,vl,dl) , checkMimeType)

#


text:
"""
	BafflingLongitudinalTracker <: LaneFollowingDriver

Extended Intelligent Driver Model. A rule based driving model that is governed by parameter
settings. The output is an longitudinal acceleration. The desired velocity, however, is purturbed within an offset range for practicality.

# Fields
- `a::Float64 = NaN` the predicted acceleration i.e. the output of the model
- `σ::Float64 = NaN` allows errorable IDM, optional stdev on top of the model, set to zero or NaN for deterministic behavior
- `k_spd::Float64 = 1.0` proportional constant for speed tracking when in freeflow [s⁻¹]
- `δ::Float64 = 4.0` acceleration exponent
- `T::Float64  = 1.5` desired time headway [s]
- `v_des::Float64 = 29.0` desired speed [m/s]
- `s_min::Float64 = 5.0` minimum acceptable gap [m]
- `a_max::Float64 = 3.0` maximum acceleration ability [m/s²]
- `d_cmf::Float64 = 2.0` comfortable deceleration [m/s²] (positive)
- `d_max::Float64 = 9.0` maximum deceleration [m/s²] (positive)
- `v_offset:Float64 = 5.0` offset from the desired speed [m/s]
"""
@with_kw mutable struct BafflingLongitudinalTracker <: LaneFollowingDriver
    a::Float64 = NaN # predicted acceleration
    σ::Float64 = NaN # optional stdev on top of the model, set to zero or NaN for deterministic behavior

    k_spd::Float64 = 1.0 # proportional constant for speed tracking when in freeflow [s⁻¹]

    δ::Float64 = 4.0 # acceleration exponent [-]
    T::Float64  = 1.5 # desired time headway [s]
    v_des::Float64 = 29.0 # desired speed [m/s]
    s_min::Float64 = 5.0 # minimum acceptable gap [m]
    a_max::Float64 = 3.0 # maximum acceleration ability [m/s²]
    d_cmf::Float64 = 2.0 # comfortable deceleration [m/s²] (positive)
    d_max::Float64 = 9.0 # maximum deceleration [m/s²] (positive)
    ΔT::Float64 = 0.2 # timestep to simulate [s]
    v_offset::Float64 = 5.0 # offset from the desired speed [m/s]
end
get_name(::BafflingLongitudinalTracker) = "BafflingLongitudinalTracker"
function AutomotiveDrivingModels.set_desired_speed!(model::BafflingLongitudinalTracker, v_des::Float64)
    model.v_des = v_des
    model
end
function AutomotiveDrivingModels.track_longitudinal!(model::BafflingLongitudinalTracker, v_ego::Float64, v_oth::Float64, headway::Float64)
    if !isnan(v_oth)
        @assert !isnan(headway)
        if headway < 0.0
            @debug("BafflingLongitudinalTracker Warning: BDM received a negative headway $headway"*
                  ", a collision may have occured.")
            model.a = max(-model.d_max, -v_ego/model.ΔT)
        else

            Δv = v_oth - v_ego
            s_des = model.s_min + v_ego*model.T - v_ego*Δv / (2*sqrt(model.a_max*model.d_cmf))
            v_ratio = model.v_des > 0.0 ? (v_ego/model.v_des) : 1.0
            model.a = model.a_max * (1.0 - v_ratio^model.δ - (s_des/headway)^2)

            if v_ego + model.ΔT * model.a < 0
                model.a = max(-model.d_max, -v_ego/model.ΔT)
            end
        end
    else
        # no lead vehicle, just drive to match desired speed
        Δv = (model.v_des + Base.rand(Float64(model.v_offset):Float64(model.v_offset)))- v_ego
        model.a = Δv*model.k_spd # predicted accel to match target speed
    end

    @assert !isnan(model.a)

    model.a = clamp(model.a, -model.d_max, model.a_max)

    return model
end
function Base.rand(model::BafflingLongitudinalTracker)
    if isnan(model.σ) || model.σ ≤ 0.0
        LaneFollowingAccel(model.a)
    else
        LaneFollowingAccel(rand(Normal(model.a, model.σ)))
    end
end
function Distributions.pdf(model::BafflingLongitudinalTracker, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a.a)
    end
end
function Distributions.logpdf(model::BafflingLongitudinalTracker, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a.a)
    end
end


text:
# Code to get PiecewiseDeterministicMarkovProcesses working with ForwardDiff
# Makes use of the trick in http://docs.juliadiffeq.org/latest/basics/faq.html#Are-the-native-Julia-solvers-compatible-with-autodifferentiation?-1
# https://diffeq.sciml.ai/stable/basics/faq/
using ForwardDiff

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
	rate::T
	dual_rate::S
end

function DiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size}}) where {T, chunk_size}
	DiffCache(u, zeros(ForwardDiff.Dual{nothing,T,chunk_size}, siz...))
end

dualcache(u::AbstractArray, N=Val{ForwardDiff.pickchunksize(length(u))}) = DiffCache(u, size(u), N)

# this is from the trick above. It fails here because x and dc.rate do not have the same dimension
get_rate(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual = reinterpret(T, dc.dual_rate)
# get_rate(dc::DiffCache, u::AbstractArray{T}) where T <: ForwardDiff.Dual = (dc.dual_rate)
get_rate(dc::DiffCache, u::AbstractArray) = dc.rate


text:
module Bfs

    export visit

    function visit(start::Number, nodes::Dict{Int,Array{Int,1}}, fn::Function)
        frontier = Number[start]
        visited  = Number[]

        function next_visit(node::Number)
            filter(n -> !in(n, visited), nodes[node])
        end

        while 0 < length(frontier)
            next = Number[]
            for node = frontier
                push!(visited, node)
                append!(next, next_visit(node))
                fn(node)
            end
            frontier = next
        end
    end

end

max_stars_repo_name: adlawson/search-algorithms

max_issues_repo_name: adlawson/search-algorithms

max_forks_repo_name: adlawson/search-algorithms


text:
__precompile__()
module ReduceAlgebra

#   This file is part of ReduceAlgebra.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

using Reduce
using ReduceLinAlg

__init__() = nothing

end # module


text:
struct RN <: UnmarkedModel
  data::UmData
  opt::UnmarkedOpt
  submodels::NamedTuple
  K::Int
end

"""
    rn(λ_formula::FormulaTerm, p_formula::FormulaTerm, data::UmData; K::Int)

Fit the mixture model of Royle and Nichols (2003), which estimates latent site
abundance using presence-absence data. Probability a species is detected at a
site is assumed to be positively related to species abundance at the site.
Covariates on abundance and detection are specified with `λ_formula` and
`p_formula`, respectively. The `UmData` object should contain `site_covs` and
`obs_covs` `DataFrames` containing columns with names matching the formulas.
The likelihood is marginalized over possible latent abundance values `0:K` for
each site. `K` should be set high enough that the likelihood of true abundance
`K` for any site is ≈ 0 (defaults to maximum observed count + 20). Larger `K`
will increase runtime.
"""
function rn(λ_formula::FormulaTerm, p_formula::FormulaTerm,
            data::UmData; K::Union{Int,Nothing}=nothing)

  abun = UmDesign(:Abundance, λ_formula, LogLink(), data.site_covs)
  det = UmDesign(:Detection, p_formula, LogitLink(), data.obs_covs)
  add_idx!([abun, det])
  np = get_np([abun, det])

  y = data.y
  N, J = size(y)
  K = check_K(K, y)
  Kmin = get_Kmin(y)

  function loglik(β::Array)

    λ = transform(abun, β)
    r = transform(det, β)

    ll = zeros(eltype(β), N)

    Threads.@threads for i = 1:N
      rsub = r[(i*J-J+1):(i*J)]
      fg = zeros(eltype(β),K+1)
      for k = Kmin[i]:K
        fk = pdf(Poisson(λ[i]), k)
        gk = 1.0
        for j = 1:J
          if ismissing(y[i,j]) continue end
          p = 1 - (1 - rsub[j])^k
          gk *= p^y[i,j] * (1-p)^(1-y[i,j])
        end
        fg[k+1] = fk * gk
      end
      ll[i] = log(sum(fg))
    end

    return -sum(ll)
  end

  opt = optimize_loglik(loglik, np)
  RN(data, opt, (abun=UnmarkedSubmodel(abun, opt),
                 det=UnmarkedSubmodel(det, opt)), K)
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
struct Model{T <: Real, N <: Integer,A,W, M <:AbstractMeltRate} <: AbstractModel{T,N,M}
    grid::Grid{T,N}
    params::Params{T,A,W}
    solver_params::SolverParams{T,N}
    initial_conditions::InitialConditions{T}
    fields::Fields{T,N}
    melt_rate::M
end

"""
    Model(;
        grid = nothing,
        bed_elevation = nothing,
        params = Params(),
        solver_params = SolverParams(),
        initial_conditions = InitialConditions(),
        melt_rate = UniformMeltRate())

Construct a WAVI.jl model object.

Keyword arguments
=================

    - `grid`: (required) an instance of a `Grid` object, which defines the computational grid
    - `bed_elevation`: (required) an array of size `grid.nx` x `grid.ny` which defines the bed elevation
    - `params`: a `Params` object that defines physical parameters
    - `solver_params`: a `SolverParams` object that defines parameters relating to the numerical scheme
    - `initial_conditions`: an `InitialConditions` object that (optionally) defines the initial ice thickness, temperature, viscosity, and damage
    - `melt_rate`: a melt rate model, responsible for controlling/setting the basal melt rate
"""
function Model(;
    grid = nothing,
    bed_elevation = nothing,
    params = Params(),
    solver_params = SolverParams(),
    initial_conditions = InitialConditions(),
    melt_rate = UniformMeltRate())

    #check that a grid and bed has been inputted
    ~(grid === nothing) || throw(ArgumentError("You must specify an input grid"))
    ~(bed_elevation === nothing) || throw(ArgumentError("You must input a bed elevation"))

    #check types
    #if a functional bed has been specified, convert to an array
    bed_array = zeros(grid.nx,grid.nx) #initialize a bed array
    try
    bed_array = get_bed_elevation(bed_elevation, grid)
    catch
    throw(ArgumentError("bed elevation must be of type function or array"))
    end
            #check the size of the bed
    #(Base.size(bed_array) = (grid.nx, grid.ny)) || throw(ArgumentError("Bed and grid sizes must be identical"))

    #Check initial conditions, and revert to default values if not
    initial_conditions = check_initial_conditions(initial_conditions, params, grid)

    ## Parameter fields checks
    #if weertman c passed as a scalar, replace weertman_c parameters with matrix of this value
    if isa(params.weertman_c, Number)
        params = @set params.weertman_c = params.weertman_c*ones(grid.nx,grid.ny)
    end
    #check size compatibility of resulting weertman C
    (size(params.weertman_c)==(grid.nx,grid.ny)) || throw(DimensionMismatch("Size of input weertman c must match grid size (i.e. $(grid.nx) x $(grid.ny))"))

    #if accumulation is passed as a scalar, replace accumulation parameters with matrix of this value
    if isa(params.accumulation_rate, Number)
        params = @set params.accumulation_rate = params.accumulation_rate*ones(grid.nx,grid.ny)
    end
    #check size compatibility of resulting weertman C
    (size(params.accumulation_rate)==(grid.nx,grid.ny)) || throw(DimensionMismatch("Size of input weertman c must match grid size (i.e. $(grid.nx) x $(grid.ny))"))

    #Setup the fields
    fields = setup_fields(grid, initial_conditions, solver_params, params, bed_array)

    #Use type constructor to build initial state with no extra physics
    model=Model(grid,params,solver_params,initial_conditions,fields,melt_rate)

    return model
end

include("model_utilities.jl")
include("update_state.jl")
include("update_velocities.jl")


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>7.0, :alg=>TMJets()),
    "X0"   => Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45]),
    "info" => Dict(
        "name" => "vanderpol",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
"""
    BinomialExpert(n, p)

PMF:

```math
P(X = k) = {n \\choose k}p^k(1-p)^{n-k},  \\quad \\text{ for } k = 0,1,2, \\ldots, n.
```

See also: [Binomial Distribution](https://en.wikipedia.org/wiki/Binomial_distribution) (Wikipedia)

"""
struct BinomialExpert{T<:Real} <: NonZIDiscreteExpert
    n::Int
    p::T
    BinomialExpert{T}(n, p) where {T<:Real} = new{T}(n, p)
end

function BinomialExpert(n::Integer, p::T; check_args=true) where {T <: Real}
    check_args && @check_args(BinomialExpert, 0 <= p <= 1 && isa(n, Integer))
    return BinomialExpert{T}(n, p)
end

## Outer constructors
# BinomialExpert(n::Integer, p::Real) = BinomialExpert(n, float(p))
BinomialExpert(n::Integer, p::Integer) = BinomialExpert(n, float(p))
BinomialExpert() = BinomialExpert(2, 0.50)

## Conversion
function convert(::Type{BinomialExpert{T}}, n::Int, p::S) where {T <: Real, S <: Real}
    BinomialExpert(n, T(p))
end
function convert(::Type{BinomialExpert{T}}, d::BinomialExpert{S}) where {T <: Real, S <: Real}
    BinomialExpert(d.n, T(d.p), check_args=false)
end
copy(d::BinomialExpert) = BinomialExpert(d.n, d.p, check_args=false)

## Loglikelihood of Expoert
logpdf(d::BinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.logpdf.(Distributions.Binomial(d.n, d.p), x...)
pdf(d::BinomialExpert, x...) = isinf(x...) ? -Inf : Distributions.pdf.(Distributions.Binomial(d.n, d.p), x...)
logcdf(d::BinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.logcdf.(Distributions.Binomial(d.n, d.p), x...)
cdf(d::BinomialExpert, x...) = isinf(x...) ? 1.0 : Distributions.cdf.(Distributions.Binomial(d.n, d.p), x...)

## expert_ll, etc
expert_ll_exact(d::BinomialExpert, x::Real) = LRMoE.logpdf(d, x)
function expert_ll(d::BinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, ceil.(yl) .- 1) - logcdf.(d, yu))
    return expert_ll
end
function expert_tn(d::BinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu))
    return expert_tn
end
function expert_tn_bar(d::BinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar = (tl == tu) ? log1mexp.(logpdf.(d, tl)) : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::BinomialExpert; exposure = 1) = d

## Parameters
params(d::BinomialExpert) = (d.n, d.p)
function params_init(y, d::BinomialExpert)
    n_init = Int(maximum(vec(y))) + 2
    p_init = mean(y) / n_init
    try
        return BinomialExpert(n_init, p_init)
    catch;
        return BinomialExpert()
    end
end

## Simululation
sim_expert(d::BinomialExpert) = Distributions.rand(Distributions.Binomial(d.n, d.p), 1)[1]

## penalty
penalty_init(d::BinomialExpert) = []
no_penalty_init(d::BinomialExpert) = []
penalize(d::BinomialExpert, p) = 0.0

## statistics
mean(d::BinomialExpert) = mean(Distributions.Binomial(d.n, d.p))
var(d::BinomialExpert) = var(Distributions.Binomial(d.n, d.p))
quantile(d::BinomialExpert, p) = quantile(Distributions.Binomial(d.n, d.p), p)

## Misc functions for E-Step

function _sum_densy_series(d::BinomialExpert, yl, yu)
    if isinf(yu)
        series = 0:(max(yl-1, 0))
        return d.n*d.p - sum(pdf.(d, series) .* series)[1]
    else
        series = yl:(min(yu, d.n))
        return sum(pdf.(d, series) .* series)[1]
    end
end

function _int_obs_Y_raw(d::BinomialExpert, yl, yu)
    return _sum_densy_series(d, yl, yu)
end

function _int_lat_Y_raw(d::BinomialExpert, tl, tu)
    return max(d.n*d.p - _sum_densy_series(d, tl, tu), 0)
end

## EM: M-Step
function EM_M_expert(d::BinomialExpert,
                    tl, yl, yu, tu,
                    exposure,
                    z_e_obs, z_e_lat, k_e;
                    penalty = true, pen_pararms_jk = [2.0 1.0])

    # Not affected by Exposurize
    expert_ll_pos = expert_ll.(d, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(d, tl, yl, yu, tu)

    # Further E-Step
    yl_yu_unique = unique_bounds(yl, yu)
    int_obs_Y_tmp = _int_obs_Y_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    Y_e_obs = exp.(-expert_ll_pos) .* int_obs_Y_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(Y_e_obs, 0.0) # get rid of NaN
    inf2num(Y_e_obs, 0.0) # get rid of inf

    tl_tu_unique = unique_bounds(tl, tu)
    int_lat_Y_tmp = _int_lat_Y_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    Y_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_Y_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(Y_e_lat, 0.0) # get rid of NaN
    inf2num(Y_e_lat, 0.0) # get rid of inf

    # Update parameters
    term_zkz_Y = (z_e_obs .* Y_e_obs) .+ (z_e_lat .* k_e .* Y_e_lat)
    term_zkz_n_Y = (z_e_obs .* (d.n .- Y_e_obs)) .+ (z_e_lat .* k_e .* (d.n .- Y_e_lat))

    p_new = sum(term_zkz_Y)[1] / (sum(term_zkz_Y)[1] + sum(term_zkz_n_Y)[1])

    # Deal with zero mass
    if ((1-p_new)^d.n > 0.999999) || (isnan((1-p_new)^d.n))
        p_new = d.p
    end

    return BinomialExpert(d.n, p_new)

end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::BinomialExpert,
                    ye, exposure,
                    z_e_obs;
                    penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Update parameters
    term_zkz_Y = z_e_obs .* ye
    term_zkz_n_Y = z_e_obs .* (d.n .- ye)

    p_new = sum(term_zkz_Y)[1] / (sum(term_zkz_Y)[1] + sum(term_zkz_n_Y)[1])

    # Deal with zero mass
    if (1-p_new)^d.n > 0.999999 || (isnan((1-p_new)^d.n))
        p_new = d.p
    end

    return BinomialExpert(d.n, p_new)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
#### Custom Filter
using ..Mesh.Filters: AbstractFilter
import ..Mesh.Filters: apply!
export AbstractCustomFilter

"""
    AbstractCustomFilter <: AbstractFilter

Dispatch type for a pointwise custom filter; see [`custom_filter!`](@ref)

!!! warning
    Modifying the prognostic state with this filter
    does not guarantee conservation of the modified
    state.
"""
abstract type AbstractCustomFilter <: AbstractFilter end

"""
    custom_filter!(
        ::AbstractCustomFilter,
        balance_law,
        state,
        aux,
    )

Apply the custom filter `AbstractCustomFilter`,
to the prognostic state for a given balance law.
"""
function custom_filter!(::AbstractCustomFilter, balance_law, state, aux) end

function apply!(
    custom_filter::AbstractCustomFilter,
    grid::DiscontinuousSpectralElementGrid,
    m::BalanceLaw,
    state_prognostic::MPIStateArray,
    state_auxiliary::MPIStateArray,
) where {F!}
    elems = grid.topology.realelems
    device = array_device(state_prognostic)
    Np = dofs_per_element(grid)
    N = polynomialorders(grid)
    knl_custom_filter! = kernel_custom_filter!(device, min(Np, 1024))
    event = Event(device)
    event = knl_custom_filter!(
        m,
        Val(dimensionality(grid)),
        Val(N),
        custom_filter,
        state_prognostic.data,
        state_auxiliary.data,
        elems,
        grid.activedofs;
        ndrange = Np * length(elems),
        dependencies = (event,),
    )
    wait(device, event)
end

@kernel function kernel_custom_filter!(
    bl::BalanceLaw,
    ::Val{dim},
    ::Val{N},
    custom_filter,
    state_prognostic,
    state_auxiliary,
    elems,
    activedofs,
) where {dim, N}
    @uniform begin
        FT = eltype(state_prognostic)
        num_state_prognostic = number_states(bl, Prognostic())
        num_state_auxiliary = number_states(bl, Auxiliary())
        vs_p = Vars{vars_state(bl, Prognostic(), FT)}
        vs_a = Vars{vars_state(bl, Auxiliary(), FT)}
        Nq = N .+ 1

        @inbounds Nqk = dim == 2 ? 1 : Nq[dim]
        @inbounds Np = Nq[1] * Nq[2] * Nqk

        local_state_prognostic = MArray{Tuple{num_state_prognostic}, FT}(undef)
        local_state_auxiliary = MArray{Tuple{num_state_auxiliary}, FT}(undef)
    end

    I = @index(Global, Linear)
    eI = (I - 1) ÷ Np + 1
    n = (I - 1) % Np + 1

    @inbounds begin
        e = elems[eI]

        active = activedofs[n + (e - 1) * Np]

        if active
            @unroll for s in 1:num_state_prognostic
                local_state_prognostic[s] = state_prognostic[n, s, e]
            end

            @unroll for s in 1:num_state_auxiliary
                local_state_auxiliary[s] = state_auxiliary[n, s, e]
            end

            custom_filter!(
                custom_filter,
                bl,
                vs_p(local_state_prognostic),
                vs_a(local_state_auxiliary),
            )

            @unroll for s in 1:num_state_prognostic
                state_prognostic[n, s, e] = local_state_prognostic[s]
            end
        end
    end
end