text:
__precompile__()

module DocOpt

export docopt

using Printf, Dates

# port of str.partition in Python
function partition(str::AbstractString, delim::Union{AbstractString, AbstractChar})
    indices = findfirst(delim, str)
    if isnothing(indices)
        empty_substring = SubString("")
        return SubString(str), empty_substring, empty_substring
    else
        return SubString(str, 1, prevind(str, first(indices))),
            SubString(delim isa AbstractChar ? string(delim) : delim),
            SubString(str, nextind(str, last(indices)))
    end
end

struct DocOptLanguageError <: Exception
    msg::AbstractString
end

struct DocOptExit <: Exception
    usage::AbstractString
end

abstract type Pattern end
abstract type LeafPattern <: Pattern end
abstract type BranchPattern <: Pattern end

mutable struct Argument <: LeafPattern
    name
    value
    Argument(name, value=nothing) = new(name, value)
end

mutable struct Command <: LeafPattern
    name
    value
    Command(name, value=false) = new(name, value)
end

mutable struct Option <: LeafPattern
    short
    long
    argcount::Int
    value

    function Option(short=nothing, long=nothing, argcount=0, value=false)
        value = value === false && argcount > 0 ? nothing : value
        new(short, long, argcount, value)
    end

    function Option(option_description::AbstractString)
        short, long, argcount, value = nothing, nothing, 0, false
        options, _, description = partition(strip(option_description), "  ")
        options = replace(options, ',' => ' ')
        options = replace(options, '=' => ' ')
        for s in split(options)
            if startswith(s, "--")
                long = s
            elseif startswith(s, '-')
                short = s
            else
                argcount = 1
            end
        end
        if argcount > 0
            matched = match(r"\[default: (.*)\]"i, description)
            value = matched == nothing ? nothing : matched.captures[1]
        end
        new(short, long, argcount, value)
    end
end

const Children = Vector{Pattern}

mutable struct Required <: BranchPattern
    children::Children
end

mutable struct Optional <: BranchPattern
    children::Children
end

mutable struct OptionsShortcut <: BranchPattern
    children::Children
    OptionsShortcut() = new(Array[])
end

mutable struct OneOrMore <: BranchPattern
    children::Children
end

mutable struct Either <: BranchPattern
    children::Children
end

mutable struct Tokens
    tokens::Vector{String}
    error::DataType
    function Tokens(source::Array, error=DocOptExit)
        new(source, error)
    end
    function Tokens(source::AbstractString, error=DocOptLanguageError)
        source = replace(source, r"([\[\]\(\)\|]|\.\.\.)" => s -> " " * s * " ")
        source = collect((m.match for m in eachmatch(r"\S*<.*?>|\S+", source)))
        new(source, error)
    end
end

function IteratorSize(::Type{Tokens})
    return Base.SizeUnknown()
end

name(pattern::LeafPattern) = pattern.name
name(o::Option) = o.long !== nothing ? o.long : o.short

function single_match(pattern::Command, left)
    for (n, pat) in enumerate(left)
        if isa(pat, Argument)
            if pat.value == name(pattern)
                return n, Command(name(pattern), true)
            else
                break
            end
        end
    end
    return nothing, nothing
end

function single_match(pattern::Argument, left)
    for (n, pat) in enumerate(left)
        if isa(pat, Argument)
            return n, Argument(name(pattern), pat.value)
        end
    end
    return nothing, nothing
end

function single_match(pattern::Option, left)
    for (n, pat) in enumerate(left)
        if name(pattern) == name(pat)
            return n, pat
        end
    end
    return nothing, nothing
end

Base.:(==)(x::Argument, y::Argument) = x.name == y.name && x.value == y.value
Base.:(==)(x::Command, y::Command) = x.name == y.name && x.value == y.value
Base.:(==)(x::Option, y::Option) = x.short == y.short && x.long == y.long && x.argcount == y.argcount && x.value == y.value
Base.:(==)(x::BranchPattern, y::BranchPattern) = x.children == y.children

function patternmatch(pattern::LeafPattern, left, collected=Pattern[])
    pos, match = single_match(pattern, left)
    if match === nothing
        return false, left, collected
    end
    # drop the pos-th match
    left_ = vcat(left[1:pos-1], left[pos+1:end])
    samename = filter(a -> name(a) == name(pattern), collected)
    if isa(pattern.value, Int) || isa(pattern.value, Array)
        if isa(pattern.value, Int)
            increment = 1
        else
            increment = isa(match.value, AbstractString) ? [match.value] : match.value
        end
        if isempty(samename)
            match.value = increment
            return true, left_, vcat(collected, [match])
        end
        if isa(samename[1].value, Int)
            samename[1].value += increment
        elseif isa(samename[1].value, Array)
            append!(samename[1].value, increment)
        else
            @assert false
        end
        return true, left_, collected
    end
    return true, left_, vcat(collected, [match])
end

function patternmatch(pattern::Union{Optional,OptionsShortcut}, left, collected=Pattern[])
    for pat in pattern.children
        m, left, collected = patternmatch(pat, left, collected)
    end
    return true, left, collected
end

function patternmatch(pattern::Required, left, collected=Pattern[])
    l = left
    c = collected
    for pat in pattern.children
        matched, l, c = patternmatch(pat, l, c)
        if !matched
            return false, left, collected
        end
    end
    return true, l, c
end

function patternmatch(pattern::Either, left, collected=Pattern[])
    outcomes = Any[]
    for pat in pattern.children
        matched, _, _ = outcome = patternmatch(pat, left, collected)
        if matched
            push!(outcomes, outcome)
        end
    end
    if !isempty(outcomes)
        m = first(outcomes)
        for outcome in outcomes
            if length(outcome[2]) < length(m[2])
                m = outcome
            end
        end
        return m
    end
    return false, left, collected
end

function patternmatch(pattern::OneOrMore, left, collected=Pattern[])
    @assert length(pattern.children) == 1
    l = left
    c = collected
    l_ = nothing
    matched = true
    times = 0
    while matched
        matched, l, c = patternmatch(pattern.children[1], l, c)
        times += matched ? 1 : 0
        if l_ == l
            break
        end
        l_ = l
    end
    if times >= 1
        return true, l, c
    end
    return false, left, collected
end

function flat(pattern::LeafPattern, types=[])
    if isempty(types) || (typeof(pattern) in types)
        return [pattern]
    else
        return Pattern[]
    end
end

function flat(pattern::BranchPattern, types=[])
    if typeof(pattern) in types
        return [pattern]
    else
        return reduce(vcat, [flat(child, types) for child in pattern.children], init=Pattern[])
    end
end

function fix(pattern::Pattern)
    fix_identities(pattern)
    fix_repeating_arguments(pattern)
end

function fix_identities(pattern::Pattern, uniq=nothing)
    if !isa(pattern, BranchPattern)
        return pattern
    end
    uniq = uniq === nothing ? unique(flat(pattern)) : uniq
    for (i, child) in enumerate(pattern.children)
        if !isa(child, BranchPattern)
            pattern.children[i] = uniq[something(findfirst(isequal(child), uniq), 0)]
        else
            fix_identities(child, uniq)
        end
    end
end

function fix_repeating_arguments(pattern::Pattern)
    either = [child.children for child in transform(pattern).children]
    for case in either
        for el in filter(child -> count(c -> c == child, case) > 1, case)
            if isa(el, Argument) || isa(el, Option) && el.argcount > 0
                if el.value === nothing
                    el.value = []
                elseif !isa(el.value, Array)
                    el.value = split(el.value)
                end
            end
            if isa(el, Command) || isa(el, Option) && el.argcount == 0
                el.value = 0
            end
        end
    end
    return pattern
end

function transform(pattern::Pattern)
    result = Any[]
    groups = Any[Pattern[pattern]]
    while !isempty(groups)
        children = popfirst!(groups)
        parents = [Required, Optional, OptionsShortcut, Either, OneOrMore]
        if any(map(t -> t in map(typeof, children), parents))
            child = first(filter(c -> typeof(c) in parents, children))
            splice!(children, something(findfirst(isequal(child), children), 0))
            if isa(child, Either)
                for c in child.children
                    push!(groups, vcat([c], children))
                end
            elseif isa(child, OneOrMore)
                push!(groups, vcat(child.children, child.children, children))
            else
                push!(groups, vcat(child.children, children))
            end
        else
            push!(result, children)
        end
    end
    return Either([Required(r) for r in result])
end

Base.hash(pattern::Pattern) = pattern |> string |> hash

Base.getindex(tokens::Tokens, i::Integer) = tokens.tokens[i]
Base.iterate(tokens::Tokens) = isempty(tokens.tokens) ? nothing : (tokens.tokens[1], 2)
Base.iterate(tokens::Tokens, i::Int) = (i > lastindex(tokens.tokens)) ? nothing : (tokens.tokens[i], i+1)
Base.length(tokens::Tokens) = length(tokens.tokens)

move!(tokens::Tokens)   = isempty(tokens.tokens) ? nothing : popfirst!(tokens.tokens)
current(tokens::Tokens) = isempty(tokens.tokens) ? nothing : tokens[1]

# parsers
function parse_long(tokens::Tokens, options)
    # long ::= '--' chars [ ( ' ' | '=' ) chars ] ;
    long, eq, value = partition(move!(tokens), '=')
    @assert startswith(long, "--")
    value = eq == value == "" ? nothing : value
    similar = filter(o -> o.long == long, options)
    if tokens.error === DocOptExit && isempty(similar)  # if no exact match
        similar = filter(o -> o.long !== nothing && startswith(o.long, long), options)
    end
    if length(similar) > 1  # might be simply specified ambiguously 2+ times?
        throw(tokens.error("$long is not a unique prefix: $(join(map(s -> s.long, similar), ","))"))
    elseif length(similar) < 1
        argcount = eq == "=" ? 1 : 0
        o = Option(nothing, long, argcount, false)
        push!(options, o)
        if tokens.error === DocOptExit
            o = Option(nothing, long, argcount, argcount > 0 ? value : true)
        end
    else
        o = Option(similar[1].short, similar[1].long,
                   similar[1].argcount, similar[1].value)
        if o.argcount == 0
            if value !== nothing
                throw(tokens.error("$long must not have an argument"))
            end
        else
            if value === nothing
                if current(tokens) in [nothing, "--"]
                    throw(tokens.error("$long requires argument"))
                end
                value = move!(tokens)
            end
        end
        if tokens.error === DocOptExit
            o.value = value !== nothing ? value : true
        end
    end
    return [o]
end

function parse_shorts(tokens, options)
    # shorts ::= '-' ( chars )* [ [ ' ' ] chars ] ;
    token = move!(tokens)
    @assert startswith(token, '-') && !startswith(token, "--")
    left = lstrip(token, '-')
    parsed = Option[]
    while !isempty(left)
        short = string('-', left[1])
        left = left[2:end]
        similar = filter(o -> o.short == short, options)
        if length(similar) > 1
            throw(tokens.error("$short is specified ambiguously $(length(similar)) times"))
        elseif length(similar) < 1
            o = Option(short, nothing, 0, false)
            push!(options, o)
            if tokens.error === DocOptExit
                o = Option(short, nothing, 0, true)
            end
        else
            o = Option(short, similar[1].long,
                       similar[1].argcount, similar[1].value)
            value = nothing
            if o.argcount != 0
                if isempty(left)
                    if current(tokens) in [nothing, "--"]
                        throw(tokens.error("$short requires argument"))
                    end
                    value = move!(tokens)
                else
                    value = left
                    left = ""
                end
            end
            if tokens.error === DocOptExit
                o.value = value !== nothing ? value : true
            end
       end
        push!(parsed, o)
    end
    return parsed
end

function parse_expr(tokens, options)
    # expr ::= seq ( '|' seq )* ;
    seq = parse_seq(tokens, options)
    if current(tokens) != "|"
        return seq
    end
    result = length(seq) > 1 ? Pattern[Required(seq)] : seq
    while current(tokens) == "|"
        move!(tokens)
        seq = parse_seq(tokens, options)
        append!(result, length(seq) > 1 ? [Required(seq)] : seq)
    end
    return length(result) > 1 ? [Either(collect(result))] : result
end

function parse_seq(tokens, options)
    # seq ::= ( atom [ '...' ] )* ;
    result = Pattern[]
    while !(current(tokens) in [nothing, "]", ")", "|"])
        atom = parse_atom(tokens, options)
        if current(tokens) == "..."
            atom = [OneOrMore(atom)]
            move!(tokens)
        end
        append!(result, atom)
    end
    return result
