text:
module Sch

using Distributed
import MemPool: DRef

import ..Dagger: Context, Thunk, Chunk, OSProc, order, free!, dependents, noffspring, istask, inputs, affinity, tochunk, @dbg, @logmsg, timespan_start, timespan_end, unrelease, procs

include("fault-handler.jl")

const OneToMany = Dict{Thunk, Set{Thunk}}

"""
    ComputeState

The internal state-holding struct of the scheduler.

Fields
- dependents::OneToMany - The result of calling `dependents` on the DAG
- finished::Set{Thunk} - The set of completed `Thunk`s
- waiting::OneToMany - Map from parent `Thunk` to child `Thunk`s that still need to execute
- waiting_data::OneToMany - Map from child `Thunk` to all parent `Thunk`s, accumulating over time
- ready::Vector{Thunk} - The list of `Thunk`s that are ready to execute
- cache::Dict{Thunk, Any} - Maps from a finished `Thunk` to it's cached result, often a DRef
- running::Set{Thunk} - The set of currently-running `Thunk`s
- thunk_dict::Dict{Int, Any} - Maps from thunk IDs to a `Thunk`
"""
struct ComputeState
    dependents::OneToMany
    finished::Set{Thunk}
    waiting::OneToMany
    waiting_data::OneToMany
    ready::Vector{Thunk}
    cache::Dict{Thunk, Any}
    running::Set{Thunk}
    thunk_dict::Dict{Int, Any}
end

"""
    SchedulerOptions

Stores DAG-global options to be passed to the Dagger.Sch scheduler. Options:
- single::Int   # Force all work onto worker with specified id. `0` disables this option.
- threads::Bool # Use multithreading if available
"""
Base.@kwdef struct SchedulerOptions
    single::Int = 0
    threads::Bool = false
end

"""
    ThunkOptions

Stores Thunk-local options to be passed to the Dagger.Sch scheduler. Options:
- single::Int   # Force thunk onto worker with specified id. `0` disables this option.
- threads::Bool # Use multithreading if available
"""
Base.@kwdef struct ThunkOptions
    single::Int = 0
    threads::Bool = false
end

function cleanup(ctx)
end

function compute_dag(ctx, d::Thunk; options=SchedulerOptions())
    if options === nothing
        options = SchedulerOptions()
    end
    ctx.options = options
    master = OSProc(myid())
    @dbg timespan_start(ctx, :scheduler_init, 0, master)

    if options.single !== 0
        @assert options.single in vcat(1, workers()) "Sch option 'single' must specify an active worker id"
        ps = OSProc[OSProc(options.single)]
    else
        ps = procs(ctx)
    end
    chan = Channel{Any}(32)
    deps = dependents(d)
    ord = order(d, noffspring(deps))

    node_order = x -> -get(ord, x, 0)
    state = start_state(deps, node_order)
    # start off some tasks
    for p in ps
        isempty(state.ready) && break
        task = pop_with_affinity!(ctx, state.ready, p, false)
        if task !== nothing
            fire_task!(ctx, task, p, state, chan, node_order)
        end
    end
    @dbg timespan_end(ctx, :scheduler_init, 0, master)

    # Loop while we still have thunks to execute
    while !isempty(state.ready) || !isempty(state.running)
        if isempty(state.running) && !isempty(state.ready)
            # Nothing running, so schedule up to N thunks, 1 per N workers
            for p in ps
                isempty(state.ready) && break
                task = pop_with_affinity!(ctx, state.ready, p, false)
                if task !== nothing
                    fire_task!(ctx, task, p, state, chan, node_order)
                end
            end
        end

        if isempty(state.running)
            # the block above fired only meta tasks
            continue
        end

        proc, thunk_id, res = take!(chan) # get result of completed thunk
        if isa(res, CapturedException) || isa(res, RemoteException)
            if check_exited_exception(res)
                @warn "Worker $(proc.pid) died on thunk $thunk_id, rescheduling work"

                # Remove dead worker from procs list
                filter!(p->p.pid!=proc.pid, ctx.procs)
                ps = procs(ctx)

                handle_fault(ctx, state, state.thunk_dict[thunk_id], proc, chan, node_order)
                continue
            else
                throw(res)
            end
        end
        node = state.thunk_dict[thunk_id]
        @logmsg("WORKER $(proc.pid) - $node ($(node.f)) input:$(node.inputs)")
        state.cache[node] = res

        @dbg timespan_start(ctx, :scheduler, thunk_id, master)
        immediate_next = finish_task!(state, node, node_order)
        if !isempty(state.ready)
            thunk = pop_with_affinity!(Context(ps), state.ready, proc, immediate_next)
            if thunk !== nothing
                fire_task!(ctx, thunk, proc, state, chan, node_order)
            end
        end
        @dbg timespan_end(ctx, :scheduler, thunk_id, master)
    end
    state.cache[d]
end

function pop_with_affinity!(ctx, tasks, proc, immediate_next)
    # allow JIT specialization on Pairs
    mapfirst(c) = first.(c)

    if immediate_next
        # fast path
        if proc in mapfirst(affinity(tasks[end]))
            return pop!(tasks)
        end
    end

    # TODO: use the size
    parent_affinity_procs = Vector(undef, length(tasks))
    # parent_affinity_sizes = Vector(undef, length(tasks))
    for i=length(tasks):-1:1
        t = tasks[i]
        aff = affinity(t)
        aff_procs = mapfirst(aff)
        if proc in aff_procs
            deleteat!(tasks, i)
            return t
        end
        parent_affinity_procs[i] = aff_procs
    end
    for i=length(tasks):-1:1
        # use up tasks without affinities
        # let the procs with the respective affinities pick up
        # other tasks
        aff_procs = parent_affinity_procs[i]
        if isempty(aff_procs)
            t = tasks[i]
            deleteat!(tasks, i)
            return t
        end
        if all(!(p in aff_procs) for p in procs(ctx))
            # no proc is ever going to ask for it
            t = tasks[i]
            deleteat!(tasks, i)
            return t
        end
    end
    return nothing
end

function fire_task!(ctx, thunk, proc, state, chan, node_order)
    @logmsg("W$(proc.pid) + $thunk ($(showloc(thunk.f, length(thunk.inputs)))) input:$(thunk.inputs) cache:$(thunk.cache) $(thunk.cache_ref)")
    push!(state.running, thunk)
    if thunk.cache && thunk.cache_ref !== nothing
        # the result might be already cached
        data = unrelease(thunk.cache_ref) # ask worker to keep the data around
                                          # till this compute cycle frees it
        if data !== nothing
            @logmsg("cache hit: $(thunk.cache_ref)")
            state.cache[thunk] = data
            immediate_next = finish_task!(state, thunk, node_order; free=false)
            if !isempty(state.ready)
                thunk = pop_with_affinity!(ctx, state.ready, proc, immediate_next)
                if thunk !== nothing
                    fire_task!(ctx, thunk, proc, state, chan, node_order)
                end
            end
            return
        else
            thunk.cache_ref = nothing
            @logmsg("cache miss: $(thunk.cache_ref) recomputing $(thunk)")
        end
    end

    if thunk.meta
        # Run it on the parent node
        # do not _move data.
        p = OSProc(myid())
        @dbg timespan_start(ctx, :comm, thunk.id, p)
        fetched = map(thunk.inputs) do x
            istask(x) ? state.cache[x] : x
        end
        @dbg timespan_end(ctx, :comm, thunk.id, p)

        @dbg timespan_start(ctx, :compute, thunk.id, p)
        res = thunk.f(fetched...)
        @dbg timespan_end(ctx, :compute, thunk.id, p)

        #push!(state.running, thunk)
        state.cache[thunk] = res
        immediate_next = finish_task!(state, thunk, node_order; free=false)
        if !isempty(state.ready)
            if immediate_next
                thunk = pop!(state.ready)
            else
                thunk = pop_with_affinity!(ctx, state.ready, proc, immediate_next)
            end
            if thunk !== nothing
                fire_task!(ctx, thunk, proc, state, chan, node_order)
            end
        end
        return
    end

    data = map(thunk.inputs) do x
        istask(x) ? state.cache[x] : x
    end
    state.thunk_dict[thunk.id] = thunk
    if thunk.options !== nothing && thunk.options.single > 0
        proc = OSProc(thunk.options.single)
    end
    async_apply(ctx, proc, thunk.id, thunk.f, data, chan, thunk.get_result, thunk.persist, thunk.cache, thunk.options)
end

function finish_task!(state, node, node_order; free=true)
    if istask(node) && node.cache
        node.cache_ref = state.cache[node]
    end
    immediate_next = false
    for dep in sort!(collect(state.dependents[node]), by=node_order)
        set = state.waiting[dep]
        pop!(set, node)
        if isempty(set)
            pop!(state.waiting, dep)
            push!(state.ready, dep)
            immediate_next = true
        end
        # todo: free data
    end
    for inp in inputs(node)
        if inp in keys(state.waiting_data)
            s = state.waiting_data[inp]
            if node in s
                pop!(s, node)
            end
            if free && isempty(s)
                if haskey(state.cache, inp)
                    _node = state.cache[inp]
                    free!(_node, force=false, cache=(istask(inp) && inp.cache))
                    pop!(state.cache, inp)
                end
            end
        end
    end
    push!(state.finished, node)
    pop!(state.running, node)
    immediate_next
end

function start_state(deps::Dict, node_order)
    state = ComputeState(
                  deps,
                  Set{Thunk}(),
                  OneToMany(),
                  OneToMany(),
                  Vector{Thunk}(undef, 0),
                  Dict{Thunk, Any}(),
                  Set{Thunk}(),
                  Dict{Int, Thunk}()
                 )

    nodes = sort(collect(keys(deps)), by=node_order)
    merge!(state.waiting_data, deps)
    for k in nodes
        if istask(k)
            waiting = Set{Thunk}(Iterators.filter(istask,
                                                  inputs(k)))
            if isempty(waiting)
                push!(state.ready, k)
            else
                state.waiting[k] = waiting
            end
        end
    end
    state
end

_move(ctx, to_proc, x) = x
_move(ctx, to_proc::OSProc, x::Union{Chunk, Thunk}) = collect(ctx, x)

@noinline function do_task(ctx, proc, thunk_id, f, data, send_result, persist, cache, options)
    @dbg timespan_start(ctx, :comm, thunk_id, proc)
    time_cost = @elapsed fetched = map(x->_move(ctx, proc, x), data)
    @dbg timespan_end(ctx, :comm, thunk_id, proc)

    @dbg timespan_start(ctx, :compute, thunk_id, proc)
    result_meta = try
        use_threads = (ctx.options !== nothing && ctx.options.threads) ||
                      (options !== nothing && options.threads)
        if use_threads
            @static if VERSION >= v"1.3.0-DEV.573"
                res = fetch(Threads.@spawn f(fetched...))
            else
                res = f(fetched...)
            end
        else
            res = f(fetched...)
        end
        (proc, thunk_id, send_result ? res : tochunk(res, persist=persist, cache=persist ? true : cache)) #todo: add more metadata
    catch ex
        bt = catch_backtrace()
        (proc, thunk_id, RemoteException(myid(), CapturedException(ex, bt)))
    end
    @dbg timespan_end(ctx, :compute, thunk_id, proc)
    result_meta
end

@noinline function async_apply(ctx, p::OSProc, thunk_id, f, data, chan, send_res, persist, cache, options)
    @async begin
        try
            put!(chan, remotecall_fetch(do_task, p.pid, ctx, p, thunk_id, f, data, send_res, persist, cache, options))
        catch ex
            bt = catch_backtrace()
            put!(chan, (p, thunk_id, CapturedException(ex, bt)))
        end
        nothing
    end
end

end # module Sch


text:
module ATN_Coupled_SIRD
export time_series_par
export time_series_eta
export time_series
using Statistics

