text:
using Test, Flux, StochasticDiffEq
println("Kolmogorov Tests")
using DiffEqDevTools , NeuralNetDiffEq
using Distributions

#Using SDEProblem for the Algorithm.
# For a diract delta take u0 = Normal(0 , sigma) where sigma --> 0
u0 = Normal(1.00 , 1.00)
xspan = (-2.0 , 6.0)
tspan = (0.0 , 1.0)
g(u , p , t) = 2.00
f(u , p , t) = -2.00
d = 1
sdealg = EM()
prob = SDEProblem(f , g , u0 , (0.0 , 1.0) ; xspan = xspan , d = d)
opt = Flux.ADAM(0.01)
m = Chain(Dense(1, 5, elu),Dense(5, 5, elu) , Dense(5 , 5 , elu) , Dense(5 , 1))
ensemblealg = EnsembleThreads()
sol = solve(prob, NNKolmogorov(m,opt , sdealg,ensemblealg) , verbose = true, dt = 0.01,
            abstol=1e-10, dx = 0.0001 , trajectories = 100000 ,  maxiters = 500)
# using Plots
#

# x_val = collect(xs)
# x_val= reshape(x_val , 1 , size(x_val)[1])
# y_val = m(x_val)
# y_val = reshape(y_val , 800001 , 1)
# x_val = collect(xs)
# plot(x_val , y_val,linewidth=3,title="Solution to the linear ODE with a thick line",
#      xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!")
# # plot(x_val , y_val)
# plot!(x_val , analytical(x_val),linewidth=3,title="Solution to the linear ODE with a thick line",
#      xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!")
#

## The solution is obtained taking the Fourier Transform.
analytical(xi) = pdf.(Normal(3 , sqrt(1.0 + 5.00)) , xi)
##Validation
xs = -5:0.00001:5
x_1 = rand(xs , 1 , 1000)
err_l2 = Flux.mse(analytical(x_1) , m(x_1))
@test err_l2 < 0.01
##

##Using the KolmogorovPDEProblem.
function phi(xi)
    y = Float64[]
    for x in xi
        y = push!(y , 1.77*x -0.015*x^3)
    end
    y = reshape(y , size(xi)[1] , size(xi)[2] )
    return y
end

xspan2 = (-6.0 , 6.0)
tspan2 = (0.0 , 1.0)
#f = mu and g  = sigma
g2(u , p , t) = 0.5*u
f2(u , p , t) = 0.5*0.25*u
d2 = 1
sdealg2 = EM()
prob2 = KolmogorovPDEProblem(f2 , g2, phi , xspan2 , tspan2, d2)
opt2 = Flux.ADAM(0.01)
m2 = Chain(Dense(1, 16, elu) , Dense(16 , 32 , elu),Dense(32 , 16 , elu), Dense(16 , 1))
sol = solve(prob2, NeuralNetDiffEq.NNKolmogorov(m2,opt2 , sdealg2, ensemblealg), verbose = true, dt = 0.01,
            dx = 0.0001 , trajectories = 1000 , abstol=1e-6, maxiters = 300)

function analytical2(xi)
    y = Float64[]
    a = 1.77*exp(0.5*(0.5)^2*1.0)
    b = -0.015*exp(0.5*(0.5*3)^2*1.0)
    for x in xi
        y = push!(y , a*x + b*x^3)
    end
    y = reshape(y , size(xi)[1] , size(xi)[2] )
    return y
end
xs2 = -5.00:0.01:5.00
x_val2 = rand(xs2 , d2 , 50)
errorl2 = Flux.mse(analytical2(x_val2) , m2(x_val2))
println("error_l2 = ", errorl2, "\n")
@test errorl2 < 0.4