end

isdash(token) = token == "-" || token == "--"

function parse_atom(tokens, options)
    # atom ::= '(' expr ')' | '[' expr ']' | 'options'
    #        | long | shorts | argument | command ;
    token = current(tokens)
    result = Pattern[]
    closing = Dict("(" => (")", Required), "[" => ("]", Optional))
    if token == "(" || token == "["
        move!(tokens)  # discard '(' or '[' token
        matching, pattern = closing[token]
        result = pattern(parse_expr(tokens, options))
        if move!(tokens) != matching
            throw(tokens.error("unmatched '$token'"))
        end
        return [result]
    elseif token == "options"
        move!(tokens)
        return [OptionsShortcut()]
    elseif startswith(token, "--") && !isdash(token)
        return parse_long(tokens, options)
    elseif startswith(token, '-') && !isdash(token)
        return parse_shorts(tokens, options)
    elseif startswith(token, '<') && endswith(token, '>') || all(isuppercase, token)
        return [Argument(move!(tokens))]
    else
        return [Command(move!(tokens))]
    end
end

function parse_argv(tokens::Tokens, options, options_first=false)
    parsed = Pattern[]
    while !isempty(tokens)
        if current(tokens) == "--"
            return append!(parsed, map(v -> Argument(nothing, v), tokens))
        elseif startswith(current(tokens), "--")
            append!(parsed, parse_long(tokens, options))
        elseif startswith(current(tokens), "-") && !isdash(current(tokens))
            append!(parsed, parse_shorts(tokens, options))
        elseif options_first
            return append!(parsed, map(v -> Argument(nothing, v), tokens))
        else
            push!(parsed, Argument(nothing, move!(tokens)))
        end
    end
    return parsed
end

function parse_pattern(source, options)
    tokens = Tokens(source)
    result = parse_expr(tokens, options)
    return Required(result)
end

function parse_section(name, source)
    pattern = Regex("^([^\\n]*$name[^\\n]*\\n?(?:[ \\t].*?(?:\\n|\$))*)", "im")
    map(strip, collect((m.match for m in eachmatch(pattern, source))))
end

function parse_defaults(doc)
    defaults = Option[]
    for s in parse_section("options:", doc)
        _, _, s = partition(s, ':')
        sp = split(s, "\n")
        sp = map(strip, sp)
        sp = filter!(s -> startswith(s, '-'), sp)
        options = map(Option, sp)
        append!(defaults, options)
    end
    return defaults
end

function formal_usage(section)
    _, _, section = partition(section, ':')
    words = split(strip(section))
    program = popfirst!(words)
    patterns = AbstractString[]
    for w in words
        if w == program
            push!(patterns, ") | (")
        else
            push!(patterns, w)
        end
    end
    return string("( ", join(patterns, ' '), " )")
end

function extras(help, version, options, doc)
    if help && any(o -> name(o) in ["-h", "--help"] && o.value, options)
        println(rstrip(doc, '\n'))
        isinteractive() || exit(0)
    end
    if version !== nothing && any(o -> name(o) == "--version" && o.value , options)
        println(version)
        isinteractive() || exit(0)
    end
end

"""
    docopt(doc::AbstractString,
           args=ARGS;
           version=nothing,
           help::Bool=true,
           options_first::Bool=false,
           exit_on_error::Bool=true)

Parse command-line arguments according to a help message.

Parsed command-line arguments are retuned as a dictionary of arguments of type
`Dict{AbstractString,Any}`; keys are argument names or flag names, and values
are argument values passed to the command-line arguments.

See http://docopt.org/ for the description language of help.

# Arguments
* `doc`: description of your command-line interface.
* `args=ARGS`: argument vector to be parsed.
* `version=nothing`: version of your command-line tool (e.g. `v"1.0.2"`).
* `help=true`: show the help when '-h' or '--help' is passed.
* `options_first=false`: force options to precede positional arguments.
* `exit_on_error=true`: print the usage and exit when parsing error happens.
"""
function docopt(doc::AbstractString,
                args=ARGS;
                version=nothing,
                help::Bool=true,
                options_first::Bool=false,
                exit_on_error::Bool=true)
    usage_sections = parse_section("usage:", doc)
    if isempty(usage_sections)
        throw(DocOptLanguageError("\"usage:\" (case-insensitive) not found."))
    elseif length(usage_sections) > 1
        throw(DocOptLanguageError("More than one \"usage:\" (case-insensitive)."))
    end
    docoptexit = DocOptExit(usage_sections[1])
    options = parse_defaults(doc)
    pattern = parse_pattern(formal_usage(docoptexit.usage), options)
    args = parse_argv(Tokens(args, DocOptExit), options, options_first)
    pattern_options = Set(flat(pattern, [Option]))
    for options_shortcut in flat(pattern, [OptionsShortcut])
        doc_options = parse_defaults(doc)
        options_shortcut.children = Pattern[x for x in setdiff(Set(doc_options), pattern_options)]
    end
    extras(help, version, args, doc)
    matched, left, collected = patternmatch(fix(pattern), args)
    if matched && isempty(left)
        ret = Dict{String,Any}()
        for a in vcat(flat(pattern), collected)
            ret[name(a)] = a.value
        end
        return ret
    end
    if exit_on_error
        @printf(stderr, "%s\n", docoptexit.usage)
        isinteractive() || exit(1)
    else
        throw(docoptexit)
    end
end

end  # DocOpt


text:
using Flux
using Flux: mse, back!
using Flux.Optimise: _update_params!
using Statistics
using Formatting
using Random

"""
DeepQAgent

    - `Q`: the action-value function, i.e., Q(s) -> Q(s, a) for all actions a ∈ A.
    - `memory`: buffer storing transitions `(s, a, r, s′)`.
    - `γ`: reward discount rate.

    # Example
    ```julia-repl
    Q = Chain(
        Conv((3, 3), 1=>32, relu),
        Conv((3, 3), 32=>32, relu),
        x -> maxpool(x, (2,2)),
        Conv((3, 3), 32=>16, relu),
        x -> maxpool(x, (2,2)),
        Conv((3, 3), 16=>10, relu),
        x -> reshape(x, :, size(x, 4)),
        Dense(90, 10),
        softmax)
    loss(x, y) = mse(x, y)
    opt = ADAM()
    agent = DeepQAgent(Q, loss, opt, 1000, 10, 0.01)
    ```
"""
mutable struct DeepQAgent
    Q  # action-value
    Q_target
    loss
    opt
    memory
    max_memory
    batch_size
    γ  # discount rate
    update_episodes
end

function DeepQAgent(Q, lr)
    Q_target = deepcopy(Q)
    loss = mse
    opt = RMSProp(lr)
    memory = []
    max_memory = 1000
    batch_size = 32
    γ = 1.0
    update_episodes = 5
    return DeepQAgent(Q,
                      Q_target,
                      loss,
                      opt,
                      memory,
                      max_memory,
                      batch_size,
                      γ,
                      update_episodes)
end
DeepQAgent(Q, loss, opt, mm) = DeepQAgent(Q, deepcopy(Q), loss, opt, [], mm, 10, 0.0, 1000)
DeepQAgent(Q, loss, opt, mm, bs) = DeepQAgent(Q, deepcopy(Q), loss, opt, [], mm, bs, 0.0, 1000)
DeepQAgent(Q, loss, opt, mm, bs, γ, us) = DeepQAgent(Q, deepcopy(Q), loss, opt, [], mm, bs, γ, us)

function action(agent::DeepQAgent, env::AbstractEnvironment, ϵ = 0.0)
    if rand() < ϵ
        return rand(env.action_space)  # explore
    else
        return argmax(agent.Q(env.state))  # exploit
    end
end

function update!(agent::DeepQAgent)
    mem_batch = rand(agent.memory, agent.batch_size)
    p = map(sars -> estimate(sars, agent), mem_batch)
    y = map(sars -> target(sars, agent), mem_batch)
    l = agent.loss(p, y)
    back!(l)
    _update_params!(agent.opt, Flux.params(agent.Q))
end

function update_target!(agent::DeepQAgent)
    agent.Q_target = deepcopy(agent.Q)
end

function estimate(sars, agent)
    s, a, _, _, done = sars
    q = agent.Q(s)
    return X = q[a, :][1]
end

function target(sars, agent)
    _, _, r, s′, done = sars
    if done
        return r
    else
        return r + agent.γ * maximum(agent.Q_target(s′).data)
    end
end

# function update!(agent::DeepQAgent)
#     mem_batch = rand(agent.memory, agent.batch_size)
#     S = hcat([m[1] for m in mem_batch]...)
#     A = hcat([m[2] for m in mem_batch]...)
#     R = hcat([m[3] for m in mem_batch]...)
#     S_ = hcat([m[4] for m in mem_batch]...)
#     done = hcat([m[5] for m in mem_batch]...)
#     target = R .+ agent.γ .* .!done .* maximum(agent.Q_target(S_).data, dims = 1)
#     estimate = agent.Q(S)[CartesianIndex.(A', axes(S, 2))]
#     loss = agent.loss(estimate, target)
#     back!(loss)
#     _update_params!(agent.opt, Flux.params(agent.Q))
# end

function init_memory!(agent::DeepQAgent, env::AbstractEnvironment)
    @info "Initializing agent memory..."
    if length(agent.memory) > 0
        agent.memory = []
        @info "> cleared existing memory"
    end
    while length(agent.memory) < agent.max_memory
        s′, done = reset!(env)
        while !done
            a = action(agent, env, 1.0)  # random action
            s = deepcopy(s′)
            s′, r, done, info = step!(env, a)
            length(agent.memory) == agent.max_memory && return
            pushfirst!(agent.memory, [s, a, r, s′, done])
        end
    end
end

"""
run_episode!(agent, env; <kwargs>)

    Train `agent` for one episode in `env` following ϵ-greedy policy.

    # Arguments
    - `ϵ`: probability of performing a non-greedy action.
"""
function run_episode!(agent::DeepQAgent, env::AbstractEnvironment; ϵ = 0.05, rendered = false, verbose = false)
    total_reward = 0.0
    steps = 0
    s′, done = reset!(env)
    env.rendered = rendered
    while !done
        a = action(agent, env, ϵ)
        s = deepcopy(s′)
        s′, r, done, info = step!(env, a)
        verbose && @info "s = $s, a = $a, r = $r, s′ = $s′, done = $done"
        length(agent.memory) == agent.max_memory && pop!(agent.memory)
        pushfirst!(agent.memory, [s, a, r, s′, done])
        total_reward += r
        steps += 1
        update!(agent)
    end
    env.rendered = false
    verbose && @info "episode finished in $steps timesteps"
    return total_reward, steps
end

"""
run_episode(agent, env; <kwargs>)

    Evaluate `agent` for one epsiode in `env` following ϵ-greedy policy.

    **NOTE** This method does **not** update the agent.

    # Arguments
    - `ϵ`: probability of performing a non-greedy action.
"""
function run_episode(agent::DeepQAgent, env::AbstractEnvironment; ϵ = 0.05, rendered = false, verbose = false)
    total_reward = 0.0
    steps = 0
    s′, done = reset!(env)
    env.rendered = rendered
    while !done
        a = action(agent, env, ϵ)
        s = deepcopy(s′)
        s′, r, done, info = step!(env, a)
        total_reward += r
        steps += 1
        verbose && @info "s = $s, a = $a, r = $r, s′ = $s′, done = $done"
        if steps == env.maxsteps
            @info "> agent reached end of episode!"
            break
        end
    end
    env.rendered = false
    verbose && @info "episode finished in $steps timesteps"
    return total_reward, steps
end

function evaluate(agent::DeepQAgent, env::AbstractEnvironment, n = 100; ϵ = 0.05, rendered = false, verbose = false)
    history = []
    @progress for i = 1:n
        total_reward, steps = run_episode(agent, env, ϵ = ϵ, rendered = rendered, verbose = verbose)
        push!(history, (total_reward = total_reward, steps = steps))
    end
    total_reward = mean([e.total_reward for e in history])
    steps = mean([e.steps for e in history])
    return total_reward, steps
end

