text:
module bingingstests

using GridapPardiso
using Test
using SparseArrays

# Define linear system

mtype = GridapPardiso.MTYPE_REAL_NON_SYMMETRIC

A = sparse([
  0. -2  3 0
  -2  4 -4 1
  -3  5  1 1
  1 -3 0 2 ])

n = A.n

b = Float64[1, 3, 2, 5]

x = zeros(Float64,n)

# Create the pardiso internal handler

pt = new_pardiso_handle()

# pardisoinit!

iparm = Vector{Int32}(new_iparm())

pardisoinit!(pt,mtype,iparm)

# pardiso! (solving the transpose of the system above)

maxfct = 1
mnum = 1
phase = GridapPardiso.PHASE_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT
a = A.nzval
ia = Vector{Int32}(A.colptr)
ja = Vector{Int32}(A.rowval)
perm = zeros(Int32,n)
nrhs = 1
msglvl = 0

err = pardiso!(
  pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x)

tol = 1.0e-13

@test err == 0

@test maximum(abs.(A'*x-b)) < tol

# pardiso! (solving the system above)

iparm[12] = 2

err = pardiso!(
  pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x)

@test err == 0

@test maximum(abs.(A*x-b)) < tol

@test pardiso_data_type(mtype,iparm) == Float64

# pardiso_getdiag!

iparm[56] = 1

pt = new_pardiso_handle()

err = pardiso!(
  pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x)

df = zeros(Float64,n)
da = zeros(Float64,n)

err = pardiso_getdiag!(pt,df,da,mnum,mtype,iparm)

@test err == 0

err = pardiso_getdiag!(pt,df,da,mnum)

@test err == 0

if Int == Int64
    # pardiso_64! (solving the transpose of the system above)
    pt = new_pardiso_handle()

    a = A.nzval
    ia = A.colptr
    ja = A.rowval
    perm = zeros(Int64,n)
    iparm = Vector{Int64}(new_iparm())

    err = pardiso_64!(
      pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x)

    @test err == 0

    # pardiso_64! (solving the transpose of the system above)

    @test maximum(abs.(A'*x-b)) < tol
end

end


text:
# TODO maybe move to MOI.Utilities alongside `MOIU.get_bounds`
function set_bounds(model::MOI.ModelLike, vi::MOI.VariableIndex, lower::T, upper::T) where T
    xval = vi.value
    c_lt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.LessThan{T}}(xval)
    c_gt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{T}}(xval)
    c_int = MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{T}}(xval)
    c_eq = MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{T}}(xval)
    if MOI.is_valid(model, c_int)
        MOI.set(model, MOI.ConstraintSet(), c_int, MOI.Interval(lower, upper))
        # It is assumed that none of the other ConstraintIndexs are valid
        return
    end
    if MOI.is_valid(model, c_eq)
        # Maybe delete it and add an interval constraint if this case is useful
        lower == upper || error("Cannot set different bounds $lower and $upper to a variables with `EqualTo` constraint.")
        MOI.set(model, MOI.ConstraintSet(), c_eq, MOI.EqualTo(lower))
        # It is assumed that none of the other `ConstraintIndex`s are valid
        return
    end
    lt = MOI.LessThan(upper)
    if MOI.is_valid(model, c_lt)
        if upper == typemin(upper)
            MOI.delete(model, c_lt)
        else
            MOI.set(model, MOI.ConstraintSet(), c_lt, lt)
        end
    elseif upper != typemax(upper)
        MOI.add_constraint(model, MOI.SingleVariable(vi), lt)
    end
    gt = MOI.GreaterThan(lower)
    if MOI.is_valid(model, c_gt)
        if lower == typemin(lower)
            MOI.delete(model, c_gt)
        else
            MOI.set(model, MOI.ConstraintSet(), c_gt, gt)
        end
    elseif lower != typemin(lower)
        MOI.add_constraint(model, MOI.SingleVariable(vi), gt)
    end
    return
end
function set_lower_bound(model::MOI.ModelLike, vi::MOI.VariableIndex, lower::T) where T
    xval = vi.value
    c_gt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{T}}(xval)
    c_int = MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{T}}(xval)
    c_eq = MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{T}}(xval)
    if MOI.is_valid(model, c_int)
        set = MOI.get(model, MOI.ConstraintSet(), c_int)
        MOI.set(model, MOI.ConstraintSet(), c_int, MOI.Interval(lower, set.upper))
        # It is assumed that none of the other ConstraintIndexs are valid
        return
    end
    if MOI.is_valid(model, c_eq)
        error("Cannot set lower bound to fixed variable.")
    end
    gt = MOI.GreaterThan(lower)
    if MOI.is_valid(model, c_gt)
        MOI.set(model, MOI.ConstraintSet(), c_gt, gt)
    else
        MOI.add_constraint(model, MOI.SingleVariable(vi), gt)
    end
    return
end
function set_upper_bound(model::MOI.ModelLike, vi::MOI.VariableIndex, upper::T) where T
    xval = vi.value
    c_lt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.LessThan{T}}(xval)
    c_int = MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{T}}(xval)
    c_eq = MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{T}}(xval)
    if MOI.is_valid(model, c_int)
        set = MOI.get(model, MOI.ConstraintSet(), c_int)
        MOI.set(model, MOI.ConstraintSet(), c_int, MOI.Interval(set.lower, upper))
        # It is assumed that none of the other ConstraintIndexs are valid
        return
    end
    if MOI.is_valid(model, c_eq)
        error("Cannot set lower bound to fixed variable.")
    end
    lt = MOI.LessThan(upper)
    if MOI.is_valid(model, c_lt)
        MOI.set(model, MOI.ConstraintSet(), c_lt, lt)
    else
        MOI.add_constraint(model, MOI.SingleVariable(vi), lt)
    end
    return
end

function generate_random_restart(m; cont=true)
    values = []
    for i=1:m.num_var
        lbi_def = true
        ubi_def = true
        if m.l_var[i] > typemin(Int64)
            lbi = m.l_var[i]
        else
            lbi = typemin(Int64)
            lbi_def = false
        end

        if m.u_var[i] < typemax(Int64)
            ubi = m.u_var[i]
        else
            ubi = typemin(Int64)
            ubi_def = false
        end

        if !ubi_def && !lbi_def
            ubi = 10
            lbi = -10
        elseif !ubi_def
            ubi = lbi+20
        elseif !lbi_def
            lbi = ubi-20
        end

        if m.var_type[i] == :Cont || cont
            push!(values,(ubi-lbi)*rand(JUNIPER_RNG;)+lbi)
        else
            lbi = Int(round(lbi))
            ubi = Int(round(ubi))
            push!(values, rand(JUNIPER_RNG,lbi:ubi))
        end
    end
    return values
end

"""
    get_reasonable_disc_vars(node, var_type, disc_vars, disc2var_idx, atol)

Get all discrete variables which aren't close to discrete yet based on atol
"""
function get_reasonable_disc_vars(node, var_type, disc_vars, disc2var_idx, atol)
    reasonable_disc_vars = zeros(Int64,0)
    for i=1:disc_vars
        idx = disc2var_idx[i]
        u_b = node.u_var[idx]
        l_b = node.l_var[idx]
        if isapprox(u_b,l_b; atol=atol) || is_type_correct(node.solution[idx],var_type[idx],atol)
            continue
        end
        push!(reasonable_disc_vars,i)
    end
    return reasonable_disc_vars
end

"""
    is_global_status(state::MOI.TerminationStatusCode)

Returns true if either ALMOST_OPTIMAL, OPTIMAL or INFEASIBLE and false otherwise
"""
function is_global_status(state::MOI.TerminationStatusCode)
    return state == MOI.ALMOST_OPTIMAL || state == MOI.OPTIMAL || state == MOI.INFEASIBLE
end

"""
    only_almost_solved(state::MOI.TerminationStatusCode)

Returns true if either ALMOST_OPTIMAL or ALMOST_LOCALLY_SOLVED
"""
function only_almost_solved(state::MOI.TerminationStatusCode)
    return state == MOI.ALMOST_OPTIMAL || state == MOI.ALMOST_LOCALLY_SOLVED
end

"""
    state_is_optimal(state::MOI.TerminationStatusCode; allow_almost=false)

Returns true if either optimal or locally solved. If allow_almost then check for `ALMOST_LOCALLY_SOLVED` and `ALMOST_OPTIMAL`
"""
function state_is_optimal(state::MOI.TerminationStatusCode; allow_almost=false)
    return state == MOI.OPTIMAL || state == MOI.LOCALLY_SOLVED ||
            (allow_almost && state == MOI.ALMOST_LOCALLY_SOLVED) || (allow_almost && state == MOI.ALMOST_OPTIMAL)
end

"""
    state_is_infeasible(state::MOI.TerminationStatusCode)

Returns true if either infeasible or locally infeasible
"""
function state_is_infeasible(state::MOI.TerminationStatusCode)
    return state == MOI.INFEASIBLE || state == MOI.LOCALLY_INFEASIBLE
end

struct ObjectiveConstraint{E<:MOI.AbstractNLPEvaluator} <: MOI.AbstractNLPEvaluator
    evaluator::E
    num_variables::Int
end
MOI.features_available(oc::ObjectiveConstraint) = MOI.features_available(oc.evaluator)
MOI.initialize(oc::ObjectiveConstraint, features) = MOI.initialize(oc.evaluator, features)
MOI.eval_objective(oc::ObjectiveConstraint, x) = MOI.eval_objective(oc.evaluator, x)
MOI.eval_objective_gradient(oc::ObjectiveConstraint, g, x) = MOI.eval_objective_gradient(oc.evaluator, g, c)
MOI.hessian_lagrangian_structure(oc::ObjectiveConstraint) = MOI.hessian_lagrangian_structure(oc.evaluator)
MOI.eval_hessian_lagrangian(oc::ObjectiveConstraint, H, x, σ, μ) = MOI.eval_hessian_lagrangian(oc.evaluator, H, x, σ, μ)
function MOI.eval_constraint(oc::ObjectiveConstraint, g, x)
    g[1] = MOI.eval_objective(oc.evaluator, x)
    MOI.eval_constraint(oc.evaluator, view(g, 2:length(g)), x)
    return
end
function MOI.jacobian_structure(oc::ObjectiveConstraint)
    sparsity = [(1, i) for i in 1:oc.num_variables]
    append!(sparsity, MOI.jacobian_structure(oc.evaluator))
    return sparsity
end
function MOI.eval_constraint_jacobian(oc::ObjectiveConstraint, J, x)
    MOI.eval_objective_gradient(oc.evaluator, view(J, 1:oc.num_variables), x)
    MOI.eval_constraint_jacobian(oc.evaluator, view(J, (oc.num_variables + 1):length(J)), x)
    return
end

"""
    add_obj_constraint(jp::JuniperProblem, rhs::Float64)

Add a constraint for the objective based on whether the objective is linear/quadratic or non linear.
If the objective sense is :MIN than add objective <= rhs else objective >= rhs
"""
function add_obj_constraint(jp::JuniperProblem, rhs::Float64)
    if jp.nlp_data.has_objective
        if jp.obj_sense == :Min
            lb, ub = -Inf, rhs
        else
            lb, ub = rhs, Inf
        end
        MOI.set(jp.model, MOI.NLPBlock(), MOI.NLPBlockData(
            [jp.nlp_data.constraint_bounds; MOI.NLPBoundsPair(lb, ub)],
            ObjectiveConstraint(jp.nlp_data.evaluator, MOI.get(jp.model, MOI.NumberOfVariables())),
            jp.nlp_data.has_objective,
        ))
    else # linear or quadratic
        if isa(jp.objective, MOI.SingleVariable)
            i = findfirst(isequal(jp.objective.variable), jp.x)
            vi = jp.objective.variable
            if jp.obj_sense == :Min
                set_upper_bound(jp.model, vi, rhs)
            else
                set_lower_bound(jp.model, vi, rhs)
            end
        else
            if jp.obj_sense == :Min
                MOI.add_constraint(jp.model, jp.objective, MOI.LessThan(rhs))
            else
                MOI.add_constraint(jp.model, jp.objective, MOI.GreaterThan(rhs))
            end
        end
    end
end

"""
    evaluate_objective(optimizer::MOI.AbstractOptimizer, jp::JuniperProblem, xs::Vector{Float64})

If no objective exists => return 0
Evaluate the objective whether it is non linear, linear or quadratic
"""
function evaluate_objective(optimizer::MOI.AbstractOptimizer, jp::JuniperProblem, xs::Vector{Float64})
    # NLP objective has priority if both are set
    if jp.nlp_data !== nothing && jp.nlp_data.has_objective
        return MOI.eval_objective(jp.nlp_data.evaluator, xs)
    elseif jp.objective !== nothing
        return MOIU.eval_variables(vi -> xs[vi.value], jp.objective)
    else
        return 0.0
    end
end