##Non-Diagonal Test
f_noise = (du,u,p,t) -> du.=1.01u
g_noise = function (du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
Σ = [1.0 0.3 ; 0.3  1.0]
uo3 = MvNormal([0.0 ; 0.0], Σ)
sdealg3= EM()
xspan3 = (-10.0 , 10.0)
tspan3 = (0.0 , 1.0)
d3 = 2
prob = SDEProblem(f_noise , g_noise , uo3 , (0.0 , 1.0) ; xspan = xspan3 , d = d3 , noise_rate_prototype=zeros(2,4))
opt = Flux.ADAM(0.01)
m3 = Chain(Dense(d3, 32, elu) ,Dense(32 , 64 , elu), Dense(64 , 1))
sol3 = solve(prob, NeuralNetDiffEq.NNKolmogorov(m3,opt , sdealg3 , EnsembleThreads()), verbose = true, dt = 0.001,
            abstol=1e-6, dx = 0.001, trajectories = 1000,maxiters = 200)
println("Non-Diagonal test working.")


text:
# Useful for LinearAlgebra tests
function generate_well_conditioned_matrix(rng::AbstractRNG, T, N)
    A = randn(rng, T, N, N)
    return A * A' + I
end

function generate_well_conditioned_matrix(rng::AbstractRNG, N)
    return generate_well_conditioned_matrix(rng, Float64, N)
end

generate_well_conditioned_matrix(N) = generate_well_conditioned_matrix(Random.GLOBAL_RNG, N)

function generate_well_conditioned_matrix(T, N)
    return generate_well_conditioned_matrix(Random.GLOBAL_RNG, T, N)
end

ext: jl


text:
export LogicCircuit, GateType, InnerGate, LeafGate,
    LiteralGate, ConstantGate, isinnergate, isleafgate, ⋁Gate, ⋀Gate,
    isliteralgate, isconstantgate, is⋁gate, is⋀gate, literal_nodes,
    literal, constant, conjoin, disjoin, op, neutral,
    variable, ispositive, isnegative, istrue, isfalse,
    conjoin, disjoin, copy, compile, pos_literals, neg_literals, literals,
    fully_factorized_circuit,
    ⋁_nodes, ⋀_nodes, or_nodes, and_nodes,
    canonical_literals, canonical_constants, tree_formula_string,
    isflat, iscnf, isdnf, has_vars_contiguous

@reexport using DirectedAcyclicGraphs

import ..Utils: variables #extend

#####################
# Abstract infrastructure for logic circuit nodes
#####################

"Root of the logic circuit node hierarchy"
abstract type LogicCircuit <: DAG end

"""
A trait hierarchy denoting types of nodes
`GateType` defines an orthogonal type hierarchy of node types, not circuit types, so we can dispatch on node type regardless of circuit type.
See @ref{https://docs.julialang.org/en/v1/manual/methods/#Trait-based-dispatch-1}
"""
abstract type GateType end

"A logical gate that is a leaf node"
abstract type LeafGate <: GateType end

"A logical gate that is an inner node"
abstract type InnerGate <: GateType end

"A trait denoting literal leaf nodes of any type"
struct LiteralGate <: LeafGate end

"A trait denoting constant leaf nodes of any type"
struct ConstantGate <: LeafGate end

"A trait denoting conjuction nodes of any type"
struct ⋀Gate <: InnerGate end

"A trait denoting disjunction nodes of any type"
struct ⋁Gate <: InnerGate end

"Get the gate type trait of the given `LogicCircuit`"
@inline GateType(instance::LogicCircuit) = GateType(typeof(instance))

# map gate type traits to graph node traits
import DirectedAcyclicGraphs.NodeType # make available for extension
@inline NodeType(::Type{N}) where {N<:LogicCircuit} = NodeType(GateType(N))
@inline NodeType(::LeafGate) = Leaf()
@inline NodeType(::InnerGate) = Inner()

#####################
# node functions that need to be implemented for each type of circuit
#####################

import DirectedAcyclicGraphs.children # make available for extension by concrete types

"Get the logical literal in a given literal leaf node"
@inline literal(n::LogicCircuit)::Lit = n.literal # override when needed

"Conjoin nodes into a single circuit"
function conjoin end

"Disjoin nodes into a single circuit"
function disjoin end

"Create new circuit nodes in the given context."
function compile end

#####################
# derived node functions
#####################

"Is the node an inner gate?"
@inline isinnergate(n) = GateType(n) isa InnerGate
"Is the node a leaf gate?"
@inline isleafgate(n) = GateType(n) isa LeafGate
"Is the node an And gate?"
@inline is⋀gate(n) = GateType(n) isa ⋀Gate
"Is the node an Or gate?"
@inline is⋁gate(n) = GateType(n) isa ⋁Gate
"Is the node a literal gate?"
@inline isliteralgate(n) = GateType(n) isa LiteralGate
"Is the node a constant gate?"
@inline isconstantgate(n) = GateType(n) isa ConstantGate

"Get the logical constant in a given constant leaf node"
@inline constant(n::LogicCircuit) = constant(GateType(n), n)
@inline constant(::ConstantGate, n::LogicCircuit) = n.constant::Bool

"Get the logical variable in a given literal leaf node"
@inline variable(n::LogicCircuit) = variable(GateType(n), n)
@inline variable(::LiteralGate, n::LogicCircuit)::Var = lit2var(literal(n))::Var

"Get the sign of the literal leaf node"
@inline ispositive(n::LogicCircuit) = ispositive(GateType(n), n)
@inline ispositive(::LiteralGate, n::LogicCircuit)::Bool = literal(n) >= 0

@inline isnegative(n::LogicCircuit)::Bool = !ispositive(n)

"Is the circuit syntactically equal to true?"
@inline istrue(n::LogicCircuit)::Bool = istrue(GateType(n), n)
@inline istrue(::GateType, n::LogicCircuit)::Bool = false
@inline istrue(::ConstantGate, n::LogicCircuit)::Bool = (constant(n) == true)

"Is the circuit syntactically equal to false?"
@inline isfalse(n::LogicCircuit)::Bool = isfalse(GateType(n), n)
@inline isfalse(::GateType, n::LogicCircuit)::Bool = false
@inline isfalse(::ConstantGate, n::LogicCircuit)::Bool = (constant(n) == false)

#####################
# methods to easily construct circuits
#####################

@inline conjoin(xs::LogicCircuit...) = conjoin(collect(xs))
@inline disjoin(xs::LogicCircuit...) = disjoin(collect(xs))

@inline Base.:&(x::LogicCircuit, y::LogicCircuit) = conjoin(x,y)
@inline Base.:&(xs::LogicCircuit...) = conjoin(xs...)
@inline Base.:|(x::LogicCircuit, y::LogicCircuit) = disjoin(x,y)
@inline Base.:|(xs::LogicCircuit...) = disjoin(xs...)

@inline Base.:-(x::LogicCircuit) = begin
    @assert isliteralgate(x) "Negation is only supported for literal gates, not arbitrary circuits."
    compile(x, -literal(x))
end

# Get the function corresponding to the gate type
@inline op(::⋀Gate)::Function = conjoin
@inline op(::⋁Gate)::Function = disjoin

# Get the neural element corresponding to the gate type
@inline neutral(::⋀Gate)::Bool = true
@inline neutral(::⋁Gate)::Bool = false

# Syntactic sugar for compile
(T::Type{<:LogicCircuit})(args...) = compile(T, args...)

compile(n::LogicCircuit, args...) = compile(typeof(n), args...)

"Get a sequence of positive literals"
pos_literals(::Type{T}, num_lits::Int) where {T<:LogicCircuit} =
    map(l -> compile(T, Lit(l)), 1:num_lits)

"Get a sequence of negative literals"
neg_literals(::Type{T}, num_lits::Int) where {T<:LogicCircuit} =
    map(l -> compile(T, -Lit(l)), 1:num_lits)

"Get a sequence of positive and negative literals"
literals(::Type{T}, num_lits::Int) where {T<:LogicCircuit} =
    zip(pos_literals(T,num_lits), neg_literals(T,num_lits))

"Generate a fully factorized circuit over the given range of variables"
function fully_factorized_circuit end

#####################
# circuit inspection
#####################

"Get the list of conjunction nodes in a given circuit"
⋀_nodes(c::LogicCircuit) = filter(is⋀gate, c)

"Get the list of And nodes in a given circuit"
@inline and_nodes(c::LogicCircuit) = ⋀_nodes(c)

"Get the list of disjunction nodes in a given circuit"
⋁_nodes(c::LogicCircuit) = filter(is⋁gate, c)

"Get the list of or nodes in a given circuit"
@inline or_nodes(c::LogicCircuit) = ⋁_nodes(c)

"Get the list of literal nodes in a given circuit"
literal_nodes(c::LogicCircuit) = filter(isliteralgate, c)

"Construct a mapping from literals to their canonical node representation"
function canonical_literals(circuit::LogicCircuit)::Dict{Lit,LogicCircuit}
    lit_dict = Dict{Lit,LogicCircuit}()
    f_lit(n)= begin
        !haskey(lit_dict, literal(n)) || error("Circuit has multiple representations of literal $(literal(n))")
        lit_dict[literal(n)] = n
    end
    foreach(circuit, noop, f_lit, noop, noop)
    lit_dict
end

"Construct a mapping from constants to their canonical node representation"
function canonical_constants(circuit::LogicCircuit)::Tuple{Union{Nothing, LogicCircuit},Union{Nothing, LogicCircuit}}
    true_node::Union{Nothing, LogicCircuit} = nothing
    false_node::Union{Nothing, LogicCircuit} = nothing
    f_con(n)= begin
        if istrue(n)
            isnothing(true_node) || error("Circuit has multiple representations of true")
            true_node = n
        else
            @assert isfalse(n)
            isnothing(false_node) || error("Circuit has multiple representations of false")
            false_node = n
        end
    end
    foreach(circuit, f_con, noop, noop, noop)
    (false_node, true_node)
end

"""
Get the formula of a given circuit as a string, expanding the formula into a tree
"""
function tree_formula_string(n::LogicCircuit)
    if isliteralgate(n)
        "$(literal(n))"
    elseif isconstantgate(n)
        "$(constant(n))"
    elseif is⋀gate(n)
        s = ""
        for (i,c) in enumerate(children(n))
            if i < length(children(n))
                s = string(s, tree_formula_string(c), " ⋀ ")
            else
                s = string(s, tree_formula_string(c))
            end
        end
        s = string("(", s, ")")
        s
    else
        @assert is⋁gate(n)
        s = ""
        for (i,c) in enumerate(children(n))
            if i < length(children(n))
                s = string(s, tree_formula_string(c), " ⋁ ")
            else
                s = string(s, tree_formula_string(c))
            end
        end
        s = string("(", s, ")")
        s
    end
end

"Is the circuit a flat DNF or CNF structure?"
isflat(circuit) = iscnf(circuit) || isdnf(circuit)

"Is the circuit a conjunction of disjunctive clauses?"
iscnf(circuit) =
    is⋀gate(circuit) && all(children(circuit)) do clause
        is⋁gate(clause) && all(isliteralgate, children(clause))
    end

"Is the circuit a disjunction of conjunctive clauses?"
isdnf(circuit) =
    is⋁gate(circuit) && all(children(circuit)) do clause
        is⋀gate(clause) && all(isliteralgate, children(clause))
    end

"Does the circuit have a contiguously indexed set of variables"
has_vars_contiguous(circuit) = begin
    vars = variables(circuit)
    (maximum(vars) == length(vars))
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using FinancialMonteCarlo,Zygote;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=0.02;
T=1.0;
d=0.01;
D=90.0;

Nsim=100;
Nstep=3;
sigma=0.2
mc=MonteCarloConfiguration(Nsim,Nstep);
toll=0.8

rfCurve=ZeroRate(r);

FwdData=Forward(T)
EUData=EuropeanOption(T,K)
AMData=AmericanOption(T,K)
BarrierData=BarrierOptionDownOut(T,K,D)
AsianFloatingStrikeData=AsianFloatingStrikeOption(T)
AsianFixedStrikeData=AsianFixedStrikeOption(T,K)
Model=BlackScholesProcess(sigma,Underlying(S0,d));

#f__(x) = pricer(BlackScholesProcess(x[1]),ZeroRate(x[2],x[3],d),mc,EUData);
#f__(x) = pricer(BlackScholesProcess(x[1],Underlying(x[2],d)),ZeroRate(x[3]),mc,FwdData);
f__(x) = sum(simulate(BlackScholesProcess(x[1],Underlying(x[2],d)),ZeroRate(x[3]),mc,T));
#x=Float64[sigma,S0]
x=Float64[sigma,S0,r]
gradient(x -> f__(x), x)


text:
if basename(pwd()) == "aoc"
    cd("2015/17")
end
using Combinatorics

function part1(filename::AbstractString, target::Integer)
    count(==(target) ∘ sum, combinations(parse.(Int, eachline(filename))))
end

@assert part1("example.txt", 25) == 4
part1("input.txt", 150)

function part2(filename::AbstractString, target::Integer)
    lengths = map(length,
        Iterators.filter(
            ==(target) ∘ sum,
            combinations(parse.(Int, eachline(filename)))))
    count(==(first(lengths)), lengths)
end

@assert part2("example.txt", 25) == 3
part2("input.txt", 150)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using MEngProject
using Test
using Images
img = convert(Array{Float32,2}, load(datadir("img", "kan_sq_cont_l.png"))
p = Laminart.kernels(img, Parameters.parameters)
p = Laminart.add_I_u_p(img, p)

@testset "MEngProject.jl" begin

    out = zeros(p.dim_i, p.dim_j)
    b = reshape(ones(p.dim_i, p.dim_j * p.K), p.dim_i, p.dim_j, p.K)
    c = 2 .* ones(p.dim_i, p.dim_j)

    @test LaminartEqConv.conv!()
    @test LaminartEqConv.fun_f!()
    @test LaminartEqConv.I_u!()
    @test LaminartEqConv.fun_dv!()
    @test LaminartEqConv.fun_x_lgn()
    @test LaminartEqConv.fun_v_C!()
    @test LaminartEqConv.fun_dx_v1!()
    @test LaminartEqConv.fun_dy!()
    @test LaminartEqConv.fun_dm!()
    @test LaminartEqConv.fun_dz!()
    @test LaminartEqConv.fun_ds!()
    @test LaminartEqConv.fun_H_z!()
    @test LaminartEqConv.fun_dx_v2()
    @test LaminartEqConv.fun_dy_v2!()
    @test LaminartEqConv.fun_x_equ!()
    @test LaminartEqConv.fun_y_equ!()

    @test LaminartEqImfilter.conv!()
    @test LaminartEqImfilter.fun_f!()
    @test LaminartEqImfilter.I_u!()
    @test LaminartEqImfilter.fun_dv!()
    @test LaminartEqImfilter.fun_x_lgn()
    @test LaminartEqImfilter.fun_v_C!()
    @test LaminartEqImfilter.fun_dx_v1!()
    @test LaminartEqImfilter.fun_dy!()
    @test LaminartEqImfilter.fun_dm!()
    @test LaminartEqImfilter.fun_dz!()
    @test LaminartEqImfilter.fun_ds!()
    @test LaminartEqImfilter.fun_H_z!()
    @test LaminartEqImfilter.fun_dx_v2()
    @test LaminartEqImfilter.fun_dy_v2!()
    @test LaminartEqImfilter.fun_x_equ!()
    @test LaminartEqImfilter.fun_y_equ!()

    @test LaminartEqImfilterGPU_FFT.conv!()
    @test LaminartEqImfilterGPU_FFT.fun_f!()
    @test LaminartEqImfilterGPU_FFT.I_u!()
    @test LaminartEqImfilterGPU_FFT.fun_dv!()
    @test LaminartEqImfilterGPU_FFT.fun_x_lgn()
    @test LaminartEqImfilterGPU_FFT.fun_v_C!()
    @test LaminartEqImfilterGPU_FFT.fun_dx_v1!()
    @test LaminartEqImfilterGPU_FFT.fun_dy!()
    @test LaminartEqImfilterGPU_FFT.fun_dm!()
    @test LaminartEqImfilterGPU_FFT.fun_dz!()
    @test LaminartEqImfilterGPU_FFT.fun_ds!()
    @test LaminartEqImfilterGPU_FFT.fun_H_z!()
    @test LaminartEqImfilterGPU_FFT.fun_dx_v2()
    @test LaminartEqImfilterGPU_FFT.fun_dy_v2!()
    @test LaminartEqImfilterGPU_FFT.fun_x_equ!()
    @test LaminartEqImfilterGPU_FFT.fun_y_equ!()

    @test LaminartEqImfilterGPU_FIR.conv!()
    @test LaminartEqImfilterGPU_FIR.fun_f!()
    @test LaminartEqImfilterGPU_FIR.I_u!()
    @test LaminartEqImfilterGPU_FIR.fun_dv!()
    @test LaminartEqImfilterGPU_FIR.fun_x_lgn()
    @test LaminartEqImfilterGPU_FIR.fun_v_C!()
    @test LaminartEqImfilterGPU_FIR.fun_dx_v1!()
    @test LaminartEqImfilterGPU_FIR.fun_dy!()
    @test LaminartEqImfilterGPU_FIR.fun_dm!()
    @test LaminartEqImfilterGPU_FIR.fun_dz!()
    @test LaminartEqImfilterGPU_FIR.fun_ds!()
    @test LaminartEqImfilterGPU_FIR.fun_H_z!()
    @test LaminartEqImfilterGPU_FIR.fun_dx_v2()
    @test LaminartEqImfilterGPU_FIR.fun_dy_v2!()
    @test LaminartEqImfilterGPU_FIR.fun_x_equ!()
    @test LaminartEqImfilterGPU_FIR.fun_y_equ!()

    @test LaminartEqImfilterGPU_IIR.conv!()
    @test LaminartEqImfilterGPU_IIR.fun_f!()
    @test LaminartEqImfilterGPU_IIR.I_u!()
    @test LaminartEqImfilterGPU_IIR.fun_dv!()
    @test LaminartEqImfilterGPU_IIR.fun_x_lgn()
    @test LaminartEqImfilterGPU_IIR.fun_v_C!()
    @test LaminartEqImfilterGPU_IIR.fun_dx_v1!()
    @test LaminartEqImfilterGPU_IIR.fun_dy!()
    @test LaminartEqImfilterGPU_IIR.fun_dm!()
    @test LaminartEqImfilterGPU_IIR.fun_dz!()
    @test LaminartEqImfilterGPU_IIR.fun_ds!()
    @test LaminartEqImfilterGPU_IIR.fun_H_z!()
    @test LaminartEqImfilterGPU_IIR.fun_dx_v2()
    @test LaminartEqImfilterGPU_IIR.fun_dy_v2!()
    @test LaminartEqImfilterGPU_IIR.fun_x_equ!()
    @test LaminartEqImfilterGPU_IIR.fun_y_equ!()

    @test LaminartFunc.LamFunction()
    @test LaminartFunc.LamFunction_allStruct()
    @test LaminartFunc.LamFunction_gpu_reuse()
    @test LaminartFunc.LamFunction_equ()
    @test LaminartFunc.LamFunction_all_struct_reuse()
    @test LaminartFunc.LamFunction_imfil_cpu()
    @test LaminartFunc.LamFunction_imfil_gpu_iir()
    @test LaminartFunc.LamFunction_imfil_gpu_fir()
    @test LaminartFunc.LamFunction_imfil_gpu_fft()

    @test LaminartInitFunc.parameterInit_conv_gpu()
    @test LaminartInitFunc.parameterInit_conv_gpu_noise()
    @test LaminartInitFunc.parameterInit_conv_cpu()
    @test LaminartInitFunc.parameterInit_imfil_cpu()

    @test LaminartInitFunc.reshape2d_4d()
    @test LaminartInitFunc.add_I_u_p()

    @test LaminartInitFunc.kernels_conv_gpu()
    @test LaminartInitFunc.kernels_conv_cpu()
    @test LaminartInitFunc.kernels_imfil_cpu()

    @test kern_d_ph()
    @test kern_d_pv()
    @test kern_d_mh()
    @test kern_d_mv()
    @test kern_d_p()
    @test kern_d_m()
    @test kern_d()
    @test kern_b()
    @test gaussian_rot()
    @test fun_R()

    @test Parameters.para_var()
    @test Parameters.para_var_k()
    @test Parameters.parameters_f32()
    @test Parameters.parameters_f64()
end


text:
"""
---------------------------------------------------
Types and Constructors
---------------------------------------------------
"""

mutable struct NOMINAL_TAKEOFF_SAMPLER <: SAMPLER
	z_final::Sampleable
	vx_final::Sampleable
	θ::Sampleable
	t_extra::Sampleable
end

function nominal_takeoff_sampler(;z_final = Normal(500,20),
						vx_final = Uniform(30,45),
						θ = Uniform(3,5),
						t_extra = Uniform(0.1,0.3)) #Normal(10,4))
	return NOMINAL_TAKEOFF_SAMPLER(z_final, vx_final, θ, t_extra)
end

mutable struct NOMINAL_TAKEOFF <: UAM_TRAJECTORY
	z_final::Float64
	vx_final::Float64
	θ::Float64
	t_extra::Float64
	sampler::SAMPLER
	dt::Float64
	vz_max::Float64
	az_max::Float64
	ax_max::Float64
	p::Array{Float64,2}
	v::Array{Float64,2}
end

function nominal_takeoff(;z_final = 0.0,
						vx_final = 0.0,
						θ = 0.0,
						t_extra = 0.0,
						sampler = nominal_takeoff_sampler(),
						dt = 1.0,
						vz_max = 550fpm2fps,
						az_max = 0.3g_ft,
						ax_max = 0.1g,
						p = zeros(1,3),
						v = zeros(1,3))
	return NOMINAL_TAKEOFF(z_final, vx_final, θ, t_extra, sampler, dt,
							vz_max, az_max, ax_max, p, v)
end

"""
---------------------------------------------------
Functions
---------------------------------------------------
"""

function solve_trajectory!(τ::NOMINAL_TAKEOFF)
	# Get some needed info
	#ttot = abs(τ.z_final/500fpm2fps) + τ.t_extra*abs(τ.z_final/500fpm2fps)
	ttot = abs(τ.z_final/τ.vz_max) + τ.t_extra*abs(τ.z_final/τ.vz_max)
	N = convert(Int64, floor(ttot/τ.dt))
	# Difference matrix
	D = zeros(N-1,N)
	for i = 1:(N-1)
	    D[i,i] = -1
	    D[i,i+1] = 1
	end
	D = D./τ.dt

	# Create and solve the optimization problem
	px = Variable(N)
	pz = Variable(N)
	vx = Variable(N)
	vz = Variable(N)

	objec = norm(D*vx) + λ*norm(vz*ft2m + vx*tand(τ.θ))

	constraints = [[px[i+1] == px[i] + vx[i]*τ.dt for i = 1:N-1];
				   [pz[i+1] == pz[i] + vz[i]*τ.dt for i = 1:N-1];
				   vx[1] == 0; px[1] == 0;
				   vz[1] == 0; pz[1] == 0;
				   vx[N] == τ.vx_final; pz[N] == τ.z_final; vz[N] == 0;
				   vx ≥ 0; pz ≥ 0; vz ≥ 0; vz ≤ τ.vz_max;
				   D*vx ≤ τ.ax_max; D*vz*ft2m ≤ τ.az_max;
				   D*vx ≥ -τ.ax_max; D*vz*ft2m ≥ -τ.az_max]

	prob = minimize(objec, constraints)
	solve!(prob, ECOS.Optimizer(verbose=0))

	# Update the trajectory
	p = zeros(N,3)
	p[:,1] = px.value[:,1]
	p[:,3] = pz.value[:,1]
	v = zeros(N,3)
	v[:,1] = vx.value[:,1]
	v[:,3] = vz.value[:,1]

	τ.p = copy(p)
	τ.v = copy(v)

	return Int(prob.status) == 1
end

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


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

JLLWrappers.@generate_main_file_header("Ipopt")
JLLWrappers.@generate_main_file("Ipopt", UUID("9cc047cb-c261-5740-88fc-0cf96f7bdcc7"))
end  # module Ipopt_jll


text:
##
cd(@__DIR__)
using DelimitedFiles
input = readdlm("input.txt", Int) |> vec;

##
count(>(0),diff(input))

## PARTE 2
using DataStructures: nil, cons

##
res = foldl(1:length(input)-2; init=nil()) do l, w
    cons(sum(input[w:w+2]), l)
end;

##
count(>(0), diff(collect(reverse(res))))


text:
# Try this out by copy pasting each line into the Julia REPL
# That helps understand how an Objective-C class with the ObjC runtime system
# using Julia
using ObjCWrap

if !classexists(:Foobar)
    Foobar = createclass(:Foobar, NSObject)
end

# Create an instance
foo = @objc [[Foobar alloc] init]

function multiply(x::Float64, y::Float64)
   x * y
end

imp = @cfunction(multiply, Cdouble, (Cdouble, Cdouble))

# Call the C function we just made
result = ccall(imp, Cdouble, (Cdouble, Cdouble), 3, 4)

# Register method on Foobar class
# Notice that the signature includes: return value, self, selector, arg1, arg2
typestr = ObjCWrap.encodetype(Cdouble, Object, Selector, Cdouble, Cdouble)
ObjCWrap.setmethod(Foobar, sel"multiply:by:", imp, typestr)

# Check that the method is actually there
m = method(Foobar, sel"multiply:by:")
types = signature(m)
ctypes = ctype(types)

# Without DSL
result = message(foo, sel"multiply:by:", 3, 4)

# With DSL
result = @objc [foo multiply:3 by:4]


text:
module Packing
using GeometryBasics

include("rectangle.jl")
include("guillotine.jl")

export RectanglePacker, GuillotinePacker

end # module


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

mutable struct P4estElementContainer{NDIMS, RealT<:Real, uEltype<:Real, NDIMSP1, NDIMSP2, NDIMSP3} <: AbstractContainer
  # Physical coordinates at each node
  node_coordinates      ::Array{RealT, NDIMSP2}   # [orientation, node_i, node_j, node_k, element]
  # Jacobian matrix of the transformation
  # [jacobian_i, jacobian_j, node_i, node_j, node_k, element] where jacobian_i is the first index of the Jacobian matrix,...
  jacobian_matrix       ::Array{RealT, NDIMSP3}
  # Contravariant vectors, scaled by J, in Kopriva's blue book called Ja^i_n (i index, n dimension)
  contravariant_vectors ::Array{RealT, NDIMSP3}   # [dimension, index, node_i, node_j, node_k, element]
  # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix)
  inverse_jacobian      ::Array{RealT, NDIMSP1}   # [node_i, node_j, node_k, element]
  # Buffer for calculated surface flux
  surface_flux_values   ::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element]

  # internal `resize!`able storage
  _node_coordinates     ::Vector{RealT}
  _jacobian_matrix      ::Vector{RealT}
  _contravariant_vectors::Vector{RealT}
  _inverse_jacobian     ::Vector{RealT}
  _surface_flux_values  ::Vector{uEltype}
end

@inline nelements(elements::P4estElementContainer) = size(elements.node_coordinates, ndims(elements) + 2)
@inline Base.ndims(::P4estElementContainer{NDIMS}) where NDIMS = NDIMS
@inline Base.eltype(::P4estElementContainer{NDIMS, RealT, uEltype}) where {NDIMS, RealT, uEltype} = uEltype

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(elements::P4estElementContainer, capacity)
  @unpack _node_coordinates, _jacobian_matrix, _contravariant_vectors,
    _inverse_jacobian, _surface_flux_values = elements

  n_dims = ndims(elements)
  n_nodes = size(elements.node_coordinates, 2)
  n_variables = size(elements.surface_flux_values, 1)

  resize!(_node_coordinates, n_dims * n_nodes^n_dims * capacity)
  elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
    (n_dims, ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_jacobian_matrix, n_dims^2 * n_nodes^n_dims * capacity)
  elements.jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix),
    (n_dims, n_dims, ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_contravariant_vectors, length(_jacobian_matrix))
  elements.contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors),
    size(elements.jacobian_matrix))

  resize!(_inverse_jacobian, n_nodes^n_dims * capacity)
  elements.inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian),
    (ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_surface_flux_values,
    n_variables * n_nodes^(n_dims-1) * (n_dims*2) * capacity)
  elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
    (n_variables, ntuple(_ -> n_nodes, n_dims-1)..., n_dims*2, capacity))

  return nothing