"""
train!(agent, env; <kwargs>)

    Train `agent` on `env`.

    # Arguments
    - `max_episodes`: maximum number of episodes to train the agent.
    - `eval_freq`: number of episodes between fitness evaluations.
    - `ϵ_schedule`: determines ϵ-annealing based on number of epsiodes.

    # Example
    ```julia-repl
    ϵ_schedule = StepDecay(1.0, 0.05, 0.1, 50)
    agent = DeepQAgent(...)
    env = CartPole(...)
    train!(agent, env, eval_freq = 20, ϵ_schedule = ϵ_schedule)
    ```
"""
function train!(agent::DeepQAgent, env::AbstractEnvironment; max_episodes = Inf, eval_freq = 20, eval_episodes = 100, eval_ϵ = 0.05, ϵ_schedule = nothing, verbose = false, patience = 10)

    init_memory!(agent, env)
    @info "training agent..."
    episodes = 0
    recorder = ExperimentRecorder("./log-$(randstring(5))")
    current_best_score = -Inf
    n_evals_since_best_score = 0
    improving = true
    ϵ = ϵ_schedule(0)
    while improving && episodes < max_episodes
        # training episode
        total_reward, steps = run_episode!(agent, env, ϵ = ϵ)
        episodes += 1
        ϵ = ϵ_schedule(episodes)
        verbose && @info "<training status message>"
        message = (episode=episodes, score=total_reward, steps=steps, epsilon=ϵ)
        stamp!(recorder, message, :train)

        agent.opt.eta = 0.99 * agent.opt.eta

        # *target* network update
        verbose && @info "updating target network..."
        if episodes % agent.update_episodes == 0
            update_target!(agent)
        end

        # fitness evaluation
        if episodes % eval_freq == 0
            score, steps = evaluate(agent, env, eval_episodes, ϵ = eval_ϵ)
            @info "episode: $episodes, score: $score, lr: $(agent.opt.eta), ϵ: $ϵ"
            message = (episode=episodes, score=score, steps=steps, epsilon=eval_ϵ)
            stamp!(recorder, message, :test)
            # evaluate(agent, env, 1, ϵ = eval_ϵ, rendered = true)

            # snapshots
            if score > current_best_score
                @info "> new best score achieved!"
                current_best_score = score
                n_evals_since_best_score = 0
                # save a snapshot of agent...
            else
                n_evals_since_best_score += 1
                n_evals_since_best_score == patience && break
            end
        end
    end
    verbose && @info "agent finished training in $episodes episodes"
    return recorder
end

ext: jl


text:
module LaplaceApproximations

using Zygote
using LinearAlgebra

export
    laplace_approx,
    GPML

include("GPML.jl")

"""
    mnll, u = laplace_approx(nll, u0; grtol = sqrt(eps(eltype(u0))))

NOT ROBUST! Uses standard Newton updates with no checking or scaling.

For negative log-likelihood function `nll`, calculate the approximate marginal
likelihood after integrating out the parameters in `u` using the Laplace
approximation. Returns the marginal negative log-likelihood and the values of
`u` at the minimum. The function `nll` should accept only the argument `u`.
Newton steps are taken until the maximum gradient component is less than `grtol`
or `maxit`, whichever comes first.
"""
function laplace_approx(nll::Function, u0;
                        grtol = sqrt(eps(eltype(u0))),
                        maxit = 100)
    u = copy(u0)
    gr = Zygote.gradient(nll, u)[1]
    he = Zygote.hessian(nll, u)
    nit = 0
    while(true)
        nit += 1
        u .-= he \ gr
        gr .= Zygote.gradient(nll, u)[1]
        he .= Zygote.hessian(nll, u)
        maximum(gr) < grtol && break
        if nit ≥ maxit
            @warn "Gradient tolerance not reached in $nit Newton steps" gr = gr
            break
        end
    end
    mnll = -length(u) / 2 * log(2π) + logdet(he) / 2 + nll(u)
    return (mnll, u)
end

end


text:
# Tracking object in real time
# similar: http://pointclouds.org/documentation/tutorials/tracking.php

using PCLCommon
using PCLIO
using PCLVisualization
using PCLFilters
using PCLTracking
using Libfreenect2
using Cxx

RT = PointXYZRGB
PT = ParticleXYZRPY

# Tracking target
target_cloud = PointCloud{RT}(Pkg.dir("PCL", "test", "data",
    "kumamon_head.pcd"))

downsampling_grid_size =  0.010;

function filterPassThrough{T}(cloud_in::PointCloud{T},
        cloud_out::PointCloud{T})
    pass = PassThrough{T}()
    setFilterFieldName(pass, "z")
    setFilterLimits(pass, 0.0, 1.4)
    setKeepOrganized(pass, false)
    setInputCloud(pass, cloud_in)
    filter(pass, cloud_out)
end

function gridSampleApprox{T}(cloud_in::PointCloud{T},
        cloud_out::PointCloud{T}, leaf_size)
    grid = ApproximateVoxelGrid{T}()
    setLeafSize(grid, leaf_size, leaf_size, leaf_size)
    setInputCloud(grid, cloud_in)
    filter(grid, cloud_out)
end

default_step_covariance = icxx"std::vector<double>(6, 0.015 * 0.015);"
icxx"$(default_step_covariance)[3] *= 40.0;"
icxx"$(default_step_covariance)[4] *= 40.0;"
icxx"$(default_step_covariance)[5] *= 40.0;"

initial_noise_covariance = icxx"std::vector<double>(6, 0.00001);"
default_initial_mean = icxx"std::vector<double>(6, 0.0);"

# Initialize tracker
global tracker = KLDAdaptiveParticleFilterOMPTracker{RT,PT}(8)

bin_size = PT()
icxx"$bin_size.x = 0.1f;"
icxx"$bin_size.y = 0.1f;"
icxx"$bin_size.z = 0.1f;"
icxx"$bin_size.roll = 0.1f;"
icxx"$bin_size.pitch = 0.1f;"
icxx"$bin_size.yaw = 0.1f;"

setMaximumParticleNum(tracker, 500)
setDelta(tracker, 0.99)
setEpsilon(tracker, 0.2)
setBinSize(tracker, bin_size)

setTrans(tracker, icxx"Eigen::Affine3f::Identity();")
setStepNoiseCovariance(tracker, default_step_covariance)
setInitialNoiseCovariance(tracker, initial_noise_covariance)
setInitialNoiseMean(tracker, default_initial_mean);
setIterationNum(tracker, 1);
setParticleNum(tracker, 400);
setResampleLikelihoodThr(tracker, 0.00);
setUseNormal(tracker, false);

# Setup coherence object for tracking
coherence = ApproxNearestPairPointCloudCoherence{RT}()

distance_coherence = DistanceCoherence{RT}()
color_coherence = HSVColorCoherence{RT}()
setWeight(color_coherence, 0.1)

addPointCoherence(coherence, distance_coherence)
addPointCoherence(coherence, color_coherence)

search = Octree{RT}(0.01)
setSearchMethod(coherence, search)
setMaximumDistance(coherence, 0.01)

setCloudCoherence(tracker, coherence)

# prepare the model of tracker's target
c = icxx"Eigen::Vector4f();"
trans = icxx"Eigen::Affine3f trans = Eigen::Affine3f::Identity(); return trans;"
transed_ref = PointCloud{RT}()
transed_ref_downsampled = PointCloud{RT}()

compute3DCentroid(target_cloud, c)
vec = icxx"Eigen::Vector3f($(c)[0], $(c)[1], $(c)[2]);"
icxx"$(trans).translation().matrix() = $vec;"
transformPointCloud(target_cloud, transed_ref, icxx"$trans.inverse();")
gridSampleApprox(transed_ref, transed_ref_downsampled, downsampling_grid_size)
@show length(target_cloud)
@show length(transed_ref)
@show length(transed_ref_downsampled)

setReferenceCloud(tracker, transed_ref_downsampled)
setTrans(tracker, trans)

global counter = 0

function update(tracker, cloud)
    cloud_pass = PointCloud{RT}()
    cloud_pass_downsampled = PointCloud{RT}()
    filterPassThrough(cloud, cloud_pass)
    gridSampleApprox(cloud_pass, cloud_pass_downsampled, downsampling_grid_size)
    if counter > 10
        setInputCloud(tracker, cloud_pass_downsampled)
        t = @elapsed compute(tracker)
        println("compute: $(1/t) Hz")
    end
    return cloud_pass_downsampled
end

function draw(tracker, viewer)
    particles = icxx"$(tracker.handle)->getParticles();"
    icxx"$particles == nullptr;" && return

    particle_cloud = PointCloud{PointXYZ}()
    icxx"""
    for (size_t i = 0; i < $particles->points.size(); i++) {
      pcl::PointXYZ point;
      point.x = $particles->points[i].x;
      point.y = $particles->points[i].y;
      point.z = $particles->points[i].z;
      $(particle_cloud.handle)->points.push_back(point);
    }
    """
    red_color = PointCloudColorHandlerCustom{PointXYZ}(particle_cloud,
        255, 0, 0)
    if !updatePointCloud(viewer, particle_cloud, red_color, id="particle cloud")
        addPointCloud(viewer, particle_cloud, red_color, id="particle cloud")
    end

    result = getResult(tracker)
    transformation = toEigenMatrix(tracker, result)
    icxx"$transformation.translation() += Eigen::Vector3f(0.0f, 0.0f, -0.005f);"
    result_cloud = PointCloud{RT}()
    icxx"""
        pcl::transformPointCloud<$RT>(*($(tracker.handle)->getReferenceCloud()),
            *$(result_cloud.handle), $transformation);
    """
    blue_color = PointCloudColorHandlerCustom{RT}(result_cloud, 0, 0, 255)
    if !updatePointCloud(viewer, result_cloud, blue_color, id="result cloud")
        addPointCloud(viewer, result_cloud, blue_color, id="result cloud")
    end
end

const w = 512
const h = 424
save_pcd = false

genfilename(ext=".pcd") =
    joinpath(dirname(@__FILE__), string(now(), "_", time_ns(), ext))

"""Get point cloud from undistored depth and reigstered color"""
function getPointCloudXYZRGB(registration::Registration, undistorted, registered)
    w = width(undistorted)
    h = height(undistorted)
    cloud = PointCloud{RT}(w, h)
    icxx"$(cloud.handle)->is_dense = false;"
    pointsptr = icxx"&$(cloud.handle)->points[0];"
    icxx"""
    for (size_t ri = 0; ri < $h; ++ri) {
        for (size_t ci = 0; ci < $w; ++ci) {
            auto p = $(pointsptr) + $w * ri + ci;
            $(registration.handle)->getPointXYZRGB($(undistorted.handle),
                $(registered.handle), ri, ci, p->x, p->y, p->z, p->rgb);
        }
    }
    """
    cloud
end

f = Freenect2()
device = openDefaultDevice(f, OpenGLPacketPipeline())
listener = SyncMultiFrameListenerPtr()
setIrAndDepthFrameListener(device, listener)
setColorFrameListener(device, listener)

start(device)

# NOTE: must be called after start(device)
registration = Registration(getIrCameraParams(device),
    getColorCameraParams(device))
undistorted = FramePtr(w, h, 4, key=Libfreenect2.FRAME_DEPTH)
registered = FramePtr(w, h, 4, key=Libfreenect2.FRAME_COLOR)

info("Prepare PCL visualizer...")
global viewer = PCLVisualizer("pcl visualizer")

global should_save = false
if !isdefined(:viewer_cb_defined) && VERSION < v"0.5.0-dev+2396"
    cxx"""
    void viewer_cb(const pcl::visualization::KeyboardEvent &event) {
        std::cout << "key event:" << event.getKeyCode() << std::endl;
        if (event.getKeyCode() == 's') {
            $:(global should_save = true; nothing);
        }
    }
    """
const viewer_cb_defined = true
end
if VERSION < v"0.5.0-dev+2396"
    icxx"$(viewer.handle)->registerKeyboardCallback(viewer_cb);"
end