"""
    set_subsolver_option!(model::MOI.ModelLike, attr::MOI.AbstractOptimizerAttribute, change::Pair)

Set the optimizer of the model if the subsolver_name is part of the name of the optimizer.
Change the option `param` of the sub_solver. i.e `change=0.1 => 1e-5` means the default
normally is 0.1 and will be changed to 1e-5.
Normally reset_subsolver_option! should be run to return to reset this change after `optimize!`
Return the previous value of the param option or if not previously set return change.first
"""
function set_subsolver_option!(model::MOI.ModelLike, attr::MOI.AbstractOptimizerAttribute, change::Pair)
    prev = try
        MOI.get(model, attr)
    catch
        change.first
    end
    MOI.set(model, attr, change.second)
    return prev
end
function set_subsolver_option!(model::MOI.ModelLike, attr::MOI.AbstractOptimizerAttribute, change)
    MOI.set(model, attr, change)
    return
end
function set_subsolver_option!(model::MOI.ModelLike, subsolver_name::String, attr::MOI.AbstractOptimizerAttribute, change)
    if occursin(subsolver_name, MOI.get(model, MOI.SolverName()))
        return set_subsolver_option!(model, attr, change)
    end
end
function set_subsolver_option!(model::MOI.ModelLike, subsolver_name::String, param::String, change)
    set_subsolver_option!(model, subsolver_name, MOI.RawParameter(param), change)
end

function set_time_limit!(optimizer, time_limit::Union{Nothing,Float64})
    old_time_limit = Inf
    if MOI.supports(optimizer, MOI.TimeLimitSec())
        old_time_limit = MOI.get(optimizer, MOI.TimeLimitSec())
        MOI.set(optimizer, MOI.TimeLimitSec(), time_limit)
    end
    return old_time_limit
end

ext: jl


text:
export Cost,
        NoLoss, NoPenalty,
        AtomicLoss, AtomicPenalty,
        ScaledLoss, ScaledPenalty,
        CompositeLoss, CompositePenalty

abstract type Cost end

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Loss: (x, y) -> L(x, y)
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

abstract type Loss <: Cost end

struct NoLoss <: Loss end

abstract type AtomicLoss <: Loss end

mutable struct ScaledLoss{AL} <: Loss where AL <: AtomicLoss
    loss::AL
    scale::Float64
end

mutable struct CompositeLoss <: Loss
    losses::Vector{ScaledLoss}
end

(sl::ScaledLoss)(x::AVR, y::AVR)    = sl.scale * sl.loss(x, y)
(cl::CompositeLoss)(x::AVR, y::AVR) = sum(loss(x, y) for loss ∈ cl.losses)

getscale(n::NoLoss)     = 0.0
getscale(l::AtomicLoss) = 1.0
getscale(l::ScaledLoss) = l.scale

# Convenient extension for classification
abstract type MultiClassLoss{c} <: AtomicLoss where c end

getc(m) = 0
getc(m, y) = 0
getc(m::MultiClassLoss{c}) where c = c
getc(m::MultiClassLoss{0}, y) = maximum(y)
getc(m::MultiClassLoss{c}, y) where c = c

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Penalty: θ -> P(θ)
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

abstract type Penalty <: Cost end

struct NoPenalty <: Penalty end

abstract type AtomicPenalty <: Penalty end

mutable struct ScaledPenalty{AP} <: Penalty where AP <: AtomicPenalty
    penalty::AP
    scale::Float64
end

mutable struct CompositePenalty <: Penalty
    penalties::Vector{ScaledPenalty}
end

(sp::ScaledPenalty)(θ::AVR)    = sp.scale * sp.penalty(θ)
(cl::CompositePenalty)(θ::AVR) = sum(penalty(θ) for penalty ∈ cl.penalties)

getscale(n::NoPenalty)     = 0.0
getscale(p::AtomicPenalty) = 1.0
getscale(p::ScaledPenalty) = p.scale

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Objective function: (x, y, θ) -> L(x, y) + P(θ)
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

mutable struct ObjectiveFunction{L<:Loss,P<:Penalty} <: Cost
    loss::L
    penalty::P
end

(J::ObjectiveFunction)(y, ŷ, θ) = J.loss(y, ŷ) + J.penalty(θ)

# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Composition of Loss & Penalty functions
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

const AL = AtomicLoss
const AP = AtomicPenalty
const SL = ScaledLoss
const CL = CompositeLoss
const SP = ScaledPenalty
const CP = CompositePenalty
const NL = NoLoss
const NP = NoPenalty
const OF = ObjectiveFunction

scale1(a::AL) = ScaledLoss(a, 1.0)
scale1(a::AP) = ScaledPenalty(a, 1.0)

# Combinations with NoLoss (NL)
*(::NL, ::Real)  = NoLoss()
+(::NL, l::Loss) = l
+(l::Loss, ::NL) = l

# Combinations with NoPenalty (NP)
*(::NP, ::Real)     = NoPenalty()
+(::NP, p::Penalty) = p
+(p::Penalty, ::NP) = p

# Combinations with AtomicLoss (AL)
+(a::AL, b::AL)           = scale1(a) + scale1(b)
+(a::AL, b::Union{SL,CL}) = scale1(a) + b
+(b::Union{SL,CL}, a::AL) = a + b
*(a::AL, c::Real)         = ScaledLoss(a, float(c))

# Combinations with AtomicPenalty (AP)
+(a::AP, b::AP)           = scale1(a) + scale1(b)
+(a::AP, b::Union{SP,CP}) = scale1(a) + b
+(b::Union{SP,CP}, a::AP) = a + b
*(a::AP, c::Real)         = ScaledPenalty(a, float(c))

# Combinations with Scaled Losses and Combined Losses
+(a::SL{T},  b::SL{T})  where {T}     = ScaledLoss(a.loss, a.scale + b.scale)
+(a::SL{T1}, b::SL{T2}) where {T1,T2} = CL([a, b])
+(a::CL, b::CL) = begin
    a_  = a.losses
    a_T = typeof.(a_)
    c_  = copy(a_)
    rem   = ones(Bool, length(b.losses))
    for (i, L) in enumerate(b.losses)
        m = findfirst(typeof(L) .== a_T)
        if m !== nothing
            c_[m] = c_[m] + L # will be SL{T} + SL{T}
            rem[i] = false
        end
    end
    CL(vcat(c_, b.losses[rem]))
end
+(a::SL, c::CL)   = CL([a]) + c
+(c::CL, a::SL)   = a + c
*(a::SL, c::Real) = SL(a.loss, c * a.scale)
*(a::CL, c::Real) = CL(a.losses .* c)

# Combinations with Scaled Penalties and Combined Penalties
+(a::SP{T},  b::SP{T})  where {T}     = ScaledPenalty(a.penalty, a.scale + b.scale)
+(a::SP{T1}, b::SP{T2}) where {T1,T2} = CP([a, b])
+(a::CP, b::CP) = begin
    a_  = a.penalties
    a_T = typeof.(a_)
    c_  = copy(a_)
    rem = ones(Bool, length(b.penalties))
    for (i, P) in enumerate(b.penalties)
        m = findfirst(typeof(P) .== a_T)
        if m !== nothing
            c_[m] = c_[m] + P # will be SP{T} + SP{T}
            rem[i] = false
        end
    end
    CP(vcat(c_, b.penalties[rem]))
end
+(a::SP, c::CP)   = CP([a]) + c
+(c::CP, a::SP)   = a + c
*(a::SP, c::Real) = SP(a.penalty, c * a.scale)
*(a::CP, c::Real) = CP(a.penalties .* c)

# higher combinations ==> OF
+(l::Loss, p::Penalty) = OF(l, p)
+(p::Penalty, l::Loss) = l + p
+(o::OF, l::Loss)      = OF(o.loss + l, o.penalty)
+(l::Loss, o::OF)      = o + l
+(o::OF, p::Penalty)   = OF(o.loss, o.penalty + p)
+(p::Penalty, o::OF)   = o + p
+(a::OF, b::OF)        = OF(a.loss+b.loss, a.penalty+b.penalty)
*(o::OF, a::Real)      = OF(a * o.loss, a * o.penalty)

# Symetric relations
*(a::Real, c::Cost) = c * a

# - and / operations with Objective Functions (just use + and *)
-(a::Cost, b::Cost) = a + (-1 * b)
/(a::Cost, c::Real) = a * (1 / c)


text:
using DotaBot
using HTTP
using ArgParse
using Sockets
using JSON

s = ArgParseSettings()
@add_arg_table! s begin
    "--breezyIP"
    help = "breezy server IP adress"
    default = "127.0.0.1"
    "--breezyPort"
    help = "breezy server port number"
    default = "8085"
    "--agentIP"
    help = "agent server IP adress"
    default = "127.0.0.1"
    "--agentPort"
    help = "agent server port number"
    default = "8086"
end

args = parse_args(ARGS, s)

function get_action(features::Array{Float64})
    println(features[1])
    rand(0:29)
end

@async HTTP.serve(parse(IPAddr, args["agentIP"]), parse(Int, args["agentPort"])) do request::HTTP.Request
   @show request
   @show request.method
   @show HTTP.header(request, "Content-Type")
   @show HTTP.payload(request)
   try
       response = server_handler(request, args, get_action)
       println("responding")
       @show response
       return response
   catch e
       println(e)
       return HTTP.Response(200, JSON.json(Dict("done"=>false)))
   end
end

sleep(1)
start_game(args)
while true
    sleep(0.1)
end


text:
# using Revise

using IncrementalInference
using Test

@testset "Test priors" begin

N=100
graphinits = [false, true]

# TEMP value when not using for loop
graphinit = true

for  graphinit = graphinits

fg = initfg()
fg.solverParams.N = N

addVariable!(fg, :x0, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:x0], Prior(Normal(-1.0, 1.0)))

addVariable!(fg, :x1, ContinuousScalar, autoinit = graphinit, N=N)

addVariable!(fg, :x2, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:x2], Prior(Normal(+1.0, 1.0)))

addFactor!(fg, [:x0; :x1], LinearConditional(Normal(0.0, 0.01)))
addFactor!(fg, [:x1; :x2], LinearConditional(Normal(0.0, 0.01)))

#solve
tree, smt, hist = solveTree!(fg)
x0_m = getKDEMean(getKDE(getVariable(fg, :x0)))[1]
x1_m = getKDEMean(getKDE(getVariable(fg, :x1)))[1]
x2_m = getKDEMean(getKDE(getVariable(fg, :x2)))[1]

@info ("Testing means = 0 with 2 priors:\ngraphinit=$graphinit\nMeans: x0: $(x0_m), x1: $x1_m, x2: $x2_m")

@test isapprox(x0_m, 0.0, atol = 0.1)
@test isapprox(x1_m, 0.0, atol = 0.1)
@test isapprox(x2_m, 0.0, atol = 0.1)

end

for  graphinit = graphinits

fg = initfg()
fg.solverParams.N = N

addVariable!(fg, :x0, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:x0], Prior(Normal(-1.0, 1.0)))

addVariable!(fg, :l0, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:l0], Prior(Normal(+1.0, 1.0)))

addVariable!(fg, :l1, ContinuousScalar, autoinit = graphinit, N=N)

addFactor!(fg, [:x0; :l0], LinearConditional(Normal(0, 0.01)))
addFactor!(fg, [:x0; :l1], LinearConditional(Normal(0, 0.01)))

addVariable!(fg, :x1, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:x0; :x1], LinearConditional(Normal(0, 0.01)))

addVariable!(fg, :x2, ContinuousScalar, autoinit = graphinit, N=N)
addFactor!(fg, [:x1; :x2], LinearConditional(Normal(0, 0.01)))

addFactor!(fg, [:x2; :l0], LinearConditional(Normal(0, 0.01)))
addFactor!(fg, [:x2; :l1], LinearConditional(Normal(0, 0.01)))

#solve
tree, smt, hist = solveTree!(fg)

x0_m = getKDEMean(getKDE(getVariable(fg, :x0)))[1]
x1_m = getKDEMean(getKDE(getVariable(fg, :x1)))[1]
x2_m = getKDEMean(getKDE(getVariable(fg, :x2)))[1]
l0_m = getKDEMean(getKDE(getVariable(fg, :l0)))[1]
l1_m = getKDEMean(getKDE(getVariable(fg, :l1)))[1]

@info ("Testing means = 0 with 2 priors:\ngraphinit=$graphinit\nMeans: x0: $(x0_m), x1: $x1_m, x2: $x2_m, l0: $l0_m, l1: $l1_m")

@test isapprox(x0_m, 0.0, atol = 0.1)
@test isapprox(x1_m, 0.0, atol = 0.1)
@test isapprox(x2_m, 0.0, atol = 0.1)
@test isapprox(l0_m, 0.0, atol = 0.1)
@test isapprox(l1_m, 0.0, atol = 0.1)

end

end


text:
function RenderObject(
        data::Dict{Symbol}, program, pre,
        bbs = Node(AABB{Float32}(Vec3f0(0),Vec3f0(1))),
        main = nothing
    )
    RenderObject(convert(Dict{Symbol,Any}, data), program, pre, bbs, main)
end

function Base.show(io::IO, obj::RenderObject)
    println(io, "RenderObject with ID: ", obj.id)
end

Base.getindex(obj::RenderObject, symbol::Symbol)         = obj.uniforms[symbol]
Base.setindex!(obj::RenderObject, value, symbol::Symbol) = obj.uniforms[symbol] = value