end

# Create element container and initialize element data
function init_elements(mesh::P4estMesh{NDIMS, RealT}, equations,
                       basis, ::Type{uEltype}) where {NDIMS, RealT<:Real, uEltype<:Real}
  nelements = ncells(mesh)

  _node_coordinates = Vector{RealT}(undef, NDIMS * nnodes(basis)^NDIMS * nelements)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
    (NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _jacobian_matrix = Vector{RealT}(undef, NDIMS^2 * nnodes(basis)^NDIMS * nelements)
  jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix),
    (NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _contravariant_vectors = similar(_jacobian_matrix)
  contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors),
    size(jacobian_matrix))

  _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements)
  inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian),
    (ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _surface_flux_values = Vector{uEltype}(undef,
    nvariables(equations) * nnodes(basis)^(NDIMS-1) * (NDIMS*2) * nelements)
  surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
    (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., NDIMS*2, nelements))

  elements = P4estElementContainer{NDIMS, RealT, uEltype, NDIMS+1, NDIMS+2, NDIMS+3}(
    node_coordinates, jacobian_matrix, contravariant_vectors,
    inverse_jacobian, surface_flux_values,
    _node_coordinates, _jacobian_matrix, _contravariant_vectors,
    _inverse_jacobian, _surface_flux_values)

  init_elements!(elements, mesh, basis)
  return elements
end

mutable struct P4estInterfaceContainer{NDIMS, uEltype<:Real, NDIMSP2} <: AbstractContainer
  u             ::Array{uEltype, NDIMSP2}       # [primary/secondary, variable, i, j, interface]
  element_ids   ::Matrix{Int}                   # [primary/secondary, interface]
  node_indices  ::Matrix{NTuple{NDIMS, Symbol}} # [primary/secondary, interface]

  # internal `resize!`able storage
  _u            ::Vector{uEltype}
  _element_ids  ::Vector{Int}
  _node_indices ::Vector{NTuple{NDIMS, Symbol}}
end

@inline ninterfaces(interfaces::P4estInterfaceContainer) = size(interfaces.element_ids, 2)
@inline Base.ndims(::P4estInterfaceContainer{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(interfaces::P4estInterfaceContainer, capacity)
  @unpack _u, _element_ids, _node_indices = interfaces

  n_dims = ndims(interfaces)
  n_nodes = size(interfaces.u, 3)
  n_variables = size(interfaces.u, 2)

  resize!(_u, 2 * n_variables * n_nodes^(n_dims-1) * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
    (2, n_variables, ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(_element_ids, 2 * capacity)
  interfaces.element_ids = unsafe_wrap(Array, pointer(_element_ids), (2, capacity))

  resize!(_node_indices, 2 * capacity)
  interfaces.node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, capacity))

  return nothing
end

# Create interface container and initialize interface data.
function init_interfaces(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_interfaces = count_required_surfaces(mesh).interfaces

  _u = Vector{uEltype}(undef, 2 * nvariables(equations) * nnodes(basis)^(NDIMS-1) * n_interfaces)
  u = unsafe_wrap(Array, pointer(_u),
    (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_interfaces))

  _element_ids = Vector{Int}(undef, 2 * n_interfaces)
  element_ids = unsafe_wrap(Array, pointer(_element_ids), (2, n_interfaces))

  _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces)
  node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces))

  interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS+2}(u, element_ids, node_indices,
                                                                _u, _element_ids, _node_indices)

  init_interfaces!(interfaces, mesh)

  return interfaces
end

function init_interfaces!(interfaces, mesh::P4estMesh)
  init_surfaces!(interfaces, nothing, nothing, mesh)

  return interfaces
end

mutable struct P4estBoundaryContainer{NDIMS, uEltype<:Real, NDIMSP1} <: AbstractContainer
  u           ::Array{uEltype, NDIMSP1}       # [variables, i, j, boundary]
  element_ids ::Vector{Int}                   # [boundary]
  node_indices::Vector{NTuple{NDIMS, Symbol}} # [boundary]
  name        ::Vector{Symbol}                # [boundary]

  # internal `resize!`able storage
  _u          ::Vector{uEltype}
end