function infected_c(
    AgentState_1,
    AgentState_2,
    AgentAct_1,
    AgentAct_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
    gamma_1,
    gamma_2,
)

    N = length(AgentAct_1)
    AgentStateNext_1 = zeros(Int, N)
    AgentStateNext_2 = zeros(Int, N)
    for k = 1:N
        AgentStateNext_1[k] = AgentState_1[k]
        AgentStateNext_2[k] = AgentState_2[k]
    end
    for i = 1:N
        #Change someone sick in group 1 recovers, mu_1
        if AgentState_1[i] == 1
            if rand() < mu_1
                AgentStateNext_1[i] = 2
            end
        end
        #Change someone sick in group 2 recovers, mu_2
        if AgentState_2[i] == 1
            if rand() < mu_2
                AgentStateNext_2[i] = 2
            end
        end

        # if AgentState_2[i] == 2
        #     if rand() < 0.01
        #         AgentStateNext_2[i] = 0
        #     end
        # end
        if AgentState_1[i] == 1
            if rand() < gamma_1
                AgentStateNext_1[i] = 3
            end
        end
        #Change someone sick in group 2 recovers, mu_2
        if AgentState_2[i] == 1
            if rand() < gamma_2
                AgentStateNext_2[i] = 3
            end
        end

        #If node is active
        if rand() < 5 * AgentAct_1[i]
            #it makes m connections

            for j = 1:m
                #choose a random index
                rInt = rand(1:N)
                #probability of staying is beta_1

                if AgentState_1[i] == 1
                    #the active node gets the new neigbor sick in POP_1
                    if ((AgentState_1[rInt] == 0) && (rand() < lambda_1))
                        AgentStateNext_1[rInt] = 1
                    end
                else
                    #the active node gets sick from the new neighbor node
                    if ((AgentState_1[rInt] == 1) && (rand() < lambda_1))
                        AgentStateNext_1[i] = 1
                    end
                end

                if rand() < beta_1
                    #the active node gets the new neigbor sick in POP_2
                    if AgentState_1[i] == 1
                        if ((AgentState_1[rInt] == 0) && (rand() < lambda_2))
                            AgentStateNext_2[rInt] = 1
                        end
                    else
                        if ((AgentState_2[rInt] == 1) && (rand() < lambda_1))
                            AgentStateNext_1[i] = 1
                        end
                    end
                end
            end
        end

        if rand() < 5 * AgentAct_2[i]
            #it makes m connections
            for j = 1:m
                #choose a random index
                rInt = rand(1:N)
                #probability of staying is beta_2

                if AgentState_2[i] == 1
                    #the active node gets the new neigbor sick in POP_1
                    if ((AgentState_2[rInt] == 0) && (rand() < lambda_2))
                        AgentStateNext_2[rInt] = 1
                    end
                else
                    #the active node gets sick from the new neighbor node
                    if ((AgentState_2[rInt] == 1) && (rand() < lambda_2))
                        AgentStateNext_2[i] = 1
                    end
                end

                if rand() < beta_2

                #the active node gets the new neigbor sick in POP_2
                    if AgentState_2[i] == 1
                        if ((AgentState_2[rInt] == 0) && (rand() < lambda_1))
                            AgentStateNext_1[rInt] = 1
                        end
                    else
                        if ((AgentState_1[rInt] == 1) && (rand() < lambda_2))
                            AgentStateNext_2[i] = 1
                        end
                    end

                end
            end

        end
    end
    return AgentStateNext_1, AgentStateNext_2

end

function infected_totals(
    tmax,
    agent_act_1,
    agent_act_2,
    m,
    mu_1,
    mu_2,
    lambda_1,
    lambda_2,
    beta_1,
    beta_2,
    gamma_1,
    gamma_2,
)
    N = length(agent_act_1)
    agent_state_1 = zeros(Int, N)
    agent_state_2 = zeros(Int, N)

    # Iend_1 = zeros(tmax)
    Iend_2 = zeros(tmax)
    for k = 1:100
        # agent_state_1[k] = 1
        agent_state_2[k] = 1
    end
    t_window = 500

    for i_time = 1:tmax
        agent_state_1, agent_state_2 = infected_c(
            agent_state_1,
            agent_state_2,
            agent_act_1,
            agent_act_2,
            m,
            mu_1,
            mu_2,
            lambda_1,
            lambda_2,
            beta_1,
            beta_2,
            gamma_1,
            gamma_2,
        )
    end
    dead_1 = 0.0
    dead_2 = 0.0
    for k = 1:N
        if agent_state_1[k] == 3
            dead_1 = dead_1 + 1
        end
        # if agent_state_2[k] == 3
        #     dead_2 = dead_2 + 1
        # end
    end
    return dead_1
end

function parameter_sweep(n_trials, agent_act_1, agent_act_2, grid)
    mu_1 = .1
    mu_2 = .1
    N = length(agent_act_1)
    beta_2 = grid[1]
    beta_1 = beta_2
    lambda_2 =  grid[2]
    lambda_1 = lambda_2
    m = 3
    tmax = 1000
    total_1 = 0.0
    total_2 = 0.0
    gamma_1 = .1
    gamma_2 = .01
    for k = 1:n_trials
        infected_1 = infected_totals(
            tmax,
            agent_act_1,
            agent_act_2,
            m,
            mu_1,
            mu_2,
            lambda_1,
            lambda_2,
            beta_1,
            beta_2,
            gamma_1,
            gamma_2,
        )
        total_1 = total_1 + infected_1 / (N)
        # total_2 = total_2 + infected_2 / (N)

    end
    return total_1 / n_trials
end
end

ext: jl


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbtestcrv()::NURBS1D

Construct a simple test curve.

# Output:
- `crv` NURBS data structure of a curve

# Examples:
```julia
julia> crv = nrbtestcrv()
```
"""
function nrbtestcrv()::NURBS1D
pnts = [0.5 1.5 4.5 3.0 7.5 6.0 8.5;
        3.0 5.5 5.5 1.5 1.5 4.0 4.5;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0];
return nrbmak(pnts,vec([0 0 0 1/4 1/2 3/4 3/4 1 1 1]));
end # nrbtestcrv


text:
abstract type ConjugatePrior end

struct OrderedNormalMeanPrior{T} <: ConjugatePrior
  priors::T
  function OrderedNormalMeanPrior(priors::T) where T
    # Assert that the first element is valid.
    priors[1] isa Normal ||
    priors[1] isa Truncated{<:Normal} ||
    error("Invalid prior. First element must be a Normal or Truncated{<:Normal}!")

    # Assert that all other elements are Truncated{<:Normal}
    if length(priors) > 1
      for k in 2:length(priors)
        priors[k] isa Truncated{<:Normal} ||
        error("Invalid prior. Element $(k) must be a Truncated{<:Normal}!")
      end
    end

    return new{T}(priors)
  end
end

function Distributions.rand(rng::Random.AbstractRNG, d::OrderedNormalMeanPrior)
  return rand.(rng, d.priors)
end

function OrderedNormalMeanPrior(K::Int, comp1::Union{Normal, Truncated{<:Normal}},
                                compk::Truncated{<:Normal})
  K >= 2 || error("`compk` was provided, but provided `K` was < 2!")
  return OrderedNormalMeanPrior([comp1; fill(compk, K - 1)])
end

isreal(p::OrderedNormalMeanPrior) = p.priors[1] isa Normal
ncomponents(p::OrderedNormalMeanPrior) = length(p.priors)

get_untruncated_mean(d::Normal) = d.μ
get_untruncated_mean(d::Truncated{<:Normal}) = get_untruncated_mean(d.untruncated)
get_untruncated_std(d::Normal) = d.σ
get_untruncated_std(d::Truncated{<:Normal}) = get_untruncated_std(d.untruncated)

function get_untruncated_mean_std(d::Union{Normal,Truncated{<:Normal}})
  return (get_untruncated_mean(d), get_untruncated_std(d))
end

function update(p::OrderedNormalMeanPrior, delta::AbstractVector{<:Real},
                y::AbstractVector{<:Real}, sigmasq::AbstractVector{<:Real},
                lambda::AbstractVector{<:Integer})
  next = copy(delta)
  for k in eachindex(next)
    next[k] = update(p, k, next, y, sigmasq, lambda)
  end
  return next
end

function update(p::OrderedNormalMeanPrior, k::Int, delta::AbstractVector{<:Real},
                y::AbstractVector{<:Real}, sigmasq::AbstractVector{<:Real},
                lambda::AbstractVector{<:Integer})
  # New value.
  next = copy(delta)

  # Set of indices such that lambda[i] > k
  Ck = findall(lambda .>= k)

  # Prior (mean, sd) of untruncated Normal
  mk, sk = get_untruncated_mean_std(p.priors[k])

  # Compute new variance for untruncated Normal
  sigmasq_k = sigmasq[Ck]
  vnew = 1 / (1/sk^2 +  sum(1 ./ sigmasq_k))

  # Compute new mean for untruncated Normal
  next[k] = 0
  mu = to_mu(next)
  lam = lambda[Ck]
  mnew = vnew * (mk/sk^2 + sum((y[Ck] - mu[lam]) ./ sigmasq_k))

  d_untruncated = Normal(mnew, sqrt(vnew))

  if p.priors[k] isa Normal
    d = d_untruncated
  else
    d = truncated(d_untruncated, p.priors[k].lower, p.priors[k].upper)
  end

  return rand(d)
end

to_mu(delta) = cumsum(delta)
to_delta(mu) = [mu[1]; diff(mu)]

function update_mu(p::OrderedNormalMeanPrior, mu::AbstractVector{<:Real},
                   y::AbstractVector{<:Real}, sigmasq::AbstractVector{<:Real},
                   lambda::AbstractVector{<:Integer})
  delta = to_delta(mu)
  new_delta = update(p, delta, y, sigmasq, lambda)
  return to_mu(new_delta)
end

"""
Make ordered prior for mu (implied by delta).
"""
function make_ordered_prior(y, K; s=nothing, upper=.9, lower=.1)
  s === nothing && (s = std(y) / K)
  comp1_mean = quantile(y, lower)
  compk_mean = (quantile(y, upper) - quantile(y, lower)) / (K - 1)
  return OrderedNormalMeanPrior(K,
                                Normal(comp1_mean, s),
                                truncated(Normal(compk_mean, s), 0, Inf))
end


text:
@testset "Testing Concat" begin
    # x, y, z = randn.(Float32, (10, 20, 30), 10)
    # X, Y, Z = OneDNN.Memory.((x, y, z))

    # # vcat
    # expected = vcat(x, y, z)
    # arg = [X, Y, Z]
    # op = OneDNN.Concat(arg, 1)
    # result = op(arg)
    # @test isa(result, OneDNN.Memory)
    # @test result == expected

    # # hcat
    # x, y, z = randn.(Float32, 10, (10, 20, 30))
    # X, Y, Z = OneDNN.Memory.((x, y, z))

    # # vcat
    # expected = hcat(x, y, z)
    # arg = [X, Y, Z]
    # op = OneDNN.Concat(arg, 2)
    # result = op(arg)
    # @test result == expected

    #####
    ##### Pullbacks
    #####

    function tester(f, fargs, g, gargs; fkw...)
        expected, back_expected = Zygote._pullback(f, fargs...; fkw...)

        Δ = Float32(0.125) * randn(Float32, size(expected))
        grads_expected = back_expected(Δ)

        result, back = Zygote._pullback(g, gargs...)
        @test isa(result, OneDNN.Memory)
        @test OneDNN.materialize(result) == expected

        # two flavors, one with a `Memory` passed in, and one with just a normal julia
        # array.
        #
        # Note that the structures of the gradients are different because `vcat` is variadic,
        # and thus gets a tuple entry for each argument.
        #
        # OneDNN.concat, hoever, takes an array, so the sensitivity to the array input comes
        # back as an array.
        grads_result = back(OneDNN.Memory(Δ))
        @test grads_result[2][1] == grads_expected[2]
        @test grads_result[2][2] == grads_expected[3]
        @test grads_result[2][3] == grads_expected[4]
        @test grads_result[3] === nothing

        grads_result = back(Δ)
        @test grads_result[2][1] == grads_expected[2]
        @test grads_result[2][2] == grads_expected[3]
        @test grads_result[2][3] == grads_expected[4]
        @test grads_result[3] === nothing

        # Test inference.
        @inferred back(OneDNN.Memory(Δ))
        @inferred back(Δ)
    end

    # vcat
    x, y, z = randn.(Float32, (10, 20, 30), 10)
    X, Y, Z = OneDNN.Memory.((x, y, z))
    tester(vcat, (x, y, z), OneDNN.concat, ([X, Y, Z], 1))
    tester(vcat, (x, y, z), OneDNN.concat, ((X, Y, Z), 1))

    # hcat
    x, y, z = randn.(Float32, 10, (10, 20, 30))
    X, Y, Z = OneDNN.Memory.((x, y, z))
    tester(hcat, (x, y, z), OneDNN.concat, ([X, Y, Z], 2))
    tester(hcat, (x, y, z), OneDNN.concat, ((X, Y, Z), 2))

    ### 3d

    # vcat
    x, y, z = randn.(Float32, (10, 20, 30), 10, 20)
    X, Y, Z = OneDNN.Memory.((x, y, z))
    tester(vcat, (x, y, z), OneDNN.concat, ([X, Y, Z], 1))
    tester(vcat, (x, y, z), OneDNN.concat, ((X, Y, Z), 1))

    # hcat
    x, y, z = randn.(Float32, 10, (10, 20, 30), 20)
    X, Y, Z = OneDNN.Memory.((x, y, z))
    tester(hcat, (x, y, z), OneDNN.concat, ([X, Y, Z], 2))
    tester(hcat, (x, y, z), OneDNN.concat, ((X, Y, Z), 2))
end


text:
# Note. `ModelProxy` is the type of a model's metadata entry (a named
# tuple). So, `info("PCA")` has this type, for example.

## BASIC IDEA

if false

    matching(model::MLJModels.ModelProxy, X) =
        !(model.is_supervised) && scitype(X) <: model.input_scitype

    matching(model::MLJModels.ModelProxy, X, y) =
        model.is_supervised &&
        scitype(X) <: model.input_scitype &&
        scitype(y) <: model.target_scitype

    matching(model::MLJModels.ModelProxy, X, y, w::AbstractVector{<:Real}) =
        model.is_supervised &&
        model.supports_weights &&
        scitype(X) <: model.input_scitype &&
        scitype(y) <: model.target_scitype

end

## IMPLEMENTATION

matching(X)       = Checker{false,false,scitype(X),missing}()
matching(X, y)    = Checker{true,false,scitype(X),scitype(y)}()
matching(X, y, w) = Checker{true,true,scitype(X),scitype(y)}()

(f::Checker{false,false,XS,missing})(model::MLJModels.ModelProxy) where XS =
    !(model.is_supervised) &&
    XS <: model.input_scitype

(f::Checker{true,false,XS,yS})(model::MLJModels.ModelProxy) where {XS,yS} =
    model.is_supervised &&
    XS <: model.input_scitype &&
    yS <: model.target_scitype

(f::Checker{true,true,XS,yS})(model::MLJModels.ModelProxy) where {XS,yS} =
    model.is_supervised &&
    model.supports_weights &&
    XS <: model.input_scitype &&
    yS <: model.target_scitype

(f::Checker)(name::String; pkg=nothing) = f(info(name, pkg=pkg))
(f::Checker)(realmodel::Model) = f(info(realmodel))

matching(model::MLJModels.ModelProxy, args...) = matching(args...)(model)
matching(name::String, args...; pkg=nothing) =
    matching(info(name, pkg=pkg), args...)
matching(realmodel::Model, args...) = matching(info(realmodel), args...)

## DUAL NOTION

struct DataChecker
    model::MLJModels.ModelProxy
end

matching(model::MLJModels.ModelProxy) = DataChecker(model)
matching(name::String; pkg=nothing) = DataChecker(info(name, pkg=pkg))
matching(realmodel::Model) = matching(info(realmodel))

(f::DataChecker)(args...) = matching(f.model, args...)


text:
# This is pretty naive and probably needs to be thought about a bit more, esp. in the periodic context.
# Right now this is sufficient to reproduce uniform fields for isolated systems.

@doc raw"""
Magnetic term ``A⋅(-i∇)``. It is assumed (but not checked) that ``∇⋅A = 0``.
"""
struct Magnetic
    Afunction::Function  # A(x,y,z) returns [Ax,Ay,Az]
                         # both [x,y,z] and [Ax,Ay,Az] are in *cartesian* coordinates
end
(M::Magnetic)(basis) = TermMagnetic(basis, M.Afunction)

struct TermMagnetic <: Term
    # Apotential[α] is an array of size fft_size for α=1:3
    Apotential::Vector{<:AbstractArray}
end
function TermMagnetic(basis::PlaneWaveBasis{T}, Afunction::Function) where T
    Apotential = [zeros(T, basis.fft_size) for α = 1:3]
    N1, N2, N3 = basis.fft_size
    rvecs = collect(r_vectors_cart(basis))
    for i = 1:N1
        for j = 1:N2
            for k = 1:N3
                Apotential[1][i, j, k],
                Apotential[2][i, j, k],
                Apotential[3][i, j, k] = Afunction(rvecs[i, j, k])
            end
        end
    end
    TermMagnetic(Apotential)
end

function ene_ops(term::TermMagnetic, basis::PlaneWaveBasis{T}, ψ, occ; kwargs...) where {T}
    ops = [MagneticFieldOperator(basis, kpoint, term.Apotential)
           for (ik, kpoint) in enumerate(basis.kpoints)]
    isnothing(ψ) && return (E=T(Inf), ops=ops)

    E = zero(T)
    for (ik, k) in enumerate(basis.kpoints)
        for iband = 1:size(ψ[1], 2)
            ψnk = @views ψ[ik][:, iband]
            # TODO optimize this
            E += basis.kweights[ik] * occ[ik][iband] * real(dot(ψnk, ops[ik] * ψnk))
        end
    end
    E = mpi_sum(E, basis.comm_kpts)

    (E=E, ops=ops)
end


text:
include("src/GEE.jl")

using .GEE
using GLM
using RCall

R"""
library(gee)
data(warpbreaks)