while !wasStopped(viewer)
    frames = waitForNewFrame(listener)
    color = frames[FrameType.COLOR]
    ir = frames[FrameType.IR]
    depth = frames[FrameType.DEPTH]

    # Depth and color registration
    t = @elapsed begin
        apply(registration, color, depth, undistorted, registered)
        cloud = getPointCloudXYZRGB(registration, undistorted, registered)
    end
    println("Registration and getPointCloudXYZRGB: $(1/t) Hz")

    # tracking
    t = @elapsed begin
        cloud_pass_downsampled = update(tracker, cloud)
    end
    println("Update tracker: $(1/t) Hz")

    # draw particles
    draw(tracker, viewer)

    if should_save
        info("save pcd file...")
        savePCDFile(genfilename(), cloud; binary_mode=true)
        should_save = false
    end

    color_handler = PointCloudColorHandlerRGBField(cloud_pass_downsampled)
    updated = updatePointCloud(viewer, cloud_pass_downsampled, color_handler,
        id="libfreenect2")
    if !updated
        addPointCloud(viewer, cloud_pass_downsampled, color_handler, id="libfreenect2")
    end
    spinOnce(viewer, 1)

    release(listener, frames)

    rand() > 0.95 && gc(false)
    global counter += 1
end

stop(device)
close(device)
close(viewer)


text:
s_2,s2_gradient= s2_grad_func(r,b)
s2_big,s2_grad_numeric = s2_grad_num(r,b)
s2_grad_ana = zeros(2)
diff1 = s2_grad_numeric-s2_gradient
s_2=s2!(r,b,s2_grad_ana)
diff2 = s2_grad_ana-s2_gradient
println("b : ",b," r: ",r," flux diff: ",s_2-s2_big," grad diff(num-auto): ",diff1," diff(ana-auto): ",diff2 )
println("gradient: ",s2_gradient," expected: ",-[2,-2./3.])
@test isapprox(s2_grad_numeric[1],s2_gradient[1],atol=1e-8)
@test isapprox(s2_grad_numeric[2],s2_gradient[2],atol=1e-8)
@test isapprox(s2_grad_ana[1],s2_gradient[1],atol=1e-8)
@test isapprox(s2_grad_ana[2],s2_gradient[2],atol=1e-8)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
register(DataDep(
    "pudl-raw-eia860",
    """
    Dataset: PUDL Raw EIA Form 860
    Website: https://doi.org/10.5281/zenodo.4127027
    Author: Gosnell, Christina, Selvans, Zane
    Date of Publication: October 24, 2020
    License:

    Raw US Energy Information Administration (EIA) Form 860 data, archived from https://www.eia.gov/electricity/data/eia860/

    EIA Form 860 data for electric power plants with 1 megawatt or greater combined nameplate capacity.

    $PUDL_DESCRIPTION
    """,
    [
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/datapackage.json",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2001.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2002.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2003.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2004.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2005.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2006.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2007.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2008.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2009.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2010.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2011.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2012.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2013.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2014.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2015.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2016.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2017.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2018.zip",
        "https://zenodo.org/api/files/64dc3ec3-5f0d-446c-966e-9793cb478960/eia860-2019.zip",
    ],
    "12c0d64081802e1df4a4bb8a960004d53dcc4204d043be9b121c9cf936c2ae88",
))


text:
using SthenoCompat
using Test

@testset "SthenoCompat.jl" begin
    # Write your tests here.
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "sgiter" begin

    nep = nep_gallery("real_quadratic")
#    nep = nep_gallery("dep_distributed")
#    nep = nep_gallery("pep0_sym")
    λ,v = sgiter(Float64,
                 nep,
                 1,
                 λ_min = -10,
                 λ_max = 0,
                 λ = -10,
                 logger = displaylevel,
                 maxit = 100, tol=1e-12)
   @test norm(compute_Mlincomb(nep,λ,v)) < 1e-9

   λ,v = sgiter(Float64,
                nep,
                2,
                logger = displaylevel,
                tol = 1e-9,
                maxit = 100,
                errmeasure = ResidualErrmeasure(nep))
  @test norm(compute_Mlincomb(nep,λ,v)) < 1e-9

end


text:
import MittagLefflerFunctions.asymptotic_start
using PyPlot

β = 1.0
m = 201
α = range(0, 1, length=m)
tol = 1e-15
min_x = 35

x = Vector{Float64}(undef, m)
N = Vector{Int64}(undef, m)
for k = 1:m
    x[k], N[k] = asymptotic_start(α[k], β, tol, min_x)
end

figure(1)
subplot(2, 1, 1)
plot(α, x, "o", markersize=2)
ylabel(L"$x$")
grid(true)
subplot(2, 1, 2)
plot(α, N, "o", markersize=2)
xlabel(L"$\alpha$")
ylabel(L"$N$")
grid(true)


text:
using Test, bad

# get curtailed simons design (implied stopping for efficacy!)
function simon(r1::Int, n1::Int, r::Int, n::Int)
    nn2 = zeros(Int, n1 + 1)
    nn2[(r1 + 2):end] .= n - n1
    cc = repeat([Inf], n1 + 1)
    cc[(r1 + 2):end] .= r
    cc2 = cc .- collect(0:n1)
    cc2[cc2 .< 0]     .= -Inf
    cc2[cc2 .>= nn2]  .= Inf
    nn2[cc2 .== -Inf] .= 0
    Design(nn2, cc2)
end

function optimal_design(pnull, palt)
    essnull = SampleSize(p | pnull)
    problem = Problem(
        minimise(essnull),
        (Power(p | pnull) <= α) |>
            cnstr -> conditional(cnstr, (-.1, 1.1) ),
        (Power(p | palt) >= 1 - β) |>
            cnstr -> conditional(cnstr, (-.1, 1.1) ),
        nmax                     = convert(Int, ceil(1.5*maximum(n(simonsdesign)))),
        n1values                 = convert(Vector{Int}, floor(.5*n1(simonsdesign)):ceil(1.5*n1(simonsdesign))),
        # relax all other heuristics
        n2ton1fctrs              = (0.0, 10.0),
        n2mincontinueabs         = 5,
        curtail_stage_one_buffer = 5,
        type                     = :GroupSequential
    )
    design = optimise(problem; verbosity = 0)
    return design, problem, size(problem)
end

α, β  = 0.05, .2
p     = Beta(1, 1) # uniform prior, only used to make conditioning valid for all p ∈ [0, 1]

# group sequential designs are notoriously hard to fit due to the global
# constraint, only consider the two cornering cases (low/high p0)
# and one in the middle, values taken from:
#
# > Simon, R. (1989). Optimal two-stage designs for phase II clinical trials.
# > Controlled clinical trials, 10(1), 1-10.

pnull, palt  = .05, .25; simonsdesign = simon(0, 9, 2, 17)
design, problem, nvar = optimal_design(pnull, palt)
@test all( problem.toer.score.((design, simonsdesign)) .<= α )
@test all( problem.power.score.((design, simonsdesign)) .>= 1 - β )
@test problem.objective(design) <= problem.objective(simonsdesign)

pnull, palt  = .1, .3; simonsdesign = simon(1, 10, 5, 29)
design, problem, nvar = optimal_design(pnull, palt)
@test all( problem.toer.score.((design, simonsdesign)) .<= α )
@test all( problem.power.score.((design, simonsdesign)) .>= 1 - β )
@test problem.objective(design) <= problem.objective(simonsdesign)

pnull, palt  = .7, .9; simonsdesign = simon(4, 6, 22, 27)
design, problem, nvar = optimal_design(pnull, palt)
@test all( problem.toer.score.((design, simonsdesign)) .<= α )
@test all( problem.power.score.((design, simonsdesign)) .>= 1 - β )
@test problem.objective(design) <= problem.objective(simonsdesign)


text:
using Fermionic
using SparseArrays
using Test

@testset "mixed.jl" begin
    @test eigensp_mixed([0.25,0.75],[State([0,0,0,im/sqrt(2),0,0,0,0,0,0,0,0,1/sqrt(2),0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],Op(4))])[1] == 0.875
    @test_throws ArgumentError rhosp_mixed([0.25,0.8],[State([0,0,0,im/sqrt(2),0,0,0,0,0,0,0,0,1/sqrt(2),0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],Op(4))])
    @test_throws ArgumentError rhosp_mixed([0.25,0.75],[State([0,0,0,im/sqrt(2),0,0,0,0,0,0,0,0,1/sqrt(2),0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0],Op(4))])
    @test_throws ArgumentError eigensp_mixed([0.25,0.8],[State([0,0,0,im/sqrt(2),0,0,0,0,0,0,0,0,1/sqrt(2),0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],Op(4))])
    @test_throws ArgumentError eigensp_mixed([0.25,0.75],[State([0,0,0,im/sqrt(2),0,0,0,0,0,0,0,0,1/sqrt(2),0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],Op(4)),State([0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0],Op(4))])
end


text:
#=
   Girl births

   From Gelman et.al "Regression and Other Stories", 69f

   Part 1 (page 69)
   """
   The probability that a baby is a girl or boy is approximately 48.8% or 51.2%, respectively,
   and these of births do not vary much across the world. Suppose that 400 babies are born in a
   hospital in a given year. How many will be girls?
   """

   Part 2 (page 70)
   """
   Accounting for twins
   We can complicate the model in various ways. For example, there is a 1/125 chance that a
   birth event results in fraternal twins, of which each has an approximate 49.5% chance of being a girl,
   and a 1/300 chance of identical twins, which have an approximate 49.5% chance of being a pair of girls.
   """

Model 1:
Chains MCMC chain (10000×3×1 Array{Float64, 3}):

Log evidence      = 0.0
Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 1.79 seconds
Compute duration  = 1.79 seconds
parameters        = girl
internals         = lp, logevidence

Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

        girl   195.2312    9.9564     0.0996    0.1216   6737.7203    0.9999     3759.8886

Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5%
      Symbol    Float64    Float64    Float64    Float64    Float64

        girl   176.0000   189.0000   195.0000   202.0000   215.0000

Credible interval for girl with mass 0.9: (178.000000..210.000000)
Credible interval for girl with mass 0.5: (186.000000..199.000000)

Model 2:
Chains MCMC chain (10000×6×1 Array{Float64, 3}):

Log evidence      = -4.462871026284194e-5
Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 5.06 seconds
Compute duration  = 5.06 seconds
parameters        = p, girl, identical_twins, fraternal_twins
internals         = lp, logevidence

Summary Statistics
       parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
           Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

  identical_twins     0.0027    0.0519     0.0005    0.0006   5793.9766    1.0000     1144.6023
  fraternal_twins     0.0109    0.1038     0.0010    0.0013   5777.4073    1.0000     1141.3290
                p     0.4881    0.0008     0.0000    0.0000   5843.3985    1.0000     1154.3656
             girl   195.6097   10.0797     0.1008    0.1158   6778.0922    0.9999     1339.0146

Quantiles
       parameters       2.5%      25.0%      50.0%      75.0%      97.5%
           Symbol    Float64    Float64    Float64    Float64    Float64

  identical_twins     0.0000     0.0000     0.0000     0.0000     0.0000
  fraternal_twins     0.0000     0.0000     0.0000     0.0000     0.0000
                p     0.4880     0.4880     0.4880     0.4880     0.4880
             girl   176.0000   189.0000   195.0000   202.0000   216.0000

Credible interval for girl with mass 0.9: (178.000000..211.000000)
Credible interval for girl with mass 0.5: (187.000000..200.000000)

=#

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

# Part 1
# """
# The probability that a baby is a girl or boy is approximately 48.8% or 51.2%, respectively,
# and these of births do not vary much across the world. Suppose that 400 babies are born in a
# hospital in a given year. How many will be girls?
# """
@model function girl_births1(n=400)
    girl ~ Binomial(n,0.488)
end

println("Model 1:")
n=400
model = girl_births1(n)

num_chains = 4

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 1_000)

display(chns)
# display(plot(chns))

credible_interval(chns,:girl,0.9)
credible_interval(chns,:girl,0.5)
println()

# Part 2
# """
# We can complicate the model in various ways. For example, there is a 1/125 chance that a
# birth event results in fraternal twins, of which each has an approximate 49.5% chance of being a girl,
# and a 1/300 chance of identical twins, which have an approximate 49.5% chance of being a pair of girls.
# """
@model function girl_births2(n=400)
    identical_twins ~ flip(1/300)
    fraternal_twins ~ flip(1/125)
    # Cannot be both identical_twins and fraternal_twins
    true ~ Dirac(identical_twins + fraternal_twins <= 1)
    # The probability of a girl:
    p ~ Dirac((identical_twins || fraternal_twins) ? 0.495 : 0.488)
    # Number of girls
    girl ~ Binomial(n,p)

end

println("\nModel 2:")
n=400
model = girl_births2(n)