@inline nboundaries(boundaries::P4estBoundaryContainer) = length(boundaries.element_ids)
@inline Base.ndims(::P4estBoundaryContainer{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(boundaries::P4estBoundaryContainer, capacity)
  @unpack _u, element_ids, node_indices, name = boundaries

  n_dims = ndims(boundaries)
  n_nodes = size(boundaries.u, 2)
  n_variables = size(boundaries.u, 1)

  resize!(_u, n_variables * n_nodes^(n_dims-1) * capacity)
  boundaries.u = unsafe_wrap(Array, pointer(_u),
    (n_variables, ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(element_ids, capacity)

  resize!(node_indices, capacity)

  resize!(name, capacity)

  return nothing
end

# Create interface container and initialize interface data in `elements`.
function init_boundaries(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_boundaries = count_required_surfaces(mesh).boundaries

  _u = Vector{uEltype}(undef, nvariables(equations) * nnodes(basis)^(NDIMS-1) * n_boundaries)
  u = unsafe_wrap(Array, pointer(_u),
    (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_boundaries))

  element_ids  = Vector{Int}(undef, n_boundaries)
  node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, n_boundaries)
  names        = Vector{Symbol}(undef, n_boundaries)

  boundaries = P4estBoundaryContainer{NDIMS, uEltype, NDIMS+1}(u, element_ids,
                                                               node_indices, names, _u)

  if n_boundaries > 0
    init_boundaries!(boundaries, mesh)
  end

  return boundaries
end

function init_boundaries!(boundaries, mesh::P4estMesh)
  init_surfaces!(nothing, nothing, boundaries, mesh)

  return boundaries
end

# Function barrier for type stability
function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh)
  # Extract boundary data
  side = unsafe_load_side(info)
  # Get local tree, one-based indexing
  tree = unsafe_load_tree(mesh.p4est, side.treeid + 1)
  # Quadrant numbering offset of this quadrant
  offset = tree.quadrants_offset

  # Verify before accessing is.full, but this should never happen
  @assert side.is_hanging == false

  local_quad_id = side.is.full.quadid
  # Global ID of this quad
  quad_id = offset + local_quad_id

  # Write data to boundaries container
  # p4est uses zero-based indexing; convert to one-based indexing
  boundaries.element_ids[boundary_id] = quad_id + 1

  # Face at which the boundary lies
  face = side.face

  # Save boundaries.node_indices dimension specific in containers_[23]d.jl
  init_boundary_node_indices!(boundaries, face, boundary_id)

  # One-based indexing
  boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1]

  return nothing
end

# Container data structure (structure-of-arrays style) for DG L2 mortars
#
# The positions used in `element_ids` are 1:3 (in 2D) or 1:5 (in 3D), where 1:2 (in 2D)
# or 1:4 (in 3D) are the small elements numbered in z-order and 3 or 5 is the large element.
# The solution values on the mortar element are saved in `u`, where `position` is the number
# of the small element that corresponds to the respective part of the mortar element.
# The first dimension `small/large side` takes 1 for small side and 2 for large side.
#
# Illustration of the positions in `element_ids` in 3D, where ξ and η are the local coordinates
# of the mortar element, which are precisely the local coordinates that span
# the surface of the smaller side.
# Note that the orientation in the physical space is completely irrelevant here.
#   ┌─────────────┬─────────────┐  ┌───────────────────────────┐
#   │             │             │  │                           │
#   │    small    │    small    │  │                           │
#   │      3      │      4      │  │                           │
#   │             │             │  │           large           │
#   ├─────────────┼─────────────┤  │             5             │
# η │             │             │  │                           │
#   │    small    │    small    │  │                           │
# ↑ │      1      │      2      │  │                           │
# │ │             │             │  │                           │
# │ └─────────────┴─────────────┘  └───────────────────────────┘
# │
# ⋅────> ξ
mutable struct P4estMortarContainer{NDIMS, uEltype<:Real, NDIMSP1, NDIMSP3} <: AbstractContainer
  u             ::Array{uEltype, NDIMSP3} # [small/large side, variable, position, i, j, mortar]
  element_ids   ::Matrix{Int}             # [position, mortar]
  node_indices  ::Matrix{NTuple{NDIMS, Symbol}} # [small/large, mortar]

  # internal `resize!`able storage
  _u            ::Vector{uEltype}
  _element_ids  ::Vector{Int}
  _node_indices ::Vector{NTuple{NDIMS, Symbol}}
end

@inline nmortars(mortars::P4estMortarContainer) = size(mortars.element_ids, 2)
@inline Base.ndims(::P4estMortarContainer{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(mortars::P4estMortarContainer, capacity)
  @unpack _u, _element_ids, _node_indices = mortars

  n_dims = ndims(mortars)
  n_nodes = size(mortars.u, 4)
  n_variables = size(mortars.u, 2)

  resize!(_u, 2 * n_variables * 2^(n_dims-1) * n_nodes^(n_dims-1) * capacity)
  mortars.u = unsafe_wrap(Array, pointer(_u),
    (2, n_variables, 2^(n_dims-1), ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(_element_ids, (2^(n_dims-1) + 1) * capacity)
  mortars.element_ids = unsafe_wrap(Array, pointer(_element_ids),
    (2^(n_dims-1) + 1, capacity))

  resize!(_node_indices, 2 * capacity)
  mortars.node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, capacity))

  return nothing
end

# Create mortar container and initialize mortar data.
function init_mortars(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_mortars = count_required_surfaces(mesh).mortars

  _u = Vector{uEltype}(undef,
    2 * nvariables(equations) * 2^(NDIMS-1) * nnodes(basis)^(NDIMS-1) * n_mortars)
  u = unsafe_wrap(Array, pointer(_u),
    (2, nvariables(equations), 2^(NDIMS-1), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_mortars))

  _element_ids = Vector{Int}(undef, (2^(NDIMS-1) + 1) * n_mortars)
  element_ids = unsafe_wrap(Array, pointer(_element_ids), (2^(NDIMS-1) + 1, n_mortars))

  _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_mortars)
  node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_mortars))

  mortars = P4estMortarContainer{NDIMS, uEltype, NDIMS+1, NDIMS+3}(u, element_ids, node_indices,
                                                                   _u, _element_ids, _node_indices)

  if n_mortars > 0
    init_mortars!(mortars, mesh)
  end

  return mortars
end

function init_mortars!(mortars, mesh::P4estMesh)
  init_surfaces!(nothing, mortars, nothing, mesh)

  return mortars
end

function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache)
  # Re-initialize elements container
  @unpack elements = cache
  resize!(elements, ncells(mesh))
  init_elements!(elements, mesh, dg.basis)

  required = count_required_surfaces(mesh)

  # resize interfaces container
  @unpack interfaces = cache
  resize!(interfaces, required.interfaces)

  # resize boundaries container
  @unpack boundaries = cache
  resize!(boundaries, required.boundaries)

  # resize mortars container
  @unpack mortars = cache
  resize!(mortars, required.mortars)

  # re-initialize containers together to reduce
  # the number of iterations over the mesh in p4est
  init_surfaces!(interfaces, mortars, boundaries, mesh)
end

# A helper struct used in initialization methods below
mutable struct InitSurfacesIterFaceUserData{Interfaces, Mortars, Boundaries, Mesh}
  interfaces  ::Interfaces
  interface_id::Int
  mortars     ::Mortars
  mortar_id   ::Int
  boundaries  ::Boundaries
  boundary_id ::Int
  mesh        ::Mesh
end

function InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh)
  return InitSurfacesIterFaceUserData{
    typeof(interfaces), typeof(mortars), typeof(boundaries), typeof(mesh)}(
      interfaces, 1, mortars, 1, boundaries, 1, mesh)
end

function init_surfaces_iter_face(info, user_data)
  # Unpack user_data
  data = unsafe_pointer_to_objref(Ptr{InitSurfacesIterFaceUserData}(user_data))

  # Function barrier because the unpacked user_data above is type-unstable
  init_surfaces_iter_face_inner(info, data)
end

# 2D
cfunction(::typeof(init_surfaces_iter_face), ::Val{2}) = @cfunction(init_surfaces_iter_face, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))
# 3D
cfunction(::typeof(init_surfaces_iter_face), ::Val{3}) = @cfunction(init_surfaces_iter_face, Cvoid, (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid}))

# Function barrier for type stability
function init_surfaces_iter_face_inner(info, user_data)
  @unpack interfaces, mortars, boundaries = user_data

  if info.sides.elem_count == 2
    # Two neighboring elements => Interface or mortar

    # Extract surface data
    sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0
      # No hanging nodes => normal interface
      if interfaces !== nothing
        init_interfaces_iter_face_inner(info, sides, user_data)
      end
    else
      # Hanging nodes => mortar
      if mortars !== nothing
        init_mortars_iter_face_inner(info, sides, user_data)
      end
    end
  elseif info.sides.elem_count == 1
    # One neighboring elements => boundary
    if boundaries !== nothing
      init_boundaries_iter_face_inner(info, user_data)
    end
  end

  return nothing
end

function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh)
  # Let p4est iterate over all interfaces and call init_surfaces_iter_face
  iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh)))
  user_data = InitSurfacesIterFaceUserData(
    interfaces, mortars, boundaries, mesh)

  iterate_p4est(mesh.p4est, user_data; iter_face_c=iter_face_c)

  return interfaces
end

# Initialization of interfaces after the function barrier
function init_interfaces_iter_face_inner(info, sides, user_data)
  @unpack interfaces, interface_id, mesh = user_data
  user_data.interface_id += 1

  # Get Tuple of local trees, one-based indexing
  trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
           unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
  # Quadrant numbering offsets of the quadrants at this interface
  offsets = SVector(trees[1].quadrants_offset,
                    trees[2].quadrants_offset)

  local_quad_ids = SVector(sides[1].is.full.quadid, sides[2].is.full.quadid)
  # Global IDs of the neighboring quads
  quad_ids = offsets + local_quad_ids

  # Write data to interfaces container
  # p4est uses zero-based indexing; convert to one-based indexing
  interfaces.element_ids[1, interface_id] = quad_ids[1] + 1
  interfaces.element_ids[2, interface_id] = quad_ids[2] + 1

  # Face at which the interface lies
  faces = (sides[1].face, sides[2].face)

  # Save interfaces.node_indices dimension specific in containers_[23]d.jl
  init_interface_node_indices!(interfaces, faces, info.orientation, interface_id)

  return nothing
end

# Initialization of boundaries after the function barrier
function init_boundaries_iter_face_inner(info, user_data)
  @unpack boundaries, boundary_id, mesh = user_data
  user_data.boundary_id += 1

  # Extract boundary data
  side = unsafe_load_side(info)
  # Get local tree, one-based indexing
  tree = unsafe_load_tree(mesh.p4est, side.treeid + 1)
  # Quadrant numbering offset of this quadrant
  offset = tree.quadrants_offset

  # Verify before accessing is.full, but this should never happen
  @assert side.is_hanging == false

  local_quad_id = side.is.full.quadid
  # Global ID of this quad
  quad_id = offset + local_quad_id

  # Write data to boundaries container
  # p4est uses zero-based indexing; convert to one-based indexing
  boundaries.element_ids[boundary_id] = quad_id + 1

  # Face at which the boundary lies
  face = side.face

  # Save boundaries.node_indices dimension specific in containers_[23]d.jl
  init_boundary_node_indices!(boundaries, face, boundary_id)

  # One-based indexing
  boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1]

  return nothing
end

# Initialization of mortars after the function barrier
function init_mortars_iter_face_inner(info, sides, user_data)
  @unpack mortars, mortar_id, mesh = user_data
  user_data.mortar_id += 1

  # Get Tuple of local trees, one-based indexing
  trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
           unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
  # Quadrant numbering offsets of the quadrants at this interface
  offsets = SVector(trees[1].quadrants_offset,
                    trees[2].quadrants_offset)

  if sides[1].is_hanging == true
    # Left is small, right is large
    faces = (sides[1].face, sides[2].face)

    local_small_quad_ids = sides[1].is.hanging.quadid
    # Global IDs of the two small quads
    small_quad_ids = offsets[1] .+ local_small_quad_ids

    # Just be sure before accessing is.full
    @assert sides[2].is_hanging == false
    large_quad_id = offsets[2] + sides[2].is.full.quadid
  else # sides[2].is_hanging == true
    # Right is small, left is large.
    # init_mortar_node_indices! below expects side 1 to contain the small elements.
    faces = (sides[2].face, sides[1].face)

    local_small_quad_ids = sides[2].is.hanging.quadid
    # Global IDs of the two small quads
    small_quad_ids = offsets[2] .+ local_small_quad_ids

    # Just be sure before accessing is.full
    @assert sides[1].is_hanging == false
    large_quad_id = offsets[1] + sides[1].is.full.quadid
  end

  # Write data to mortar container, 1 and 2 are the small elements
  # p4est uses zero-based indexing; convert to one-based indexing
  mortars.element_ids[1:end-1, mortar_id] .= small_quad_ids[:] .+ 1
  # Last entry is the large element
  mortars.element_ids[end, mortar_id] = large_quad_id + 1

  init_mortar_node_indices!(mortars, faces, info.orientation, mortar_id)

  return nothing
end

# Iterate over all interfaces and count
# - (inner) interfaces
# - mortars
# - boundaries
# and collect the numbers in `user_data` in this order.
function count_surfaces_iter_face(info, user_data)
  if info.sides.elem_count == 2
    # Two neighboring elements => Interface or mortar

    # Extract surface data
    sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0
      # No hanging nodes => normal interface
      # Unpack user_data = [interface_count] and increment interface_count
      ptr = Ptr{Int}(user_data)
      id = unsafe_load(ptr, 1)
      unsafe_store!(ptr, id + 1, 1)
    else
      # Hanging nodes => mortar
      # Unpack user_data = [mortar_count] and increment mortar_count
      ptr = Ptr{Int}(user_data)
      id = unsafe_load(ptr, 2)
      unsafe_store!(ptr, id + 1, 2)
    end
  elseif info.sides.elem_count == 1
    # One neighboring elements => boundary

    # Unpack user_data = [boundary_count] and increment boundary_count
    ptr = Ptr{Int}(user_data)
    id = unsafe_load(ptr, 3)
    unsafe_store!(ptr, id + 1, 3)
  end

  return nothing
end

# 2D
cfunction(::typeof(count_surfaces_iter_face), ::Val{2}) = @cfunction(count_surfaces_iter_face, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))
# 3D
cfunction(::typeof(count_surfaces_iter_face), ::Val{3}) = @cfunction(count_surfaces_iter_face, Cvoid, (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid}))

function count_required_surfaces(mesh::P4estMesh)
  # Let p4est iterate over all interfaces and call count_surfaces_iter_face
  iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh)))

  # interfaces, mortars, boundaries
  user_data = [0, 0, 0]

  iterate_p4est(mesh.p4est, user_data; iter_face_c=iter_face_c)

  # Return counters
  return (interfaces = user_data[1],
          mortars    = user_data[2],
          boundaries = user_data[3])
end

# Return direction of the face, which is indexed by node_indices
@inline function indices2direction(indices)
  if indices[1] === :begin
    return 1
  elseif indices[1] === :end
    return 2
  elseif indices[2] === :begin
    return 3
  elseif indices[2] === :end
    return 4
  elseif indices[3] === :begin
    return 5
  else # if indices[3] === :end
    return 6
  end
end

include("containers_2d.jl")
include("containers_3d.jl")

end # @muladd


text:
using Documenter, GRUtils

include("api_page.jl")

makedocs(sitename="GRUtils.jl",
    pages = [
        "Introduction" => "index.md",
        "Multiple plots" => "multiple.md",
        "Animations" => "animations.md",
        "Color management" => "color.md",
        "API reference" => [
            "api/plotting.md",
            "api/attributes.md",
            "api/control.md"
        ],
        "Internals" => [
            "internals/structure.md",
            "internals/drawing.md",
            "internals/createplots.md",
            "internals/extending.md"
        ]
    ],
    expandfirst = ["index.md"]
)


text:
function kv!(r::Robot)# - главная функция
    for side in (Nord, West, Sud, Ost) # - перебор всех возможных направлений
        move_till_border!(r, side)
        putmarkers_border!(r, side)
    end
    show!(r)
end

function move_till_border!(r::Robot, side::HorizonSide)
    while isborder(r, Nord) == false
        move!(r, Nord)
    end
end

function putmarkers_border!(r::Robot, side::HorizonSide)
    for side in (Ost, Sud, West, Nord)
        while isborder(r, side) == false
            putmarker!(r)
            move!(r, side)
        end
    end
end


text:
# in this attempt, I shall play a set of games with the policy and move the policy using
# the data
using Game2048
using Flux
using Flux: logitcrossentropy, crossentropy
using StatsBase: wsample, sample, countmap
using Statistics: mean, std
using Base.Threads: @threads