fit <- summary(gee(breaks ~ tension, id=wool, data=warpbreaks, corstr="exchangeable"))

library(geeM)
summary(geem(breaks ~ tension, id=wool, data=warpbreaks, corstr="exchangeable", useP = FALSE))
"""

warpbreaks = rcopy(R"warpbreaks")

fitted = GEE.initialfit(
    GEE.GeneralizedEstimatingEquation,
    @formula(breaks ~ tension),
    warpbreaks,
    Normal(),
    IdentityLink()
)

id = warpbreaks.wool
Yᵢ = [warpbreaks.breaks[id .== i, :] for i in unique(id)]
Xᵢ = [fitted.pp.X[id .== i, :] for i in unique(id)]
μ̂ᵢ = [predict(fitted)[id .== i] for i in unique(id)]
residᵢ = [warpbreaks.breaks[id .== i] .- predict(fitted)[id .== i] for i in unique(id)]

# this should be called Pearsons
e = GEE.StandardizedResiduals.(Normal(), Yᵢ, μ̂ᵢ)

N = length(id)
n = length(unique(id))
m = length(Yᵢ[1])
r = zeros(m, m)

phi = GEE.ϕ(e)
# phi = 141.1481 # when useP = TRUE in geeM

# this will create the unstructured correlation matrix
function alpha(n, m, e, phi)
    R = zeros(m, m)
    for i in 1:n
        for j in 1:m
            for k in 1:m
                if j == k
                    R[j, k] = 1
                else
                    R[j, k] += e[i][j]*e[i][k]
                end
            end
        end
    end
    R = R / (n * phi)
    R[diagind(R)] .= 1
    return R
end

unstr = alpha(n, m, e, phi)

# for i in 1:n
#     for j in 1:27
#         for k in 1:27
#             if j == k
#                 r[j, k] = 1
#             else
#                 r[j, k] += e[i][j]*e[i][k]
#             end
#         end
#     end
# end

# unstr = r / (n * phi)

function offdiag2(A::AbstractMatrix)
    [A[ι] for ι in CartesianIndices(A) if ι[1] ≠ ι[2]]
end

# we can calculate the exchangeable alpha from the unstructured correlation
a = mean(offdiag2(r / 2phi)) # I am successfully calculating exchangeable alpha when useP = FALSE in geeM

# R should theoretically be the correlation matrix, R()
# R = GEE.Independence(m)
R = GEE.Exchangeable(m, a)

Vᵢ = [GEE.V(R, mu) for mu in μ̂ᵢ]
Dᵢ = GEE.D.(Normal(), Xᵢ, μ̂ᵢ)

gee = GEE.GeeComponents(n, fitted.pp.beta0, Dᵢ, Vᵢ, residᵢ, 0, 0, 0)
GEE.updateβ(gee)

A = GEE.bread.(Dᵢ, Vᵢ)
B = GEE.meat.(Dᵢ, Vᵢ, residᵢ)
GEE.sandwich(A, B)
GEE.robustse(A, B)


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

# the following functions use the Plan object to create simulation data

export zeroarray, zerodict, zerodata, steadystatearray, steadystatedict, steadystatedata

"""
    zeroarray(model, plan)
    zeroarray(model, range)

Create a matrix of the proper dimension for a simulation with the given model
with the given plan or over the given range. If a range is given, the data is prepared for the
default plan. This means that appropriate number of periods are added before and
after the range to account for initial and final conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function zeroarray end
@inline zeroarray(m::Model, rng::AbstractUnitRange) = zeroarray(m, Plan(m, rng))
@inline zeroarray(m::Model, p::Plan) = inverse_transform(zeros(Float64, size(p.exogenous)), m)

"""
    zerodict(model, plan)
    zerodict(model, range)

Create a dictionary containing a [`TSeries`](@ref) of the appropriate range for
each variable in the model for a simulation with the given plan or over the
given range. If a range is given rather than a plan, the data is prepared for
the default plan over that range. This means that appropriate number of periods
are added before and after the range to account for initial and final
conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function zerodict end
@inline zerodict(m::Model, rng::AbstractUnitRange) = zerodict(m, Plan(m, rng))
function zerodict(m::Model, p::Plan)
    data = Dict{String,Any}()
    for v in m.varshks
        push!(data, string(v.name) => TSeries(p.range, inverse_transform(0.0, v)))
    end
    return data
end

"""
    zerodata(model, plan)
    zerodata(model, range)

Create a `NamedTuple` containing a [`TSeries`](@ref) of the appropriate range for each
variable in the model for a simulation with the given plan or over the given
range. If a range is given rather than a plan, the data is prepared for the
default plan over that range. This means that appropriate number of periods are
added before and after the range to account for initial and final conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function zerodata end
@inline zerodata(m::Model, rng::AbstractUnitRange) = zerodata(m, Plan(m, rng))
@inline zerodata(m::Model, p::Plan) = SimData(firstdate(p), m.varshks, zeroarray(m, p))

##################

"""
    steadystatearray(model, plan)
    steadystatearray(model, range)

Create a matrix of the proper dimensions for a simulation with the given model
with the given plan or over the given range. The matrix is initialized with the
steady state level of each variable. If a range is given rather than a plan, it
is augmented with periods before and after the given range in order to
accommodate initial and final conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function steadystatearray end
@inline steadystatearray(m::Model, rng::AbstractUnitRange; ref=first(rng)) = steadystatearray(m, Plan(m, rng), ref=ref)
function steadystatearray(m::Model, p::Plan; ref=firstdate(p) + m.maxlag)
    vs = keys(p.varshks)
    data = Matrix{Float64}(undef, length(p), length(vs))
    for (vi, v) ∈ enumerate(vs)
        data[:, vi] = m.sstate.:($v)[p.range, ref=ref]
    end
    return data
end

"""
    steadystatearray(model, plan)
    steadystatearray(model, range)

Create a dictionary containing a [`TSeries`](@ref) of the appropriate range for each
variable in the model for a simulation with the given plan or over the given
range. The matrix is initialized with the steady state level of each variable.
If a range is given rather than a plan, it is augmented with periods before and
after the given range in order to accommodate initial and final conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function steadystatedict end
@inline steadystatedict(m::Model, rng::AbstractUnitRange; ref=first(rng)) = steadystatedict(m, Plan(m, rng), ref=ref)
function steadystatedict(m::Model, p::Plan; ref=firstdate(p) + m.maxlag)
    data = Dict{String,Any}()
    for v ∈ keys(p.varshks)
        push!(data, string(v) => TSeries(p.range, m.sstate.:($v)[p.range, ref=ref]))
    end
    return data
end

"""
    steadystatedata(model, plan)
    steadystatedata(model, range)

Create a [`SimData`](@ref) containing a [`TSeries`](@ref) of the appropriate range for each
variable in the model for a simulation with the given plan or over the given
range. The matrix is initialized with the steady state level of each variable.
If a range is given rather than a plan, it is augmented with periods before and
after the given range in order to accommodate initial and final conditions.

See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
[`steadystatedict`](@ref)

"""
function steadystatedata end
@inline steadystatedata(m::Model, rng::AbstractUnitRange; ref=first(rng)) = steadystatedata(m, Plan(m, rng), ref=ref)
function steadystatedata(m::Model, p::Plan; ref=firstdate(p) + m.maxlag)
    return SimData(firstdate(p), m.varshks, steadystatearray(m, p, ref=ref))
end


text:
@testset "Angular momenta" begin
    @testset "Total angular momentum" begin
        𝐉 = TotalAngularMomentum()
        for ℓ = 0:10
            for j = ℓ .+ (-half(1):half(1))
                for m = -j:j
                    b = SpinOrbital(RelativisticOrbital(:m, ℓ, j), m)
                    for (±,∓) in [(+,-),(-,+)]
                        a = SpinOrbital(RelativisticOrbital(:k, ℓ, j), m)
                        # (V13.2.41,2)
                        @test dot(a, TensorComponent(𝐉, 0), b) ≈ m
                        # (V13.2.42)
                        @test dot(a, cartesian_tensor_component(𝐉, :z), b) ≈ m

                        m′ = m ± 1
                        abs(m′) > j && continue
                        a = SpinOrbital(RelativisticOrbital(:k, ℓ, j), m′)

                        # (V13.2.41)
                        @test dot(a, TensorComponent(𝐉, 0 ± 1), b) ≈ 0 ∓ √((j ± m + 1)*(j ∓ m)/2)
                        # (V13.2.42)
                        @test dot(a, cartesian_tensor_component(𝐉, :x), b) ≈ 0.5*√((j ± m + 1)*(j ∓ m))
                        @test dot(a, cartesian_tensor_component(𝐉, :y), b) ≈ 0 ∓ 0.5im*√((j ± m + 1)*(j ∓ m))
                    end
                end
            end
        end
    end