num_chains = 4

# chns = sample(model, Prior(), 100_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 1_000)

display(chns)
# display(plot(chns))

credible_interval(chns,:girl,0.9)
credible_interval(chns,:girl,0.5)
println()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using DifferentialEquations
using Plots
plotlyjs()

#constants
g=9.81       #acceleration due to gravity [ms^-2]
m=3000       #mass of capsule [Kg]

#initial conditions
u₀=[3,0.1,0.1,0.5]
r₀=6
v₀=0

#function to compute lagrangian
function lagrangian!(du,u,p,t)
    θ,ω,ϕ,ψ=u
    v=v₀
    r=r₀+v*t
    du[1]=dθ=ω
    du[2]=dω=sin(θ)*cos(θ)*ψ^2-((2*v*ω)/r)+((g*sin(θ))/r)
    du[3]=dϕ=ψ
    du[4]=dψ=-((2*v*ψ)/r)-(2*cot(θ)*ω*ψ)
end

#time interval
tspan=(0.0,100.0)

#defining the problem and solving
prob=ODEProblem(lagrangian!,u₀,tspan)
sol=solve(prob,saveat=0.1)

#extracting solution
U=sol[1:end,:]
(θ,ω,ϕ,ψ)=[U[x,:] for x in 1:size(U,1)]
t=sol.t
r=r₀ .+ v₀*t
v=v₀*ones(size(t))
x=r.*sin.(θ).*cos.(ϕ)
y=r.*sin.(θ).*sin.(ϕ)
z=r.*cos.(θ)
uₜ=U[:,end]
rₜ=r[end]
vₜ=v[end]
FC=[rₜ,vₜ,uₜ]

#plotting the solution
display(plot(x,y,z,label="trajectory",xaxis="x (m)",
    yaxis="y (m)",zaxis="z (m)",grid=(:on,:black)))
savefig("./plots/trajectory-sp.png")

p1=plot(sol.t,r,label="r",color="black",
    xaxis="t (s)", yaxis="r (m)")
p2=plot(sol.t,v,label="v",color="red",
    xaxis="t (s)", yaxis="v (m/s)")
p3=plot(sol,vars=(0,1),label="θ",
    color="green", xaxis="t (s)", yaxis="θ (rad)")
p4=plot(sol,vars=(0,2),label="ω",
    color="orange", xaxis="t (s)", yaxis="ω (rad/s)")
p5=plot(sol,vars=(0,3),label="ϕ",
    color="purple", xaxis="t (s)", yaxis="ϕ (rad)")
p6=plot(sol,vars=(0,4),label="ψ",
    color="blue", xaxis="t (s)", yaxis="ψ (rad/s)")
display(plot(p1,p2,p3,p4,p5,p6,layout=(3,2)))
savefig("./plots/plots-sp.png")

#energy plot
KE=0.5*m*(v.^2 .+ (r.*ω).^2 .+ (r.*ψ.*sin.(θ)).^2)
PE=m*g.*z
TE=KE.+PE
p7=plot(t,KE,label="KE",
    color="black",xaxis="t (s)", yaxis="KE (J)")
p8=plot(t,PE,label="PE",
    color="green",xaxis="t (s)", yaxis="PE (J)")
p9=plot(t,TE,label="TE",
    color="red",xaxis="t (s)", yaxis="TE (J)")
display(plot!(p7,p8,p9))
savefig("./plots/energyplot-sp.png")

max_stars_repo_path: Modeling the Dragon Module as a Pendulum/Scripts/spherical pendulum.jl


text:
# some common functions
# Query operations are located in the derived.jl file.

abstract type AbstractIndexableBitVector <: AbstractVector{Bool} end
const AbstractBitVector = Union{BitVector,AbstractIndexableBitVector}

Base.size(b::AbstractIndexableBitVector) = (length(b),)


text:
#=
epigraph of perspective of squared Euclidean norm (AKA rotated second-order cone)
(u, v, w) : v ≥ 0, u ≥ v / 2 * ‖w / v‖² = ‖w‖² / 2v
self-dual

extended formulation
∃ λ, u ≥ Σᵢ λᵢ, (λᵢ, v, wᵢ) ∈ EpiPerSquare
i.e. λᵢ ≥ 0, 2 λᵢ v ≥ wᵢ²
=#

mutable struct EpiPerSquare{E <: NatExt} <: Cache
    oa_s::Vector{AE}
    d::Int
    λ::Vector{VR}
    EpiPerSquare{E}() where {E <: NatExt} = new{E}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.EpiPerSquareCone{RealF},
    opt::Optimizer,
)
    dim = MOI.dimension(cone)
    @assert dim == length(oa_s)
    d = dim - 2
    E = nat_or_ext(opt, d)
    cache = EpiPerSquare{E}()
    cache.oa_s = oa_s
    cache.d = d
    return cache
end

function per_square(v::RealF, w::AbstractVector{RealF})
    den = 2 * max(v, 1e-9)
    return sum(w_i / den * w_i for w_i in w)
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::EpiPerSquare,
    opt::Optimizer,
)
    return _get_cuts(z[2], z[3:end], cache, opt)
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::EpiPerSquare,
    opt::Optimizer,
)
    us = s[1]
    vs = s[2]
    @views ws = s[3:end]
    if us - per_square(vs, ws) > -1e-7 # TODO option
        return AE[]
    end

    # gradient cut is (1, ‖w‖² / (2 * vs²), -ws / vs)
    r = ws / -max(vs, 1e-9)
    q = 0.5 * sum(abs2, r)
    return _get_cuts(q, r, cache, opt)
end

# unextended formulation

function MOIPajarito.Cones.add_init_cuts(cache::EpiPerSquare{Nat}, opt::Optimizer)
    # add variable bounds
    u = cache.oa_s[1]
    v = cache.oa_s[2]
    JuMP.@constraint(opt.oa_model, u >= 0)
    JuMP.@constraint(opt.oa_model, v >= 0)
    opt.use_init_fixed_oa || return

    # add cuts (1, 2, ±2eᵢ) on (u, v, wᵢ)
    @views w = cache.oa_s[3:end]
    d = cache.d
    JuMP.@constraints(opt.oa_model, begin
        [i in 1:d], u + 2v + 2w[i] >= 0
        [i in 1:d], u + 2v - 2w[i] >= 0
    end)
    return
end

function _get_cuts(q::RealF, r::Vector{RealF}, cache::EpiPerSquare{Nat}, opt::Optimizer)
    clean_array!(r) && return AE[]
    # strengthened cut is (‖r‖² / 2q, q, r)
    p = per_square(q, r)
    u = cache.oa_s[1]
    v = cache.oa_s[2]
    @views w = cache.oa_s[3:end]
    cut = JuMP.@expression(opt.oa_model, p * u + q * v + JuMP.dot(r, w))
    return [cut]
end

# extended formulation

MOIPajarito.Cones.num_ext_variables(cache::EpiPerSquare{Ext}) = cache.d

function MOIPajarito.Cones.extend_start(
    cache::EpiPerSquare{Ext},
    s_start::Vector{RealF},
    opt::Optimizer,
)
    u_start = s_start[1]
    v_start = s_start[2]
    w_start = s_start[3:end]
    if max(u_start, v_start) < 1e-8
        return zeros(cache.d)
    end
    return [w_i / 2u_start * w_i for w_i in w_start]
end

function MOIPajarito.Cones.setup_auxiliary(cache::EpiPerSquare{Ext}, opt::Optimizer)
    @assert cache.d >= 2
    λ = cache.λ = JuMP.@variable(opt.oa_model, [1:(cache.d)], lower_bound = 0)
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= sum(λ))
    return λ
end

function MOIPajarito.Cones.add_init_cuts(cache::EpiPerSquare{Ext}, opt::Optimizer)
    # add variable bounds
    u = cache.oa_s[1]
    v = cache.oa_s[2]
    JuMP.@constraint(opt.oa_model, u >= 0)
    JuMP.@constraint(opt.oa_model, v >= 0)
    opt.use_init_fixed_oa || return

    # add disaggregated cuts (1, 2, ±2) on (λᵢ, v, wᵢ)
    @views w = cache.oa_s[3:end]
    d = cache.d
    λ = cache.λ
    JuMP.@constraints(opt.oa_model, begin
        [i in 1:d], λ[i] + 2v + 2w[i] >= 0
        [i in 1:d], λ[i] + 2v - 2w[i] >= 0
    end)
    return
end

function _get_cuts(q::RealF, r::Vector{RealF}, cache::EpiPerSquare{Ext}, opt::Optimizer)
    clean_array!(r) && return AE[]
    p = per_square(q, r)
    v = cache.oa_s[2]
    @views w = cache.oa_s[3:end]
    λ = cache.λ
    cuts = AE[]
    for i in 1:(cache.d)
        r_i = r[i]
        iszero(r_i) && continue
        # p = ‖r‖² / 2q, strengthened disaggregated cut on (λᵢ, v, wᵢ) is (p, rᵢ² / 2p, rᵢ)
        # TODO check math
        q_i = r_i / 2p * r_i
        cut = JuMP.@expression(opt.oa_model, p * λ[i] + q_i * v + r_i * w[i])
        push!(cuts, cut)
    end
    return cuts
end


text:
sense_parity_states(n) = RecordArrays.fill((sense = true, parity = 0), n; align = 64)

struct DisseminationBarrier <: SyncBarriers.DisseminationBarrier
    n::Int
    flags::Array{OneWayObservable{Bool},3}
    locals::typeof(sense_parity_states(1))
end

DisseminationBarrier(n::Integer) = DisseminationBarrier(
    n,
    [OneWayObservable{Bool}(false) for _ in 1:ceillog2(n), _ in 1:n, _ in 1:2],
    sense_parity_states(n),
)

function inc_sense_parity!(state)
    parity = state.parity[]
    sense = state.sense[]
    if parity == 1
        state.sense[] = !sense
    end
    state.parity[] = 1 - parity
    return sense, parity
end

function prev_sense_parity(state)
    parity = 1 - state.parity[]
    sense = state.sense[]
    return ((parity == 1 ? !sense : sense), parity)
end

function SyncBarriers.cycle!(
    handle::BarrierHandle{DisseminationBarrier},
    spin::Union{Nothing,Integer} = nothing,
)
    i = handle.i
    barrier = handle.barrier

    sense, parity = inc_sense_parity!(view(barrier.locals, handle.i))

    shift = 1
    for flags in eachrow(@view barrier.flags[:, :, parity+1])
        @_assert flags[i][] != sense
        flags[i][] = sense

        j = i - shift
        if j <= 0
            j += lastindex(flags)
        end
        shift *= 2

        waitif(x -> x != sense, flags[j], spin)
        @_assert flags[j][] == sense
    end
end


text:
"""
This module provides an implementation of NetMSA algorithm in Julia, which
can be used for multiple sequence alignment.
"""
module NetMSA

using StatsBase

export createPeerMatrix, matrixalignment

@doc raw"""
Store the position of a given particle. Position ``x_{s_{i}}(r)`` of the particle
``p_{s_{i}}`` is defined by using the row ``r`` that contains the symbol ``s_i`` as
well as locations of the symbol ``s_i`` in the different columns (indexes of the
columns that contain ``s_i`` ) in the row ``r``.
"""
mutable struct Position
    row::Int64
    indexes::Vector{Int64}
end

"""
A particle that is used for creating swarms.

# Fields
- value::Char : Value of the particle, e.g. 'b' or 'c'
- updated::Int64 : Number of turns till last updated
- pos::Position : The original position of the particle
- best::Position : The best local position of the particle
- bestvalue::Float64 : Best local score
"""
mutable struct Particle
    value::Char
    updated::Int64
    pos::Position
    best::Position
    bestvalue::Float64

    """
        Particle(value::Char, pos::Position)

    Initializes the Particle, where the best position is intialized
    to be the current position.
    """
    function Particle(value::Char, pos::Position)
        return new(value, 0, pos, pos, 0.0)
    end
end

"""
    createPeerMatrix(inputStrings::Vector{String})::Matrix{Union{Missing,Char}}

Create and return a Peer matrix, containing charachters as elements, where each input sequence,
provided in the inputStrings, is represented as a column. Missing values are represented in the
matrix by the `missing` keyword.

# Examples
```jldoctest
julia> NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'
```
"""
function createPeerMatrix(inputStrings::Vector{String})::Matrix{Union{Missing,Char}}
		col_size = size(inputStrings, 1);
		row_size = max([length(s) for s in inputStrings]...);
		peer_matrix = Matrix{Union{Missing,Char}}(missing, row_size, col_size);
		for (i, s) in enumerate(inputStrings)
			peer_matrix[1:length(s), i] = collect(s);
		end
		return peer_matrix;