# based on testing gpu is slower than cpu, so forget that idea for now
# see benchmarks/gpu_slower_than_cpu.jl
# using CUDA

const γ = 0.9 # the reward discount rate
const α = 0.001 # learning rate

function play_game_with_policy(policy; verbose=false, epsilon=0, greedy=false, min_remarkable=0.3)
    board = initboard()

    histories = Array{Int8, 2}[]
    dir_histories = Dirs[]
    possibles = Vector{Bool}[]
    max_prob = 0
    remarkable=false
    rewards = Int[]
    while true
        res = move.(Ref(board), DIRS)
        possible = [p for (_, (_, p)) in res]
        #possible::Vector{Bool} = [p for (_, (_, p)) in res]

        if all(!, possible)
            return histories, dir_histories, possibles, rewards, board
        end

        choices = (1:4)[possible]

        if greedy
            epsed = false
            prob = policy(board)[possible]
            remarkable = maximum(prob) > min_remarkable
            if maximum(prob) > max_prob
                max_prob = maximum(prob)
                # remarkable = true
            end
            if any(isnan, prob)
                println(prob)
                error("wtf")
            end
            # sample a direction from policy
            _, ii = findmax(prob)
            i = choices[ii]
        elseif rand() < epsilon
            epsed = true
            i = rand(choices)
        else
            epsed = false
            prob = policy(board)[possible]
            remarkable = maximum(prob) > min_remarkable
            if maximum(prob) > max_prob
                max_prob = maximum(prob)
                # remarkable = true
            end
            if any(isnan, prob)
                println(prob)
                error("wtf")
            end
            # sample a direction from policy
            i = wsample(choices, prob)
        end

        dir = DIRS[i]
        if verbose && remarkable
            println("")
            println("")
            println("****************************begin: one move")
            display(board)
            println(dir)
            println(DIRS[choices])
            if epsed
                "epsilon'ed"
            else
                println(prob)
            end
        end
        push!(histories, board)
        push!(dir_histories, dir)
        push!(possibles, possible)
        push!(rewards, res[i][2][1])

        board = res[i][1]

        if verbose && remarkable
            display(board)
            println("--------------------------------end: one move")
            println("")
            println("")
        end
        remarkable = false
        gen_new_tile!(board)
    end
end

function normalise_scale_rewards!(rewards)
    mean_lvl = mean(rewards)
    rewards .= (rewards .- mean_lvl) ./ std(rewards)
    rewards
end

function many_games(policy; n=20, kwargs...)
    @assert n > 3
    rewards = Float32[]
    X = Vector{Array{Int8, 4}}()
    dirs = []
    l = Int[]
    p = Vector{Bool}[]
    for i in 1:n
        # global X, dirs, l, p, rewards
        # if iseven(i)
            histories, dir_histories, possibles, r, last_board = play_game_with_policy(policy; kwargs...)
        # else
        #     histories, dir_histories, possibles, r, last_board = play_game_with_policy(policy; greedy = true, kwargs...)
        # end
        @assert length(histories) == length(r)

        X = vcat(X, histories)
        dirs = vcat(dirs, dir_histories)
        push!(l, length(dir_histories))
        p = vcat(p, possibles)

        # compute the discounted rewards
        tmp_r = Vector{Float32}(undef, length(r))
        tmp_r[end] = log(1+r[end])
        for i in length(r)-1:-1:1
            tmp_r[i] = log(1+r[i]) + tmp_r[i+1]*γ
        end
        rewards = vcat(rewards, tmp_r)
    end

    # rewards .= rewards ./ 10

    # mean(rewards[vcat([1], cumsum(l)[1:end-1] .+ 1)] .- mean(rewards))/std(rewards)
    # normalise the rewards
    # the rewards have to be at least rewards
    normalise_scale_rewards!(rewards)
    # println(maximum(rewards))

    @assert length(X) == length(dirs) == length(rewards)
    actual_history = zip(X, dirs, rewards) |> collect

    # for histories with no possible directions
    @assert length(X) == length(p)
    data_impossible = []
    for (h, p1) in zip(X, p)
        for (dir, p1) in zip(DIRS, p1)
            if !p1
                # make that direction impossible!
                push!(data_impossible, (h, dir, -3))
            end
        end
    end

    v1 = vcat(actual_history, data_impossible)

    # now we need to rotate the game
    v2 = eltype(v1)[]
    for (board, dir, reward) in v1
        for (new_board, new_dir) in rotate_mirror(board, dir)
            push!(v2, (new_board, new_dir, reward))
        end
    end
    println("ok")
    v2
end

function showmaxtile(policy; n=100, kwargs...)
    res = Vector{Int}(undef, n)
    @threads for i in 1:n
        res[i] = play_game_with_policy(policy; verbose=false, kwargs...)[end] |> maximum
    end
    countmap(res) |> sort
end

function update_param!(ps, board, dir, reward)
    grad = gradient(()-> reward*log(1+policy(board)[Int(dir)+1]), ps)
    for p1 in ps
        Flux.update!(p1, -α*grad[p1])
    end
end

policy = Chain(
    board -> Float32.(board),
    board -> board ./ maximum(board),
    board -> reshape(board, 16),
    Dense(16, 20, relu),
    IdentitySkip(Dense(20, 20, relu), identity),
    IdentitySkip(Dense(20, 20, relu), identity),
    Dense(20, 4),
    softmax
)

ps = Flux.params(policy)

if false
    # show change in policy
    board = initboard()
    policy(board)
    update_param!(ps, board, left, 3)
    policy(board)
    update_param!(ps, board, left, -3)
    policy(board)

    board = initboard()
    board
    function doit()
        display(board)
        pb=policy(board)
        dir = wsample(DIRS, pb)
        move!(board, dir);
        gen_new_tile!(board)
        println(dir)
        println(pb)
        display(board)
    end

    doit()
end

# display(showmaxtile(policy; n=10000))
# random
# OrderedCollections.OrderedDict{Int64,Int64} with 6 entries:
#   4 => 31
#   5 => 666
#   6 => 3651
#   7 => 4817
#   8 => 832
#   9 => 3
display(showmaxtile(policy; n=100))

if false
    i= 1

    while !any(isnan, policy(board))
        global i
        board, dir, reward = Xy[i]
        println(policy(board))
        # display(board)
        dir
        update_param!(ps, board, dir, reward)
        println(policy(board))
        i+=1
    end
end

function ok()
    batch=1
    while true
        println("batch $batch")
        @time Xy = many_games(policy; n=1000, epsilon=0.1, greedy = false);
        Xy = sample(Xy, length(Xy), replace=false)
        for (board, dir, reward) in Xy
            policy(board)
            update_param!(ps, board, dir, reward)
        end
        play_game_with_policy(policy, verbose=true, min_remarkable=0.3)
        display(showmaxtile(policy; greedy=true))
        display(showmaxtile(policy; greedy=false))

        batch += 1
    end
end

ok()


text:
export simplify!

function simulated_annealing(f, x0, temp)
    l_old = f(x0)
    loss = l_old
    @show loss
    candidate = copy(x0)
    for epoch in 1:1000
        # local update
        idx = rand(1:length(candidate))
        candidate_idx_old = candidate[idx]
        candidate[idx] = rand(rewrite_rules)

        l_new = f(candidate)
        t = temp / epoch
        ratio = exp((l_old - l_new)/t)
        if l_new < l_old || rand() < ratio
            # pick candidate
        else
            candidate[idx] = candidate_idx_old
        end
        loss = l_old
        @show loss
        l_old = l_new
    end
    return candidate
end

function simplify_instance(tait) # -> cost
    return function instance(rules::Vector{Symbol}) # -> cost
        for rule in rules
            matches = match(Rule(rule), tait)
            isempty(matches) && continue
            m = first(matches)
            # @show rule, m
            rewrite!(tait, m)
        end
        return tait_loss(tait)
    end
end

function tait_loss(q::Tait)
    return length(q.quon_params)/2 + length(q.genuses) * 10
end

function simplify!(tait::Tait)
    x0 = rand(rewrite_rules, 10)
    instance = simplify_instance(tait)
    rules = simulated_annealing(instance, x0, 1.0)
    return tait
end

function simplify!(q::Tait, r::Rule{R}) where {R}
    ms = match(r, q)
    while length(ms) > 0
        for m in ms
            if check(q, m)
                rewrite!(q, m)
            end
        end
        ms = match(r, q)
    end
    return q
end


text:
# Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
# Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
# pp. 829–858 (2017).

using Printf
using ProximalCore: Zero
using LinearAlgebra
using Printf

"""
    DavisYinIteration(; <keyword-arguments>)

Iterator implementing the Davis-Yin splitting algorithm [1].

This iterator solves convex optimization problems of the form

    minimize f(x) + g(x) + h(x),

where `f` is smooth.

See also [`DavisYin`](@ref).

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `g=Zero()`: proximable objective term.
- `h=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of h.
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).

# References
1. Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization Applications", Set-Valued and Variational Analysis, vol. 25, no. 4, pp. 829-858 (2017).
"""
Base.@kwdef struct DavisYinIteration{R,C<:Union{R,Complex{R}},T<:AbstractArray{C},Tf,Tg,Th}
    f::Tf = Zero()
    g::Tg = Zero()
    h::Th = Zero()
    x0::T
    lambda::R = real(eltype(x0))(1)
    Lf::Maybe{R} = nothing
    gamma::Maybe{R} = Lf !== nothing ? (1 / Lf) : error("You must specify either Lf or gamma")
end

Base.IteratorSize(::Type{<:DavisYinIteration}) = Base.IsInfinite()

struct DavisYinState{T}
    z::T
    xg::T
    grad_f_xg::T
    z_half::T
    xh::T
    res::T
end

function Base.iterate(iter::DavisYinIteration)
    z = copy(iter.x0)
    xg, = prox(iter.g, z, iter.gamma)
    grad_f_xg, = gradient(iter.f, xg)
    z_half = 2 .* xg .- z .- iter.gamma .* grad_f_xg
    xh, = prox(iter.h, z_half, iter.gamma)
    res = xh - xg
    z .+= iter.lambda .* res
    state = DavisYinState(z, xg, grad_f_xg, z_half, xh, res)
    return state, state
end

function Base.iterate(iter::DavisYinIteration, state::DavisYinState)
    prox!(state.xg, iter.g, state.z, iter.gamma)
    gradient!(state.grad_f_xg, iter.f, state.xg)
    state.z_half .= 2 .* state.xg .- state.z .- iter.gamma .* state.grad_f_xg
    prox!(state.xh, iter.h, state.z_half, iter.gamma)
    state.res .= state.xh .- state.xg
    state.z .+= iter.lambda .* state.res
    return state, state
end

default_stopping_criterion(tol, ::DavisYinIteration, state::DavisYinState) = norm(state.res, Inf) <= tol
default_solution(::DavisYinIteration, state::DavisYinState) = state.xh
default_display(it, ::DavisYinIteration, state::DavisYinState) = @printf("%5d | %.3e\n", it, norm(state.res, Inf))

"""
    DavisYin(; <keyword-arguments>)

Constructs the Davis-Yin splitting algorithm [1].

This algorithm solves convex optimization problems of the form

    minimize f(x) + g(x) + h(x),

where `f` is smooth.

The returned object has type `IterativeAlgorithm{DavisYinIteration}`,
and can be called with the problem's arguments to trigger its solution.

See also: [`DavisYinIteration`](@ref), [`IterativeAlgorithm`](@ref).

# Arguments
- `maxit::Int=10_000`: maximum number of iteration
- `tol::1e-8`: tolerance for the default stopping criterion
- `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
- `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
- `verbose::Bool=false`: whether the algorithm state should be displayed
- `freq::Int=100`: every how many iterations to display the algorithm state
- `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
- `kwargs...`: additional keyword arguments to pass on to the `DavisYinIteration` constructor upon call

# References
1. Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization Applications", Set-Valued and Variational Analysis, vol. 25, no. 4, pp. 829–858 (2017).
"""
DavisYin(;
    maxit=10_000,
    tol=1e-8,
    stop=(iter, state) -> default_stopping_criterion(tol, iter, state),
    solution=default_solution,
    verbose=false,
    freq=100,
    display=default_display,
    kwargs...
) = IterativeAlgorithm(DavisYinIteration; maxit, stop, solution, verbose, freq, display, kwargs...)


text:
function zlb_rule()
    AltPolicy(:zlb_rule, zlb_rule_eqcond, zlb_rule_solve,
              forecast_init = zlb_rule_forecast_init,
              color = RGB(0., 0., 0.5430)) # dark blue
end