end


text:
testOnAllContexts("CrossEntropy") do ctx
    println("\tsolves a simple problem")
    testSimpleProblem(
        comparer ->
        begin
            ce = CrossEntropy(ctx, 60, comparer, dna -> dna)
            @test ce.dnaSize == 60
            @test ce.populationSize == 100
            @test ce.mutationChance == 0.05f0
            @test ce.mutationStrength == 0.02f0
            @test ce.selectionStdDev == 0.3f0
            @test ce.keepElitesRate == 0.05f0
            ce
        end)

    println("\tsolves a problem with synthesizer")
    testSimpleProblem(
        ctx,
        (comparer, synth) ->
        begin
            ce = CrossEntropy(
                ctx,
                dnaSize!(synth),
                comparer,
                dna -> decode(synth, dna),
                mutationChance = 0.05f0,
                selectionStdDev = 0.3f0)
            @test ce.dnaSize == dnaSize!(synth)
            @test ce.populationSize == 100
            @test ce.mutationChance == 0.05f0
            @test ce.mutationStrength == 0.02f0
            @test ce.selectionStdDev == 0.3f0
            @test ce.keepElitesRate == 0.05f0
            ce
        end)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module Combinatorics

export AbstractCombinatorics
export Combinations, DulCombinations, Permutations, DulPermutations

"""
    AbstractCombinatorics{M, C}

Abstract combinatorial algorithms.
"""
abstract type AbstractCombinatorics{M, C} end
@inline Base.eltype(::Type{<:AbstractCombinatorics{M, C}}) where {M, C} = NTuple{M, eltype(C)}
@inline @generated gettuple(contents, indices, ::Val{M}) where {M} = Expr(:tuple, [:(contents[indices[$i]]) for i = 1:M]...)

"""
    Combinations{M}(contents::C) where {M, C}

Combinations of M elements from contents. Duplicates are not allowed.
"""
struct Combinations{M, C} <: AbstractCombinatorics{M, C}
    contents::C
    N::Int
    Combinations{M}(contents::C) where {M, C} = new{M, C}(contents, length(contents))
end
@inline Base.length(c::Combinations{M}) where {M} = binomial(c.N, M)
Base.iterate(c::Combinations{M}) where {M} = (M > c.N) ? nothing : (M == 0) ? ((), [c.N+2]) : (gettuple(c.contents, 1:M, Val(M)), nextmstate!(collect(1:M), c.N, M))
Base.iterate(c::Combinations{M}, state) where M = (state[1] > c.N-M+1) ? nothing : (gettuple(c.contents, state, Val(M)), nextmstate!(state, c.N, M))
function nextmstate!(state::Vector{Int}, N::Int, M::Int)
    for i = M:-1:1
        state[i] += 1
        (state[i] > N-(M-i)) && continue
        for j = i+1:M
            state[j] = state[j-1] + 1
        end
        break
    end
    state
end

"""
    DulCombinations{M}(contents::C) where {M, C}

Combinations of M elements from contents. Duplicates are allowed.
"""
struct DulCombinations{M, C} <: AbstractCombinatorics{M, C}
    contents::C
    N::Int
    DulCombinations{M}(contents::C) where {M, C} = new{M, C}(contents, length(contents))
end
@inline Base.length(c::DulCombinations{M}) where {M} = binomial(c.N+M-1, c.N-1)
function Base.iterate(c::DulCombinations{M}) where M
    (M == 0) ? ((), [c.N+1]) : (gettuple(c.contents, ntuple(i->1, M), Val(M)), nextdmstate!(collect(ntuple(i->1, M)), c.N, M))
end
Base.iterate(c::DulCombinations{M}, state) where {M} = (state[1] > c.N) ? nothing : (gettuple(c.contents, state, Val(M)), nextdmstate!(state, c.N, M))
function nextdmstate!(state::Vector{Int}, N::Int, M::Int)
    for i = M:-1:1
        state[i] += 1
        (state[i] > N) && continue
        for j = i+1:M
            state[j] = state[j-1]
        end
        break
    end
    state
end

"""
    Permutations{M}(contents::C) where {M, C}

Permutations of M elements from contents. Duplicates are not allowed.
"""
struct Permutations{M, C} <: AbstractCombinatorics{M, C}
    contents::C
    N::Int
    Permutations{M}(contents::C) where {M, C} = new{M, C}(contents, length(contents))
end
@inline Base.length(p::Permutations{M}) where {M} = (0 <= M <= p.N) ? prod((p.N-M+1):p.N) : 0
function Base.iterate(p::Permutations{M}) where M
    ((p.N == 0) && (M > 0) || (0 < p.N < M)) ? nothing : (state = collect(1:p.N); (gettuple(p.contents, state, Val(M)), nextpstate!(state, p.N, M)))
end
function Base.iterate(p::Permutations{M}, state) where M
    ((p.N == 0) && (M > 0) || (0 < p.N < max(state[1], M))) ? nothing : (gettuple(p.contents, state, Val(M)), nextpstate!(state, p.N, M))
end
function nextpstate!(state::Vector{Int}, N::Int, M::Int)
    (M <= 0) && return [N+1]
    (M < N) && (j = M + 1; while (j <= N) && (state[M] >= state[j]); j += 1; end)
    if (M < N) && (j <= N)
        state[M], state[j] = state[j], state[M]
    else
        (M < N) && reverse!(state, M+1)
        i = M - 1
        while (i >= 1) && (state[i] >= state[i+1]); i -= 1; end
        if i > 0
            j = N
            while (j > i) && (state[i] >= state[j]); j -= 1; end
            state[i], state[j] = state[j], state[i]
            reverse!(state, i+1)
        else
            state[1] = N + 1
        end
    end
    return state
end

"""
    DulPermutations{M}(contents::C) where {M, C}

Permutations of M elements from contents. Duplicates are allowed.
"""
struct DulPermutations{M, C} <: AbstractCombinatorics{M, C}
    contents::C
    N::Int
    DulPermutations{M}(contents::C) where {M, C} = new{M, C}(contents, length(contents))
end
@inline Base.length(p::DulPermutations{M}) where M = p.N^M
function Base.iterate(p::DulPermutations{M}) where M
    indices = CartesianIndices(ntuple(i->p.N, M|>Val))
    index = iterate(indices)
    isnothing(index) && return nothing
    return gettuple(p.contents, reverse(index[1].I), M|>Val), (indices, index[2])
end
function Base.iterate(p::DulPermutations{M}, state) where M
    index = iterate(state[1], state[2])
    isnothing(index) && return nothing
    return gettuple(p.contents, reverse(index[1].I), M|>Val), (state[1], index[2])
end

end # module


text:
"""
    sequencingRRB(setting::RRBSetting, node::RestrictedSequencingNode, packet)

Call the correct rrb method based on the settings

# Arguments
- `setting::RRBSetting`: Which function to call
- `node::RestrictedSequencingNode`: The current node
- `packet`: Data with the necesary info
"""
function sequencingRRB(setting::RRBSetting, node::RestrictedSequencingNode, packet)
    if setting == sop
        return sopRRB(node.visited, getState(node), node.value, packet[1],packet[2])
    end
end


text:
@with_kw struct FastGrid
    resolution::Float64 = 1.0
    #tight::Bool         = false
end

# This is the main function
function solve(solver::FastGrid, problem::Problem) #original
    result = true
    delta = solver.resolution

    center = problem.input.center
    radius = problem.input.radius[1]

    (W, b) = (problem.network.layers[1].weights, problem.network.layers[1].bias)

    input = forward_affine_map(solver, W, b, problem.input)

    lower, upper = low(input), high(input)
    n_hypers_per_dim = BigInt.(max.(ceil.(Int, (upper-lower) / delta), 1))

    # preallocate work arrays
    local_lower, local_upper, CI = similar(lower), similar(lower), similar(lower)
    for i in 1:prod(n_hypers_per_dim)
        n = i
        for j in firstindex(CI):lastindex(CI)
            n, CI[j] = fldmod1(n, n_hypers_per_dim[j])
        end
        @. local_lower = lower + delta * (CI - 1)
        @. local_upper = min(local_lower + delta, upper)
        hyper = Hyperrectangle(low = local_lower, high = local_upper)

        k_1 = size(W, 1) #length(lower)
        #I = zeros(k_1, k_1) #Identity matrix
        #IN = zeros(k_1, k_1) #negative Identity matrix
        #for k = 1:k_1
            #I[k, k] = 1.0
            #IN[k, k] = -1.0
        #end
        #Hi
        #C_1 = vcat(I, IN)
        C_1 = vcat(Diagonal(ones(k_1)), -Diagonal(ones(k_1)))
        d_1 = vcat(local_upper, -local_lower)
        #Pi
        C = C_1 * W #size(2k_1, k_0)
        d = d_1 - C_1 * b #size(2k_1)

        inter = true #if the intersection of Pi and the input set is empty
        for l = 1:length(d)
            sum = 0.0
            for m = 1:size(C, 2)
                if C[l,m] > 0.0
                    sum += low(problem.input, m) * C[l,m]
                elseif C[l,m] < 0.0
                    sum += high(problem.input, m) * C[l,m]
                end
            end
            if sum > d[l]
                inter = false
                break
            end
        end

        if inter
            hull = false #if Hi is a hull of Z, or if Pi and center of the input set are close
            for l = 1:length(d)
                if distance(center, C[l,:], d[l]) > radius
                    hull = true
                    break
                end
            end

            if hull
                reach = forward_network(solver, problem.network, hyper)
                if !issubset(reach, problem.output)
                    result = false
                end
            end
        end
    end
    if result
        return BasicResult(:holds)
    end
    return BasicResult(:violated)
end

function forward_network(solver::FastGrid, nnet::Network, input::Hyperrectangle)
    layers = nnet.layers
    act = layers[1].activation
    reach = Hyperrectangle(low = act.(low(input)), high = act.(high(input)))

    for i in 2:length(layers)
        reach = forward_layer(solver, layers[i], reach)
    end
    return reach
end

# This function is called by forward_network
function forward_layer(solver::FastGrid, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_affine_map(solver::FastGrid, W::Matrix, b::Vector, input::Hyperrectangle)
    center = W * input.center + b
    radius = abs.(W) * input.radius
    return Hyperrectangle(center, radius)
end

function forward_node(solver::FastGrid, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    β    = node.act(output)  # TODO expert suggestion for variable name. beta? β? O? x?
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    return ((βmax + βmin)/2, (βmax - βmin)/2)
    #if solver.tight
        #return ((βmax + βmin)/2, (βmax - βmin)/2)
    #else
        #return (β, max(abs(βmax - β), abs(βmin - β)))
    #end
end

#distance of point to hyperPlane cx = d
function distance(point::Vector, c::Vector, d::Real)
    if length(point) == length(c)
        return abs(sum(point .* c) - d)/sqrt(sum(abs2.(c)))
    else
        error("Dimesion dismatch for point and constraint")
    end
end


text:
import TexasHoldem
const TH = TexasHoldem
using TexasHoldem
using BenchmarkTools
using Logging

struct BotCheckCall <: AbstractAI end

TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CheckRaiseFold) = check!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallRaiseFold) = call!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallAllInFold) = call!(game, player)
TH.player_option!(game::Game, player::Player{BotCheckCall}, ::CallFold) = call!(game, player)

players() = ntuple(i->(Player(BotCheckCall(), i)), 4)

@benchmark with_logger(NullLogger()) do
    play!(Game(players()))
end


text:
# ============================================================
# Type Hierarchy
# ============================================================

export AbstractSolverSet
export AbstractSetup, AbstractCondition
export AbstractProperty, AbstractGas, AbstractPlasma
export AbstractPhysicalSpace, AbstractPhysicalSpace1D, AbstractPhysicalSpace2D
export AbstractStructPhysicalSpace, AbstractUnstructPhysicalSpace
export AbstractVelocitySpace, AbstractVelocitySpace1D, AbstractVelocitySpace2D
export AbstractControlVolume, AbstractControlVolume1D, AbstractControlVolume2D
export AbstractUnstructControlVolume
export AbstractInterface, AbstractInterface1D, AbstractInterface2D
export AbstractSolution, AbstractSolution1D, AbstractSolution2D
export AbstractFlux, AbstractFlux1D, AbstractFlux2D
export AbstractParticle, AbstractParticle1D, AbstractParticle2D

export Setup
export Scalar, Radiation, Gas, DiatomicGas, Mixture, Plasma1D, Plasma2D
export IB, IB1F, IB2F, IB3F, IB4F
export ControlVolume1D, ControlVolume1D1F, ControlVolume1D2F, ControlVolume1D3F, ControlVolume1D4F
export ControlVolume2D, ControlVolume2D1F, ControlVolume2D2F, ControlVolume2D3F
export ControlVolumeUS, ControlVolumeUS1F, ControlVolumeUS2F
export Interface1D, Interface1D1F, Interface1D2F, Interface1D3F, Interface1D4F
export Interface2D, Interface2D1F, Interface2D2F
export Solution1D, Solution1D1F, Solution1D2F
export Solution2D, Solution2D1F, Solution2D2F
export Flux1D, Flux1D1F, Flux1D2F
export Flux2D, Flux2D1F, Flux2D2F
export Particle, Particle1D, Particle2D
export ControlVolumeParticle1D, ControlVolumeParticle2D