Base.getindex(obj::RenderObject, symbol::Symbol, x::Function)     = getindex(obj, Val(symbol), x)
Base.getindex(obj::RenderObject, ::Val{:prerender}, x::Function)  = obj.prerenderfunctions[x]
Base.getindex(obj::RenderObject, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x]

Base.setindex!(obj::RenderObject, value, symbol::Symbol, x::Function)     = setindex!(obj, value, Val(symbol), x)
Base.setindex!(obj::RenderObject, value, ::Val{:prerender}, x::Function)  = obj.prerenderfunctions[x] = value
Base.setindex!(obj::RenderObject, value, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x] = value

const empty_signal = Node(false)
post_empty() = push!(empty_signal, false)

"""
Function which sets an argument of a Context/RenderObject.
If multiple RenderObjects are supplied, it'll try to set the same argument in all
of them.
"""
function set_arg!(robj::RenderObject, sym, value)
    current_val = robj[sym]
    set_arg!(robj, sym, current_val, value)
    # GLVisualize relies on reactives event system no for rendering
    # so if a change should be visible there must be an event to indicate change
    post_empty()
    nothing
end
function set_arg!(robj::Context, sym, value)
    set_arg!(robj.children, sym, value)
    nothing
end
function set_arg!(robj::Vector, sym, value)
    for elem in robj
        set_arg!(elem, sym, value)
    end
    nothing
end

function set_arg!(robj::RenderObject, sym, to_update::GPUArray, value)
    update!(to_update, value)
end
function set_arg!(robj::RenderObject, sym, to_update, value)
    robj[sym] = value
end
function set_arg!(robj::RenderObject, sym, to_update::Node, value::Node)
    robj[sym] = value
end
function set_arg!(robj::RenderObject, sym, to_update::Node, value)
    push!(to_update, value)
end

"""
Represents standard sets of function applied before rendering
"""
struct StandardPrerender
    transparency::Node{Bool}
    overdraw::Node{Bool}
end

function (sp::StandardPrerender)()
    if sp.overdraw[]
        # Disable depth testing if overdrawing
        glDisable(GL_DEPTH_TEST)
    else
        glEnable(GL_DEPTH_TEST)
        glDepthFunc(GL_LEQUAL)
    end
    # Disable depth write for transparent objects
    glDepthMask(sp.transparency[] ? GL_FALSE : GL_TRUE)
    # Disable cullface for now, untill all rendering code is corrected!
    glDisable(GL_CULL_FACE)
    # glCullFace(GL_BACK)
    enabletransparency()
end

struct StandardPostrender
    vao::GLVertexArray
    primitive::GLenum
end
function (sp::StandardPostrender)()
    render(sp.vao, sp.primitive)
end
struct StandardPostrenderInstanced{T}
    main::T
    vao::GLVertexArray
    primitive::GLenum
end
function (sp::StandardPostrenderInstanced)()
    renderinstanced(sp.vao, to_value(sp.main), sp.primitive)
end

struct EmptyPrerender
end
function (sp::EmptyPrerender)()
end
export EmptyPrerender
export prerendertype