function zlb_rule_replace_eq_entries(m::AbstractDSGEModel,
                                      Γ0::Matrix{Float64}, Γ1::Matrix{Float64},
                                      C::Vector{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64}, reg::Int)

    eq             = m.equilibrium_conditions
    endo           = m.endogenous_states
    exo            = m.exogenous_shocks

    # Zero out old MP rule
    Γ0[eq[:eq_mp], :] .= 0.
    Γ1[eq[:eq_mp], :] .= 0.

    # Implement zero rate
    Γ0[eq[:eq_mp], endo[:R_t]] = 1.0
    # NOTE: For now, let's do it via settings rather than using regime-switching parameters
    #       (to avoid adding too many extra dimensions to parameter output)
    C[eq[:eq_mp]] = (haskey(get_settings(m), :zlb_rule_values) ? get_setting(m, :zlb_rule_values)[reg] / 4. :
                     (haskey(get_settings(m), :zlb_rule_value) ? get_setting(m, :zlb_rule_value) / 4. :
                      forecast_zlb_value(m))) - m[:Rstarn]

    # Default behavior of ZLB rule is to remove anticipated monetary policy shocks,
    # but this setting can turn it off.
    if !haskey(get_settings(m), :zlb_rule_remove_mon_anticipated_shocks) ||
        get_setting(m, :zlb_rule_remove_mon_anticipated_shocks)
        Γ1[eq[:eq_rm], endo[:rm_tl1]]   = 0.
        Γ0[eq[:eq_rml1], endo[:rm_tl1]] = 1.
        Ψ[eq[:eq_rml1], exo[:rm_shl1]]  = 0.

        if n_mon_anticipated_shocks(m) > 1
            for i in 2:n_mon_anticipated_shocks(m)
                Γ1[eq[Symbol("eq_rml$(i-1)")], endo[Symbol("rm_tl$i")]] = 0.
                Γ0[eq[Symbol("eq_rml$i")], endo[Symbol("rm_tl$i")]]     = 1.
                Ψ[eq[Symbol("eq_rml$i")], exo[Symbol("rm_shl$i")]]      = 0.
            end
        end
    end

    return Γ0, Γ1, C, Ψ, Π
end

"""
```
zlb_rule_eqcond(m::AbstractDSGEModel)
```

Solves for the transition equation of `m` under a price level
targeting rule (implemented by adding a price-gap state)
"""
function zlb_rule_eqcond(m::AbstractDSGEModel, reg::Int = 1)

    # Get equilibrium condition matrices
    Γ0, Γ1, C, Ψ, Π = eqcond(m, reg)
    Γ0, Γ1, C, Ψ, Π = zlb_rule_replace_eq_entries(m, Γ0, Γ1, C, Ψ, Π, reg)

    return Γ0, Γ1, C, Ψ, Π
end

"""
```
zlb_rule_solve(m::AbstractDSGEModel)
```

Solves for the transition equation of `m` under a price level
targeting rule (implemented by adding a price-gap state)
"""
function zlb_rule_solve(m::AbstractDSGEModel; regime_switching::Bool = false, regimes::Vector{Int} = Int[1])

    # Get equilibrium condition matrices
    if length(regimes) == 1
        Γ0, Γ1, C, Ψ, Π  = zlb_rule_eqcond(m, regimes[1])

        TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6, verbose = :low)

        # Check for LAPACK exception, existence and uniqueness
        if eu[1] != 1 || eu[2] != 1
            throw(GensysError())
        end

        TTT_gensys = real(TTT_gensys)
        RRR_gensys = real(RRR_gensys)
        CCC_gensys = real(CCC_gensys)

        # Augment states
        TTT, RRR, CCC = DSGE.augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys; regime_switching = regime_switching,
                                            reg = regimes[1])
        return TTT, RRR, CCC
    else
        Γ0s = Vector{Matrix{Float64}}(undef, length(regimes))
        Γ1s = Vector{Matrix{Float64}}(undef, length(regimes))
        Cs = Vector{Vector{Float64}}(undef, length(regimes))
        Ψs = Vector{Matrix{Float64}}(undef, length(regimes))
        Πs = Vector{Matrix{Float64}}(undef, length(regimes))
        for reg in regimes
            Γ0s[reg], Γ1s[reg], Cs[reg], Ψs[reg], Πs[reg]  = zlb_rule_eqcond(m, reg)
        end

        n_regimes = length(regimes)
        TTTs= Vector{Matrix{Float64}}(undef, n_regimes)
        RRRs = Vector{Matrix{Float64}}(undef, n_regimes)
        CCCs = Vector{Vector{Float64}}(undef, n_regimes)

        # Solve model
        for reg in regimes
            TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0s[reg], Γ1s[reg], Cs[reg], Ψs[reg], Πs[reg], 1+1e-6)

            if !((eu[1] == 1) & (eu[2] == 1))
                throw(GensysError("Gensys does not give existence"))
            end
            TTT_gensys = real(TTT_gensys)
            RRR_gensys = real(RRR_gensys)
            CCC_gensys = reshape(CCC_gensys, size(CCC_gensys, 1))

            # Augment states
            TTTs[reg], RRRs[reg], CCCs[reg] = DSGE.augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys;
                                                                  regime_switching = regime_switching,
                                                                  reg = reg)
        end
    end

    return TTTs, RRRs, CCCs
end

"""
```
init_zlb_rule_forecast(m::AbstractDSGEModel, shocks::Matrix{T}, final_state::Vector{T})
```

Adjust shocks matrix and final state vector for forecasting under the ZLB_RULE rule
"""
function zlb_rule_forecast_init(m::AbstractDSGEModel, shocks::Matrix{T}, final_state::Vector{T};
                                 cond_type::Symbol = :none) where {T<:AbstractFloat}
    return shocks, final_state
end


text:
import ClassImbalance
import Random
import Test

Random.seed!(999) # seed the global random number generator

include("install_r.jl")

include("install_dmwr.jl")

Test.@testset "miscellaneous.jl" begin
    include("miscellaneous.jl")
end

Test.@testset "smote_example.jl" begin
    include("smote_example.jl")
end

Test.@testset "test_rose.jl" begin
    include("test_rose.jl")
end

Test.@testset "test_smote_exs.jl" begin
    include("test_smote_exs.jl")
end

Test.@testset "test_utils.jl" begin
    include("test_utils.jl")
end


text:
using HDF5

#=

=#

function to_dict(hd5, key)
    list = read(hd5, key)
    Dict(d => i for (i, d) in enumerate(list))
end

struct Data
    file
    day_key
    stat_key
end

function load_data(file_name)
    file = h5open(file_name)
    day_keys = to_dict(file, "day-order")
    stat_keys = to_dict(file, "stat-order")
    Data(file, day_keys, stat_keys)
end

function get_data(data::Data, line, day, stat)
    d_i = data.day_key[day]
    s_i = data.stat_key[stat]
    mat = read(data.file, "traffic/$line")
    mat[:, d_i, s_i]
end

function close_data(data::Data)
    close(data.file)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# # src/eago_semiinfinite/types.jl
# Defines intermediate types used by SIP routines.
#############################################################################

abstract type AbstractSIPAlgo end

abstract type AbstractSubproblemType end
struct LowerLevel1 <: AbstractSubproblemType end
struct LowerLevel2 <: AbstractSubproblemType end
struct LowerLevel3 <: AbstractSubproblemType end
struct LowerProblem <: AbstractSubproblemType end
struct UpperProblem <: AbstractSubproblemType end
struct ResProblem <: AbstractSubproblemType end

for ptype in (LowerLevel1, LowerLevel2, LowerLevel3,
                  LowerProblem, UpperProblem, ResProblem)
    pstring = String(Symbol(ptype))[6:end]
    @eval Base.show(io::IO, x::$ptype) = print(io, $pstring)
    @eval Base.show(io::IO, m::MIME"text/plain", x::$ptype) = print(io, $pstring)
end

"""
    SIPResult

Structure storing the results of the SIPres algorithm.
"""
mutable struct SIPResult
    iteration_number::Int64
    res_iteration_number::Int64
    upper_bound::Float64
    lower_bound::Float64
    feasibility::Bool
    xsol::Vector{Float64}
    psol::Vector{Float64}
    solution_time::Float64
end
SIPResult() = SIPResult(1, 1, Inf, -Inf, true, Float64[], Float64[], 0.0)
SIPResult(nx::Int, np::Int) = SIPResult(1, 1, Inf, -Inf, true, zeros(nx), zeros(np), 0.0)

"""
      SIPProblem

Structure storing problem information for the solution routine.
"""
Base.@kwdef mutable struct SIPProblem
    x_l::Vector{Float64}            = Float64[]
    x_u::Vector{Float64}            = Float64[]
    p_l::Vector{Float64}            = Float64[]
    p_u::Vector{Float64}            = Float64[]
    np::Int                         = 0
    nSIP::Int                       = 0
    nx::Int                         = 0

    abs_tolerance::Float64     = 1E-3
    cons_tolerance::Float64   = 1E-3
    iteration_limit::Int            = 100
    res_iteration_limit::Int        = 3
    initial_eps_g::Float64          = 1.0
    initial_r::Float64              = 2.0
    local_solver::Bool              = false

    return_hist::Bool               = false
    header_interval::Int            = 20
    print_interval::Int             = 1
    verbosity::Int                  = 1

    kwargs_llp1::Dict{String,Any}   = Dict{String,Any}()
    kwargs_llp2::Dict{String,Any}   = Dict{String,Any}()
    kwargs_llp3::Dict{String,Any}   = Dict{String,Any}()
    kwargs_lbd::Dict{String,Any}    = Dict{String,Any}()
    kwargs_ubd::Dict{String,Any}    = Dict{String,Any}()
    kwargs_res::Dict{String,Any}    = Dict{String,Any}()
end

get_sip_kwargs(s::LowerLevel1, p::SIPProblem) = p.kwargs_llp1
get_sip_kwargs(s::LowerLevel2, p::SIPProblem) = p.kwargs_llp2
get_sip_kwargs(s::LowerLevel3, p::SIPProblem) = p.kwargs_llp3
get_sip_kwargs(s::LowerProblem, p::SIPProblem) = p.kwargs_lbd
get_sip_kwargs(s::UpperProblem, p::SIPProblem) = p.kwargs_ubd
get_sip_kwargs(s::ResProblem, p::SIPProblem) = p.kwargs_res