include("abstract.jl")
include("struct_general.jl")
include("struct_ib.jl")
include("struct_ctr.jl")
include("struct_face.jl")
include("struct_sol.jl")
include("struct_flux.jl")
include("struct_ptc.jl")

function copy_ctr!(
    ctr::T,
    ctr0::T,
) where {T<:Union{ControlVolume1D,ControlVolume2D,ControlVolumeUS}}
    ctr.w = ctr0.w
    ctr.prim .= ctr0.prim
    ctr.sw = ctr0.sw

    return nothing
end

function copy_ctr!(
    ctr::T,
    ctr0::T,
) where {T<:Union{ControlVolume1D1F,ControlVolume2D1F,ControlVolumeUS1F}}
    ctr.w .= ctr0.w
    ctr.prim .= ctr0.prim
    ctr.sw .= ctr0.sw
    ctr.f .= ctr0.f

    return nothing
end

function copy_ctr!(
    ctr::T,
    ctr0::T,
) where {T<:Union{ControlVolume1D2F,ControlVolume2D2F,ControlVolumeUS2F}}
    ctr.w .= ctr0.w
    ctr.prim .= ctr0.prim
    ctr.sw .= ctr0.sw
    ctr.h .= ctr0.h
    ctr.b .= ctr0.b

    return nothing
end

function copy_ctr!(
    ctr::T,
    ctr0::T,
) where {T<:Union{ControlVolume1D3F,ControlVolume2D3F}}
    ctr.w .= ctr0.w
    ctr.prim .= ctr0.prim
    ctr.sw .= ctr0.sw
    ctr.h0 .= ctr0.h0
    ctr.h1 .= ctr0.h1
    ctr.h2 .= ctr0.h2
    ctr.E .= ctr0.E
    ctr.B .= ctr0.B
    ctr.ϕ = ctr0.ϕ
    ctr.ψ = ctr0.ψ
    ctr.lorenz .= ctr0.lorenz

    return nothing
end

function copy_ctr!(
    ctr::T,
    ctr0::T,
) where {T<:ControlVolume1D4F}
    ctr.w .= ctr0.w
    ctr.prim .= ctr0.prim
    ctr.sw .= ctr0.sw
    ctr.h0 .= ctr0.h0
    ctr.h1 .= ctr0.h1
    ctr.h2 .= ctr0.h2
    ctr.h3 .= ctr0.h3
    ctr.E .= ctr0.E
    ctr.B .= ctr0.B
    ctr.ϕ = ctr0.ϕ
    ctr.ψ = ctr0.ψ
    ctr.lorenz .= ctr0.lorenz

    return nothing
end


text:
function priors_onegroup_alpha_delta_variant_cities(θ)
	@unpack β₀,β_school,β_home,β_other,β_work,ϵ,χ,p_test,E₀,inc_R_αvar,time_scale_αvar,mid_point_αvar,inc_R_δvar,time_scale_δvar,mid_point_δvar,init_scale= θ #,S₁,S₂,S₃,S₄,S₅,S₆

    LP = 0.
    LP += logpdf(Beta(50,50),ϵ)
    LP += logpdf(Gamma(3,4.5/3),χ)
    LP += logpdf(Gamma(3,5/3),p_test) ##units relative to 1% detection chance

	LP += logpdf(Gamma(3,1000/3),E₀)

    LP += logpdf(Gamma(10,1.5/10),β₀)
    LP += logpdf(Gamma(10,1.5/10),β_other)
    LP += logpdf(Gamma(10,1.5/10),β_work)
    LP += logpdf(Gamma(10,1.5/10),β_school)
    LP += logpdf(Gamma(10,1.5/10),β_home)

	LP += logpdf(Gamma(15,0.4/15),inc_R_αvar)
	LP += logpdf(Gamma(15,0.15/15),time_scale_αvar)
	LP += logpdf(Gamma(5,0.5/5),mid_point_αvar) #In units of 30 days after Feb 1st
	LP += logpdf(Gamma(15,0.6/15),inc_R_δvar)
	LP += logpdf(Gamma(15,0.2/15),time_scale_δvar)
	LP += logpdf(Gamma(5,0.5/5),mid_point_δvar)
	LP += logpdf(Gamma(5,1.5/5),init_scale)

    return LP
end

function priors_onegroup_alpha_delta_variant_outside_cities(θ)
	@unpack β₀,β_school,β_home,β_other,β_work,ϵ,χ,p_test,E₀,inc_R_αvar,time_scale_αvar,mid_point_αvar,inc_R_δvar,time_scale_δvar,mid_point_δvar,init_scale = θ
	LP = 0.
    LP += logpdf(Beta(50,50),ϵ)
    LP += logpdf(Gamma(3,4.5/3),χ)
    LP += logpdf(Gamma(3,2/3),p_test) #units relative to 1% detection chance

	LP += logpdf(Gamma(3,1000/3),E₀)

    LP += logpdf(Gamma(10,1.5/10),β₀)
    LP += logpdf(Gamma(10,1.5/10),β_other)
    LP += logpdf(Gamma(10,1.5/10),β_work)
    LP += logpdf(Gamma(10,1.5/10),β_school)
    LP += logpdf(Gamma(10,1.5/10),β_home)

	LP += logpdf(Gamma(15,0.4/15),inc_R_αvar)
	LP += logpdf(Gamma(15,0.15/15),time_scale_αvar)
	LP += logpdf(Gamma(5,0.5/5),mid_point_αvar) #In units of 30 days after Feb 1st
	LP += logpdf(Gamma(15,0.6/15),inc_R_δvar)
	LP += logpdf(Gamma(15,0.2/15),time_scale_δvar)
	LP += logpdf(Gamma(5,0.5/5),mid_point_δvar)
	LP += logpdf(Gamma(5,1.5/5),init_scale)

    return LP
end

function priors_onegroup_alpha_delta_variant_noncities_fitted(θ)
    @unpack β₀, β_school, β_home, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ

    LP = 0.0
    LP += logpdf(Beta(10, 10), ϵ)

    LP += logpdf(fitted_priors.d_χ, χ)  #
    LP += logpdf(fitted_priors.d_ptest, p_test) # units relative to 1% detection chance

    LP += logpdf(Gamma(3, 1000 / 3), E₀)

    LP += logpdf(Gamma(5, 1.5 / 5), β₀)
    LP += logpdf(Gamma(5, 1.5 / 5), β_other)
    LP += logpdf(Gamma(5, 1.5 / 5), β_work)
    LP += logpdf(Gamma(5, 1.5 / 5), β_school)
    LP += logpdf(Gamma(5, 1.5 / 5), β_home)

    LP += logpdf(Gamma(15, 0.4 / 15), inc_R_αvar)
    LP += logpdf(Gamma(15, 0.15 / 15), time_scale_αvar)
    LP += logpdf(Gamma(5, 0.5 / 5), mid_point_αvar) #In units of 30 days after Feb 1st
    LP += logpdf(Gamma(15, 0.6 / 15), inc_R_δvar)
    LP += logpdf(Gamma(15, 0.2 / 15), time_scale_δvar)
    LP += logpdf(Gamma(5, 0.5 / 5), mid_point_δvar)

	LP += logpdf(Gamma(3, 1.2 / 3), init_scale)

    return LP
end


text:
module GeneticBitArrays

import Base.show,
       Base.length,
       Base.getindex,
       Base.setindex!,
       Base.==,
       Base.rand

import StatsBase.Weights,
       StatsBase.sample

abstract type GeneticSeq end

include("DNA.jl")
include("RNA.jl")
include("Common.jl")

export GeneticSeq, DNASeq, RNASeq, Weights

end # module


text:
using Test

using StaticArrays

using Polyhedra

include("solvers.jl")