function instanced_renderobject(data, program, bb = Node(AABB(Vec3f0(0), Vec3f0(1))), primitive::GLenum=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(convert(Dict{Symbol,Any}, data), program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrenderInstanced(main, robj.vertexarray, primitive)
    robj
end

function std_renderobject(data, program, bb = Node(AABB(Vec3f0(0), Vec3f0(1))), primitive=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(convert(Dict{Symbol,Any}, data), program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrender(robj.vertexarray, primitive)
    robj
end

prerendertype(::Type{RenderObject{Pre}}) where {Pre} = Pre
prerendertype(::RenderObject{Pre}) where {Pre} = Pre

extract_renderable(context::Vector{RenderObject}) = context
extract_renderable(context::RenderObject) = RenderObject[context]
extract_renderable(context::Vector{T}) where {T <: Composable} = map(extract_renderable, context)
function extract_renderable(context::Context)
    result = extract_renderable(context.children[1])
    for elem in context.children[2:end]
        push!(result, extract_renderable(elem)...)
    end
    result
end
transformation(c::RenderObject) = c[:model]
function transformation(c::RenderObject, model)
    c[:model] = const_lift(*, model, c[:model])
end
function transform!(c::RenderObject, model)
    c[:model] = const_lift(*, model, c[:model])
end

function _translate!(c::RenderObject, trans::TOrSignal{Mat4f0})
    c[:model] = const_lift(*, trans, c[:model])
end
function _translate!(c::Context, m::TOrSignal{Mat4f0})
    for elem in c.children
        _translate!(elem, m)
    end
end

function translate!(c::Composable, vec::TOrSignal{T}) where T <: Vec{3}
     _translate!(c, const_lift(translationmatrix, vec))
end
function _boundingbox(c::RenderObject)
    bb = to_value(c[:boundingbox])
    bb == nothing && return AABB()
    to_value(c[:model]) * bb
end
function _boundingbox(c::Composable)
    robjs = extract_renderable(c)
    isempty(robjs) && return AABB()
    mapreduce(_boundingbox, union, robjs)
end
"""
Copy function for a context. We only need to copy the children
"""
function Base.copy(c::Context{T}) where T
    new_children = [copy(child) for child in c.children]
    Context{T}(new_children, c.boundingbox, c.transformation)
end

"""
Copy function for a RenderObject. We only copy the uniform dict
"""
function Base.copy(robj::RenderObject{Pre}) where Pre
    uniforms = Dict{Symbol, Any}([(k,v) for (k,v) in robj.uniforms])
    robj = RenderObject{Pre}(
        robj.main,
        uniforms,
        robj.vertexarray,
        robj.prerenderfunction,
        robj.postrenderfunction,
        robj.boundingbox,
    )
    Context(robj)
end

# """
# If you have an array of OptimizedPrograms, you only need to put PreRender in front.
# """
# type OptimizedProgram{PreRender}
#     program::GLProgram
#     uniforms::FixedDict
#     vertexarray::GLVertexArray
#     gl_parameters::PreRender
#     renderfunc::Callable
#     visible::Boolean
# end


text:
using Clang.Generators
using Catalyst_jll

cd(@__DIR__) do
    include_dir = joinpath(Catalyst_jll.artifact_dir, "include", "catalyst-2.0") |> normpath
    catalyst_h = joinpath(include_dir, "catalyst_api.h") |> normpath
    @assert isfile(catalyst_h)
    blueprint_mesh_h = joinpath(include_dir, "conduit_blueprint.h") |> normpath
    @assert isfile(blueprint_mesh_h)

    options = load_options(joinpath(@__DIR__, "generator.toml"))
    args = push!(get_default_args(), "-I$include_dir")
    ctx = create_context([catalyst_h, blueprint_mesh_h], args, options)
    build!(ctx)
end


text:
# This file contains a function to
# read generator parameters from a file

"""
    @skip()

Skip two lines and check for EOF each time.

# Usage

```julia
using Scats.Internal.Gen: @skip

@skip() #= ===
begin
    readline(io)
    if eof(io)
        throw(ScatsGenEOF(file))
    end
    readline(io)
    if eof(io)
        throw(ScatsGenEOF(file))
    end
end =#
```
"""
macro skip()
    return esc(quote
        readline(io)
        if eof(io)
            throw(ScatsGenEOF(file))
        end
        readline(io)
        if eof(io)
            throw(ScatsGenEOF(file))
        end
    end)
end

"""
    @read(element::AbstractString, type::AbstractString)

Read a scalar element.

# Usage

```julia
using Scats.Internal.Gen: @read

@read("N", "IT") #= ===
begin
    try
        Gen.N = parse(IT, split(readline(io))[1])
    catch
        throw(ScatsGenWR_N(file))
    end
end =#
```
"""
macro read(element::AbstractString, type::AbstractString)
    return esc(
        Meta.parse("""
        try
            Gen.$element = parse($type, split(readline(io))[1])
        catch
            throw(ScatsGenWR_$element(file))
        end
        """)
    )
end

"""
    @read_array(element::AbstractString, type::AbstractString)

Read an array element.

# Usage

```julia
using Scats.Internal.Gen: @read_array

@read_array("A", "RT") #= ===
begin
    try
        Gen.A = parse.(RT, split(readline(io))[1:Gen.r])
    catch
        throw(ScatsGenWR_A(file))
    end
end =#
```
"""
macro read_array(element::AbstractString, type::AbstractString)
    return esc(
        Meta.parse("""
        try
            Gen.$element = parse.($type, split(readline(io))[1:Gen.r])
        catch
            throw(ScatsGenWR_$element(file))
        end
        """)
    )
end

"""
    read!(Gen::GenStruct, file::AbstractString)

Read generator parameters from a file to an instance of [`GenStruct`](@ref).

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Gen.example(file)
s.Gen.read!(file)

# output

```
"""
function read!(Gen::GenStruct, file::AbstractString)

    # Strip the line
    file = strip(file)

    # Check if the file exists
    if !isfile(file)
        throw(ScatsGenNotAFile(file))
    end

    # Open a file for reading
    open(file, "r") do io

        # Check for EOF
        if eof(io)
            throw(ScatsGenEOF(file))
        end

        readline(io)

        # Check for EOF
        if eof(io)
            throw(ScatsGenEOF(file))
        end

        @read "N" "IT"
        @skip

        @read "Δt" "RT"
        @skip

        @read "q" "RT"
        @skip

        @read "α" "RT"
        @skip

        @read "β" "RT"
        @skip

        @read "r" "IT"
        @skip

        @read_array "A" "RT"
        @skip

        @read_array "ν" "RT"
        @skip

        @read_array "ϕ" "RT"
        @skip

        @read "γ" "RT"

    end

    return nothing

end


text:
"""
    training_loop!(::updateInfo, ::lossInfo, ::trainingInfo, recorder::NamedTuple)

Train model defined in lossInfo, using update rule defined in updateInfo, and with training parameters defined in trainingInfo
"""
function training_loop!(upd::updateInfo, li::lossInfo, tr::trainingInfo, rec::NamedTuple; cb = () -> ())
    store = make_temp_storage(li, upd, tr, rec)
    store = initialise_updates!(store, upd)
    store = set_not_updated!(store, rec)
    rec, store = add_records!(rec,store,0)

    for i in 1:tr.iters
        store = do_updates!(store, upd)
        rec, store = add_records!(rec,store,i)
    end
    println("$(tr.iters) iterations")
    return store, rec
end


text:
import Base: iterate

struct APiterator{TA, TB, TX} # where {TA <:FeasibleSet, TB <:FeasibleSet, TX}
    setA::TA
    setB::TB
    x⁰::TX
end

mutable struct APstate{TX}
    xᵏ⁻¹::TX
    yᵏ::TX
    xᵏ::TX
end

function iterate(iter::APiterator)
    # state = APstate(copy(iter.x⁰), project(iter.x⁰,iter.setB), project(project(iter.x⁰,iter.setB),iter.setA))
    if iter.x⁰ === nothing x0 = complex(one(iter.setA.amp))
    else x0 = copy(iter.x⁰)
    end
    y1 = project(x0, iter.setB)
    x1 = project(y1, iter.setA)
    state = APstate(x0, y1, x1)
    return state, state
end

function iterate(iter::APiterator, state::APstate)
    state.xᵏ⁻¹ .= state.xᵏ
    project!(state.yᵏ, state.xᵏ⁻¹, iter.setB)
    project!(state.xᵏ, state.yᵏ, iter.setA)
    return state, state
end

Iterators.IteratorSize(::APiterator) = Iterators.IsInfinite()

# wrappers (see https://lostella.github.io/2018/07/25/iterative-methods-done-right.html)

struct HaltingIterable{I, F}
    iter::I
    fun::F
end

function iterate(iter::HaltingIterable)
    next = iterate(iter.iter)
    return dispatch(iter, next)
end

function iterate(iter::HaltingIterable, (instruction, state))
    if instruction == :halt return nothing end
    next = iterate(iter.iter, state)
    return dispatch(iter, next)
end

function dispatch(iter::HaltingIterable, next)
    if next === nothing return nothing end
    return next[1], (iter.fun(next[1]) ? :halt : :continue, next[2])
end

halt(iter::I, fun::F) where {I, F} = HaltingIterable{I, F}(iter, fun)

struct TeeIterable{I, F}
    iter::I
    fun::F
end

function iterate(iter::TeeIterable, args...)
    next = iterate(iter.iter, args...)
    if next !== nothing iter.fun(next[1]) end
    return next
end

tee(iter::I, fun::F) where {I, F} = TeeIterable{I, F}(iter, fun)

struct SamplingIterable{I}
    iter::I
    period::UInt
end

function iterate(iter::SamplingIterable, state=iter.iter)
    current = iterate(state)
    if current === nothing return nothing end
    for i = 1:iter.period-1
        next = iterate(state, current[2])
        if next === nothing return current[1], Iterators.rest(state, current[2]) end
        current = next
    end
    return current[1], Iterators.rest(state, current[2])
end

sample(iter::I, period) where I = SamplingIterable{I}(iter, period)

struct StopwatchIterable{I}
    iter::I
end

function iterate(iter::StopwatchIterable)
    t0 = time_ns()
    next = iterate(iter.iter)
    return dispatch(iter, t0, next)
end

function iterate(iter::StopwatchIterable, (t0, state))
    next = iterate(iter.iter, state)
    return dispatch(iter, t0, next)
end

function dispatch(iter::StopwatchIterable, t0, next)
    if next === nothing return nothing end
    return (time_ns()-t0, next[1]), (t0, next[2])
end

stopwatch(iter::I) where I = StopwatchIterable{I}(iter)

function loop(iter)
    x = nothing
    for y in iter x = y end
    return x
end


text:
###############################################################################
###          TISTING JULIA
###############################################################################
var_a=2+2
println(var_a)

a=Ref{Cfloat}(3)
b=Ref{Cfloat}(0)
c=Ref{Cfloat}(0)

ccall(("testf","./example"),Cvoid,(Ref{Cfloat},Ptr{Cfloat},Ptr{Cfloat}),a,b,c)
#b=ccall(("testf","./example"),Cfloat,(Ref{Cfloat},),a)
b=getindex(b)
e=getindex(c)
println(b)
#println(d)
println(c)
println(e)

max_stars_repo_name: jdlaubrie/shell-elem

max_issues_repo_name: jdlaubrie/shell-elem

max_forks_repo_name: jdlaubrie/shell-elem


text:
import Pkg
Pkg.activate(".")
Pkg.instantiate()

using Rebugger

# This doesn't work with Revise, so it's pretty useless.
# push!(LOAD_PATH, pwd() * "/lib")


text:
# Modelo de gerador de vapor modificado
type Gerador_VaporMod
	Gerador_VaporMod()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Brief=>"Taxa de calor gerado na caldeira",
				:PosX=>1,
				:PosY=>0.5
			)),
			Eficiencia(Dict{Symbol,Any}(
				:Brief=>"Eficiencia do gerador de vapor"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Taxa de calor no pre aquecedor de ar"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Taxa de calor nos superaquecedores"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Taxa de calor no evaporador"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Taxa de calor nos economizadores"
			)),
			Corrente(Dict{Symbol,Any}(
				:Brief=>"Agua de alimentacao",
				:Symbol=>"_{in_a}",
				:PosX=>0.5,
				:PosY=>1
			)),
			Corrente(Dict{Symbol,Any}(
				:Brief=>"Vapor Superaquecido",
				:Symbol=>"_{out_sa}",
				:PosX=>0.5,
				:PosY=>0
			)),
			Corrente(Dict{Symbol,Any}(
				:Brief=>"Evaporador"
			)),
			Corrente(Dict{Symbol,Any}(
				:Brief=>"Economizadores"
			)),
			[
				:([Feco.S,Feco.H] = PP2.propPTv(Feco.P,Feco.T)),
				:(Qec = Feco.F * (Feco.H - Fin_a.H)),
				:(Fvap.F = Feco.F),
				:([Fvap.S,Fvap.H] = PP2.propPTv(Fvap.P,Fvap.T)),
				:(Qca = Fvap.F * (Fvap.H - Feco.H)),
				:(Fout_sa.F = Fvap.F),
				:([Fout_sa.S,Fout_sa.H] = PP2.propPTv(Fout_sa.P,Fout_sa.T)),
				:(Qsa = Fout_sa.F * (Fout_sa.H - Fvap.H)),
				:(Q_GV * EF_GV = Qec + Qca + Qsa + Qpre),
			],
			[
				"Economizadores ECO1 + ECO1","","Evaporador - Camisa dagua","","","Superaquecedores BT + AT","","","Caldeira",
			],
			[:PP2,],
			[:Q_GV,:EF_GV,:Qpre,:Qsa,:Qca,:Qec,:Fin_a,:Fout_sa,:Fvap,:Feco,]
		)
	end
	PP2::DanaPlugin
	Q_GV::Pot_sinal
	EF_GV::Eficiencia
	Qpre::Potencia
	Qsa::Potencia
	Qca::Potencia
	Qec::Potencia
	Fin_a::Corrente
	Fout_sa::Corrente
	Fvap::Corrente
	Feco::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Gerador_VaporMod
function setEquationFlow(in::Gerador_VaporMod)
	#	[Fin_a.S,Fin_a.H] = PP2.propPTl(Fin_a.P,Fin_a.T); # Reduntante no ciclo fechado
	#	Feco.F = Fin_a.F; # Reduntante no ciclo fechado
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
end
function atributes(in::Gerador_VaporMod,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/caldeira"
	drive!(fields,_)
	return fields
end
Gerador_VaporMod(_::Dict{Symbol,Any})=begin
	newModel=Gerador_VaporMod()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(
        scene, 1, 2;
        rowsizes = Relative(1),
        colsizes = [Auto(), Fixed(200)],
        addedcolgaps = Fixed(30),
        alignmode = Outside(30, 30, 30, 30)
    )

    gl = maingl[1, 1] = GridLayout(
        3, 3;
        rowsizes = Relative(1/3),
        colsizes = Auto(),
        addedcolgaps = Fixed(20),
        addedrowgaps = Fixed(20),
        alignmode = Outside()
    )

    las = []
    for i in 1:3, j in 1:3
        # la =
        # al = ProtrusionLayout(gl, la.protrusions, la.bboxnode)
        la = gl[i, j] = LAxis(scene)
        sc = scatter!(
            la,
            rand(100, 2) .* 90 .+ 5,
            color=RGBf0(rand(3)...),
            raw=true, markersize=5)
        push!(las, la)
    end

    # buttons need change in abstractplotting to correctly update frame position

    glside = maingl[1, 2] = GridLayout(7, 1, alignmode=Outside())

    but = glside[1, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Titles")
    on(but.clicks) do c
        for la in las
            la.titlevisible[] = !la.titlevisible[]
        end
    end

    but2 = glside[2, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Labels")
    on(but2.clicks) do c
        for la in las
            la.xlabelvisible[] = !la.xlabelvisible[]
            la.ylabelvisible[] = !la.ylabelvisible[]
        end
    end

    but3 = glside[3, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Ticklabels")
    on(but3.clicks) do c
        for la in las
            la.xticklabelsvisible[] = !la.xticklabelsvisible[]
            la.yticklabelsvisible[] = !la.yticklabelsvisible[]
        end
    end

    but4 = glside[4, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Ticks")
    on(but4.clicks) do c
        with_updates_suspended(maingl) do
            for la in las
                la.xticksvisible[] = !la.xticksvisible[]
                la.yticksvisible[] = !la.yticksvisible[]
            end
        end
    end

    but5 = glside[5, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Tick Align")
    on(but5.clicks) do c
        maingl.block_updates = true
        for la in las
            la.xtickalign[] = la.xtickalign[] == 1 ? 0 : 1
            la.ytickalign[] = la.ytickalign[] == 1 ? 0 : 1
        end
        maingl.block_updates = true
        maingl.needs_update[] = true
    end

    but6 = glside[6, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Grids")
    on(but6.clicks) do c
        for la in las
            la.xgridvisible[] = !la.xgridvisible[]
            la.ygridvisible[] = !la.ygridvisible[]
        end
    end

    but7 = glside[7, 1] = LButton(scene, width = 200, height = 50, label = "Toggle Spines")
    on(but7.clicks) do c
        for la in las
            la.xspinevisible[] = !la.xspinevisible[]
            la.yspinevisible[] = !la.yspinevisible[]
            la.xoppositespinevisible[] = !la.xoppositespinevisible[]
            la.yoppositespinevisible[] = !la.yoppositespinevisible[]
        end
    end
end

map(la -> la.xticklabelpad = 5, las)
map(la -> la.xlabelpadding = 0, las)
map(la -> la.spinewidth = 1, las)
map(la -> la.xticklabelrotation = pi/2, las)
map(la -> la.xticklabelalign = (:right, :center), las)
map(la -> la.xticklabelspace = 50, las)

begin
    begin
        for i in 1:9
            las[i].ylabelvisible[] = false
            las[i].xlabelvisible[] = false
            sleep(0.05)
        end

        for i in 1:9
            las[i].yticklabelsvisible[] = false
            las[i].xticklabelsvisible[] = false
            sleep(0.05)
        end

        for i in 1:9
            las[i].titlevisible[] = false
            sleep(0.05)
        end

        for i in 1:9
            las[i].ylabelvisible[] = true
            las[i].xlabelvisible[] = true
            sleep(0.05)
        end

        for i in 1:9
            las[i].yticklabelsvisible[] = true
            las[i].xticklabelsvisible[] = true
            sleep(0.05)
        end

        for i in 1:9
            las[i].titlevisible[] = true
            sleep(0.05)
        end

        for i in 1:9
            las[i].title[] = "Big\nTitle"
            las[i].ylabel[] = "Big\ny label"
            las[i].xlabel[] = "Big\nx label"
            sleep(0.05)
        end

        for i in 1:9
            las[i].title[] = "Title"
            las[i].ylabel[] = "y label"
            las[i].xlabel[] = "x label"
            sleep(0.05)
        end
    end
    begin
        for i in 1:9
            las[i].ylabelsize[] = 30
            las[i].xlabelsize[] = 30
            las[i].yticklabelsize[] = 30
            las[i].xticklabelsize[] = 30
            sleep(0.05)
        end

        for i in 1:9
            las[i].ylabelsize[] = 20
            las[i].xlabelsize[] = 20
            las[i].yticklabelsize[] = 20
            las[i].xticklabelsize[] = 20
            sleep(0.05)
        end
    end
end


text:
function hess_approx(h, L, D, η; ϵ = sqrt(eps()))
    n = size(D)[1]
    θ = Vector{Float64}(undef, n)
    ω = Vector{Float64}(undef, n)
    ν = Vector{Float64}(undef, n)

    θ[n] = 0
    ν[1] = h[1, 1]

    # ξ = maximum(abs.(LDLt.L))
    ξ = maximum(abs.(UnitLowerTriangular(h)- Matrix(Diagonal(ones(n)))))
    γ = maximum(abs.(diag(h)))
    β² = max(ϵ, γ, ξ/n)
    β = sqrt(β²)
    # println("ξ/n = $(ξ/n) γ = $γ β = $β")

    ω[1] = max(β, abs(ν[1]), (θ[1]/β)^2)

    good_approx = true; j = 1

    # println("avant le while")

    while (good_approx) && (j <= n)
        # println("dans le while j = $j")
        if j > 1
            ν[j] = h[j, j] - sum(L[j, s]^2 * ω[s] for s = 1:j-1)
        end
        # println("ν[j] = $(ν[j])")

        for i=j+1:n
            # println("dans le for i = $i")
            # println(" θ[j] = $(θ[j])")
            # println(" maximum( abs.( h[i,j])) = $(maximum( abs.( h[i,j])))")
            if j > 1
                # println("on est dans le if qui fait la somme")
                sum_ljs_lis = sum(L[j, s] * L[i, s] * ω[s] for s=1:j-1)
                # println("sum_ljs_lis = $sum_ljs_lis")
            else
                # println("on est dans le if qui ne fait pas la somme")
                sum_ljs_lis = 0.0
                # println("sum_ljs_lis = $sum_ljs_lis")
            end
            θ[j] = max(θ[j],maximum( abs.( h[i,j] .- sum_ljs_lis)))
        end

        if (θ[j]/β)^2 >= η * ν[j]
            good_approx = false
        else
            good_approx = true
        end
        # println("j = $j good_approx = $good_approx")
        # println("ν[j] = $(ν[j]) θ[j] = $(θ[j])")
        j += 1
    end

    return good_approx
end


text:
using SpecialFunctions

"""
    G(x::Real)

Computes the second synchrotron function at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.

``G(x) = x K_{\\frac{2}{3}}(x)``
"""
G(x::Real) = x * besselk(2/3, x)

"""
    𝒢(x::Real)

Computes the second synchrotron function at a given frequency ratio ``x = \\frac{\\nu}{\\nu_0}``.

``G(x) = x K_{\\frac{2}{3}}(x)``

Wrapper for [`G`](@ref).
"""
𝒢(x::Real) = G(x)


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@inline coolDown(method::Symbol, args...) = coolDown(Val(method), args...)

"""
    coolDown(::Val{:linear}, (start, stop), steps)

Return a linear cooldown function.

The range `start - stop` is devided into equidistant `steps`.
"""
function coolDown(::Val{:linear}, (start, stop), steps)
    return step -> start - ((start - stop) / steps) * (step - 1)
end

"""
    coolDown(::Val{:leadInOut}, (start, stop), steps; leadIn = 0.1, leadOut = 0.95)

Return a linear cooldown function with lead-in and lead-out.

The range `start - stop` is devided into equidistant `steps`. If the given
`step` is below the `leadIn` quantile, `start` is returned; if it's above the
`leadOut` quantile, `stop` is returned.
"""
function coolDown(::Val{:leadInOut}, (start, stop), steps; leadIn = 0.1, leadOut = 0.95)
    @assert leadIn < leadOut
    @assert leadOut - leadIn >= 0
    return step -> if (step/steps <= leadIn) start
                   elseif (step/steps >= leadOut) stop
                   else start - ((start - stop) / steps) * (step - 1)
                   end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RandomPath

Traverse a spatial object with `N` points in a random
permutation of `1:N`.
"""
struct RandomPath <: AbstractPath end

traverse(object, path::RandomPath) = randperm(nelms(object))


text:
# LaleOperator abstract type is a subtype
# of the AutoMLPipeline Machine supertype
#
# abstract type Machine end
# abstract type Computer     <: Machine  end
# abstract type Workflow     <: Machine  end
# abstract type Learner      <: Computer end
# abstract type Transformer  <: Computer end
#
# abstract type LaleOperator <: Learner end
#
# typeof(LaleLearner)        <: LaleOperator
# typeof(LalePreprocessor)   <: LaleOperator
# typeof(LaleOptimizer)      <: LaleOperator
# ----------------------------------------

using Lale
using InteractiveUtils

using Random
using Statistics
using Test
using DataFrames: DataFrame
using AutoMLPipeline: Utils

iris = getiris()
Xreg = iris[:,1:3] |> DataFrame
Yreg = iris[:,4]   |> Vector
Xcl  = iris[:,1:4] |> DataFrame
Ycl  = iris[:,5]   |> Vector

# lale ops
pca     = laleoperator("PCA")
rb      = laleoperator("RobustScaler","autogen")
noop    = laleoperator("NoOp","lale")
rfr     = laleoperator("RandomForestRegressor")
rfc     = laleoperator("RandomForestClassifier")
treereg = laleoperator("DecisionTreeRegressor")
# amlp ops
ohe  = OneHotEncoder()
catf = CatFeatureSelector()
numf = NumFeatureSelector()

# Lale regression
lalepipe    = (pca + noop) >>  (rfr | treereg )
lale_hopt   = LalePipeOptimizer(lalepipe,max_evals = 10,cv = 3)
laletrained = fit(lale_hopt,Xreg,Yreg)
lalepred    = transform(laletrained,Xreg)
lalermse    = score(:rmse,lalepred,Yreg)

# AutoMLPipeline regression
amlpipe = @pipeline  (pca + noop) |> (rfr * treereg)
amlpred = fit_transform!(amlpipe,Xreg,Yreg)
crossvalidate(amlpipe,Xreg,Yreg,"mean_squared_error")
amlprmse = score(:rmse,amlpred,Yreg)

# Lale classification
lalepipe    = (rb + pca) |> rfc
lale_hopt   = LalePipeOptimizer(lalepipe,max_evals = 10,cv = 3)
laletrained = fit(lale_hopt,Xcl,Ycl)
lalepred    = transform(laletrained,Xcl)
laleacc     = score(:accuracy,lalepred,Ycl)

# AutoMLPipeline classification
amlpipe = @pipeline  (pca + rb) |> rfc
amlpred = fit_transform!(amlpipe,Xcl,Ycl)
crossvalidate(amlpipe,Xcl,Ycl,"accuracy_score")
amlpacc = score(:accuracy,amlpred,Ycl)

# AutoMLPipeline regression
plr = @pipeline (catf |> ohe) + (numf |> rb |> pca) |> rfr;
crossvalidate(plr,Xreg,Yreg,"mean_absolute_error",10,false)

#AutoMLPipeline classification
plc = @pipeline (catf |> ohe) + (numf |> rb |> pca) |> rfc;
crossvalidate(plc,Xcl,Ycl,"accuracy_score",10,false)


text:
vint = Vector{Int}
function make_frame(width::Int, height::Int, window::Int)
    frame = fill(-1f0, width+(2*window), height+(2*window))
    frame[window+1:end-window,window+1:end-window] .= 0
    frame
end

small_dd_builder() = DefaultDict{Int, Float32}(0f0)
big_dd_builder() = DefaultDict{Int, DefaultDict{Int, Float32}}(()->small_dd_builder())
bigger_dd_builder() = DefaultDict{NTuple{2, Float32}, DefaultDict{Int, DefaultDict{Int, Float32}}}(()->big_dd_builder())

rand_pos(env::NatureEnv) = rand(2) .* size(env) .|> x->ceil(Int,x)
outofbounds(env::NatureEnv, pos::Tuple{Int, Int}) = !all((1,1) .<= pos .<= size(env))

($)(a::Real, b::Real) = a:sign(b):a+b-sign(b)

function onehot(idx::Int, len::Int)
    z = zeros(len)
    z[idx] = 1
    z
end

function ecat(args...)
    cat(args..., dims=ndims(args[1])+1)
end

# Clips an n-dimensional vec/tuple between bounds at each dim
function clip(x::Union{Tuple,Vector},mins::Tuple,maxes::Tuple)
    Tuple(clamp(args...) for args in zip(x, mins, maxes))
end


text:
#=
Given a function that generates perfectly random numbers between 1 and k (inclusive), where k is an input, write a function that shuffles a deck of cards represented as an array using only swaps.

It should run in O(N) time.

Hint: Make sure each one of the 52! permutations of the deck is equally likely.
=#

using Test
include("Solutions/problem51_shuffledeck.jl")

@test sort(shuffle(collect(1:52))) == collect(1:52)

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Distributions

include("ParseUtils.jl")

function bearing(x1::Array{Float64,1}, x2::Array{Float64,1})::Float64
    return atan2(x2[2] - x1[2], x2[1] - x1[1])
end

function relative_bearing(pose::Array{Float64, 1}, point::Array{Float64, 1})::Float64
    pose2pointbearing = atan2(point[2] - pose[2], point[1] - pose[1])
    pose_bearing = pose[3]
    b = pose2pointbearing - pose_bearing
    if b < -pi
        b += 2.0*pi
    end
    if b >= pi
        b -= 2.0*pi
    end
    return b
end

function range_meas(x1::Array{Float64,1}, x2::Array{Float64,1})::Float64
    return norm(x2 - x1)
end

function rb_to_xy(meas::Array{Float64}, curr_pose::Array{Float64})::Array{Float64}
    x, y, w = curr_pose
    r, b = meas

    theta = w + b
    lm_x = x + r * cos(theta)
    lm_y = y + r * sin(theta)
    return [lm_x; lm_y]
end

function landmarks_in_range(pose::Array{Float64,1}, landmarks::Array{Float64}, sensor_max_range::Float64)::Array{Float64}
    lm_in_range = Array{Float64,2}(size(landmarks,1),0)
    for i in 1:size(landmarks,2)
        lm = landmarks[:,i]
        if range_meas(pose[1:2], lm[1:2]) < sensor_max_range
            lm_in_range = [lm_in_range lm]
        end
    end
    return lm_in_range
end

# Return all landmarks in field of view
# params:
#     pose::Array{Float64,1} estimated robot pose
#     landmarks::Array{Float64} array of landmarks to test
#     sensor_max_range::Float64 maximum sensor range (m)
#     sensor_fov::Float64 total sensor field of view (rad)
# return:
#     Array{Float64} array of landmarks within the field of view of the robot from the estimated pose
function landmarks_in_fov(pose::Array{Float64,1}, landmarks::Array{Float64}, sensor_max_range::Float64, sensor_fov::Float64)::Array{Float64}
    lm_in_fov = Array{Float64,2}(size(landmarks,1),0)
    for i in 1:size(landmarks,2)
        lm = landmarks[:,i]
        if range_meas(pose[1:2], lm[1:2]) < sensor_max_range && abs(relative_bearing(pose[1:3], lm[1:2])) < sensor_fov/2.0
            lm_in_fov = [lm_in_fov lm]
        end
    end
    return lm_in_fov
end

function measure(pose::Array{Float64,1},
                 landmarks::Array{Float64},
                 sensor_max_range::Float64,
                 bearing_noise::Float64,
                 range_noise::Float64)::Array{Float64}

    meas = Array{Float64,2}(2,0)
    lm_in_range = landmarks_in_range(pose, landmarks, sensor_max_range)
    if isempty(lm_in_range)
        return []
    end
    for i in 1:size(lm_in_range, 2)
        lm = lm_in_range[:,i]
        z = range_meas(pose[1:2], lm) + rand(Normal(0, range_noise))
        b = bearing(pose[1:2], lm) - pose[3] + rand(Normal(0, bearing_noise))

        # Correction for bearing?
        if b < 0.0
            b = b + 2*pi
        end
        if b > 2*pi
            b = b - 2*pi
        end

        meas = [meas [z; b]]
    end
    return meas
end

# transform coordinate p2 (originally relative to p1)
# into global coordinate frame [0.0; 0.0; 0.0]
function transform_relative_to_global(p1, p2)
    x1, y1, w1 = p1
    total_w = p2[3] + p1[3]
    R = [cos(w1) -sin(w1); sin(w1) cos(w1)]
    A = [R [x1; y1]; 0.0 0.0 1.0]
    p2 = [p2[1:2]; 1.0]
    global_p2 = A*p2
    global_p2[3] = total_w
    return global_p2
end

# transform coordinate p2 (originally in global frame)
# into relative coordinate frame p1
function transform_global_to_relative(p1, p2)
    x1, y1, w1 = p1
    dw = p2[3] - p1[3]
    R = [cos(w1) -sin(w1); sin(w1) cos(w1)]
    A = [R [x1; y1]; 0.0 0.0 1.0]
    p2 = [p2[1:2]; 1.0]
    rel_p2 = inv(A)*p2
    rel_p2[3] = dw
    return rel_p2
end

function measure_semantic(pose::Array{Float64,1},
                          landmarks::Array{Float64},
                          sensor_max_range::Float64,
                          sensor_fov::Float64,
                          bearing_noise::Float64,
                          range_noise::Float64,
                          n_classes::Int64)::Array{Float64}

    meas = Array{Float64,2}(2 + n_classes,0)
    @show lm_in_fov = landmarks_in_fov(pose, landmarks, sensor_max_range, sensor_fov)
    @show size(lm_in_fov)
    if isempty(lm_in_fov)
        return []
    end
    for i in 1:size(lm_in_fov,2)
        lm = lm_in_fov[:,i]
        z = range_meas(pose[1:2], lm[1:2]) + rand(Normal(0, range_noise))
        b = bearing(pose[1:2], lm[1:2]) - pose[3] + rand(Normal(0, bearing_noise))

        # Here the confusion matrix is set up for a classifier with 90% accuracy, independent of class
        # This can be made a parameter in future
        p = ones(n_classes) * (0.1 / (n_classes - 1))
        p[convert(Int64, lm[3])] = 0.9
        pred = rand(Categorical(p)) # Generate a random sample from confusion matrix probabilities
        s = ones(n_classes) * (0.1 / (n_classes - 1))
        s[pred] = 0.9

        # Correction for bearing
        if b < -pi
            b = b + 2*pi
        end

        if b >= pi
            b = b - 2*pi
        end

        meas = [meas [z; b; s]]
    end
    return meas
end

function run_sim(poses, landmarks, sensor_max_range, sensor_fov, range_noise, bearing_noise, odom_noise, n_classes; save_data=false)
    sbr_data = Array{Float64, 2}(1 + 2 + n_classes, 0)
    odom_data = Array{Float64, 2}(5, 0)
    for i in 1:size(poses,2)
        pose = poses[:,i]
        if i > 1
            previous_pose = poses[:,i-1]
            rel_xyz = transform_global_to_relative(previous_pose, pose + rand(MvNormal(zeros(3), odom_noise)))
            odom_meas = rel_xyz
            odom_meas = [i - 1; i; odom_meas]
            odom_data = [odom_data odom_meas]
        end
        sbr_meas = measure_semantic(pose, landmarks, sensor_max_range, sensor_fov, range_noise, bearing_noise, n_classes)
        if !isempty(sbr_meas)
          sbr_meas = [i; sbr_meas]
          sbr_data = [sbr_data sbr_meas]
        end

    end

    # Write results to files
    # Could param out a data directory
    if save_data
        stamp = Dates.format(Dates.now(), "yyyy-mm-dd-HH-MM-SS")
        @show stamp
        write_poses_landmarks(poses, landmarks, string("../data/poses_", stamp, ".txt"), string("../data/landmarks_", stamp, ".txt"))
        write_measurements(odom_data, string("../data/odom_", stamp, ".txt"))
        write_measurements(sbr_data, string("../data/sbr_", stamp, ".txt"))
    end
end

# Measurement parameters
sensor_max_range = 3.5
sensor_fov = 2.0*pi/3.0
landmark_sample_noise = 0.1
range_model_noise = 0.1
bearing_model_noise = 0.025
odom_noise = 0.001*eye(3)
n_classes = 2

# Set up landmarks; would be cool to randomly generate these
landmarks = [[1.0; 1.0; 1] [1.0; -1.0; 2] [3.0; 1.0; 2] [4.0; -1.0; 2] [6.0; 1.0; 1] [6.0; -1.0; 1] [8.0; -1.0; 2] [9.0; 1.0; 1] [11.0; -1.0; 2] [11.0; 2.0; 1] [11.0; 4.0; 2] [11.0; 5.0; 1] [11.0; 7.0; 2] [11.0; 9.0; 1] [11.0; 11.0; 2] [9.0; 3.0; 2] [9.0; 4.0; 1] [9.0; 6.0; 1] [9.0; 9.0; 2] [8.0; 9.0; 1] [8.0; 11.0; 2] [7.0; 11.0; 1] [5.0; 9.0; 2] [5.0; 11.0; 2] [4.0; 11.0; 1] [3.0; 9.0; 1] [1.0; 6.0; 1] [-1.0; 10.0; 2] [-1.0; 5.0; 1]]

# True robot poses; would be cool to randomly generate these
poses = [[0.0; 0.0; 0.0] [2.0; 0.0; 0.0] [4.0; 0.0; 0.0] [6.0; 0.0; 0.0] [8.0; 0.0; 0.0] [10.0; 0.0; pi/2] [10.0; 2.0; pi/2] [10.0; 4.0; pi/2] [10.0; 6.0; pi/2] [10.0; 8.0; pi/2] [10.0; 10.0; pi] [8.0; 10.0; pi] [6.0; 10.0; pi] [4.0; 10.0; pi] [2.0; 10.0; pi] [0.0; 10.0; 3*pi/2] [0.0; 8.0; 3*pi/2] [0.0; 6.0; 3*pi/2] [0.0; 4.0; 3*pi/2] [0.0; 2.0; 3*pi/2] [0.0; 0.0; 0.0]]

# Run the simulation and save data files
# run_sim(poses, landmarks, sensor_max_range, sensor_fov, range_model_noise, bearing_model_noise, odom_noise, n_classes, save_data=true)


text:
"""
    scatterlines(xs, ys, [zs]; kwargs...)

Plots `lines` between sets of x and y coordinates provided,
as well as plotting those points using `scatter`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(ScatterLines) do scene
    merge(default_theme(scene, Scatter), default_theme(scene, Lines))
end

function plot!(p::Combined{scatterlines, <:NTuple{N, Any}}) where N
    plot!(p, Lines, attributes_from(Lines, p), p[1:N]...)
    plot!(p, Scatter, attributes_from(Scatter, p), p[1:N]...)
end


text:
import Base.push!

export DSLBackend

# *************************************************************************** #
#                               DSL Backend
# *************************************************************************** #

"
The dsl backend will get functions that act on a tensor network circuit to
write dsl commands. It also holds filenames for files containing dsl commands,
tensor data and the result of fully contracting a tensor network circuit.
"
mutable struct DSLBackend <: AbstractBackend
    dsl_filename::String
    tensor_data_filename::String
    output_data_filename::String
    metrics::Metrics

    function DSLBackend(dsl::String="contract_network.tl",
                        tensor_data="tensor_data.h5"::String,
                        output::String="", overwrite::Bool=false)
        # Create an empty dsl file if none exists or overwrite is selected
        if !isfile(dsl) || overwrite
            close(open(dsl, "w"))
        end

        # If no output filename given by the user, use the tensor data file
        # to store output
        if output==""
            output = tensor_data
        end

        # empty the file if it exists
        if !isfile(tensor_data) || overwrite
            close(h5open(tensor_data, "w"))
        end

        new(dsl, tensor_data, output, Metrics())
    end
end

"""
    push!(dsl::DSLBackend, instruction::String)

Function to append a dsl command to the dsl script contained in a dsl backend.
"""
function push!(dsl::DSLBackend, instruction::String)
    open(dsl.dsl_filename, "a") do io
        write(io, instruction * "\n")
    end
end

"""
    save_tensor_data(backend::DSLBackend,
                              node_label::Symbol,
                              tensor_label::Symbol,
                              tensor_data::Array{<:Number})

Function to save tensor data to the hdf5 file specified by a dsl backend. It
also appends a load command to the dsl script contained in the backend to load
the saved tensor data as a tensor with the given node label.
"""
function save_tensor_data(backend::DSLBackend,
                          tensor_label::Symbol,
                          tensor_data::Array{<:Number})

    # Write the tensor data to the file specified by the dsl backend
    # and append commands to the dsl script to the tensor.
    h5open(backend.tensor_data_filename, "cw") do file
        tensor_label = string(tensor_label)
        if exists(file, tensor_label)
            o_delete(file, tensor_label)
        end
        write(file, tensor_label, tensor_data)
        push!(backend, "tensor $tensor_label $tensor_label")
    end
end

"""
    load_tensor_data(backend::DSLBackend,
                              tensor_label::Symbol)

Function to load tensor data from backend storage (if present)
"""
function load_tensor_data(backend::DSLBackend,
                          tensor_label::Symbol)

      if tensor_label == :result
          file = backend.output_data_filename
      else
          file = backend.tensor_data_filename
      end

      h5open(file, "r") do file
          tensor_label = string(tensor_label)
          if exists(file, tensor_label)
              return read(file, tensor_label)
          end
      end
end

"""
    contract_tensors(backend::DSLBackend,
                              A_label::Symbol, A_ncon_indices::Array{Int, 1},
                              B_label::Symbol, B_ncon_indices::Array{Int, 1},
                              C_label::Symbol)

Function to add dsl commands to a dsl backend's script that contract two tensors
A and B to create a third tensor C.
"""
function contract_tensors(backend::DSLBackend,
                          A_label::Symbol, A_ncon_indices::Array{Int, 1},
                          B_label::Symbol, B_ncon_indices::Array{Int, 1},
                          C_label::Symbol)

    # Create and save the dsl commands for contracting A and B.
    ncon_str = "ncon " * string(C_label) * " "
    ncon_str *= string(A_label) * " " * join(A_ncon_indices, ",") * " "
    ncon_str *= string(B_label) * " " * join(B_ncon_indices, ",")
    push!(backend, ncon_str)

    # Delete the old tensors.
    push!(backend, "del " * string(A_label))
    push!(backend, "del " * string(B_label))
end

"""
    save_output(backend::DSLBackend, node::Symbol, name::String)

Function to save the result of fully contracting a network under the given
name.
"""
function save_output(backend::DSLBackend, node::Symbol, name::String="result")
    command = "save $(node) $(backend.output_data_filename) $name"
    push!(backend, command)
end

"""
    reshape_tensor(backend::DSLBackend,
                            tensor::Symbol,
                            groups::Array{Array{Int, 1}, 1})

Function to add dsl command that reshapes a given tensor.
"""
function reshape_tensor(backend::DSLBackend,
                        tensor::Symbol,
                        groups::Array{Array{Int, 1}, 1})
    command = "reshape $tensor " * join([join(y, ",") for y in groups], ";")
    push!(backend, command)
end

"""
    permute_tensor(nothing, tensor::Symbol, axes)

Function to permute the axes of the given tensor
"""
function permute_tensor(backend::DSLBackend,
                        tensor::Symbol,
                        axes::Array{Int, 1})
    command = "permute $tensor " * join(axes, ",")
    push!(backend, command)
end

"""
    decompose_tensor!(backend::DSLBackend,
                               tensor::Symbol,
                               left_indices::Array{Int, 1},
                               right_indices::Array{Int, 1};
                               threshold::AbstractFloat=1e-13,
                               max_rank::Int=0,
                               left_label::Symbol,
                               right_label::Symbol)

Function to decompose a single tensor into two tensors. Returns 0 as an
indication that the dimension of the new virtual edge cannot be determined
until runtime.
"""
function decompose_tensor!(backend::DSLBackend,
                           tensor::Symbol,
                           left_indices::Array{Int, 1},
                           right_indices::Array{Int, 1};
                           threshold::AbstractFloat=1e-13,
                           max_rank::Int=0,
                           left_label::Symbol,
                           right_label::Symbol)
    cmd_str = "decompose $tensor "
    cmd_str *= "$left_label " * join(left_indices, ",")
    cmd_str *= " $right_label " * join(right_indices, ",")
    cmd_str *= " {\"threshold\":$threshold, \"max_rank\":$max_rank}"
    push!(backend, cmd_str)
    0
end

"""
    delete_tensor!(backend::DSLBackend, tensor_label::Symbol)

Mark tensor for deletion
"""
function delete_tensor!(backend::DSLBackend, tensor_label::Symbol)
    push!(backend, "del $tensor_label")
end

"""
    view_tensor!(backend::DSLBackend, view_node, node, bond_idx, bond_range)

Create a view on a tensor
"""
function view_tensor!(backend::DSLBackend, view_node, node, bond_idx, bond_range)
    bond_range_str = join(bond_range, ",")
    push!(backend, "view $view_node $node $bond_idx $bond_range_str")
end


text:
function generate_satisfy(domain::Domain, state::State,
                          domain_type::Symbol, state_type::Symbol)
    satisfy_def = quote
        function check(::$domain_type, state::$state_type, term::Const)
            return getfield(state, term.name)
        end
        function check(domain::$domain_type, state::$state_type, term::Compound)
            val = if term.name == :and
                all(check(domain, state, a) for a in term.args)
            elseif term.name == :or
                any(check(domain, state, a) for a in term.args)
            elseif term.name == :imply
                !check(domain, state, term.args[1]) |
                check(domain, state, term.args[2])
            elseif term.name == :not
                !check(domain, state, term.args[1])
            elseif term.name in keys(GLOBAL_PREDICATES)
                GLOBAL_PREDICATES[term.name](evaluate(domain, state, term.args[1]),
                                    evaluate(domain, state, term.args[2]))
            elseif !is_ground(term)
                missing
            else
                evaluate(domain, state, term)
            end
            return val
        end
        function satisfy(domain::$domain_type, state::$state_type, term::Term)
            val = check(domain, state, term)
            val !== missing ? val : !isempty(satisfiers(domain, state, term))
        end
        function satisfy(domain::$domain_type, state::$state_type,
                         terms::AbstractVector{<:Term})
            val = all(check(domain, state, t) for t in terms)
            val !== missing ? val : !isempty(satisfiers(domain, state, terms))
        end
    end
    return satisfy_def
end

function generate_satisfy(domain::AbstractedDomain, state::State,
                          domain_type::Symbol, state_type::Symbol)
    satisfy_def = quote
        function check(::$domain_type, state::$state_type, term::Const)
            return getfield(state, term.name)
        end
        function check(domain::$domain_type, state::$state_type, term::Compound)
            val = if term.name == :and
                (&)((check(domain, state, a) for a in term.args)...)
            elseif term.name == :or
                (|)((check(domain, state, a) for a in term.args)...)
            elseif term.name == :imply
                !check(domain, state, term.args[1]) |
                check(domain, state, term.args[2])
            elseif term.name == :not
                !check(domain, state, term.args[1])
            elseif term.name in keys(GLOBAL_PREDICATES)
                GLOBAL_PREDICATES[term.name](evaluate(domain, state, term.args[1]),
                                    evaluate(domain, state, term.args[2]))
            elseif !is_ground(term)
                missing
            else
                evaluate(domain, state, term)
            end
            return val
        end
        function satisfy(domain::$domain_type, state::$state_type, term::Term)
            val = check(domain, state, term)
            return val !== missing ? (val == true || val == both) :
                !isempty(satisfiers(domain, state, term))
        end
        function satisfy(domain::$domain_type, state::$state_type,
                         terms::AbstractVector{<:Term})
            val = (&)((check(domain, state, t) for t in terms)...)
            return val !== missing ? (val == true || val == both) :
                !isempty(satisfiers(domain, state, terms))
        end
    end
    return satisfy_def
end

function satisfiers(domain::CompiledDomain, state::CompiledState,
                    terms::AbstractVector{<:Term})
    gen_state = GenericState(state)
    return satisfiers(get_source(domain), gen_state, terms)
end

function satisfiers(domain::CompiledDomain, state::CompiledState, term::Term)
    gen_state = GenericState(state)
    return satisfiers(get_source(domain), gen_state, term)
end


text:
# Column-based operations: Some operations on rows are faster when considering columns

# In `map`, the output shouldn't alias inputs, so copies are made
Base.map(::typeof(identity), t::Union{FlexTable, Table}) = copy(t)

Base.map(::typeof(merge), t::Union{FlexTable, Table}) = copy(t)

function Base.map(::typeof(merge), t1::Table, t2::Table)
    return copy(Table(merge(columns(t1), columns(t2))))
end

function Base.map(::typeof(merge), df1::Union{Table{<:Any, N}, FlexTable{N}}, df2::Union{Table{<:Any, N}, FlexTable{N}}) where {N}
    return copy(FlexTable{N}(merge(columns(df1), columns(df2))))
end

function Base.map(::GetProperty{name}, t::Union{Table{<:Any, N}, FlexTable{N}}) where {name, N}
    return copy(getproperty(t, name::Symbol))::AbstractArray{<:Any, N}
end

# In `mapview`, the output should alias the inputs
function SplitApplyCombine.mapview(::typeof(merge), t1::Table, t2::Table)
    return Table(merge(columns(t1), columns(t2)))
end

function SplitApplyCombine.mapview(::typeof(merge), df1::Union{Table{<:Any, N}, FlexTable{N}}, df2::Union{Table{<:Any, N}, FlexTable{N}}) where {N}
    return FlexTable{N}(merge(columns(df1), columns(df2)))
end

@inline function SplitApplyCombine.mapview(f::GetProperty{name}, t::Union{Table{<:Any, N}, FlexTable{N}}) where {name,  N}
    return getproperty(t, name::Symbol)::AbstractArray{<:Any, N}
end

# broadcast

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, ::typeof(merge), ts::Table{<:Any, N}...) where {N}
	Table(merge(map(columns, ts)...))
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, ::typeof(merge), ts::Union{Table{<:Any, N},FlexTable{N}}...) where {N}
	FlexTable{N}(merge(map(columns, ts)...))
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, f::GetProperty{names}, t::Table{<:Any, N}) where {N, name}
	return getproperty(t, name::Symbol)
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, f::GetProperty{names}, t::FlexTable{N}) where {N, name}
	return getproperty(t, name::Symbol)::AbstractArray{<:Any, N}
end


text:
module MaxwellFDFD

using Reexport
@reexport using MaxwellBase
using AbbreviatedTypes

# The order of inclusion matters: if types or functions in file A are used in file B, file A
# must be included first.
include("source/source.jl")
# include("equation.jl")
include("model/model.jl")

end # module MaxwellFDFD


text:
using Multilane
using StatsBase
using MCTS
using JLD

behaviors=Dict{String,Any}()
for p in linspace(0., 3/4, 4)
    wv = (1-p)/8*ones(9)
    wv[1] = p
    behaviors[@sprintf("agents_%03d", 100*p)] =
            DiscreteBehaviorSet(Multilane.NINE_BEHAVIORS, WeightVec(wv))
end

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())
solvers = Dict{String, Any}(
    "dpw"=>dpws,
    "assume_normal"=>SingleBehaviorSolver(dpws, Multilane.NORMAL)
)

curve = TestSet(lambda=[0.1, 1.0, 2.15, 4.64, 10., 50.,], N=500, p_appear=1.0)

tests = []
for p in linspace(0., 3/4, 4)
    push!(tests, TestSet(curve, solver_key="dpw",
                         behaviors=@sprintf("agents_%03d", p*100),
                         key=@sprintf("upper_bound_%03d", p*100)))
    push!(tests, TestSet(curve, solver_key="assume_normal",
                         behaviors=@sprintf("agents_%03d", p*100),
                         key=@sprintf("assume_normal_%03d", p*100)))
end

objects = gen_initials(tests, behaviors=behaviors, generate_physical=true)

@show objects["param_table"]
objects["solvers"] = solvers
for is in values(objects["initial_states"])
    @assert !isnull(is.cars[1].behavior)
end

sbatch_spawn(tests, objects,
             batch_size=50,
             time_per_batch="1:00:00",
             submit_command="sbatch",
             template_name="sherlock.sh")

# results = evaluate(tests, objects, parallel=true)
#
# filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
# results["histories"] = nothing
# save(filename, results)
# println("results saved to $filename")

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
using Scryfall
using Test
@test getOracle("Black lotus")=="{T}, Sacrifice Black Lotus: Add three mana of any one color."


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Counts
using Nemo, NumberTheory

export ModuleCounts
export L000961, L002808, L005117, L013928, L025528, L065515
export L065855, L069637, L246547, L246655, L000720
export A007917, A151800, A257993
export PreviousPrime, NextPrime, PrimePiList
export takeFirst, Nth, Count, List, HilbertHotel

"""
* PreviousPrime, NextPrime, PrimePiList, takeFirst, Nth, Count, List, HilbertHotel
* L000961, L002808, L005117, L013928, L025528, L065515, L065855, L069637, L246547, L246655, L000720, A007917, A151800, A257993
"""
const ModuleCounts = ""

"""
Return a list of length len of integers ``≥ 0`` which are isA.
"""
function List(len, isA::Function)
    len ≤ 0 && return fmpz[]
    j, c = Int(0), Int(1)
    A = fill(ZZ(0), len)
    while c <= len
        if isA(j)
            A[c] = fmpz(j)
            c += 1
        end
        j += 1
    end
    A
end

"""
Iverson brackets.
"""
ι(b) = b ? 1 : 0

"""
Inverse Iverson brackets.
"""
ιι(n) = n == 0 ? true : false
 # The Unix way: A successful command returns a 0, while an unsuccessful one
 # returns a non-zero value that usually can be interpreted as an error code.

"""
Return a iterator of length n which has value 1 if isA(i) is true and otherwise 0.
"""
function Indicators(n, isA::Function)
    (ι(isA(i)) for i in 0:n - 1)
end

"""
Return a list of length len which gives the numbers of integers ≤ n which are isA. Integers start at ``n=0``.
```
julia> CountList(8, isPrime)
[0, 0, 1, 2, 2, 3, 3, 4]
```
"""
CountList(len::Int, isA::Function) = Accumulate(Indicators(len, isA))

# Consider two sequences A and invA. We say invA is the left inverse of A iff
# invA(A(n)) = n and A(n) is the least number m such that A(invA(m)) = A(n).
# A(invA(n)) = n if and only if isA(n) = true where 'isA' is the indicator
# function of the sequence A.
# We introduce 'Count' as the left inverse of 'Nth'. For example
# if Nth(n) = NthPrime(n) then Count(n) = PrimePi(n) (A000720).

"""
Return the numbers of integers in the range 0:n which are isA.
```
julia> Count(8, isPrime)
4
```
"""
Count(n::Int, isAb::Function) = Base.count((isAb(i) for i in 0:n))

"""
Return the numbers of integers in the range a:b which are isA.
```
julia> Count(3:8, isPrime)
3
```
"""
Count(r, isAb::Function) = Base.count((isAb(i) for i in r))

"""
Return a SeqArray listing the values satisfying the predicate isA for arguments ``0 ≤ x ≤ `` bound.
```
julia> FindUpTo(7, isPrime)
[2, 3, 5, 7]
```
"""
function FindUpTo(bound, isA::Function)
    bound < 0 && return fmpz[]
    filter(isA, 0:bound)
end

"""
Return the first ``n`` numbers satisfying the predicate isA.
"""
takeFirst(isA, n) = Iterators.take(Iterators.filter(isA, Iterators.countfrom(1)), n)

"""
Return a iterator listing the values satisfying the predicate isA for arguments in ``0 ≤ n ≤ bound .``
"""
function IterateUpTo(bound, isA::Function)
    (i for i in 0:bound if isA(i))
end

"""
Returns an integer which is the highest index in `b` for the value `a`. Whenever `a` is not a member of `b` it returns -1.
```
julia> L = List(10, isPrime); IndexIn(13, L)
5
```
"""
function IndexIn(a, b::AbstractArray)
    bdict = Dict(zip(b, 0:length(b)))
    get(bdict, fmpz(a), -1)
end

"""
Return the Nth integer which is isA. (For N ≤ 0 return 0.)
```
julia> Nth(7, isPrime)
17
```
"""
function Nth(N, isA::Function)
    N ≤ 0 && return 0
    n, c = Int(0), Int(0)
    while c < N
        i = isA(n)
        i && (c += 1)
        n += 1
    end
    n - 1
end

"""
Return the cumulative sum of an SeqArray.
"""
function Accumulate(A)
    R = fill(ZZ(0), length(A))
    i, acu = 1, 0
    for a in A
        acu += a
        R[i] = acu
        i += 1
    end
    R
end

"""
Return the smallest list of indicators of isA with ∑(A) = count.
```
julia> IndicatorsFind(7, isPrime)
[0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1]
```
"""
function IndicatorsFind(count, isA::Function)
    count ≤ 0 && return []
    n, c = Int(0), Int(0)
    A = Int[]
    while c < count
        i = isA(n)
        i && (c += 1)
        push!(A, ι(i))
        n += 1
    end
    A
end

"""
Return the first integer ``n ≥ 0`` such that isA(n) = true.
```
julia> First(isPrime)
2
```
"""
function First(isA::Function)
    n = 0
    while !isA(n)
        n += 1
    end
    n
end

First(A::Array{Int}) = A == [] ? nothing : first(A)

"""
Return the element at the end of the list A if A is not empty, 0 otherwise.
"""
Last(A) = A == [] ? "undef" : A[end]

"""
Trick described by David Hilbert in a 1924 lecture "Über das Unendliche".
"""
HilbertHotel(guest, hotel) = prepend!(hotel, guest)

"""
Return largest ``0 < k < n`` such that isA(k) = true or nothing if no such ``k`` exists.
```
julia> Previous(7, isPrime)
5
```
"""
function Previous(n, isA::Function)
    n == nothing && return First(isA)
    while true
        n -= 1
        isA(n) && break
        # This definition avoids throwing an exception. Alternatively:
        # n < 0 && throw(ArgumentError("Not defined for $n."))
        n < 0 && return First(isA)
    end
    n
end

"""
Return least ``k > n ≥ 0`` such that isA(k) = true. NOTE: It is assumed that such a ``k`` exists! (If not, the function will run forever.)
```
julia> Next(7, isPrime)
11
```
"""
function Next(n, isA::Function)
    ((n ≤ 0) || (n == nothing)) && return First(isA)
    while true
        n += 1
        isA(n) && break
    end
    n
end

"""
Return a list of composite numbers of length len. (Numbers which have more than one prime divisor.)
```
julia> L002808(8)
[4, 6, 8, 9, 10, 12, 14, 15]
```
"""
L002808(len) = List(len, isComposite)

"""
Return a list of the number of composite numbers ``≤ n``.
```
julia> L065855(8)
[0, 0, 0, 0, 1, 1, 2, 2]
```
"""
L065855(len) = CountList(len, isComposite)

"""
Return a list of squarefree numbers of length len. (Numbers which are not divisible by a square greater than 1.)
```
julia> L005117(8)
[1, 2, 3, 5, 6, 7, 10, 11]
```
"""
L005117(len) = List(len, isSquareFree)

"""
Return a list of the number of squarefree numbers ``< n``.
```
julia> L013928(8)
[0, 1, 2, 3, 3, 4, 5, 6]
```
"""
L013928(len) = CountList(len, isSquareFree)

"""
Return a list of powers of primes of length len. (Numbers of the form ``p^k`` where ``p`` is a prime and ``k ≥ 0``.)
```
julia> L000961(8)
[1, 2, 3, 4, 5, 7, 8, 9]
```
"""
L000961(len) = List(len, isPowerOfPrimes)

"""
Return the number of powers of primes ``≤ n``. (Powers of primes are numbers of the form ``p^k`` where ``p`` is a prime and ``k ≥ 0``.)
```
julia> L065515(8)
[0, 1, 2, 3, 4, 5, 5, 6]
```
"""
L065515(len) = CountList(len, isPowerOfPrimes)

"""
Return a list of prime powers of length len. (Numbers of the form ``p^k`` where ``p`` is a prime and ``k ≥ 1``.)
```
julia> L246655(8)
[2, 3, 4, 5, 7, 8, 9, 11]
```
"""
L246655(len) = List(len, isPrimePower)

"""
Return a list of the number of prime powers ``≤ n`` with exponents ``k ≥ 1``.
```
julia> L025528(8)
[0, 0, 1, 2, 3, 4, 4, 5]
```
"""
L025528(len) = CountList(len, isPrimePower)

"""
Return a list of perfect powers of length len. (Numbers of the form ``p^k`` where ``p`` is a prime and ``k ≥ 2``.
```
julia> L246547(8)
[4, 8, 9, 16, 25, 27, 32, 49]
```
"""
L246547(len) = List(len, isPerfectPower)

"""
Return a list of the number of prime powers ``≤ n`` with exponents ``k ≥ 2``.
```
julia> L069637(8)
[0, 0, 0, 0, 1, 1, 1, 1]
```
"""
L069637(len) = CountList(len, isPerfectPower)

# cf. also:
# A067535 Smallest squarefree number >= n.
# A070321 Largest squarefree number <= n.
# A025528 Number of prime powers <= n with exponents > 0.
# A000015 Smallest prime power >= n.
# A031218 Largest prime power <= n
# A167184 Smallest prime power >= n that is not prime.
# A081676 Largest perfect power <= n

"""
Return the largest prime in ``N`` (the semiring of natural numbers including zero) less than n for ``n ≥ 0``.
 (The `prev_prime` function of Mathematica, Maple, Magma and SageMath.)
"""
A007917(n::Int) = Previous(n, isPrime)

"""
Return the largest prime in ``Z`` (the ring of all integers) less than ``n`` for ``n ≥ 0`` (cf. A007917).
"""
PreviousPrime(n::Int) = n ∈ [0, 1, 2] ? -2 : Previous(n - 1, isPrime)

"""
Return least prime ``> n``. The next_prime function of Mathematica, Maple, Magma and SageMath (cf. A151800).
"""
NextPrime(n::Int) = Next(n, isPrime)

"""
Return least prime ``> n``. The `next_prime` function of Mathematica, Maple, Magma and SageMath.
"""
A151800(n::Int) = Next(n, isPrime)

"""
Return the list of number of primes ``≤ n`` for ``n ≥ 0``.

```
julia> PrimePiList(8)

[0, 0, 1, 2, 2, 3, 3, 4]
```
"""
PrimePiList(len::Int) = CountList(len, isPrime)

"""
Return the list of number of primes ``≤ n`` for ``n ≥ 0``.

```
julia> L000720(8)

[0, 0, 1, 2, 2, 3, 3, 4]
```
"""
L000720(len::Int) = PrimePiList(len)

"""
Return the index of the least prime not dividing ``n``.
"""
function A257993(n::Int)
    c, p = 1, 2
    while n % p == 0
        p = NextPrime(p)
        c += 1
    end
    c
end

#START-TEST-########################################################

using Test, SeqTests, SeqUtils

function test()

    indicators = [isPositive, isEven, isSquare, isPrime]
    indicatorNames = ["isPositive", "isEven", "isSquare", "isPrime"]

    len = 14
    @testset "Counts" begin

        @test Nth(96, isPrime) == 503
        @test Nth(97, isPrime) == 509
        @test Nth(98, isPrime) == 521

        # In other words: 97 is the 25-th prime.
        @test Count(96, isPrime) == 24
        @test Count(97, isPrime) == 25
        @test Count(98, isPrime) == 25

        @test List(24, isPrime)[end] == 89
        @test List(25, isPrime)[end] == 97

        for (i, isA) in enumerate(indicators)
            # This test shows that the logic behind 'Nth' and 'Count' is OK.
            for n in 1:len
               @test isA(n) == (Nth(Count(n, isA), isA) == n)
               @test     n  ==  Count(Nth(n, isA), isA)
            end
        end

        a = [A257993(n) for n in 1:10]
        b = [1, 2, 1, 2, 1, 3, 1, 2, 1, 2]
        @test all(a .== b)

        if is_oeis_installed()
            L = [
            L000961, L002808, L005117, L013928, L246547, L246655
            # L025528, L065515, L065855, L069637, L000720
            ]
            SeqTest(L, 'L')
        end
    end
end

function demo()

    indicators = [isNonnegative, isPositive, isEven, isComposite,
    isSquare, isSquareFree, isPrimePower, isPowerOfPrimes,
    isPerfectPower, isPrime
    ]
    indicatorNames = ["isNonnegative", "isPositive", "isEven", "isComposite",
    "isSquare", "isSquareFree", "isPrimePower", "isPowerOfPrimes",
    "isPerfectPower", "isPrime"
    ]

    println()
    len = 14
    for (i, isA) in enumerate(indicators)

        println("Predicate      ", indicatorNames[i])
        println("First          ", First(isA))
        println("---")

        println("Nth            ", [Nth(n, isA) for n in 0:len])
        print("List           ") ; Println(List(len, isA))
        println("FindUpTo       ", FindUpTo(len, isA))
        println("IterateUpTo    ", [k for k in IterateUpTo(len, isA)])
        println("---")

        println("IndicatorsFind ", IndicatorsFind(len, isA))
        println("Indicators     ", [k for k in Indicators(len, isA)])
        println("---")

        println("IndexIn (list) ", [IndexIn(fmpz(n), List(len, isA)) for n in 0:len])

        println("Count   (list) ", [Count(n, isA) for n in 0:len - 1])
        print("CountList      ") ; Println(CountList(len, isA))

        println("Previous(n)    ", [Previous(n, isA) for n in 0:len])
        println("Next(n)        ", [Next(n, isA)  for n in 0:len])
        println("Nth(Count(n))  ", [Nth(Count(n, isA), isA) for n in 0:len])
        println("Count(Nth(n))  ", [Count(Nth(n, isA), isA) for n in 0:len])

        println('='^66)
    end
end

"""
[A257993(n) for n in 1:10000]
    0.000635 seconds (9.20 k allocations: 221.859 KiB)
PrimePiList(10000)
    0.002290 seconds (35.82 k allocations: 637.891 KiB)
"""
function perf()
    @time [A257993(n) for n in 1:10000]
    GC.gc()
    @time PrimePiList(10000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
import DataFrames
import MLBase
import StatsBase

"""
"""
function singlelabelregressionytrue(
        labels::AbstractVector;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    result = float_type.(labels)
    return result
end

"""
"""
function singlelabelregressionypred(
        labels::AbstractVector;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    result = float_type.(labels)
    return result
end

"""
"""
function singlelabelregressionmetrics_resultdict(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        )
    ytrue = singlelabelregressionytrue(
        labels_df[single_label_name],
        )
    predictionsalllabels = predict(estimator, features_df)
    ypred = singlelabelregressionypred(
        predictionsalllabels[single_label_name],
        )
    results = Dict()
    results[:r2_score] = r2_score(
        ytrue,
        ypred,
        )
    results[:mean_square_error] = mean_square_error(
        ytrue,
        ypred,
        )
    results[:root_mean_square_error] = root_mean_square_error(
        ytrue,
        ypred,
        )
    results = fix_type(results)
    return results
end

"""
"""
function singlelabelregressionmetrics(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        )
    vectorofestimators = AbstractFittable[estimator]
    result = singlelabelregressionmetrics(
        vectorofestimators,
        features_df,
        labels_df,
        single_label_name,
        )
    return result
end

"""
"""
function singlelabelregressionmetrics(
        vectorofestimators::AbstractVector{AbstractFittable},
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol;
        kwargs...
        )
    metricsforeachestimator = [
        singlelabelregressionmetrics_resultdict(
            est,
            features_df,
            labels_df,
            single_label_name,
            )
            for est in vectorofestimators
        ]
    result = DataFrames.DataFrame()
    result[:metric] = [
        "R^2 (coefficient of determination)",
        "Mean squared error (MSE)",
        "Root mean square error (RMSE)",
        ]
    for i = 1:length(vectorofestimators)
        result[Symbol(vectorofestimators[i].name)] = [
            metricsforeachestimator[i][:r2_score],
            metricsforeachestimator[i][:mean_square_error],
            metricsforeachestimator[i][:root_mean_square_error],
            ]
    end
    return result
end


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools
using Optim

# Wavelength range [nm]
λ = 400:1000
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (p->1.0, s->0.0, weighted->between 0.0 and 1.0)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Indices of refraction
n1 = RIdb.air(beam.λ)
n2 = RIdb.silicon(beam.λ)
n3 = RIdb.glass(beam.λ)

# Define the RI models to use
layers = [
    LayerTMMO(n1), # 1
    ModelFit(:looyenga; N=(ninc=n1, nhost=n2)), # 2
    ModelFit(:looyenga; N=(ninc=n1, nhost=n2)), # 3
    LayerTMMO(n3), # 4
]

# Set the order of the layers (build the system) with the layers
order = [1, # incident medium
         2, 3, 2, 3, 2, 3, 2, 3, # top DBR of the MC
         2, 2, # defect of the MC
         3, 2, 3, 2, 3, 2, 3, 2, # bottom DBR of the MC
         4, # substrate
]

# Absolute reflectance spectrum to fit stored in Utils
Rexp = SpectraDB.fp_pectrum(beam.λ)
# plot(Spectrum1D(), beam.λ, Rexp)
# gui()

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=true,
);

# Seeds for each ModelFit layer defined above (No alpha)
seed = [
    [119.0, 0.8], # layers[2]
    [76.0, 0.5], # layers[3]
]

solOptim = fit_tmm_optics(
    Reflectance(Rexp), seed, beam, layers; order=order, options=options, alg=LBFGS(),
)

plot(FitSpectrum(), solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit)
gui()

# Seeds for each ModelFit layer defined above plus the alpha (with alpha)
seed2 = [
    [120, 0.8],
    [83.0, 0.5],
    [0.995], # alpha
]

solOptim2 = fit_tmm_optics(
    Reflectance(Rexp), seed2, beam, layers;
    order=order, options=options, alg=LBFGS(), alpha=true,
)

plot(FitSpectrum(), solOptim2.beam.λ, solOptim2.spectrumExp, solOptim2.spectrumFit)
gui()


text:
# This function constructs the field F_q \cong parent(h)/h and a homomorphism
# from parent(h) to F_q.
function field_extension(h::Union{ fq_nmod_poly, fq_poly })
  m = FqPolyRingToFqMor(h)
  return codomain(m), m
end

mutable struct FqPolyRingToFqMor{S, T, PolyType, MatType} <: Map{S, T, HeckeMap, FqPolyRingToFqMor}
  header::MapHeader{S, T}
  h::PolyType
  mat::MatType
  mat_inv::MatType

  function FqPolyRingToFqMor{S, T, PolyType, MatType}(h::PolyType) where {
           S <: Union{ FqNmodPolyRing, FqPolyRing },
           T <: Union{ FqNmodFiniteField, FqFiniteField },
           PolyType <: Union{ fq_nmod_poly, fq_poly },
           MatType <: Union{ gfp_mat, Generic.MatSpaceElem{Generic.ResF{fmpz}} }
    }

    z = new{S, T, PolyType, MatType}()
    z.h = h

    isnmod = (T == FqNmodFiniteField)

    Fq = base_ring(h)
    p = characteristic(Fq)
    if isnmod
      Fp = GF(Int(p); cached = false)
    else
      Fp = GF(p; cached = false)
    end
    Fpx = PolynomialRing(Fp, "x", cached = false)[1]
    g = Fpx()
    if isnmod
      pt = ccall((:fq_nmod_ctx_modulus, :libflint), Ptr{nmod_poly}, (Ref{FqNmodFiniteField}, ), Fq)
      ccall((:nmod_poly_set, :libflint), Nothing, (Ref{gfp_poly}, Ptr{nmod_poly}), g, pt)
    else
      pt = ccall((:fq_ctx_modulus, :libflint), Ptr{fmpz_mod_poly}, (Ref{FqFiniteField}, ), Fq)
      ccall((:fmpz_mod_poly_set, :libflint), Nothing, (Ref{gfp_fmpz_poly}, Ptr{fmpz_mod_poly}), g, pt)
    end
    n = degree(Fq)
    @assert n == degree(g)
    m = degree(h)
    if isnmod
      Fqm = FqNmodFiniteField(p, n*m, :$, false)
    else
      Fqm = FqFiniteField(p, n*m, :$, false)
    end
    Fqmy, y = PolynomialRing(Fqm, "y", cached = false)
    gy = Fqmy([ Fqm(coeff(g, i)) for i = 0:degree(g) ])
    a = roots(gy)[1]
    aa = Vector{typeof(a)}()
    push!(aa, Fqm(1))
    for i = 2:n
      push!(aa, aa[i - 1]*a)
    end
    @assert parent(a) == Fqm
    hy = Fqmy()
    for i = 0:m
      c = coeff(h, i)
      d = Fqm()
      for j = 0:n - 1
        d += coeff(c, j)*aa[j + 1]
      end
      setcoeff!(hy, i, d)
    end

    b = roots(hy)[1]
    bb = Vector{typeof(b)}()
    push!(bb, Fqm(1))
    for i = 2:m
      push!(bb, bb[i - 1]*b)
    end
    @assert parent(b) == Fqm

    # We have now a second basis for Fqm given by a^ib^j for 0 \leq i\leq n - 1, 0 \leq j\leq m - 1
    M = zero_matrix(Fp, n*m, n*m)
    c = Fqm()
    for i = 0:n - 1
      for j = 0:m - 1
        c = mul!(c, aa[i + 1], bb[j + 1])
        for k = 0:n*m - 1
          M[k + 1, i + j*n + 1] = Fp(coeff(c, k))
        end
      end
    end
    z.mat = M
    invM = inv(M)
    z.mat_inv = invM

    function _image(f::PolyType)
      f = mod(f, h)
      x_mat = zero_matrix(Fp, n*m, 1)
      for j = 0:m - 1
        t = coeff(f, j)
        for i = 0:n - 1
          x_mat[i + j*n + 1, 1] = Fp(coeff(t, i))
        end
      end
      x_mat = M*x_mat
      t = parent(g)([ x_mat[k, 1] for k = 1:n*m ])
      x = Fqm()
      if isnmod
        ccall((:fq_nmod_set, :libflint), Nothing, (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), x, t, Fqm)
      else
        ccall((:fq_set, :libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), x, t, Fqm)
      end
      return x
    end

    function _preimage(x::Union{ fq_nmod, fq })
      @assert typeof(x) == elem_type(T)
      x_mat = zero_matrix(Fp, n*m, 1)
      for k = 0:n*m - 1
        x_mat[k + 1, 1] = Fp(coeff(x, k))
      end
      x_mat = invM*x_mat
      f = parent(h)()
      t = Fqm()
      for j = 0:m - 1
        tt = parent(g)([ x_mat[i + n*j + 1, 1] for i = 0:n - 1 ])
        if isnmod
          ccall((:fq_nmod_set, :libflint), Nothing, (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}), t, tt, Fq)
        else
          ccall((:fq_set, :libflint), Nothing, (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}), t, tt, Fq)
        end
        setcoeff!(f, j, t)
      end
      return f
    end
    z.header = MapHeader{S, T}(parent(h), Fqm, _image, _preimage)
    return z
  end
end

function FqPolyRingToFqMor(h::fq_poly)
  return FqPolyRingToFqMor{FqPolyRing, FqFiniteField, fq_poly, Generic.MatSpaceElem{Generic.ResF{fmpz}}}(h)
end

function FqPolyRingToFqMor(h::fq_nmod_poly)
  return FqPolyRingToFqMor{FqNmodPolyRing, FqNmodFiniteField, fq_nmod_poly, gfp_mat}(h)
end