end

"""
    getposition(value, rowindex, matrix)

Return the Position (rowindex, [colindex1, colindex2, ...]) of the Particle
represented by `value`, at the `rowindex` in the `matrix`.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.getposition('b', 2, M)
 NetMSA.Position(2, [1, 3, 4])
```
"""
function getposition(value, rowindex, matrix)
  indexes = findall(i -> i == value, skipmissing(matrix[rowindex, :]))
  return Position(rowindex, indexes);
end

"""
    mostfrequent(row)

Return a tuple containing the most frequent element occuring in the `row`, along
with its frequency.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.mostfrequent(M[2, :])
 (3, 'b')
```
"""
function mostfrequent(row)
  counts = countmap(row);
  delete!(counts, '-');
  max = isempty(counts) ? 0 : findmax(counts);
  return max;
end

"""
    aligned(row)::Bool

Return whether a row is aligned or not.

A row is *aligned* if it only contains different occurrences of the same symbol.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.aligned(M[1, :])
 true

 juila> NetMSA.aligned(M[2, :])
 false
```
"""
function aligned(row)::Bool
  row = Set(row)
  return (length(row) == 1 && !(missing in row)) || (length(row) == 2 && ('-' in row || missing in row))
end

"""
    aligned(row)::Bool

Return whether a row is full or not.

An aligned row r is called full if no gaps (—) are added in the row r .
That is, the number of occurrences of the symbol in the row is equal to
the number of columns in the matrix.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.full(M[1, :])
 true

 juila> NetMSA.full(M[2, :])
 false
```
"""
function full(row)::Bool
  return length(Set(row)) == 1 && !(missing in Set(row))
end

@doc raw"""
    weight(row; w1=0.25, w2=0.5, w3=1.0)

Return the weight of the row, calculated as:

```math
w(r) = \begin{cases}
  w_1  \times  \frac{x}{c}; & \text{ if r is not aligned}   \\
  w_2  \times  \frac{n_s}{c}; & \text{ if r is aligned}   \\
  w_3; & \text{ if r is full}   \\
\end{cases}
```
where ``n_s`` is the number of occurrences of the symbol ``s`` in the aligned row ``r``,
and ``c`` is the total number of columns in the row. The value of ``x`` is equal to zero
if every symbol in the row ``r`` occurred at most once, otherwise ``x`` is equal to the
max number of occurrences (matches) of some symbol in ``r``.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.weight(M[1, :])
 1.0

 juila> NetMSA.weight(M[2, :])
 0.1875
```
"""
function weight(row; w1=0.25, w2=0.5, w3=1.0)
  if full(row)
    return w3;
    end

  c = length(row);
  max = mostfrequent(row)[1];

  if aligned(row)
    return w2 * max / c;
  else
    x = max <= 1 ? 0 : max;
    return w1 * x / c;
  end
end

@doc raw"""
    objective(M, rowindex; endindex=0)

Return the objective score of the row, calculated as follows:

```math
f(x_s(r)) = \frac{A(r) \times C(r)}{1 + Gaps(r)} \times  \sum_{j=r}^{k} w(j)
```
where ``A(r)`` is the number of aligned rows in M from ``r`` to the last row,
``C(r)`` is the maximum number of matched charachters in the current row,
``Gaps(r)`` is the number of gaps added to the matrix M from row ``r`` to the
last row, and ``w(r)`` is the weight of the row ``r``.

`endindex` is used to reduce the search area for Gaps, and if it is not provided,
it would default to `size(M)[1]`.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

 juila> NetMSA.objective(M, 2)
 2.625
```
"""
function objective(M, rowindex; endindex=0)
  weights = sum(weight.(eachrow(M[rowindex:end, :])))
  C = mostfrequent(M[rowindex, :])[1];
  A = sum(aligned.(eachrow(M))[rowindex:end])

  endindex = endindex == 0 ? size(M)[1] : endindex;
  if endindex > size(M)[1]
    throw(ArgumentError("endindex exceeds the matrix size"));
    end

  counts = countmap(M[rowindex:endindex, :]);
  Gaps = get(counts, '-', 0);

  return weights * (A * C) / (1 + Gaps)
end

"""
    createswarm(rowindex::Int64, M)

Create a swarm containing unique Particles in the current row.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

juila> NetMSA.createswarm(2, M)
2-element Array{NetMSA.Particle,1}:
 NetMSA.Particle('c', 0, NetMSA.Position(2, [2]), NetMSA.Position(2, [2]), 0.0)
 NetMSA.Particle('b', 0, NetMSA.Position(2, [1, 3, 4]), NetMSA.Position(2, [1, 3, 4]), 0.0)
```
"""
function createswarm(rowindex::Int64, M)
  unique = Set(skipmissing(M[rowindex, :]))
  swarm = Vector{Particle}(undef, length(unique))
  for (i, c) in enumerate(unique)
    swarm[i] = Particle(c, getposition(c, rowindex, M))
  end
  return swarm
end

function criteria3(p::Particle, newindex::Int, M)
  return length(p.pos.indexes) != length(getposition(p.value, newindex, M).indexes)
end

function criteria2(p::Particle; threshold::Int=5)
return p.updated > threshold;
end

"""
    stopcriteria(p::Particle, newindex, M; threshold::Int=5, debug=false)

Check whether a stopping criteria has been met. Two stopping criteria
are checked in this function:
  - Criteria 2: If a particle hasn't updated its best score in the last `threshold` turns.
  - Criteria 3: If a particle moves to a new row which already contains the same symbol as that of the particle.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

juila> p = NetMSA.Particle('b', NetMSA.getposition('b', 2, M));
NetMSA.Particle('b', 0, Main.NetMSA.Position(2, [1, 3, 4]), Main.NetMSA.Position(2, [1, 3, 4]), 0.0)

julia> NetMSA.stopcriteria(p, 3, M; debug=true)
"Terminating because of criteria 3"
true
```
"""
function stopcriteria(p::Particle, newindex, M; threshold::Int=5, debug=false)
  c3 = criteria3(p, newindex, M);
  c2 = criteria2(p; threshold=threshold);
  if debug && c3
    display("Terminating because of criteria 3")
  elseif debug && c2
    display("Terminating because of criteria 2")
  end
  return c3 || c2;
end

function remove_missing_rows(M)
  return M[[length(Set(skipmissing(r))) != 0 for r in eachrow(M)], :]
end

"""
    flydown(p, M; stride=1)

Fly down the given particle by `stride`.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

juila> p = NetMSA.Particle('b', NetMSA.getposition('b', 2, M));
NetMSA.Particle('b', 0, Main.NetMSA.Position(2, [1, 3, 4]), Main.NetMSA.Position(2, [1, 3, 4]), 0.0)

julia> NetMSA.flydown(p, M)
9×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 '-'  'c'      '-'      '-'
 'b'  'b'      'b'      'b'
 'c'  'c'      'c'      'c'
 'b'  'f'      'h'      'b'
 'c'  'g'      'i'      'c'
 'd'  missing  'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

julia> NetMSA.flydown(p, M; stride=3)
11×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 '-'  'c'      '-'      '-'
 '-'  'b'      '-'      '-'
 '-'  'c'      '-'      '-'
 'b'  'f'      'b'      'b'
 'c'  'g'      'c'      'c'
 'b'  missing  'h'      'b'
 'c'  missing  'i'      'c'
 'd'  missing  'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'
```
"""
function flydown(p, M; stride=1)
  notpcols = setdiff(collect(1:size(M, 2)), p.pos.indexes)
  colsize = size(M, 2)
  pos = p.pos
  newrows = fill('-', (stride, colsize))
  M = vcat(M[1:pos.row - 1, :], reshape(newrows, stride, colsize), M[pos.row:end, :])
  for i in collect(pos.row + stride:size(M, 1))
    M[i - stride,notpcols] = M[i, notpcols]
    M[i, notpcols] .= missing
  end
  M = remove_missing_rows(M)
end

"""
    rowalignment(rowindex, M)

Return Particle with best position that aligns (maximizes the objective score)
the given row in the matrix.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

juila> p = NetMSA.Particle('b', NetMSA.getposition('b', 2, M));
NetMSA.Particle('b', 0, Main.NetMSA.Position(2, [1, 3, 4]), Main.NetMSA.Position(2, [1, 3, 4]), 0.0)

julia> NetMSA.rowalignment(2, M)
NetMSA.Particle('c', 0, Main.NetMSA.Position(2, [2]), Main.NetMSA.Position(3, [1]), 9.0)
```
"""
function rowalignment(rowindex, M)
  row = M[rowindex, :];
  if aligned(row)
    return nothing;
  end

  swarm = createswarm(rowindex, M);

  gₒ = g = swarm[1];
  gₒvalue = gvalue = objective(M, rowindex, endindex=rowindex);
  cols = size(row, 1)

  for p in swarm

    t = rowindex;
    N = copy(M);

    p.bestvalue = objective(M, rowindex, endindex=t)

    missingp = setdiff(collect(1:size(N, 2)), p.pos.indexes)
    maxlen = maximum([length(collect(skipmissing(col))) for col in eachcol(N[:, missingp])])
    criteria1 = maxlen;

    while stopcriteria(p, t, N) != true && t < criteria1
      t += 1;
      p.updated += 1;

      N = flydown(p, N);
      score = objective(N, rowindex);
      if score > p.bestvalue
        p.bestvalue = score;
        p.updated = 0;
      end

      if score > gvalue
        gvalue = score;
        g = deepcopy(p);
        g.best = getposition(p.value, t, N);
        g.bestvalue = score;
    end

    end
  end

  if gvalue == gₒvalue
    return nothing;
  end

  return g;
end

"""
    matrixalignment(M)

Align the matrix using NetMSA algorithm.

# Examples
```jldoctest
julia> M = NetMSA.createPeerMatrix(["abcbcdem", "acbcfg", "abchimn", "abcbcjkm"])
8×4 Array{Union{Missing, Char},2}:
 'a'  'a'      'a'      'a'
 'b'  'c'      'b'      'b'
 'c'  'b'      'c'      'c'
 'b'  'c'      'h'      'b'
 'c'  'f'      'i'      'c'
 'd'  'g'      'm'      'j'
 'e'  missing  'n'      'k'
 'm'  missing  missing  'm'

juila> NetMSA.matrixalignment(M)
9×4 Array{Union{Missing, Char},2}:
 'a'  'a'  'a'  'a'
 'b'  '-'  'b'  'b'
 'c'  'c'  'c'  'c'
 'b'  'b'  '-'  'b'
 'c'  'c'  '-'  'c'
 'd'  'f'  'h'  'j'
 'e'  'g'  'i'  'k'
 'm'  '-'  'm'  'm'
 '-'  '-'  'n'  '-'
```
"""
function matrixalignment(M)
  for (index, row) in enumerate(eachrow(M))
    g = rowalignment(index, M)
    if !isnothing(g)
      M = flydown(g, M, stride=g.best.row - g.pos.row)
    end
  end
  replace!(M, missing => '-')
end

end


text:
using HCIToolbox
using LossFunctions
using Optim
using PSFModels: Gaussian

# generative model
function model(X::AbstractVector{T}) where T
    position = @view X[1:2] # x, y position
    fwhm     = X[3]         # fwhm
    amp      = exp(X[4])    # amplitude
    return amp * Gaussian{T}(position, fwhm)
end

# objective function
function loss(X::AbstractVector{T}, target::SubArray) where T
    # cheap way to enforce positivity
    all(>(0), X) || return T(Inf)
    # get generative model
    m = model(X)
    # l2-distance loss (χ² loss) (LossFunctions.jl)
    stamp = @view m[target.indices...]
    return value(L2DistLoss(), target, stamp, AggMode.Sum())
end