function SIPProblem(x_l::Vector{Float64}, x_u::Vector{Float64},
                    p_l::Vector{Float64}, p_u::Vector{Float64},
                    gSIP, kwargs)

    prob = SIPProblem()

    for key in keys(kwargs)
        string_key = String(key)
        if string_key[1:5] === "llp1_"
            prob.kwargs_llp1[string_key[6:end]] = kwargs[key]
        elseif string_key[1:5] === "llp2_"
            prob.kwargs_llp2[string_key[6:end]] = kwargs[key]
        elseif string_key[1:5] === "llp3_"
            prob.kwargs_llp3[string_key[6:end]] = kwargs[key]
        elseif string_key[1:4] === "lbd_"
            prob.kwargs_lbd[string_key[5:end]] = kwargs[key]
        elseif string_key[1:4] === "ubd_"
            prob.kwargs_ubd[string_key[5:end]] = kwargs[key]
        elseif string_key[1:4] === "res_"
            prob.kwargs_res[string_key[5:end]] = kwargs[key]
        elseif key in fieldnames(SIPProblem)
            setfield!(prob, key, kwargs[key])
        else
            error("Keyword = `$string_key` is not supported by SIPProblem. Use
                  `llp1_$string_key` to pass this argument to the optimizer used by
                  the first lower level problem, `lbd_$string_key` to pass the
                  keyword to the optimizer used in the lower bounding problem
                  and so on...")
        end
    end

    (prob.initial_r <= 1.0) && error("initial_r must be greater than 1")
    (prob.initial_eps_g <= 0.0) && error("eps_g must be greater than 0")

    prob.np = length(p_l)
    prob.nx = length(x_l)
    prob.nSIP = length(gSIP)
    append!(prob.x_l, x_l)
    append!(prob.x_u, x_u)
    append!(prob.p_l, p_l)
    append!(prob.p_u, p_u)

    return prob
end

struct SIPCallback
    f
    gSIP
end

mutable struct SubProblemInfo
    sol::Vector{Float64}
    obj_val::Float64
    obj_bnd::Float64
    feas::Bool
    tol::Vector{Float64}
end
function SubProblemInfo(nd::Int, ng::Int, tol::Float64)
    SubProblemInfo(zeros(nd), 0.0, 0.0, false, fill(tol, ng))
end

"""
    SIPBuffer

Hold objective value, solution, discretization set, and feasibility status of
each subproblem encountered by SIP algorithm.
"""
Base.@kwdef mutable struct SIPSubResult
    fRes::Float64 = Inf
    lbd::SubProblemInfo
    ubd::SubProblemInfo
    res::SubProblemInfo
    llp1::SubProblemInfo
    llp2::SubProblemInfo
    llp3::SubProblemInfo
    r_g::Float64 = 5.0
    r_l::Float64 = 5.0
    r_u::Float64 = 5.0
    eps_g::Vector{Float64} = Float64[]
    eps_l::Vector{Float64} = Float64[]
    eps_u::Vector{Float64} = Float64[]
    disc_l_buffer::Vector{Float64} = Float64[]
    disc_u_buffer::Vector{Float64} = Float64[]
    disc_l::Vector{Vector{Vector{Float64}}} = Vector{Vector{Float64}}[]
    disc_u::Vector{Vector{Vector{Float64}}} = Vector{Vector{Float64}}[]
end
function SIPSubResult(nx::Int, np::Int, ng::Int, tol::Float64)
    buffer = SIPSubResult(lbd  = SubProblemInfo(nx, 1, tol),
                          ubd  = SubProblemInfo(nx, 1, tol),
                          res  = SubProblemInfo(nx, 1, tol),
                          llp1 = SubProblemInfo(np, ng, tol),
                          llp2 = SubProblemInfo(np, ng, tol),
                          llp3 = SubProblemInfo(np, ng, tol))
    append!(buffer.eps_g, fill(1E-3, ng))
    append!(buffer.eps_l, fill(1E-3, ng))
    append!(buffer.eps_u, fill(1E-3, ng))
    append!(buffer.disc_l_buffer, zeros(np))
    append!(buffer.disc_u_buffer, zeros(np))
    for _ in 1:ng
        push!(buffer.disc_l, Vector{Float64}[])
        push!(buffer.disc_u, Vector{Float64}[])
    end
    return buffer
end

get_eps(::LowerProblem, sr::SIPSubResult, i::Int) = 0.0
get_eps(::ResProblem, sr::SIPSubResult, i::Int) = 0.0
get_eps(::UpperProblem, sr::SIPSubResult, i::Int) = sr.eps_g[i]

const SUBPROB_SYM = Dict{Symbol,Symbol}(:LowerProblem => :lbd,
                                        :UpperProblem => :ubd,
                                        :ResProblem   => :res,
                                        :LowerLevel1  => :llp1,
                                        :LowerLevel2  => :llp2,
                                        :LowerLevel3  => :llp3)

for (typ, fd) in SUBPROB_SYM
    @eval function load!(::$typ, sr::SIPSubResult, feas::Bool, val::Float64,
                                 bnd::Float64, x::Vector{Float64})

        sr.$fd.feas = feas
        sr.$fd.obj_val = val
        sr.$fd.obj_bnd = bnd
        sr.$fd.sol .= x
        return nothing
    end
end


text:
###############################################################################
#
# psofunction
#
###############################################################################
"""
    psofunction(K::FT,
                k::FT,
                Ω::FT,
                LAI::FT,
                q::FT,
                dso::FT,
                xl::FT
    ) where {FT<:AbstractFloat}

# TODO explain the variables
Return the probability of observing a sunlit leaf at depth `xl` (`pso`, see eq
    31 in vdT 2009), given
- `xl` Leaf depth in the canopy
"""
function psofunction(
            K::FT,
            k::FT,
            Ω::FT,
            LAI::FT,
            q::FT,
            dso::FT,
            xl::FT
) where {FT<:AbstractFloat}
    # [nl+1]  factor for correlation of Ps and Po
    if dso != 0
        alf = dso / q * 2 / (k + K);

        return exp( (K+k)*Ω*LAI*xl +
                    sqrt(K*k)*Ω*LAI / alf * (1 - exp(xl*alf)) )
    else
        return exp( (K + k -sqrt(K*k)) * Ω * LAI * xl )
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
"""
    perk([rng], y, T, xDists, noiseDist, signalModels, kernel, ρ)
    perk([rng], y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

Train PERK and then estimate latent parameters.

# Arguments
- `rng::AbstractRNG = Random.GLOBAL_RNG`: Random number generator to use
- `y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}}`:
  Test data points [D,N] or \\[N\\] (if D = 1) or scalar (if D = N = 1)
- `ν::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Known
  parameters [K,N] or \\[N\\] (if K = 1) or scalar (if K = N = 1); omit this
  parameter if K = 0
- `T::Integer`: Number of training points
- `xDists`: Distributions of latent parameters \\[L\\] or scalar (if L = 1);
  `xDists` can be any object such that `rand(xDists, ::Integer)` is defined (or
  a collection of such objects)
- `νDists`: Distributions of known parameters \\[K\\] or scalar (if K = 1);
  `νDists` can be any object such that `rand(νDists, ::Integer)` is defined (or
  a collection of such objects); omit this parameter if K = 0
- `noiseDist`: Distribution of noise (assumes same noise distribution for both
  real and imaginary channels in complex case); `noiseDist` can be any object
  such that `rand(noiseDist, ::Integer)` is defined
- `signalModels::Union{<:Function,<:AbstractVector{<:Function}}`: Signal models
  used to generate noiseless data \\[numSignalModels\\]; each signal model
  accepts as inputs L latent parameters (scalars) first, then K known parameters
  (scalars); user-defined parameters (e.g., scan parameters in MRI) should be
  built into the signal model
- `kernel::Kernel`: Kernel to use
- `ρ::Real`: Tikhonov regularization parameter

## Note
- D is the combined number of outputs from all signal models (e.g., number of
  scans in MRI)
- K is the number of known parameters
- L is the number of unknown or latent parameters to be estimated
- N is the number of test points (e.g., number of voxels in MRI)

# Return
- `xhat::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Estimated latent parameters [L,N] or \\[N\\] (if L = 1) or \\[L\\] (if N = 1)
  or scalar (if L = N = 1)
"""
function perk(
    rng::AbstractRNG,
    y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}},
    T::Integer,
    xDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}},
    kernel::Kernel,
    ρ::Real
)

    trainData = train(rng, T, xDists, noiseDist, signalModels, kernel, ρ)

    return perk(y, trainData, kernel)

end

function perk(
    rng::AbstractRNG,
    y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}},
    ν::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    T::Integer,
    xDists,
    νDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}},
    kernel::Kernel,
    ρ::Real
)

    trainData = train(rng, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    return perk(y, ν, trainData, kernel)

end

function perk(
    y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}},
    args...
)

    return perk(Random.GLOBAL_RNG, y, args...)

end

"""
    perk(y, trainData, kernel)
    perk(y, ν, trainData, kernel)

Estimate latent parameters using the provided training data.

# Arguments
- `y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}}`:
  Test data points [D,N] or \\[N\\] (if D = 1) or scalar (if D = N = 1)
- `ν::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Known
  parameters [K,N] or \\[N\\] (if K = 1) or scalar (if K = N = 1); omit this
  parameter if K = 0
- `trainData::TrainingData`: Training data
- `kernel::Kernel`: Kernel to use

# Return
- `xhat::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Estimated latent parameters [L,N] or \\[N\\] (if L = 1) or \\[L\\] (if N = 1)
  or scalar (if L = N = 1)
"""
function perk(
    y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}},
    trainData::TrainingData,
    kernel::Kernel
)

    eltype(y) <: Complex && (y = complex2real(y))

    return krr(y, trainData, kernel)

end

function perk(
    y::Union{<:Number,<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}},
    ν::Union{<:Real,<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    trainData::TrainingData,
    kernel::Kernel
)

    eltype(y) <: Complex && (y = complex2real(y))

    q = combine(y, ν)

    return krr(q, trainData, kernel)

end


text:
function _bioclim_score(x)
    if isnothing(x)
        return nothing
    end
    @assert 0. <= x <= 1.
    if x >= 0.5
        return 2(1-x)
    else
        return 2x
    end
end

function bioclim(layer::T, records::GBIFRecords) where {T <: SimpleSDMLayer}
    bioclim_prediction = convert(SimpleSDMResponse, copy(layer))
    qf = StatsBase.ecdf(Float32.(filter(!isnothing, layer[records])))
    qfgrid = replace(x -> isnothing(x) ? nothing : qf(x), layer.grid)
    bioclim_prediction.grid = _bioclim_score.(qfgrid)
    return bioclim_prediction
end


text:
#!/usr/bin/env julia

function fact(n::Int)::Int64
    result = 1
    for i in 2:n
        result *= i
    end
    return result
end

for n in 0:5
    println("fact($n) = $(fact(n))")
end


text:
import Chipmunk
Cp = Chipmunk
import SFML
Sf = SFML

cd(dirname(Base.source_path()))

function sf_vec(vect)
    Sf.Vector2f(vect.x, -vect.y)
end

settings = Sf.ContextSettings()
settings.antialiasing_level = 3
window = Sf.RenderWindow(Sf.VideoMode(800, 600), "Physics", settings, Sf.window_defaultstyle)
Sf.set_framerate_limit(window, 60)

event = Sf.Event()

circles = Sf.CircleShape[]
bodies = Cp.Body[]

gravity = Cp.Vect(0, -100)
space = Cp.Space()
Cp.set_gravity(space, gravity)

static = Cp.StaticBody()
groundshape = Cp.add_shape(space, Cp.SegmentShape(static, Cp.Vect(0, -600), Cp.Vect(800, -600), 0))
Cp.set_friction(groundshape, 1)
leftwallshape = Cp.add_shape(space, Cp.SegmentShape(static, Cp.Vect(0, 0), Cp.Vect(0, -600), 0))
Cp.set_friction(leftwallshape, 1)
rightwallshape = Cp.add_shape(space, Cp.SegmentShape(static, Cp.Vect(800, 0), Cp.Vect(800, -600), 0))
Cp.set_friction(rightwallshape, 1)

smiley_texture = Sf.Texture("smiley.png")
Sf.set_smooth(smiley_texture, true)

function make_ball(position, r)
    body = Cp.Body(1, 1)
    shape = Cp.CircleShape(body, r, Cp.Vect(0, 0))
    Cp.set_position(body, position)
    Cp.set_friction(shape, 0.8)

    Cp.add_body(space, body)
    Cp.add_shape(space, shape)
    push!(bodies, body)

    circle = Sf.CircleShape()
    Sf.set_radius(circle, r)
    Sf.set_texture(circle, smiley_texture)
    Sf.set_origin(circle, Sf.Vector2f(r, r))
    push!(circles, circle)
end

timestep = 1./60
while Sf.isopen(window)
    while Sf.pollevent(window, event)
        if Sf.get_type(event) == Sf.EventType.CLOSED
            Sf.close(window)
        end
    end

    if Sf.is_mouse_pressed(0)
        mousepos = Sf.get_mousepos(window)
        make_ball(Cp.Vect(mousepos.x, -mousepos.y), rand(10:30))
    end

    Cp.step(space, timestep)

    Sf.clear(window, Sf.white)
    for i = 1:length(circles)
        circle = circles[i]
        Sf.set_rotation(circle, -rad2deg(Cp.get_angle(bodies[i])))
        Sf.set_position(circle, sf_vec(Cp.get_position(bodies[i])))
        Sf.draw(window, circle)
    end
    Sf.display(window)
end


text:
"""
    ReducePartitionBy(f, rf, [init])

Reduce partitions determined by `isequal` on the output value of `f` with an
associative reducing function `rf`. Partitions are reduced on-the-fly and no
intermediate arrays are allocated.

# Examples

Consider the input `1:6` "keyed" by a function `x -> x ÷ 3`:

```jldoctest
julia> map(x -> x ÷ 3, 1:6)
6-element Vector{Int64}:
 0
 0
 1
 1
 1
 2
```

i.e., there are three partitions with the key values `0`, `1`, and `2`. We
can use `ReducePartitionBy` to compute, e.g., the length and the sum of each
partition by:

```jldoctest
julia> using Transducers

julia> 1:6 |> ReducePartitionBy(x -> x ÷ 3, Map(_ -> 1)'(+)) |> collect
3-element Vector{Int64}:
 2
 3
 1

julia> 1:6 |> ReducePartitionBy(x -> x ÷ 3, +) |> collect
3-element Vector{Int64}:
  3
 12
  6
```
"""
struct ReducePartitionBy{F,RF,Init} <: Transducer
    f::F
    rf::RF
    init::Init
end
ReducePartitionBy(f, rf) = ReducePartitionBy(f, rf, Init)

Adapt.adapt_structure(to, xf::ReducePartitionBy) where {inbounds} = ReducePartitionBy(
    Adapt.adapt(to, xf.f),
    Adapt.adapt(to, xf.rf),
    Adapt.adapt(to, xf.init),
)

struct PartitionChunk{K,V}
    kr::K
    right::V
end

struct PartitionVacant{KL, KR, Left, Right}
    kl::KL
    kr::KR
    left::Left
    right::Right
end

is_prelude(::PartitionChunk) = true

start(rf::R_{ReducePartitionBy}, acc) = wrap(rf, Unseen(), start(inner(rf), acc))

next(rf::R_{ReducePartitionBy}, acc, input) =
    wrapping(rf, acc) do state, iacc
        k = xform(rf).f(input)
        if state isa Unseen
            init = start(xform(rf).rf, xform(rf).init)
            (PartitionChunk(k, next(xform(rf).rf, init, input)), iacc)
        else
            if isequal(k, state.kr)
                acc0 = next(xform(rf).rf, state.right, input)
                ((@set state.right = acc0), iacc)
            else
                if state isa PartitionVacant
                    y = complete(xform(rf).rf, state.right)
                    iacc′ = next(inner(rf), iacc, y)
                    left = state.left
                    kl = state.kl
                else
                    iacc′ = iacc
                    left = state.right
                    kl = state.kr
                end
                init = start(xform(rf).rf, xform(rf).init)
                right = next(xform(rf).rf, init, input)
                (PartitionVacant(kl, k, left, right), iacc′)
            end
        end
    end

function complete(rf::R_{ReducePartitionBy}, acc)
    state, iacc1 = unwrap(rf, acc)
    yr = complete(xform(rf).rf, state.right)
    iacc2 = next(inner(rf), iacc1, yr)
    if state isa PartitionVacant
        yl = complete(xform(rf).rf, state.left)
        iacc0 = next(inner(rf), start(inner(rf), DefaultInit), yl)
        iacc3 = combine(inner(rf), iacc0, iacc2)
    else
        iacc3 = iacc2
    end
    return complete(inner(rf), iacc3)
end