@testset "Redundancy removal" begin
    @testset "LP-based" begin
        @testset "VRepresentation with $(typeof(x))" for (x, y, z) in [
            ([1, 0], [0, 1], [1, 1]),
            ([1.0, 0.0], [0.0, 1.0], [1.0, 1.0]),
            ((@SVector [1, 0]), (@SVector [0, 1]), (@SVector [1, 1]))]
            T = eltype(x)
            AT = typeof(x)
            d = Polyhedra.FullDim(x)
            @testset "VLinearSpace" begin
                vr = vrep(typeof(Line(x))[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa LinesHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(Line(x))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa LinesHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 1
                vr = convexhull(Line(x), Line(y))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa LinesHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 2
                vr = convexhull(Line(x), Line(y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa LinesHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 2
                vrepm1 = Polyhedra.emptyspace(vr)
                @test vrepm1 === @inferred removevredundancy(vrepm1, lp_solver)
            end
            @testset "Polyhedra.RaysHull" begin
                vr = vrep(typeof(Line(x))[], typeof(Ray(x))[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(Ray(x))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test !haslines(rm)
                for vr in [convexhull(Ray(x), Ray(y), Ray(z)),
                           convexhull(Ray(x), Ray(z), Ray(y)),
                           convexhull(Ray(z), Ray(x), Ray(y))]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.RaysHull{T, AT}
                    @test npoints(rm) == 1
                    @test nrays(rm) == 2
                    rs = collect(rays(rm))
                    @test rs[1] == Ray(x)
                    @test rs[2] == Ray(y)
                    @test !haslines(rm)
                end
                vr = convexhull(Ray(x), Ray(y), Ray(-z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 3
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test nlines(rm) == 1
                vr = convexhull(Ray(x), Ray(-y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test nlines(rm) == 1
            end
            @testset "Polyhedra.PointsHull" begin
                vr = vrep(typeof(x)[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test !haslines(rm)
                for vr in [convexhull(2x, 2y, z),
                           convexhull(2x, z, 2y),
                           convexhull(z, 2x, 2y)]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.PointsHull{T, AT}
                    @test npoints(rm) == 2
                    ps = collect(points(rm))
                    @test ps[1] == 2x
                    @test ps[2] == 2y
                    @test !hasrays(rm)
                    @test !haslines(rm)
                end
                vr = convexhull(x, y, z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y, -z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test nlines(rm) == 1
                vr = convexhull(x, -y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test first(points(rm)) == -y
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, -y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 1
            end
        end
    end
    @testset "Representation-based" begin
        @testset "isredundant" begin
            @test isredundant(HyperPlane([0, 1], 0) ∩ HalfSpace([1, 1], 1), [0, 0], nl=0)
            @test !isredundant(HalfSpace([1, 1], 1) ∩ HyperPlane([0, 1], 0) ∩ HalfSpace([-1, -1], 1), [1, 0], nl=0)
            #@test !ishredundant(convexhull(SymPoint([1, 0])), HyperPlane([0, 1], 0), d=1)
            #@test isredundant(convexhull(SymPoint([1, 0])), HalfSpace([0, 1], 1), d=1)
            @test isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0)
            @test !isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0, strongly=true)
        end
#    @testset "Non-redundant SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1]))
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1) ∩ HalfSpace([-1, 1], 1)
#        vr = @inferred Polyhedra.removevredundancy(vr, hr)
#        @test vr isa Polyhedra.PointsHull{2,Int,Vector{Int}}
#        @test collect(sympoints(vr)) == [SymPoint([-1, 0]), SymPoint([0, 1])]
#        @test !haspoints(vr)
#        @test !hasallrays(vr)
#    end
        @testset "Point" begin
            # x - y ≤  1
            # x + y ≥ -1
            hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
            # [-1, 0] is weakly redundant as it belongs to the facet x + y = -1
            @test isredundant(hr, [-1, 0], nl=0)
            @test !isredundant(hr, [-1, 0], strongly=true, nl=0)
            # [0, -1] is not redundant as it is an extreme point
            @test !isredundant(hr, [0, -1], nl=0)
            @test !isredundant(hr, [0, -1], strongly=true, nl=0)
            # [0, 0] is strongly redundant as it is in the relative interior
            @test isredundant(hr, [0, 0], nl=0)
            @test isredundant(hr, [0, 0], strongly=true, nl=0)
            vr = convexhull([-1, 0], [0, -1], [0, 0]) + conichull([1, 1], [-1, 1])
            @test collect(removevredundancy(points(vr), hr, nl=0)) == [[0, -1]]
            @test collect(removevredundancy(points(vr), hr, strongly=true, nl=0)) == [[-1, 0], [0, -1]]
            vrr = Polyhedra.removevredundancy(vr, hr)
            @test collect(points(vrr)) == [[0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            vrr = Polyhedra.removevredundancy(vr, hr, strongly=true)
            @test collect(points(vrr)) == [[-1, 0], [0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            p = polyhedron(vr)
            Polyhedra.computehrep!(p)
            removevredundancy!(p, strongly=true)
            @test collect(points(p)) == [[-1, 0], [0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
            removevredundancy!(p)
            @test collect(points(p)) == [[0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
        end
#    @testset "Split SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1])) + conichull([1, 1], [-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([1, 1]), Ray([-1, 1])]
#
#        vr = convexhull(SymPoint([1, 0]), SymPoint([0, 1])) + conichull([-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[1, 0], [0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([-1, 1])]
#    end
    end
end
@testset "Duplicate removal" begin
    h = removeduplicates(hrep([1 1; -1 -1; 1 0; 0 -1], [1, -1, 1, 0]))
    @test nhyperplanes(h) == 1
    @test first(hyperplanes(h)) == HyperPlane([1, 1], 1)
    @test first(hyperplanes(h)) == HyperPlane([-1, -1], -1)
    @test nhalfspaces(h) == 1
    h = removeduplicates(hrep([1 1; -1 -1], [2, -1]))
    @test !hashyperplanes(h)
    @test nhalfspaces(h) == 2
#    v = removeduplicates(convexhull([1, 2], SymPoint([1, 2]), [-1, 1]))
#    @test typeof(v) == Polyhedra.PointsHull{2,Int,Vector{Int}}
#    @test collect(sympoints(v)) == [SymPoint([1, 2])]
#    @test collect(points(v)) == [[-1, 1]]
#    @test !hasallrays(v)
    for v in (#removeduplicates(convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1]) + Line([0, 1])),
              #removeduplicates(Line([0, 1]) + convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1])),
              removeduplicates(Line([0, 1]) + convexhull([-1, -2], [2, 1], [1, 2], [-1, 1], [-1, -2])),)
        @test typeof(v) == Polyhedra.Hull{Int, Vector{Int}, Int}
        #@test collect(sympoints(v)) == [SymPoint([1, 2])]
        @test collect(points(v)) == [[-1, -2], [2, 1], [1, 2]]
        @test !hasrays(v)
        @test collect(lines(v)) == [Line([0, 1])]
    end
end


text:
function new_material(name::String; use_nodes=false)
    mat = data.materials.new(name=name)
    mat.use_nodes = use_nodes
    mat
end

function connections(node)
    println("Node \"$(node.name)\"")
    println(node)
    indices = 1:max(length(node.inputs),length(node.outputs))
    inames = vcat([i.name for i in node.inputs], fill("", indices[end]-length(node.inputs)))
    onames = vcat([o.name for o in node.outputs], fill("", indices[end]-length(node.outputs)))
    pretty_table(hcat(indices, inames, onames),
                 header=["#", "Input", "Output"],
                 alignment=[:r,:l,:l],
                 vlines=[1], hlines=[1])
end

function find_input(o, name)
    i = findfirst(i -> o.inputs[i].name == name, 1:length(o.inputs))
    isnothing(i) && error("No such input $(name)")
    o.inputs[i]
end

function find_output(o, name)
    i = findfirst(i -> o.outputs[i].name == name, 1:length(o.outputs))
    isnothing(i) && error("No such output $(name)")
    o.outputs[i]
end

function set_input!(o, name, value)
    find_input(o, name).default_value = value
end

export new_material, connections,
    find_input, find_output, set_input!


text:
using LinearAlgebra
using DataDrivenDiffEq
using Plots
using ModelingToolkit

# Frist we define a set of variables and parameters
@variables u[1:3]
@parameters t w[1:2]

# Now the equations which form our basis
h = [u[1]; u[2]; cos(w[1]*u[2]+w[2]*u[3])]

# Then we simply create a basis
b= Basis(h, u, parameters = w)
ODESystem(b)

# Basis are callable with out of place
x = b([1;2;3])
# And in place
dx = similar(x)
b(dx, [1;2;3])
dx
# And look at the corresponding eqs
println(b)

# Suppose we want to add another equation, say sin(u[1])
# The basis behaves like an array
push!(b, sin(u[1]))
size(b) # (4)

# Adding an equation which is already present, does not change the basis
push!(b, sin(u[1]))
size(b) # Still 4

# We can iterate over the basis
for bi in b
    println(bi)
end

# Index specific eqs
b[3]

# And of course evaluate over trajectories
X = randn(3, 40)
Y_p = b(X)
# With parameter and time
t = independent_variable(b)
push!(b, sin(t))
Y = b(X, [2;4], 0:39)

# This allows you to transform a basis simply via
@variables x[1:2]
y = [sin(x[1]); cos(x[1]); x[2]]

b2 = Basis(b(y, parameters(b), independent_variable(b)), x, parameters = w, iv = t)
println(b2)
b2([1;2;3], [1 3], 0.0)

# We can merge basis
b3 = merge(b, b2)

# Also in place
merge!(b3, b2)
println(b3)

# Get the variables or parameters
variables(b)
parameters(b)
independent_variable(b)

# We can also check if two bases are equal
b == b

# Every function which is defined over Operations can be parsed
f(u, p, t) = [u[3]; u[2]*u[1]; sin(u[1])*u[2]]
b = Basis(f, u)

using Flux
NNlib.σ(x::Operation) = 1 / (1+exp(-x))
# Build a fully
c = Chain(Dense(3,5,σ), Dense(5, 2, σ))
ps, re = Flux.destructure(c)
@parameters p[1:length(ps)]
g(u, p, t) = re(p)(u)
b = Basis(g, u, parameters = p)
println(b)


text:
using Test
include("module-3-exercise-4.jl")

@testset "$(get_module_3_excersise_4_author_name()): problem 4 tests" begin
    @test get_module_3_excersise_4_author_name() != "YOUR NAME HERE"
    @test problem_4(3) == ([0.9999999999999997, 1.9999999999999978, 2.0000000000000036, 3.9999999999999982, 4.0, 4.999999999999997], [-0.40824829046386313 -0.2886751345948132 0.5000000000000001 0.28867513459481314 0.4999999999999999 -0.40824829046386263; 0.40824829046386246 -0.28867513459481337 -0.5000000000000004 -0.2886751345948126 0.49999999999999967 -0.40824829046386346; -0.4082482904638625 0.5773502691896258 3.88578058618805e-16 -0.5773502691896252 -4.996003610813204e-16 -0.4082482904638641; 0.4082482904638628 -0.2886751345948139 0.4999999999999996 -0.28867513459481203 -0.5000000000000006 -0.408248290463863; -0.4082482904638637 -0.2886751345948127 -0.5000000000000001 0.28867513459481364 -0.49999999999999983 -0.4082482904638618; 0.4082482904638635 0.5773502691896248 0.0 0.5773502691896264 0.0 -0.4082482904638631])
    @test problem_4(10) == ([8.000000000000004, 9.0, 9.000000000000005, 11.0, 11.000000000000002, 11.999999999999998], [0.40824829046386224 -0.2886751345948138 0.49999999999999944 0.5000000000000004 0.2886751345948123 -0.4082482904638637; -0.40824829046386124 -0.288675134594812 -0.5000000000000021 0.4999999999999977 -0.28867513459481375 -0.40824829046386424; 0.4082482904638627 0.5773502691896262 2.1649348980190608e-15 -2.2759572004815654e-15 -0.5773502691896251 -0.4082482904638638; -0.408248290463863 -0.2886751345948146 0.49999999999999895 -0.500000000000001 -0.2886751345948111 -0.40824829046386235; 0.4082482904638647 -0.28867513459481187 -0.49999999999999944 -0.5000000000000004 0.28867513459481414 -0.4082482904638609; -0.408248290463864 0.5773502691896247 0.0 0.0 0.5773502691896264 -0.408248290463863])
end


text:
module SurfaceCouplingTests

using Test
using Gridap
import Gridap: ∇
using LinearAlgebra: tr, ⋅

# Analytical functions

u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
∇u(x) = TensorValue( 2*x[1], 4*x[2], -2*x[1], zero(x[1]) )
Δu(x) = VectorValue( 6, -2 )

p(x) = x[1] + 3*x[2]
∇p(x) = VectorValue(1,3)

s(x) = -Δu(x)
f(x) = -Δu(x) + ∇p(x)
g(x) = tr(∇u(x))

∇(::typeof(u)) = ∇u
∇(::typeof(p)) = ∇p

# Geometry + Integration

n = 20
mesh = (n,n)
domain = 2 .* (0,1,0,1) .- 1
order = 1
model = CartesianDiscreteModel(domain, mesh)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[6,7,8])

Ω = Triangulation(model)

const R = 0.4

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

cell_to_coords = get_cell_coordinates(Ω)
cell_to_is_solid = lazy_map(is_in,cell_to_coords)
cell_to_is_fluid = lazy_map(!,cell_to_is_solid)

Ωs = Triangulation(model,cell_to_is_solid)
Ωf = Triangulation(model,cell_to_is_fluid)

Λ = BoundaryTriangulation(model,labels,tags="neumann")
Γ = InterfaceTriangulation(Ωf,Ωs)

n_Λ = get_normal_vector(Λ)
n_Γ = get_normal_vector(Γ)

order = 2
degree = 2*order

dΩ = Measure(Ω,degree)
dΩs = Measure(Ωs,degree)
dΩf = Measure(Ωf,degree)
dΛ = Measure(Λ,degree)
dΓ = Measure(Γ,degree)

# FE Spaces

reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P)

V = TestFESpace(Ω,reffe_u,conformity=:H1,labels=labels,dirichlet_tags="dirichlet")
Q = TestFESpace(Ωf,reffe_p,conformity=:L2)
U = TrialFESpace(V,u)
P = Q

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

#uh, ph = FEFunction(X,rand(num_free_dofs(X)))
#vh, qh = FEFunction(Y,rand(num_free_dofs(Y)))
#writevtk(Ω,"trian",cellfields=["uh"=>uh,"ph"=>ph,"vh"=>vh,"qh"=>qh])

# Weak form

a((u,p),(v,q)) =
  ∫( ∇(v)⊙∇(u) )*dΩs +
  ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩf

l((v,q)) =
  ∫( v⋅s )*dΩs +
  ∫( v⋅f + q*g )*dΩf +
  ∫( v⋅(n_Λ⋅∇u) - (n_Λ⋅v)*p )*dΛ +
  ∫( - (n_Γ.⁺⋅v.⁺)*p )*dΓ

# FE problem

op = AffineFEOperator(a,l,X,Y)
uh, ph = solve(op)

# Visualization

eu = u - uh
ep = p - ph

#writevtk(Ω,"trian_Ω",cellfields=["uh"=>uh,"ph"=>ph,"eu"=>eu,"ep"=>ep])
#writevtk(Ωs,"trian_Ωs",cellfields=["uh"=>uh,"ph"=>ph,"eu"=>eu,"ep"=>ep])
#writevtk(Γ,"trian_Γ",cellfields=["uh+"=>uh.⁺,"p"=>p,"n+"=>n_Γ.⁺])
#writevtk(Λ,"trian_Λ",cellfields=["uh"=>uh,"ph"=>ph,"eu"=>eu,"ep"=>ep,"n"=>n_Λ])
#writevtk(Ωf,"trian_Ωf",cellfields=["uh"=>uh,"ph"=>ph,"eu"=>eu,"ep"=>ep])

# Errors

eu_l2 = sqrt(sum(∫( eu⋅eu )*dΩ))
eu_h1 = sqrt(sum(∫( eu⋅eu + ∇(eu)⊙∇(eu) )*dΩ))
ep_l2 = sqrt(sum(∫( ep*ep )*dΩf))

tol = 1.0e-9
@test eu_l2 < tol
@test eu_h1 < tol
@test ep_l2 < tol

end # module


text:
"""
    TransformedKernel(k::Kernel,t::Transform)

Return a kernel where inputs are pretransformed by `t` : `k(t(x),t(x'))`
Can also be called via [transform](@ref) : `transform(k, t)`
"""
struct TransformedKernel{Tk<:Kernel,Tr<:Transform} <: Kernel
    kernel::Tk
    transform::Tr
end

(k::TransformedKernel)(x, y) = k.kernel(k.transform(x), k.transform(y))

# Optimizations for scale transforms of simple kernels to save allocations:
# Instead of a multiplying every element of the inputs before evaluating the metric,
# we perform a scalar multiplcation of the distance of the original inputs, if possible.
function (k::TransformedKernel{<:SimpleKernel,<:ScaleTransform})(
    x::AbstractVector{<:Real},
    y::AbstractVector{<:Real},
)
    return kappa(k.kernel, _scale(k.transform, metric(k.kernel), x, y))
end

function _scale(t::ScaleTransform, metric::Euclidean, x, y)
    return first(t.s) * evaluate(metric, x, y)
end
function _scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y)
    return first(t.s)^2 * evaluate(metric, x, y)
end
_scale(t::ScaleTransform, metric, x, y) = evaluate(metric, t(x), t(y))

"""
```julia
    transform(k::Kernel, t::Transform) (1)
    transform(k::Kernel, ρ::Real) (2)
    transform(k::Kernel, ρ::AbstractVector) (3)
```
(1) Create a TransformedKernel with transform `t` and kernel `k`
(2) Same as (1) with a `ScaleTransform` with scale `ρ`
(3) Same as (1) with an `ARDTransform` with scales `ρ`
"""
transform

transform(k::Kernel, t::Transform) = TransformedKernel(k, t)

function transform(k::TransformedKernel, t::Transform)
    return TransformedKernel(k.kernel, t ∘ k.transform)
end

transform(k::Kernel, ρ::Real) = transform(k, ScaleTransform(ρ))

transform(k::Kernel, ρ::AbstractVector) = transform(k, ARDTransform(ρ))

kernel(κ) = κ.kernel

Base.show(io::IO, κ::TransformedKernel) = printshifted(io, κ, 0)

function printshifted(io::IO, κ::TransformedKernel, shift::Int)
    printshifted(io, κ.kernel, shift)
    print(io,"\n" * ("\t" ^ (shift + 1)) * "- $(κ.transform)")
end

# Kernel matrix operations

function kerneldiagmatrix!(K::AbstractVector, κ::TransformedKernel, x::AbstractVector)
    return kerneldiagmatrix!(K, κ.kernel, map(κ.transform, x))
end

function kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix!(K, kernel(κ), map(κ.transform, x))
end

function kernelmatrix!(
    K::AbstractMatrix,
    κ::TransformedKernel,
    x::AbstractVector,
    y::AbstractVector,
)
    return kernelmatrix!(K, kernel(κ), map(κ.transform, x), map(κ.transform, y))
end

function kerneldiagmatrix(κ::TransformedKernel, x::AbstractVector)
    return kerneldiagmatrix(κ.kernel, map(κ.transform, x))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix(kernel(κ), map(κ.transform, x))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix(kernel(κ), map(κ.transform, x), map(κ.transform, y))
end


text:
module aITaskDispatcher

using DataFrames

include("aIGenerateDataSet.jl")

export aI_generate_00

end


text:
#=
   https://www.allendowney.com/blog/2018/10/21/the-game-of-ur-problem/
   """
   Here’s a probability puzzle to ruin your week.

   In the Royal Game of Ur, players advance tokens along a track with 14 spaces.
   To determine how many spaces to advance, a player rolls 4 dice with 4 sides. Two corners
   on each die are marked; the other two are not. The total number of marked corners —
   which is 0, 1, 2, 3, or 4 — is the number of spaces to advance.

   For example, if the total on your first roll is 2, you could advance a token to space 2.
   If you roll a 3 on the next roll, you could advance the same token to space 5.

   Suppose you have a token on space 13. How many rolls did it take to get there?
   """

   See:
   https://www.allendowney.com/blog/lions-and-tigers-and-bears/

   Allen Downey's solution:
   http://nbviewer.jupyter.org/github/AllenDowney/ThinkBayes2/blob/master/solutions/game_of_ur_soln.ipynb?flush=true

   cf ~/blog/game_of_ur_problem.blog
      ~/webppl/game_of_ur_problem.wppl

=#

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

@model function game_of_ur_problem(s=13)
    numRolls ~ DiscreteUniform(3,20)

    zroll = tzeros(numRolls)
    for i in 1:numRolls
        zroll[i] ~ DiscreteUniform(4)
    end

    sumRoll ~ DiscreteUniform(0,4*numRolls)
    sumRoll ~ Dirac(sum(zroll))

    true ~ Dirac(sumRoll == s)

end

model = game_of_ur_problem(13)

num_chns = 4

# chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)

# Got weird errors for MH()
#   - BoundsError: attempt to access 8-element Array{Int64,1} at index [2:20]
#
# chns = sample(model, MH(), 10_000)
# chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)

# chns = sample(model, PG(15), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, PG(20), 1_000)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, SMC(1000), 10_000)

# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(1000,0.65), 1_000)
# chns = sample(model, HMC(0.1,5), 1_000)
# chns = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)
# chns = sample(model, Gibbs(MH(:gender),NUTS(10,0.65,:height)), 1_000)
# chns = sample(model, Gibbs(MH(:gender),HMC(0.1,5,:height)), 1_000)
# chns = sample(model, Gibbs(PG(10,:gender),HMC(0.1,5,:height)), 1_000)
# chns = sample(model, Gibbs(MH(:gender),NUTS(1_000,0.65,:height)), 1_000)

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

show_var_dist_pct(chns, :numRolls,20)
println()
show_var_dist_pct(chns, :sumRoll,20)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Type definition
struct Constant{T} <: FunctionSpace{T}
    coeffs::Base.RefValue{T}
    Constant(x::T) where T = new{T}(Ref(x))
end

# Basic properties
domain(f::Constant) = (-Inf,+Inf)
period(f::Constant) = 0

coefficients(F::Constant) = F.coeffs

truncate!(F::Constant; kwargs...) = F

Base.:(==)(F1::Constant, F2::Constant) = (F1[] == F2[])

# Indexing, getting and setting coefficients
Base.eachindex(F::Constant) = 0:0
Base.getindex(F::Constant) = F.coeffs[]
Base.getindex(F::Constant, i::Integer) = iszero(i) ? getindex(F) : throw(BoundsError(F, i))
Base.setindex!(F::Constant, v) = F.coeffs[] = v
Base.getindex(F::Constant, v, i::Integer) =
    iszero(i) ? setindex!(F, v) : throw(BoundsError(F, i))

# Use as function
(F::Constant)(x) = F[]

# special purpose constructor
Base.similar(F::Constant, ::Type{T} = scalartype(F)) where {T} =
    Constant(zero(T)*F[])

Base.zero(F::Constant) = Constant(zero(F[]))
Base.one(F::Constant) = Constant(one(F[]))

# Arithmetic (out of place)
Base.copy(F::Constant) = Constant(copy(F[]))

Base.:-(F::Constant) = Constant(-F[])

Base.:*(F::Constant, a) = Constant(F[]*a)
Base.:*(a, F::Constant) = Constant(a*F[])
Base.:/(F::Constant, a) = Constant(F[]/a)
Base.:\(a, F::Constant) = Constant(a\F[])

Base.:+(F1::Constant, F2::Constant) = Constant(F1[]+F2[])
Base.:-(F1::Constant, F2::Constant) = Constant(F1[]-F2[])
Base.:*(F1::Constant, F2::Constant) = Constant(F1[]*F2[])
Base.:/(F1::Constant, F2::Constant) = Constant(F1[]/F2[])
Base.:\(F1::Constant, F2::Constant) = Constant(F1[]\F2[])

Base.conj(F::Constant) = Constant(conj(F[]))
Base.adjoint(F::Constant) = Constant(adjoint(F[]))
Base.transpose(F::Constant) = Constant(transpose(F[]))

LinearAlgebra.tr(F::Constant) = Constant(tr(F[]))

Base.real(F::Constant) = Constant(real(F[]))
Base.imag(F::Constant) = Constant(imag(F[]))

# Arithmetic (in place / mutating methods)
function Base.copy!(Fdst::Constant, Fsrc::Constant)
    if eltype(Fdst) <: Number
        Fdst[] = Fsrc[]
    else
        copy!(Fdst[], Fsrc[])
    end
    return Fdst
end
function LinearAlgebra.rmul!(F::Constant, α)
    if eltype(F) <: Number || !(α isa Number)
        F[] = F[]*α
    else
        rmul!(F[], α)
    end
    return F
end
function LinearAlgebra.lmul!(α, F::Constant)
    if eltype(F) <: Number || !(α isa Number)
        F[] = α*F[]
    else
        lmul!(α, F[])
    end
    return F
end

function LinearAlgebra.mul!(Fdst::Constant, α, Fsrc::Constant)
    if eltype(Fdst) <: Number
        Fdst[] = α*Fsrc[]
    else
        mul!(Fdst[], α, Fsrc[])
    end
    return Fdst
end
function LinearAlgebra.mul!(Fdst::Constant, Fsrc::Constant, α)
    if eltype(Fdst) <: Number
        Fdst[] = Fsrc[]*α
    else
        mul!(Fdst[], Fsrc[], α)
    end
    return Fdst
end
function LinearAlgebra.axpy!(α, Fx::Constant, Fy::Constant)
    if eltype(Fy) <: Number
        Fy[] = α*Fx[] + Fy[]
    else
        axpy!(α, Fx[], Fy[])
    end
    return Fy
end
function LinearAlgebra.axpby!(α, Fx::Constant, β, Fy::Constant)
    if eltype(Fy) <: Number
        Fy[] = α*Fx[] + β*F[y]
    else
        axpby!(α, Fx[], β, Fy[])
    end
    return Fy
end

LinearAlgebra.mul!(F::Constant, F1::Constant, F2::Constant,
                    α = true, β = false) = truncmul!(F, F1, F2, α, β)

function truncmul!(F::Constant, F1::Constant, F2::Constant,
                    α = true, β = false;
                    Kmax::Integer = 0, tol::Real = 0)

    if eltype(F) <: Number || (eltype(F1) <: Number && eltype(F2) <: Number)
        F[] = α * F1[] + β * F[2]
    else
        if eltype(F1) <: Number
            axpby!(α*F1[], F2[], β, F[])
        elseif eltype(F2) <: Number
            axpby!(α*F2[], F1[], β, F[])
        else
            mul!(F[], F1[], F2[], α, β)
        end
    end
    return F
end

# Inner product and norm
localdot(F1::Constant, F2::Constant) = Constant(dot(F1[], F2[]))

LinearAlgebra.dot(F1::Constant, F2::Constant) = dot(F1[], F2[])
LinearAlgebra.norm(F::Constant) = norm(F[])

function LinearAlgebra.isapprox(x::Constant, y::Constant;
                                atol::Real=0,
                                rtol::Real=defaulttol(x[0]))
    return norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))
end

differentiate(F::Constant) = zero(F)
integrate(F::Constant, (a,b)::Tuple{Real,Real}) = F[]*(b-a)

# Fit constant: take average over N points
fit(f, ::Type{Constant}, (a,b)::Tuple{Real,Real}; N = 5) =
    Constant(sum(f, range(a, b; length = N))/N)

# Inverse and square root
Base.inv(F::Constant) = Constant(inv(F[]))
Base.sqrt(F::Constant) = Constant(sqrt(F[]))


text:
struct GameState
    candidates::Vector{Candidate}
    trials::Vector{Trial}
    constraints::Constraints

    function GameState(use_wordle::Bool)
        candidates = candidate_pool(use_wordle)
        trials = Trial[]
        constraints = Constraints()
        new(candidates, trials, constraints)
    end
end

struct Config{T}
    weights::T
    filtering::Bool
    scoring::Bool
    use_wordle::Bool
    verbose::Bool
    debug::Bool
end

# TODO: Restore interactive mode -- interactive candidate, interactive scoring
# function prompt(text::String)
#     print(text)
#     input = String(chomp(readline()))
#     input
# end

# function parse_input()
#     candidate = Candidate(prompt("Candidate:\t"))
#     ints_string = prompt("Result Codes:\t")
#     codes = tuple(map(c -> ResultCode(parse(Int, c)), split(ints_string, ""))...)
#     trial = Trial(candidate, Result(codes))
# end

function print_state(state::GameState, candidate, iteration, config, n_candidates)
    if config.verbose
        @printf("Iteration: %s\n", iteration)
        @printf("Candidate: %s\n", String(candidate))
        pretty_print(state.trials)
    end

    if config.debug
        pretty_print(state.constraints)
        @printf(
            "Number of Candidates: %d / %d [%.3f%%]\n",
            length(state.candidates),
            n_candidates,
            100 * length(state.candidates) / n_candidates,
        )
        if length(state.candidates) <= 10
            @printf("Candidates: %s\n\n", map(String, state.candidates))
        end
    end
end

function play(
    truth_string::String;
    weights::Vector = fill(1 / 3, 3),
    filtering::Bool = true,
    scoring::Bool = true,
    use_wordle::Bool = true,
    verbose::Bool = false,
    debug::Bool = false,
)
    truth = Candidate(truth_string)
    config = Config(weights, filtering, scoring, use_wordle, verbose, debug)
    state = GameState(config.use_wordle)
    n_candidates = length(state.candidates)

    iteration = 1
    while iteration != 11
        i, candidate = select_candidate(state.candidates, state.trials, config.weights, config.scoring)

        trial = try_candidate(candidate, truth)

        deleteat!(state.candidates, i)

        push!(state.trials, trial)

        update_constraints!(state.constraints, trial)

        if filtering
            filter_candidates!(state.candidates, state.constraints)
        end

        print_state(state, candidate, iteration, config, n_candidates)

        if candidate == truth
            if config.verbose
                println("Success!")
            end
            break
        end

        if length(state.candidates) == 0
            error("Exhausted all candidates")
        end

        iteration += 1
    end

    if iteration == 11
        -1
    else
        iteration
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 0ce5327a-393b-11eb-1c38-4d48c8b2cdb6
begin
	using DrWatson

	function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
	end
end;

# ╔═╡ 18e233a2-393b-11eb-2a6e-e77e896a6ce3
begin
	@quickactivate

	using PlaymatSimulator

	md"""
	## GAME RULES (INACTIVE)

	Provide the game's rules (things to check for every frame) by defining custom types, functions, and objects below. Some basic rules funcionality is provided.
	"""
end

# ╔═╡ 2e3ae5ee-5019-11eb-0e4f-fd8ceaa8a8f7
begin
	abstract type AbstractDice end

	abstract type AbstractCounter end

	abstract type AbstractToken end

	abstract type AbstractCard end

	abstract type AbstractAbility end

	abstract type AbstractLand <: AbstractCard end

	abstract type AbstractSpell <: AbstractCard end

	abstract type TriggeredAbility <: AbstractAbility end

	abstract type ActivatedAbility <: AbstractAbility end
end

# ╔═╡ 62ea38ec-5019-11eb-0060-c5e3ae62838c
stack = Union{AbstractSpell,AbstractAbility}[]

# ╔═╡ 09703f54-4fcf-11eb-3a14-45a4cfafbe5d
mutable struct Card <: AbstractCard
	id::String
    name::String
	owner::String
	controller::String
	faces::Vector{Actor}
	tapped::Bool
	flipped::Bool
	scale::Vector{Float32}
	data::Dict{Symbol,Any}
end

# ╔═╡ 32ee705e-558c-11eb-2e89-77f8060cf4f6
mutable struct Ability <: AbstractAbility
	source::Card
	effect::Function
end

# ╔═╡ 357a5b80-558c-11eb-20bf-2188af834a12
mutable struct Spell <: AbstractSpell
	card::Card
	ability::Ability
end

# ╔═╡ 94d2e13c-55ec-11eb-0c93-91aa4e14ceff
mutable struct Dice <: AbstractDice
	id::String
	faces::Vector{Actor}
	sides::Int
end

# ╔═╡ 49b31c24-5012-11eb-2d79-9bfc1514c3aa
mutable struct Token <: AbstractToken
	id::String
	name::String
	owner::String
	controller::String
	faces::Vector{Actor}
	tapped::Bool
	flipped::Bool
	scale::Vector{Float32}
	data::Dict{Symbol,Any}
end

# ╔═╡ 7b5b091c-500d-11eb-35c5-45af2e6c4e29
mutable struct Zone
	owner::String
	name::String
	cards::Vector{Card}
end

# ╔═╡ d61fed94-5016-11eb-1b20-edc6a8b86bb5
mutable struct Planeswalker
	card::Card
	loyalty::Int
end

# ╔═╡ fbce7bf4-4fce-11eb-3205-b31b312ef4e5
mutable struct Player
    name::String
	portrait::Actor
    health::Int
    poison::Int
    priority::Bool
    resources::Dict{Symbol,Any}
    effects::Vector{Ability}
    zones::Vector{Zone}
end

# ╔═╡ 3819f224-558c-11eb-2498-95b826027aa2
mutable struct Counter <: AbstractCounter
	id::String
	faces::Vector{Actor}
	on::Union{Card,Player,Nothing}
end

# ╔═╡ 88d425c0-5011-11eb-11f9-271c42ea3c50
mutable struct PlaymatSimulatorGame
	player_names::Vector{Player}
	global_effects::Vector{Ability}
	server_address::String
	stack::Vector{Ability}
	gs::Dict
end

# ╔═╡ a938ca5c-393a-11eb-0557-d1580caa7d9d
gr = Dict(
    :turn_order => [:untap, :upkeep, :draw, :main_phase_one, :declare_attackers, :declare_blockers, :main_phase_two, :end_phase ],
	:win_conditions => [ :all_opponents_lost, :no_opponents_left ],
	:loss_conditions => [ :zero_health, :cant_draw, :poisoned, :commander_dmg ],
	)

# ╔═╡ 26a74f36-393b-11eb-296b-5dce05b5e22c
game_step = [
    "untap",
    "upkeep",
    "draw",
    "precombat_main",
    "combat",
    "postcombat_main",
    "end",
    "cleanup",
    ]

# ╔═╡ 96f9024e-3aa7-11eb-2916-c3f03b94b5e6
begin
	resource_types = ["{W}","{U}","{B}","{R}","{G}",:E,:X,:L,:E,:P]

	md"""
	Resource types:

	Mana
	:W - white
	:U - blue
	:B - black
	:R - red
	:G - green
	:C - colorless

	Other
	:L - life
	:E - energy
	:P - poison

	"""
end

# ╔═╡ 084fbf90-393c-11eb-0e3f-e1cf9b126e31
function declare_creature_attack(target::Union{Player, Planeswalker}, creatures::Dict, opponents::Dict, planeswalkers::Dict)
end

# ╔═╡ 712f3b76-393c-11eb-0ae4-f1c8a5eac20a
function declare_blockers(p::Player, creatures::Dict)
end

# ╔═╡ 877f03ac-393c-11eb-1f6f-3bfb1f3a008a
function shuffle_deck(p::Player)
end

# ╔═╡ 9147eeda-393c-11eb-04ee-e1587dfb35c4
function draw_card(p::Player, num::Int)
end

# ╔═╡ 973e3aa4-393c-11eb-08a4-59338f595c15
function discard_card(p::Player, amount::Int, random=false)
end

# ╔═╡ 9cdcbf34-393c-11eb-15f3-c716ca844b1b
function add_mana_to_pool(p::Player, mana_type)
end

# ╔═╡ 9ffdeb8c-393c-11eb-3900-8d09fa3ce1fb
function subtract_mana_from_pool(p::Player, mana_type)
end

# ╔═╡ a3426304-393c-11eb-047a-71ebc4f070fd
function tap_card(p::Player, card)
end

# ╔═╡ a73a4918-393c-11eb-0ecc-b58671c0c518
function untap_card(p::Player, card)
end

# ╔═╡ aa969f8a-393c-11eb-3f11-5f32374e939c
function play_land(p::Player)
end

# ╔═╡ ae6acabe-393c-11eb-235f-59a50435ecff
function play_spell(p::Player)
end

# ╔═╡ b13dba94-393c-11eb-21f4-ab7d7f18db74
function take_damage(p::Player, creature=nothing)
end

# ╔═╡ b3d88946-393c-11eb-2e79-9d3c657bb6d8
function heal_damage(p::Player, creature=nothing)
end

# ╔═╡ b85ee3d4-393c-11eb-3062-139a585a91a4
function gain_life(p::Player)
end

# ╔═╡ bb80fa16-393c-11eb-0f56-e713be77eebe
function lose_life(p::Player)
end

# ╔═╡ bebd3406-393c-11eb-1210-85373fd1d8ad
function end_phase(p::Player)
end

# ╔═╡ c1ff5a6a-393c-11eb-3bb0-6d4be308bbc6
function end_turn(p::Player)
end

# ╔═╡ c4f23038-393c-11eb-3bd8-b73d5c00f8f5
function run_upkeep(p::Player)
end

# ╔═╡ Cell order:
# ╟─0ce5327a-393b-11eb-1c38-4d48c8b2cdb6
# ╟─18e233a2-393b-11eb-2a6e-e77e896a6ce3
# ╠═2e3ae5ee-5019-11eb-0e4f-fd8ceaa8a8f7
# ╠═62ea38ec-5019-11eb-0060-c5e3ae62838c
# ╠═09703f54-4fcf-11eb-3a14-45a4cfafbe5d
# ╠═32ee705e-558c-11eb-2e89-77f8060cf4f6
# ╠═357a5b80-558c-11eb-20bf-2188af834a12
# ╠═3819f224-558c-11eb-2498-95b826027aa2
# ╠═94d2e13c-55ec-11eb-0c93-91aa4e14ceff
# ╠═49b31c24-5012-11eb-2d79-9bfc1514c3aa
# ╠═7b5b091c-500d-11eb-35c5-45af2e6c4e29
# ╠═d61fed94-5016-11eb-1b20-edc6a8b86bb5
# ╠═fbce7bf4-4fce-11eb-3205-b31b312ef4e5
# ╠═88d425c0-5011-11eb-11f9-271c42ea3c50
# ╠═a938ca5c-393a-11eb-0557-d1580caa7d9d
# ╠═26a74f36-393b-11eb-296b-5dce05b5e22c
# ╟─96f9024e-3aa7-11eb-2916-c3f03b94b5e6
# ╟─084fbf90-393c-11eb-0e3f-e1cf9b126e31
# ╟─712f3b76-393c-11eb-0ae4-f1c8a5eac20a
# ╟─877f03ac-393c-11eb-1f6f-3bfb1f3a008a
# ╟─9147eeda-393c-11eb-04ee-e1587dfb35c4
# ╟─973e3aa4-393c-11eb-08a4-59338f595c15
# ╟─9cdcbf34-393c-11eb-15f3-c716ca844b1b
# ╟─9ffdeb8c-393c-11eb-3900-8d09fa3ce1fb
# ╟─a3426304-393c-11eb-047a-71ebc4f070fd
# ╟─a73a4918-393c-11eb-0ecc-b58671c0c518
# ╟─aa969f8a-393c-11eb-3f11-5f32374e939c
# ╟─ae6acabe-393c-11eb-235f-59a50435ecff
# ╟─b13dba94-393c-11eb-21f4-ab7d7f18db74
# ╟─b3d88946-393c-11eb-2e79-9d3c657bb6d8
# ╟─b85ee3d4-393c-11eb-3062-139a585a91a4
# ╟─bb80fa16-393c-11eb-0f56-e713be77eebe
# ╟─bebd3406-393c-11eb-1210-85373fd1d8ad
# ╟─c1ff5a6a-393c-11eb-3bb0-6d4be308bbc6
# ╟─c4f23038-393c-11eb-3bd8-b73d5c00f8f5


text:
using SpectralMeasures, ValidatedNumerics, ValidatedNumerics.RootFinding,
            ComplexPhasePortrait

function tripleplot(a,b,Z=linspace(-3, 3, 300).+linspace(3,-3,300)'*im)
  # Build the measure and resolvents
  μ = spectralmeasure(a,b)
  R = principalresolvent(a,b)
  r = discresolvent(a,b)

  # Create the subplots
  p1 = plot(μ,xlims=(-2,2),ylims=(0,1.5))
  p2 = plot(portrait(R(Z),PTstepmod),xlims=(-3,3),ylims=(-3,3),aspect_ratio=1)
  plot!([-1.,1.],[0.,0.],linewidth=3,color=:black)
  p3 = plot(portrait(r(Z),PTstepmod),xlims=(-3,3),ylims=(-3,3),aspect_ratio=1)
  plot!(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=3,color=:black)

  # Plot the subplots
  l = @layout [a{0.7w}; b c]
  plot(p1,p2,p3,layout=l)
end

function validatedspectrum(a,b)
  # Chop the a and b down
  a = chop!(a); b = .5+chop!(b-.5)
  n = max(2,length(a),length(b)+1)
  a = [a;zeros(n-length(a))]; b = [b;.5+zeros(n-length(b))]

  # Finds C such that J*C = C*Toeplitz([0,1/2])
  C = SpectralMeasures.connectioncoeffsoperator(a,b)
  c = Fun(Taylor,C.T.nonnegative)

  rts = find_roots(x->c(x),-1,1)
  if length(rts) > 0
     eigs=real(map(x->joukowsky(x.interval),rts))
     eigserrs = map(midpoint_radius,eigs)
     spectrum = ApproxFun.Interval(-1,1)
     maxerr = 0
     for (eig,err) in eigserrs
       maxerr = max(maxerr,err)
       spectrum = ApproxFun.Point(eig) ∪ spectrum
     end
  else
    spectrum = ApproxFun.Interval(-1,1)
    maxerr = 0
  end
  spectrum, maxerr
end


text:
module NumericalRationalFunctions

using DocStringExtensions

include("aaa.jl")
include("lightning_laplace.jl")

export aaa, laplace

end


text:
"""
"""
module FeatureEng

include("numeric/binning.jl")
include("numeric/transforming.jl")
include("numeric/scaling.jl")
include("numeric/interactions.jl")
include("datetime.jl")
include("categorical.jl")

export
# General
    fit_transform!, apply_transform,

# Numeric -- Binning
    BinFixedWidth, BinExponential, BinQuantile,
# Numeric -- Scaling
    ScaleMinMax, ScaleVariance, ScaleL2,
# Numeric -- Transforming
    transformLog, transformRoot, transformBoxCox,
# Numeric -- Interaction terms
    polynomial,

# Categorical
    encode_onehot, encode_dummy, encode_hash,

# Datetime
    strp_datetimes, extract_datetime_features, extract_date_features, extract_time_features,
    get_month, get_weekday

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule XSLT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("XSLT")
JLLWrappers.@generate_main_file("XSLT", UUID("aed1982a-8fda-507f-9586-7b0439959a61"))
end  # module XSLT_jll