"""
    gaussian_fit(frame, guess=center(frame))

Fit a 2D Gaussian PSF profile to `frame` centered around `guess` and return the best fitting parameters `[x, y, FWHM, log(Amp)]`.
"""
function gaussian_fit(frame::AbstractMatrix{T}, guess=center(frame); sub_size=100) where T
    # crop target view for speed
    half_size = sub_size ÷ 2
    indsx = round(Int, guess[1]-half_size):round(Int, guess[1]+half_size)
    indsy = round(Int, guess[2]-half_size):round(Int, guess[2]+half_size)
    target = @view frame[indsy, indsx]
    X0 = T[guess..., 7.5, log(maximum(target))]
    res = optimize(P -> loss(P, target), X0, LBFGS(); autodiff=:forward)
    # turn best-fit location into index
    return Optim.minimizer(res)
end

"""
    gaussian_fit_offset(frame, guess=center(frame))

Fit a 2D Gaussian PSF profile to `frame` centered around `guess` and return the offset from the center. Combine this with `HCIToolbox.shift_frame` for registering frames.
"""
function gaussian_fit_offset(frame::AbstractMatrix{T}, args...; kwargs...) where T
    Xbest = gaussian_fit(frame, args...; kwargs...)
    cy, cx = center(frame)
    dx = cx - Xbest[1]
    dy = cy - Xbest[2]
    return dx, dy