function combine(rf::R_{ReducePartitionBy}, a, b)
    a1, a2 = unwrap(rf, a)
    b1, b2 = unwrap(rf, b)
    a1 isa Unseen && return wrap(rf, b1, b2)
    b1 isa Unseen && return wrap(rf, a1, a2)
    # @show a1 b1
    if a1 isa PartitionChunk
        if b1 isa PartitionChunk
            if isequal(a1.kr, b1.kr)
                c1 = PartitionChunk(a1.kr, combine(xform(rf).rf, a1.right, b1.right))
            else
                c1 = PartitionVacant(a1.kr, b1.kr, a1.right, b1.right)
            end
        else
            if isequal(a1.kr, b1.kl)
                c1 = PartitionVacant(a1.kr, b1.kr, combine(xform(rf).rf, a1.right, b1.left), b1.right)
            else
                y = complete(xform(rf).rf, b1.left)
                # @show y
                a2 = next(inner(rf), a2, y)
                c1 = PartitionVacant(a1.kr, b1.kr, a1.right, b1.right)
            end
        end
    else
        if b1 isa PartitionChunk
            if isequal(a1.kr, b1.kr)
                c1 = PartitionVacant(a1.kl, a1.kr, a1.left, combine(xform(rf).rf, a1.right, b1.right))
            else
                y = complete(xform(rf).rf, a1.right)
                # @show y
                a2 = next(inner(rf), a2, y)
                c1 = PartitionVacant(a1.kl, b1.kr, a1.left, b1.right)
            end
        else
            if isequal(a1.kr, b1.kl)
                acc = combine(xform(rf).rf, a1.right, b1.left)
                y = complete(xform(rf).rf, acc)
                # @show y
                a2 = next(inner(rf), a2, y)
                c1 = PartitionVacant(a1.kl, b1.kr, a1.left, b1.right)
            else
                ya = complete(xform(rf).rf, a1.right)
                yb = complete(xform(rf).rf, b1.left)
                # @show ya yb
                a2 = next(inner(rf), a2, ya)
                a2 = next(inner(rf), a2, yb)
                c1 = PartitionVacant(a1.kl, b1.kr, a1.left, b1.right)
            end
        end
    end
    # @show c1
    c2 = combine(inner(rf), a2, b2)
    return wrap(rf, c1, c2)
end

left(l, r) = l
left(l) = l
InitialValues.@def_monoid left

function array_partitionby(f::F, xs) where {F}
    @inline getkey(i) = f(@inbounds xs[i])
    @inline getview((i, j),) = @inbounds view(xs, i:j)
    return eachindex(xs) |>
        ReducePartitionBy(getkey, TeeRF((left, right)), DefaultInit) |>
        Map(getview)
end

(f::PartitionBy)(xs::AbstractArray) = array_partitionby(f.f, xs)

if VERSION < v"1.3"
    Base.:|>(xs::AbstractArray, f::PartitionBy) = array_partitionby(f.f, xs)
end


text:
fnction mySqrt(a)
    return sqrt(a)
end

function lots_of_math(n)
    var = 0
    for i = 1:n
        var += sqrt(2*pi*i)*(((i-1)+(i+2))/((n+5)*(5*i))
    end
    return var
end

function more_math()
    i = 10
    var1
    = sqrt(2*pi*i)*(((i-1)+(i+2))/((i+5)*(5*i)))
    println(var1)

    i = 20
    var2 =
    sqrt(2*pi*i)*(((i-1)+(i+2))/
        ((i+5)*(5*i)))
    println(var2)
end

function array_init(n)
    array = Array{Float64,2}(n,n)
    for i = 1:size(array,1):
        for j = 1:size(array,2)
            array[i,j] = i*j
        end
    end
end

# Compute the first n terms of the Taylor expansion of e^x centered around x=0
function taylorExpansion(x, n)
    sum = 0
    for i = 1:n
        sum += x^n/factorial(n)
    end
    return sum
end

taylorExpansion(0.5, 25)

# Confuse the audience with a compicated math expression
function lots_of_math2(n)
    var = 0
    for i = 1:n
        var += sqrt(2*pi*i)(((i-1)+(i+2))/((n+5)*(5*i)))
    end
    return var
end

lots_of_math2(5)

# Print out a comparison between x and y
function compare(x,y)
    if x < y
        relation = "less than"
    elseif x > y
        relation = "greater than"
    end
    println(x, " is ", relation, " ", y)
end

compare(5,5)

# Initialize a 2d array to hold first m x n integer products
function array_init2(m,n)
    array = Matrix{Float64}(undef,m,n)
    for i = 1:n
        for j = 1:m
            array[i,j] = i*j
        end
    end
    return array
end

array_init2(5,4)

# Divide an input number by 0
function divide_by_zero(input)
    return input/0
end

divide_by_zero("string")

# Divide an input number by 0
function divide_by_zero(input::Number) # Only allows an input of type Number
    return input/0
end

divide_by_zero("string")

# Returns the sum of the first n integers
function array_sum(n)
    # Initialization
    array = Vector{Float64}(undef,n)
    for i = 1:n-1
        array[i] = i
    end

    # Summation
    sum = 0
    for i = 1:n
        sum += array[i]
    end
    return sum
end

# Testing cell for array_sum()

# Returns the sum of the first n integers
function array_sum(n)
    # Initialization
    array = Vector{Float64}(undef,n)
    for i = 1:n-1
        array[i] = i
    end
    println(array) # TODO: Debugging

    # Summation
    sum = 0
    println("Sum before = ", sum) # TODO: Debugging
    for i = 1:n
        sum += array[i]
    end
    println("Sum after = ", sum) # TODO: Debugging
    return sum
end

array_sum(1);
array_sum(2);
array_sum(5);


text:
# This example is basically the same as the instanced_depth peel but now we are also adding an osscilation system to it.
#%%
using Glimpse
const Gl = Glimpse
dio = Gl.Diorama(; background = RGBAf0(0.0, 0.0, 0.0, 1.0));
function spherepoint(r, theta, phi)
    return r * Point3f0(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta))
end
nspheres = 20000
radii = Iterators.cycle(range(1.0f0, 200.0f0; length = div(nspheres, 40)))
angs = Iterators.cycle(range(0, 2pi; length = div(nspheres, 40)))
angs2 = Iterators.cycle(range(-pi, pi; length = div(nspheres, 60)))

ps       = [spherepoint(r, a1, a2) for (r, a1, a2, i) in zip(radii, angs, angs2, 1:nspheres)]
progtag  = Gl.ProgramTag{Gl.InstancedPeelingProgram}()
cs       = [RGBAf0(0.1, 0.5, 0.9, 0.4), RGBAf0(0.9, 0.1, 0.5, 0.4)]
sph_geom = Gl.PolygonGeometry(Gl.Sphere(Point3f0(0.0), 1.0f0))
spring   = Gl.Spring()
for i in 1:nspheres
    color = cs[mod1(i, 2)]
    Entity(dio, Gl.Spatial(; position = ps[i]), Gl.Material(), Gl.Shape(),
           Gl.UniformColor(color), Gl.Dynamic(), Gl.Selectable(), progtag, sph_geom, spring)
end

Gl.renderloop(dio)
#%%


text:
using McCormick, EAGO, JuMP, NaNMath, CSV, DataFrames

include("worse_case_subsea_surrogate.jl")

#validity_surrogate_path = "validity_surrogate"
#performance_surrogate_path = "performance_surrogate"
#ANN_valid = surrogate_build_model(validity_surrogate_path)
#ANN_perform = surrogate_build_model(performance_surrogate_path)

API = 35
SG_G = 0.6
SG_W = 1.0
ρᵒ_W = 1000.0

P_well = 5.52*10^6 # Pa

k_GLS = 8.33*10^(-3)
L_GLS = 4.0        # Length (m) of LLS
R_GLS = 0.6        # Radius (m) of LLS
P_GLS = 4.0*10^6   # Pa

k_LLS = 1.67*10^(-4)
L_LLS = 5.0        # Length (m) of LLS
R_LLS = 0.8        # Radius (m) of LLS
P_LLS = 4*10^6     # Operating Pressure (Pa) of LLS
H_LLS = 0.6        # Liquid Level in the LLS (m)

C_v1 = 1.67*10^(-2)
C_v2 = 0.1675

SG_O = 141.5/(131/5 + API)

# Fixed relations from Control Value (V-1)
P_1 = P_well
P_2 = P_GLS
ξO1 = 0.4

# Fixed relations from Gas-liquid Separator
P_3 = P_2
ξ3 = [1.0; 0.0; 0.0]

# Fixed relations from Control Value (V-2)
P_5 = P_LLS

# Fixed relations from Liquid-liquid separator
V_LLS = L_LLS*((H_LLS - R_LLS)*sqrt(2.0*R_LLS*H_LLS - H_LLS^2) + (R_LLS^2)*acos(1.0 - H_LLS/R_LLS))
ξ6 = [0.0; 1.0; 0.0]
ξ8 = [1.0; 0.0; 0.0]
P_8 = P_LLS

ξG7_spec = 0.05

guard_tol = 1E-8

f(x) = -x[1]
function gSIP_perform(x, p)
    u1 = x[2]
    u2 = x[3]
    ξG1 = p[1]
    ξW1 = 1 - ξO1 - ξG1
    SG_mix_inv_1 = ξG1/SG_G + ξW1/SG_W + ξO1/SG_O
    sqrt_arg1 = positive((P_1 - P_2)*SG_mix_inv_1) + guard_tol
    m2 = u1*C_v1*sqrt(sqrt_arg1)
    ξG2 = ξG1
    ξW2 = ξW1
    HGLS, ξG4, P4, ρ4 = ANN_perform((m2, u2, ξG2, ξW2,))
    sqrt_arg2 = positive(ρᵒ_W*(P4 - P_LLS)/positive(ρ4)) + guard_tol
    m4 = u2*C_v2*sqrt(sqrt_arg2)
    ξG7 = ξG4*exp(-k_LLS*m4*V_LLS/(positive(ρ4)+ guard_tol))
    return ξG7 - ξG7_spec
end

function gSIP_valid(x, p)
    u1 = x[2]
    u2 = x[3]
    ξG1 = p[1]
    ξW1 = 1 - ξO1 - ξG1
    SG_mix_inv_1 = ξG1/SG_G + ξW1/SG_W + ξO1/SG_O
    sqrt_arg1 = positive((P_1 - P_2)*SG_mix_inv_1) + guard_tol
    m2 = u1*C_v1*sqrt(sqrt_arg1)
    ξG2 = ξG1
    ξW2 = ξW1
    valid_flout = ANN_valid((m2, u2, ξG2, ξW2,))
    0.01*(1.0 - 1.0*valid_flout[1])
end

function gSIP(x,p)
    g1 = gSIP_valid(x, p)
    g2 = gSIP_perform(x, p)
    return x[1] - max(g1, g2)
end
function solve_worst_case_sip()

    x_l = [-0.05, 0.35, 0.35]
    x_u = [0.05, 0.8, 0.8]
    p_l = [0.35]
    p_u = [0.5]

    opt = EAGO.Optimizer

    sip_tol = 1E-4
    abs_tol = sip_tol/10.0
    rel_tol = sip_tol/10.0

    sip_result = explicit_sip_solve(x_l, x_u, p_l, p_u, f, [gSIP],
                                    sip_absolute_tolerance = 1E-3,
                                    sip_relative_tolerance = 1E-3,
                                    sip_optimizer = opt,
                                    absolute_tolerance = abs_tol,
                                    relative_tolerance = rel_tol,
                                    #sip_verbosity = 4,
                                    sip_sense = :min,
                                    verbosity = 1,
                                    output_iterations = 10000,
                                    upper_bounding_depth = 50)
    return sip_result
end

start_time = time_ns()
worst_case_sip_result = solve_worst_case_sip()
end_time = time_ns()
elapsed_time = end_time - start_time
@show elapsed_time/(1000^3)

max_stars_repo_name: PSORLab/RobustHybridModels

max_issues_repo_name: PSORLab/RobustHybridModels

max_forks_repo_name: PSORLab/RobustHybridModels


text:
__precompile__()
module COSMO

using SparseArrays, LinearAlgebra, SuiteSparse

export  assemble!, warmStart!, empty_model!

const DefaultFloat = Float64
const DefaultInt   = Int64

include("./algebra.jl")
include("./projections.jl")
include("./settings.jl")            # TODO: unmodified - revisit
include("./types.jl")               # some types still need tidying
include("./constraint.jl")          # TODO: unmodified - revisit
include("./parameters.jl")          # TODO: unmodified - revisit
include("./residuals.jl")           # TODO: unmodified - revisit
include("./scaling.jl")             # TODO: set scaling / E scaling is broken
include("./kkt.jl")                 # TODO: unmodified - revisit.  Add lin solver type
include("./infeasibility.jl")       # TODO: stylistic fixes needed
include("./printing.jl")            # TODO: unmodified - revisit
include("./setup.jl")               # TODO: unmodified - revisit (short - consolidate?)
include("./solver.jl")              # TODO: unmodified - revisit
include("./interface.jl")           # TODO: unmodified - revisit
include("./MOIWrapper.jl")

end #end module


text:
function comment_using_revise(filename)
    lines = readlines(filename)
    idx = findfirst(lines .== "using Revise")
    if isnothing(idx)
        @warn "Could not find \"using Revise\" in file $filename, doing nothing."
    else
        lines[idx] = "# using Revise"
        open(filename, "w") do f
            for l in lines
                println(f,l)
            end
        end
    end
end

testfiles = filter(x->startswith(x,"test_"),readdir("test"))
testfiles = ["test/"*s for s in testfiles]
for f in testfiles
    comment_using_revise(f)
end