end

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
function updateterms!(system::CVSystem,n::Int64,i::Int64)
    # arteriolar compartment V, P
    system.branches.term[i].V[n+2,2] = (
        system.branches.term[i].V[n+1,2]+
        system.solverparams.h*(system.branches.term[i].Q[n+1,1]-
        system.branches.term[i].Q[n+1,2]));
    system.branches.term[i].P[n+2,2] = (
        1/system.branches.term[i].C[2]*(
        system.branches.term[i].V[n+2,2]-
        system.branches.term[i].V0[2]));
    # capillary compartment V, P
    system.branches.term[i].V[n+2,3] = (
        system.branches.term[i].V[n+1,3]+
        system.solverparams.h*(system.branches.term[i].Q[n+1,2]-
        system.branches.term[i].Q[n+1,3]));
    system.branches.term[i].P[n+2,3] = (
        1/system.branches.term[i].C[3]*(
        system.branches.term[i].V[n+2,3]-
        system.branches.term[i].V0[3]));
    # venular compartment V, P
    system.branches.term[i].V[n+2,4] = (
        system.branches.term[i].V[n+1,4]+
        system.solverparams.h*(system.branches.term[i].Q[n+1,3]-
        system.branches.term[i].Q[n+1,4]));
    system.branches.term[i].P[n+2,4] = (
        1/system.branches.term[i].C[4]*(
        system.branches.term[i].V[n+2,4]-
        system.branches.term[i].V0[4]));
    # venous compartment
    if i == 12 # inferior mesenteric flows into splenic
        system.branches.term[i].V[n+2,5] = (
            system.branches.term[i].V[n+1,5] + system.solverparams.h*
            (system.branches.term[21].Q[n+1,5] + system.branches.term[i].Q[n+1,4]-
            system.branches.term[i].Q[n+1,5]));
    else
        system.branches.term[i].V[n+2,5] = (
            system.branches.term[i].V[n+1,5]+
            system.solverparams.h*(system.branches.term[i].Q[n+1,4]-
            system.branches.term[i].Q[n+1,5]));
    end
    system.branches.term[i].P[n+2,5] = (
        1/system.branches.term[i].C[5]*(
        system.branches.term[i].V[n+2,5]-
        system.branches.term[i].V0[5]));
    if system.branches.group[i] == "lower"
        if i == 12 || i == 13 || i == 15 # divert flow to hepatic portal vein
            system.branches.term[i].Q[n+2,5] = (
                (1 - system.branches.term[i].R[5]/
                system.branches.term[i].L[5]*system.solverparams.h)*
                system.branches.term[i].Q[n+1,5]+
                system.solverparams.h/system.branches.term[i].L[5]*(
                system.branches.term[i].P[n+1,5]-
                system.liver.P[n+1,1]));
        elseif i == 14 # hepatic artery flows into hepatic veins
            system.branches.term[i].Q[n+2,5] = (
                (1 - system.branches.term[i].R[5]/
                system.branches.term[i].L[5]*system.solverparams.h)*
                system.branches.term[i].Q[n+1,5]+
                system.solverparams.h/system.branches.term[i].L[5]*(
                system.branches.term[i].P[n+1,5]-
                system.liver.P[n+1,4]));
        elseif i == 21 # inferior mesenteric joins splenic
            system.branches.term[i].Q[n+2,5] = (
                (1 - system.branches.term[i].R[5]/
                system.branches.term[i].L[5]*system.solverparams.h)*
                system.branches.term[i].Q[n+1,5]+
                system.solverparams.h/system.branches.term[i].L[5]*(
                system.branches.term[i].P[n+1,5]-
                system.branches.term[12].P[n+1,5]));
        else
            system.branches.term[i].Q[n+2,5] = (
                (1 - system.branches.term[i].R[5]/
                system.branches.term[i].L[5]*system.solverparams.h)*
                system.branches.term[i].Q[n+1,5]+
                system.solverparams.h/system.branches.term[i].L[5]*(
                system.branches.term[i].P[n+1,5]-
                system.ivc.P[n+1]));
        end
    else
        system.branches.term[i].Q[n+2,5] = (
            (1 - system.branches.term[i].R[5]/
            system.branches.term[i].L[5]*system.solverparams.h)*
            system.branches.term[i].Q[n+1,5]+
            system.solverparams.h/system.branches.term[i].L[5]*(
            system.branches.term[i].P[n+1,5]-
            system.svc.P[n+1]));
    end
    # Q for other compartments
    system.branches.term[i].Q[n+2,1] = (
        (system.branches.term[i].P[n+2,1]-
        system.branches.term[i].P[n+2,2])/
        system.branches.term[i].R[1]);
    system.branches.term[i].Q[n+2,2] = (
        (system.branches.term[i].P[n+2,2]-
        system.branches.term[i].P[n+2,3])/
        system.branches.term[i].R[2]);
    system.branches.term[i].Q[n+2,3] = (
        (system.branches.term[i].P[n+2,3]-
        system.branches.term[i].P[n+2,4])/
        system.branches.term[i].R[3]);
    system.branches.term[i].Q[n+2,4] = (
        (system.branches.term[i].P[n+2,4]-
        system.branches.term[i].P[n+2,5])/
        system.branches.term[i].R[4]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using SIMDDualNumbers, Aqua, ForwardDiff, VectorizationBase
using Test

function toaos(d::ForwardDiff.Dual{TAG,Vec{W,T},P}) where {TAG,W,T,P}
  let v = ForwardDiff.value(d),
    p = ForwardDiff.partials(d)
    ntuple(Val(W)) do w
      ForwardDiff.Dual{TAG}(v(w), ForwardDiff.Partials( ntuple(i -> p[i](w), Val(P)) ))
    end
  end
end

function test(ref, vanswer)
  answer = toaos(vanswer)
  for i ∈ eachindex(ref)
    @test ref[i] ≈ answer[i]
  end
end

@testset "SIMDDualNumbers.jl" begin

  dx = ForwardDiff.Dual(
    Vec(ntuple(_ -> randn(), VectorizationBase.pick_vector_width(Float64))...),
    Vec(ntuple(_ -> randn(), VectorizationBase.pick_vector_width(Float64))...),
    Vec(ntuple(_ -> randn(), VectorizationBase.pick_vector_width(Float64))...)
  )
  dy = ForwardDiff.Dual(
    Vec(ntuple(_ -> rand(), VectorizationBase.pick_vector_width(Float64))...),
    Vec(ntuple(_ -> rand(), VectorizationBase.pick_vector_width(Float64))...),
    Vec(ntuple(_ -> rand(), VectorizationBase.pick_vector_width(Float64))...)
  )

  dxaos = toaos(dx)
  dyaos = toaos(dy)
  for uf ∈ [SIMDDualNumbers.tanh_fast, SIMDDualNumbers.sigmoid_fast, abs, VectorizationBase.relu]
    test(map(uf, dxaos), uf(dx))
  end
  for bf ∈ [max, min]
    test(map(bf, dxaos, dyaos), bf(dx, dy))
  end

  vz = Vec(ntuple(_ -> rand(), VectorizationBase.pick_vector_width(Float64))...)
  tz = Tuple(vz)
  cmp = dx > dy
  cmpaos = dxaos .> dyaos
  test(map(ifelse, cmpaos, dxaos, dyaos), SIMDDualNumbers.ifelse(cmp, dx, dy))
  test(map(ifelse, cmpaos, tz, dyaos), SIMDDualNumbers.ifelse(cmp, vz, dy))
  test(map(ifelse, cmpaos, dxaos, tz), SIMDDualNumbers.ifelse(cmp, dx, vz))

  Aqua.test_all(SIMDDualNumbers, ambiguities=false) #TODO: test ambiguities once ForwardDiff fixes them, or once ForwardDiff is dropped
end


text:
# -----------------------------------------------------
# Precompute relevant 1-D basis elements
# -----------------------------------------------------

# Efficiency criticality: MED
# Computation done once, but is done upon import
# This could be improved

KMAX = 5
LMAX = 10

precomputed_diffs = Dict{NTuple{4, Int}, Array{Float64, 2}}()

function precompute_diffs()
	if length(precomputed_diffs) == 0
		print("Precomputing 1D derivative matrix elements... ")
		flush(STDOUT)
		for k in 1:KMAX
			for level in 0:LMAX
			    for cell in 1:(1<<max(0, level-1))
			        for mode in 1:k
			            precomputed_diffs[(k,level,cell,mode)] = diff_basis_DG(k,level,cell,mode)
			        end
			    end
			end
		end
		println("done.")
		flush(STDOUT)
		return true
	end
	return false
end


text:
#$##############################################################################
#
#    Misc/Infinite.jl: Infinity
#
# This file is part of Hecke.
#
# Copyright (c) 2015-2019: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2019 Tommy Hofmann
#
################################################################################

export PosInf, inf, IntExt

# This is a type for positive infinity for use in valuations.

struct PosInf
end

const inf = PosInf()

+(::Int, ::PosInf) = inf

+(::PosInf, ::Int) = inf

+(::PosInf, ::PosInf) = inf

max(::Int, ::PosInf) = inf

max(::PosInf, ::Int) = inf

Base.isless(::Int, ::PosInf) = true

Base.isless(::PosInf, ::Int) = false

Base.isfinite(::PosInf) = false

Base.iszero(::PosInf) = false

const IntExt = Union{Int, PosInf}


text:
module Messages
    include("../../utils/Utils.jl")
    import .Utils

    using AMQPClient

    export create_and_publish

    """Функция создания сообщения с его последующей публикацией.

    channel: канал связи с сервисом очередей,
    exchanger: имя обменника,
    route: имя маршрута,
    data: данные, которые требуется передать в сообщении,
    props: настройки для передаваемого сообщения
    """
    function create_and_publish(channel::AMQPClient.MessageChannel, exchanger::String, route::String, data::Any, props::Dict{String, Any})::Bool
        try
            data_to_send = Vector{UInt8}(data)
            message::AMQPClient.Message = AMQPClient.Message(
                data_to_send,
                content_type= get(props, "CONTENT_TYPE", nothing),
                content_encoding = get(props, "CONTENT_ENCODING", nothing),
                headers = get(props, "HEADERS", nothing),
                delivery_mode = get(props, "DELIVERY_MODE", nothing),
                priority = get(props, "PRIORITY", 0),
                correlation_id = get(props, "CORRELATION_ID", nothing),
                reply_to = get(props, "REPLY_TO", nothing),
                expiration = get(props, "EXPIRATION", nothing),
                message_id = get(props, "MESSAGE_ID", nothing),
                timestamp = get(props, "TIMESTAMP", nothing),
                message_type = get(props, "MESSAGE_TYPE", nothing),
                user_id = get(props, "USER_ID", nothing),
                app_id = get(props, "APPLICATION_ID", nothing)
            )
            AMQPClient.basic_publish(channel, message; exchange=exchanger, routing_key=route)
            return true
        catch error
            @error Utils.Chronometer.message_with_time("Создание и отправка сообщения завершилось ошибой") error
            throw(error)
        end
    end

end


text:
# ffr, 19 April 2016

# for the moment, I decided to not include L, because it can
# be calculated from A and B
struct LF1dGrid
    N::Int64
    A::Float64
    B::Float64
    h::Float64
    grid::Array{Float64,1}
    D1jl::Array{Float64,2}
    D2jl::Array{Float64,2}
end

# Periodic LF
function init_LF1d_p( N::Int64, A::Float64, B::Float64, verbose=false )
    # Check argument
    if N % 2 == 0
        error("N must be an odd number\n")
    end
    L = (B - A)
    h = L/N
    grid = zeros(Float64,N)
    for i = 1:N
        grid[i] = A + 0.5*(B-A)*(2*i-1)/N
    end
    D1jl = zeros(Float64,N,N)
    D2jl = zeros(Float64,N,N)
    #
    # Diagonal elements
    #
    Nprimed = (N-1)/2
    for j = 1:N
        # the diagonal elements of D1jl are already zero
        D2jl[j,j] = -(2.0*pi/L)^2 * Nprimed * (Nprimed+1)/3.
    end
    #
    # Off diagonal elements
    #
    for j = 1 : N
        for l = j+1 : N
            #
            nn = j - l
            #
            tt1 = pi/L * (-1.0)^nn
            tt2 = sin(pi*nn/N)
            #
            tt3 = (2.0*pi/L)^2 * (-1.0)^nn * cos(pi*nn/N)
            tt4 = 2.0*sin(pi*nn/N)^2
            #
            D1jl[j,l] =  tt1/tt2
            D1jl[l,j] = -tt1/tt2
            #
            D2jl[j,l] = -tt3/tt4
            D2jl[l,j] = -tt3/tt4
        end
    end
    #
    LF = LF1dGrid( N, A, B, h, grid, D1jl, D2jl )
    if verbose
        @printf("Allocated: periodic 1d LBF grid: N, h = %d , %f\n", N, h)
    end
    #
    return LF
end # function init_LF1d_p

# cluster Lagrange function
function init_LF1d_c( N::Int64, A::Float64, B::Float64, verbose=false )
    L = B - A
    h = L/(N + 1.)
    #
    grid = zeros(Float64,N)
    for i = 1:N
        grid[i] = A + i*(B-A)/(N+1)
    end
    #
    D1jl = zeros(Float64,N,N)  # XXX this is not yet used
    D2jl = zeros(Float64,N,N)
    #
    # Diagonal part
    #
    pre = -pi^2/(2.0*L^2)
    for i = 1:N
        t1 = ( 2.0*(N+1)^2 + 1 )/3.0
        t2 = sin( i*pi/(N+1) )^2
        D2jl[i,i] = pre*( t1 - 1.0/t2 )
    end
    #
    # Off-diagonal
    #
    for l = 1 : N
        for j = l+1 : N
            nnm = l - j
            nnp = l + j
            pre = -pi^2 / (2*L^2) * (-1.0)^nnm
            t1 = sin( pi*nnm/2.0/(N+1) )^2
            t2 = sin( pi*nnp/2.0/(N+1) )^2
            #
            D2jl[l,j] = pre*( 1.0/t1 - 1.0/t2 )
            D2jl[j,l] = pre*( 1.0/t1 - 1.0/t2 )  # XXX is it faster to just recalculate?
        end
    end
    LF = LF1dGrid( N, A, B, h, grid, D1jl, D2jl )
    if verbose
        @printf("Allocated: cluster 1d LBF grid: N, h = %5d , %10.5f\n", N, h)
    end
    return LF
end

# Lagrange-sinc function
function init_LF1d_sinc( N::Int64, h::Float64, verbose=false )
    # Choice for A and B
    A = -(N-1)*h/2.0
    B =  (N-1)*h/2.0
    #
    grid = zeros(Float64,N)
    for i = 1:N
        grid[i] = A + (i-1)*h
    end
    #
    D1jl = zeros(Float64,N,N)  # XXX this is not yet used
    D2jl = zeros(Float64,N,N)
    #
    # Diagonal part
    #
    for i = 1:N
        D2jl[i,i] = -pi^2 / 3.0 / h^2
    end
    #
    # Off-diagonal
    #
    for j = 1 : N
        for i = j+1 : N
            D2jl[i,j] = -2.0*(-1.0)^(i-j)/( grid[i] - grid[j] )^2
            D2jl[j,i] = D2jl[i,j]  # XXX is it faster to just recalculate?
        end
    end
    LF = LF1dGrid( N, A, B, h, grid, D1jl, D2jl )
    if verbose
        @printf("Allocated: 1d Lagrange-sinc grid: N, h = %5d , %10.5f\n", N,  h)
    end
    return LF
end


text:
using Pkg
Pkg.activate(".")
using Test

@testset "1-data_loading_cleanup_processing" begin
	using CSV
	using DataFrames
	path = pwd();
	pythonfile = CSV.File(path*"/../data/cp_data_cleaned.csv") |> DataFrame
	juliafile = CSV.File(path*"/../data/cp_data_cleaned_jl.csv") |> DataFrame;
	@test  isequal(pythonfile[!,1],juliafile[!,1])
	@test all(isapprox.(pythonfile[!,2],juliafile[!,2],atol=1.0e-3))
	@test all(isapprox.(pythonfile[!,3],juliafile[!,3],atol=1.0e-3))
end

@testset "2-data_splitting" begin
	 using CSV
	 using DataFrames
	 path = pwd();
	 dataset = ("train","val","test")
	 @testset "Comparing number of formulae in $(d) set" for d in dataset
	 	pythonfile = CSV.File(path*"/../data/cp_$(d).csv") |> DataFrame
		juliafile = CSV.File(path*"/../data/cp_$(d)_jl.csv") |> DataFrame;
	 	@test length(unique(pythonfile[!,:formula])) == length(unique(juliafile[!,:formula]))
	end
end

@testset "3-modeling_classical_models" begin
	 path= pwd();
	 juliavalue = parse(Float64,readline(open(path * "/../data/R2_GradientBoostingRegressor.fit")))
	 @test (juliavalue > 0.8442) && (juliavalue < 0.9900)
end


text:
module Lunettes

#Very basic struct
"""
   my_lens = Lens{A,B}()

Creates a new `Lens` for which the functions `getr` and `setr` can be dispatched on.
"""
struct Lens{T,U} end

#compose two lenses
"""
   Lens{A,B}() ∘ Lens{C,D}()

Create a new `Lens` based upon the composition of two other `Lens`es.
This composed `Lens` will automatically have defined for it both
`getr` and `setr` in such a way which is consist with the Lens composition law.
"""
function Base.:∘(::Lens{T,U},::Lens{V,W}) where {T,U,V,W}
    Lens{Lens{T,U}, Lens{V,W}}()
end

# getter law for composed lenses
"""
   getr(::Lens{A,B}, a)

Returns the value that Lens{A,B}() "sees" for input a.

Note that `getr` is automatically curried i.e.
```
   getr(Lens{A,B}(), a) == getr(Lens{A,B}())(a)
```
"""
function getr(::Lens{Lens{T,U},Lens{V,W}}, a) where {T,U,V,W}
    getr(Lens{V,W}(), getr(Lens{T,U}(), a))
end
# curried version
function getr(l::Lens)
    function out(a)
        getr(l,a)
    end
end

# setter law for composed lens
"""
   setr(::Lens{A,B}, a, c)

Takes value `a` and returns a new value similar to `a` except `c` is now in the position
that Lens{A,B}() "sees".

Note that `setr` is automatically curried i.e.
```
setr(Lens{A,B}(), a, c) == setr(Lens{A,B}(), c)(a)
```

"""
function setr(::Lens{Lens{T,U},Lens{V,W}}, a, c) where {T,U,V,W}
    setr(Lens{T,U}(), a, setr(Lens{V,W}(), getr(Lens{T,U}(), a), c))
end
# curried version
function setr(l::Lens, c)
    function out(a)
        setr(l, a, c)
    end
end

array_lens(N) = Lens{Vector,N}()
function getr(::Lens{Vector,N}, a) where {N}
    a[N]
end
function setr(::Lens{Vector,N}, a, c) where {N}
    if length(a) < N
        throw("Bounds Error")
    end
    map(1:length(a)) do index
        if index == N
            c
        else
            a[index]
        end
    end
end

include("./macro.jl")
export Lens, setr, @lens, array_lens, getr

end  #module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# Porting note:
# There is no need for a `Thing` type because it's an abstract thing that could just be
# an entity, relation, of attribute. See `instantiate` function below.

function instantiate(t::Proto.Thing)
    if t._type.encoding == Proto.Type_Encoding.ENTITY_TYPE
        return Entity(t)
    elseif t._type.encoding == Proto.Type_Encoding.RELATION_TYPE
        return Relation(t)
    elseif t._type.encoding == Proto.Type_Encoding.ATTRIBUTE_TYPE
        return Attribute(t)
    else
        throw(TypeDBClientException(CONCEPT_BAD_ENCODING))
    end
end

"""
    set_has(transaction::AbstractCoreTransaction, thing::AbstractThing, attribute::Attribute)

For a given Thing the attribute will set as associated.
"""
function set_has(transaction::AbstractCoreTransaction, thing::AbstractThing, attribute::Attribute)
    res_proto_attribute = get_proto_thing(ConceptManager(transaction), attribute.iid)
    has_req = ThingRequestBuilder.set_has_req(thing.iid, res_proto_attribute)
    execute(transaction, has_req)

    return nothing
end

"""
    unset_has(transaction::AbstractCoreTransaction, thing::AbstractThing, attribute::Attribute)

Here an given attribute will be detached from a given thing
"""
function unset_has(transaction::AbstractCoreTransaction, thing::AbstractThing, attribute::Attribute)
    res_proto_attribute = get_proto_thing(ConceptManager(transaction), attribute.iid)
    unset_has_req = ThingRequestBuilder.unset_has_req(thing.iid, res_proto_attribute)
    execute(transaction, unset_has_req)

    return nothing
end

"""
    get_has(transaction::AbstractCoreTransaction,
        thing::AbstractThing,
        attribute_type::Optional{AttributeType} = nothing,
        attribute_types::Optional{Vector{<:AbstractAttributeType}} = nothing,
        keys_only = false)

With this function it is possible to get attributes of an entity or a relation.
Optional it is possible to restrict the type(s) of attributes which will retrieved from
the database. You can decide between only one ore more than one type.
"""
function get_has(transaction::AbstractCoreTransaction,
                thing::AbstractThing,
                attribute_type::Optional{AttributeType} = nothing,
                attribute_types::Optional{Vector{<:AbstractAttributeType}} = nothing,
                keys_only = false)

    if count([attribute_type !== nothing, attribute_types !== nothing, keys_only]) > 1
        throw(TypeDBClientException(GET_HAS_WITH_MULTIPLE_FILTERS))
    end

    attribute_types_intern = attribute_type !== nothing ? [attribute_type] : []
    attribute_types_intern = attribute_types !== nothing ? attribute_types : attribute_types_intern

    all_attribute_types = get_owns(as_remote(thing.type, transaction))

    res_attr_types = !isempty(intersect(attribute_types_intern, all_attribute_types)) ?
                        intersect(attribute_types_intern, all_attribute_types) :
                        all_attribute_types

    attribute_types_proto = proto.(res_attr_types)
    has_req = keys_only ?
                ThingRequestBuilder.get_has_req(thing.iid, keys_only) :
                ThingRequestBuilder.get_has_req(thing.iid, attribute_types_proto)
    res_has = stream(transaction, has_req)

    return instantiate.(collect(Iterators.flatten(
        r.thing_res_part.thing_get_has_res_part.attributes for r in res_has)))
end

"""
    delete(r::RemoteConcept{<:AbstractThing})
To delete an Entity, Relation or Attribute this have to be packed with the transaction
to a RemoteConcept. To accomplish this see the function as_remote.
"""
function delete(r::RemoteConcept{<:AbstractThing})
    del_req = ThingRequestBuilder.delete_req(r.concept.iid)
    execute(r.transaction, del_req)

    return nothing
end

function is_deleted(r::RemoteConcept{<:AbstractThing})
    res_thing = get(ConceptManager(r.transaction), r.concept.iid)
    return res_thing === nothing
end


text:
import Base.convert

# generated datatypes
include("generated/hsa_datatypes.jl")
# custom datatypes and wrappers
include("custom/all.jl")
# generated wrappers
# these have to go last because they use the custom datatypes
include("generated/hsa_wrappers.jl")


text:
module RecyclersTests

include("utils.jl")
include("test_recyclers.jl")
include("test_global.jl")
include("test_aqua.jl")

end  # module RecyclersTests