text:
# Ordinary Maximum Likelihood Estimation

############################################################################
### Types
############################################################################

struct SemML{INV,M,M2,B,FT,GT,HT} <: SemLossFunction
    inverses::INV #preallocated inverses of imp_cov
    mult::M
    meandiff::M2
    approx_H::B

    objective::FT
    gradient::GT
    hessian::HT
end

############################################################################
### Constructors
############################################################################

function SemML(;observed, n_par, approx_H = false, parameter_type = Float64, kwargs...)
    isnothing(obs_mean(observed)) ?
        meandiff = nothing :
        meandiff = copy(obs_mean(observed))
    return SemML(
        copy(obs_cov(observed)),
        copy(obs_cov(observed)),
        meandiff,
        approx_H,

        zeros(parameter_type, 1),
        zeros(parameter_type, n_par),
        zeros(parameter_type, n_par, n_par)
        )
end

############################################################################
### functors
############################################################################

# for symbolic imply type
function (semml::SemML)(
    par,
    F,
    G,
    H,
    model::Sem{O, I, L, D}) where {O, I <: SemImplySymbolic, L, D}
    semml.inverses .= Σ(imply(model))
    a = cholesky!(Symmetric(semml.inverses); check = false)

    if !isposdef(a)
        if G semml.gradient .= 1.0 end
        if H semml.hessian .= 1.0 end
        if F semml.objective[1] = Inf end
    else
        ld = logdet(a)
        semml.inverses .= LinearAlgebra.inv!(a)

        # without means
        if isnothing(μ(imply(model)))

            if G && H
                J = (vec(semml.inverses)-vec(semml.inverses*obs_cov(observed(model))*semml.inverses))'
                gradient = J*∇Σ(imply(model))
                semml.gradient .= gradient'
                if semml.approx_H
                    hessian = 2*∇Σ(imply(model))'*kron(semml.inverses, semml.inverses)*∇Σ(imply(model))
                end
                if !semml.approx_H
                    M = semml.inverses*obs_cov(observed(model))*semml.inverses
                    H_outer =
                        2*kron(M, semml.inverses) -
                        kron(semml.inverses, semml.inverses)
                    hessian = ∇Σ(imply(model))'*H_outer*∇Σ(imply(model))
                    ∇²Σ_function(imply(model))(∇²Σ(imply(model)), J, par)
                    hessian += ∇²Σ(imply(model))
                end
                semml.hessian .= hessian
            end

            if G && !H
                gradient = (vec(semml.inverses)-vec(semml.inverses*obs_cov(observed(model))*semml.inverses))'*∇Σ(imply(model))
                semml.gradient .= gradient'
            end

            if !G && H
                J = (vec(semml.inverses)-vec(semml.inverses*obs_cov(observed(model))*semml.inverses))'
                if semml.approx_H
                    hessian = 2*∇Σ(imply(model))'*kron(semml.inverses, semml.inverses)*∇Σ(imply(model))
                end
                if !semml.approx_H
                    M = semml.inverses*obs_cov(observed(model))*semml.inverses
                    H_outer =
                        2*kron(M, semml.inverses) -
                        kron(semml.inverses, semml.inverses)
                    hessian = ∇Σ(imply(model))'*H_outer*∇Σ(imply(model))
                    ∇²Σ_function(imply(model))(∇²Σ(imply(model)), J, par)
                    hessian += ∇²Σ(imply(model))
                end
                semml.hessian .= hessian
            end

            if F
                mul!(semml.mult, semml.inverses, obs_cov(observed(model)))
                semml.objective[1] = ld + tr(semml.mult)
            end
        else
        # with means
        μ_diff = obs_mean(observed(model)) - μ(imply(model))
        diff⨉inv = μ_diff'*semml.inverses
            if H throw(DomainError(H, "hessian of ML + meanstructure is not implemented yet")) end
            if G
                gradient =
                    vec(
                        semml.inverses*(
                            LinearAlgebra.I -
                            obs_cov(observed(model))*semml.inverses -
                            μ_diff*diff⨉inv))'*∇Σ(imply(model)) -
                    2*diff⨉inv*∇μ(imply(model))
                semml.gradient .= gradient'
            end
            if F
                mul!(semml.mult, semml.inverses, obs_cov(observed(model)))
                semml.objective[1] = ld + tr(semml.mult) + diff⨉inv*μ_diff
            end
        end
    end
end

# for non-symbolic imply type
function (semml::SemML)(par, F, G, H, model::Sem{O, I, L, D}) where {O, I <: RAM, L, D}

    if H
        throw(DomainError(H, "hessian for ML estimation with non-symbolic imply type is not implemented"))
    end

    semml.inverses .= Σ(imply(model))
    a = cholesky!(Symmetric(semml.inverses); check = false)

    if !isposdef(a)
        if G semml.gradient .= 1.0 end
        if H semml.hessian .= 1.0 end
        if F semml.objective[1] = Inf end
    else
        ld = logdet(a)
        semml.inverses .= LinearAlgebra.inv!(a)

        # without means
        if isnothing(μ(imply(model)))

            if G
                gradient = SemML_gradient(
                    S(imply(model)),
                    F⨉I_A⁻¹(imply(model)),
                    semml.inverses,
                    I_A(imply(model)),
                    ∇A(imply(model)),
                    ∇S(imply(model)),
                    obs_cov(observed(model)))
                semml.gradient .= gradient'
            end

            if F
                mul!(semml.mult, semml.inverses, obs_cov(observed(model)))
                semml.objective[1] = ld + tr(semml.mult)
            end

        else # with means

            μ_diff = obs_mean(observed(model)) - μ(imply(model))
            diff⨉inv = μ_diff'*semml.inverses

            if G
                gradient = SemML_gradient(
                        S(imply(model)),
                        F⨉I_A⁻¹(imply(model)),
                        semml.inverses,
                        I_A(imply(model)),
                        ∇A(imply(model)),
                        ∇S(imply(model)),
                        obs_cov(observed(model))) +
                    SemML_gradient_meanstructure(
                        diff⨉inv,
                        F⨉I_A⁻¹(imply(model)),
                        I_A(imply(model)),
                        S(imply(model)),
                        M(imply(model)),
                        ∇M(imply(model)),
                        ∇A(imply(model)),
                        ∇S(imply(model)))
                semml.gradient .= gradient'
            end

            if F
                mul!(semml.mult, semml.inverses, obs_cov(observed(model)))
                semml.objective[1] = ld + tr(semml.mult) + diff⨉inv*μ_diff
            end

        end
    end
end

############################################################################
### recommended methods
############################################################################

objective(lossfun::SemML) = lossfun.objective
gradient(lossfun::SemML) = lossfun.gradient
hessian(lossfun::SemML) = lossfun.hessian

update_observed(lossfun::SemML, observed::SemObsMissing; kwargs...) =
    throw(ArgumentError("ML estimation does not work with missing data - use FIML instead"))

function update_observed(lossfun::SemML, observed::SemObs; kwargs...)
    if (size(lossfun.inverses) == size(obs_cov(observed))) & (isnothing(lossfun.meandiff) == isnothing(obs_mean(observed)))
        return lossfun
    else
        return SemML(;observed = observed, kwargs...)
    end
end

############################################################################
### additional functions
############################################################################

function SemML_gradient_common(F⨉I_A⁻¹, obs_cov, Σ⁻¹)
    M = transpose(F⨉I_A⁻¹)*transpose(LinearAlgebra.I-obs_cov*Σ⁻¹)*Σ⁻¹*F⨉I_A⁻¹
    return M
end

function SemML_gradient(S, F⨉I_A⁻¹, Σ⁻¹, I_A⁻¹, ∇A, ∇S, obs_cov)
    M = SemML_gradient_common(F⨉I_A⁻¹, obs_cov, Σ⁻¹)
    G = 2vec(M*S*I_A⁻¹')'∇A + vec(M)'∇S
    return G
end

function SemML_gradient_meanstructure(diff⨉inv, F⨉I_A⁻¹, I_A⁻¹, S, M, ∇M, ∇A, ∇S)
    k = diff⨉inv*F⨉I_A⁻¹
    G = -2k*∇M - 2vec(k'M'I_A⁻¹')'∇A - 2vec(k'k*S*I_A⁻¹')'∇A - vec(k'k)'∇S
    return G
end

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemML)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
using Plots, ApproxFun
# The Blasius boundary layer problem (https://en.wikipedia.org/wiki/Blasius_boundary_layer)
# describes the boundary layer on a semi-infinite plane with a head-on steady fluid flow.
# When properly non-dimensionalized, it presents itself as a relatively simple nonlinear BVP.
# The Falkner-Skan problem generalizes the Blasius boundary layer by considering a wedge at an
# angle to the steady flow.
# In both cases, no-slip conditions are in place at the wall and the boundary layer
# grows linearly (u'(∞)→1) in the far-field.
# There are more variants which include different nonlinearities which arise in magneto-
# hydrodynamics, for example.

# We begin with Blasius:
#   solves 2u''' + uu'' = 0 ,  u(0.) = u'(0.) = 0, u'[∞] = 1
x=Fun(identity,0..4π)
d=domain(x)
B=[ldirichlet(d),lneumann(d),rneumann(d)]
D=Derivative(d)
κ = 0.33205733621519630     # This is diff(u,2)(0.), due to Boyd.
u0 = (1//2) * κ * x^2      # First-order approximation
u = 0.5x^2                  # Other initial solutions may fail to converge

f = (u)->(2.0*D^3*u + u*D^2*u)
df = (u)->(2.0*D^3 + u*D^2 + D^2*u)

plot([u,u0])      # Initial estimate and first-order approximation

for k=1:10
    global u  -= [B; df(u)]\[u(0.);u'(0.);u'(rightendpoint(d))-1.;f(u)]
end
norm(f(u))             # should be zero
abs(u''(0.)-κ)         # should also be zero
plot([u,u0])           # Solution and First-order approximation
# Now for Falkner-Skan:
#     solves 2u''' + uu'' - β(1 - (u')^2) = 0 ,  u(0.) = u'(0.) = 0, u'[∞] = 1
m = 0.11                # m ∈ [-0.0905, 2]
β = 2m/(1+m)
F = (u)->(2.0*D^3*u + u*D^2*u - β*(1.0 - (D*u)*(D*u)))
dF = (u)->(2.0*D^3 + u*D^2 + D^2*u + 2β*D*u*D)
v=u
norm(F(v))               # should be non-zero, as we've perturbed the Blasius equation

for k=1:10
  global v  -= [B; dF(v)]\[v(0.);v'(0.);v'(rightendpoint(d))-1.;F(v)]
end
norm(F(v))               # should be zero

plot([u,v])    # Blasius and Falkner-Skan solutions

ext: jl

max_stars_repo_path: ODEs/Blasius and Falkner-Skan.jl

max_issues_repo_path: ODEs/Blasius and Falkner-Skan.jl

max_forks_repo_path: ODEs/Blasius and Falkner-Skan.jl


text:
#=
# Plotting Graphs with `GraphMakie.jl`
## The `graphplot` Command
Plotting your first `AbstractGraph` from [`LightGraphs.jl`](https://juliagraphs.org/LightGraphs.jl/latest/)
is as simple as
=#
using CairoMakie
CairoMakie.activate!(type="png") # hide
set_theme!(resolution=(800, 400)) #hide
CairoMakie.inline!(true) # hide
using GraphMakie
using LightGraphs
import Random; Random.seed!(2) # hide

g = wheel_graph(10)
f, ax, p = graphplot(g)
hidedecorations!(ax); hidespines!(ax)
ax.aspect = DataAspect()
f # hide

#=
The `graphplot` command is a recipe which wraps several steps
- layout the graph in space using a layout function,
- create a `scatter` plot for the nodes and
- create a `linesegments` plot for the edges.

The default layout is `Spring()` from
[`NetworkLayout.jl`](https://github.com/JuliaGraphs/NetworkLayout.jl). The
layout attribute can be any function which takes an `AbstractGraph` and returns
a list of `Point{dim,Ptype}` (see [`GeometryBasics.jl`](https://github.com/JuliaGeometry/GeometryBasics.jl)
objects where `dim` determines the dimensionality of the plot.

Besides that there are some common attributes which are forwarded to the
underlying plot commands. See [`graphplot`](@ref).
=#
using GraphMakie.NetworkLayout

g = SimpleGraph(5)
add_edge!(g, 1, 2); add_edge!(g, 2, 4);
add_edge!(g, 4, 3); add_edge!(g, 3, 2);
add_edge!(g, 2, 5); add_edge!(g, 5, 4);
add_edge!(g, 4, 1); add_edge!(g, 1, 5);

## define some edge colors
edgecolors = [:black for i in 1:ne(g)]
edgecolors[4] = edgecolors[7] = :red

f, ax, p = graphplot(g, layout=Shell(),
                     node_color=[:black, :red, :red, :red, :black],
                     edge_color=edgecolors)

hidedecorations!(ax); hidespines!(ax)
ax.aspect = DataAspect()
f #hide

#=
We can interactively change the attributes as usual with Makie.
=#

fixed_layout(_) = [(0,0), (0,1), (0.5, 1.5), (1,1), (1,0)]
## set new layout
p.layout = fixed_layout; autolimits!(ax)
## change edge width & color
p.edge_width = 5.0
p.edge_color[][3] = :green;
p.edge_color = p.edge_color[] # trigger observable
f #hide

#=
## Adding Node Labels
=#
Random.seed!(2)
g = wheel_graph(10)

colors = [:black for i in 1:nv(g)]
colors[1] = :red

f, ax, p = graphplot(g,
                     nlabels=repr.(1:nv(g)),
                     nlabels_color=colors,
                     nlabels_align=(:center,:center))
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

# This is not very nice, lets change the offsets based on the `node_positions`

offsets = 0.15 * (p[:node_positions][] .- p[:node_positions][][1])
offsets[1] = Point2f0(0, 0.3)
p.nlabels_offset[] = offsets
autolimits!(ax)
f # hide

# ## Adding Edge Labels
Random.seed!(42)
g = barabasi_albert(6, 2)

labels =  repr.(1:ne(g))

f, ax, p = graphplot(g, elabels=labels,
                     elabels_color=[:black for i in 1:ne(g)],
                     edge_color=[:black for i in 1:ne(g)])
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
The position of the edge labels is determined by several plot arguments.
All possible arguments are described in the docs of the [`graphplot`](@ref) function.

Basicially, each label is placed in the middle of the edge and rotated to match the edge rotation.
The rotaion for each label can be overwritten with the `elabels_rotation` argument.
=#
p.elabels_rotation[] = Vector{Union{Nothing, Float64}}(nothing, ne(g))
p.elabels_rotation[][5] = 0.0 # set absolute rotation angle for label 5
p.elabels_rotation[] = p.elabels_rotation[]
nothing #hide

#=
One can shift the label along the edge with the `elabels_shift` argument and determine the distance
in pixels using the `elabels_distance` argument.
=#
p.elabels_opposite[] = [1,2,8,6]

p.elabels_offset[] = [Point2f0(0.0, 0.0) for i in 1:ne(g)]
p.elabels_offset[][5] = Point2f0(-0.4,0)
p.elabels_offset[] = p.elabels_offset[]

p.elabels_shift[] = [0.5 for i in 1:ne(g)]
p.elabels_shift[][1] = 0.6
p.elabels_shift[][7] = 0.4
p.elabels_shift[] = p.elabels_shift[]

p.elabels_distance[] = zeros(ne(g))
p.elabels_distance[][8] = 15
p.elabels_distance[] = p.elabels_distance[]

f # hide

#=
## Indicate Edge Direction

It is possible to put arrows on the edges using the `arrow_show` parameter. This parameter
is `true` for `SimpleDiGraph` by default. The position and size of each arrowhead can be
change using the `arrow_shift` and `arrow_size` parameters.
=#
g = wheel_digraph(10)
arrow_size = [10+i for i in 1:ne(g)]
arrow_shift = range(0.1, 0.8, length=ne(g))
f, ax, p = graphplot(g; arrow_size, arrow_shift)
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
## Plot Graphs in 3D
If the layout returns points in 3 dimensions, the plot will be in 3D. However this is a bit
experimental. Feel free to file an issue if there are any problems.
=#
set_theme!(resolution=(800, 800)) #hide
g = smallgraph(:cubical)
elabels_shift = [0.5 for i in 1:ne(g)]
elabels_shift[[2,7,8,9]] .= 0.3
elabels_shift[10] = 0.25
graphplot(g; layout=Spring(dim=3, seed=5),
          elabels="Edge ".*repr.(1:ne(g)),
          elabels_textsize=12,
          elabels_opposite=[3,5,7,8,12],
          elabels_shift,
          elabels_distance=3,
          arrow_show=true,
          arrow_shift=0.9,
          arrow_size=15)

#=
Using [`JSServe.jl`](https://github.com/SimonDanisch/JSServe.jl) and [`WGLMakie.jl`](https://github.com/JuliaPlots/WGLMakie.jl)
we can also add some interactivity:
=#
using JSServe
Page(exportable=true, offline=true)
#
using WGLMakie
WGLMakie.activate!()
set_theme!(resolution=(800, 600))
g = smallgraph(:dodecahedral)
graphplot(g, layout=Spring(dim=3), node_size=100)


text:
using Test

# # State debug statistics
#
# Set up a basic environment
using MPI
using StaticArrays
using Random
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.MPIStateArrays
using ClimateMachine.StateCheck
using ClimateMachine.GenericCallbacks

@testset "$(@__FILE__)" begin

    ClimateMachine.init()
    FT = Float64

    # Define some dummy vector and tensor abstract variables with associated types
    # and dimensions
    F1 = @vars begin
        ν∇u::SMatrix{3, 2, FT, 6}
        κ∇θ::SVector{3, FT}
    end
    F2 = @vars begin
        u::SVector{2, FT}
        θ::SVector{1, FT}
    end

    # Create ```MPIStateArray``` variables with arrays to hold elements of the
    # vectors and tensors
    Q1 = MPIStateArray{Float32, F1}(
        MPI.COMM_WORLD,
        ClimateMachine.array_type(),
        4,
        9,
        8,
    )
    Q2 = MPIStateArray{Float64, F2}(
        MPI.COMM_WORLD,
        ClimateMachine.array_type(),
        4,
        3,
        8,
    )

    # ### Create a call-back
    cb = ClimateMachine.StateCheck.sccreate(
        [(Q1, "My gradients"), (Q2, "My fields")],
        1;
        prec = 15,
    )

    # ### Invoke the call-back
    #     Compare on local "realdata", fill via broadcast to keep GPU happy.
    Q1.data = rand(MersenneTwister(0), Float32, size(Q1.data))
    Q2.data = rand(MersenneTwister(0), Float64, size(Q2.data))
    Q1.realdata .= Q1.data
    Q2.realdata .= Q2.data
    GenericCallbacks.init!(cb, nothing, nothing, nothing, nothing)
    GenericCallbacks.call!(cb, nothing, nothing, nothing, nothing)

    # ### Check with reference values
    include("test_statecheck_refvals.jl")

    # This should return true (MPI rank != 0 always returns true)
    test_stat_01 = ClimateMachine.StateCheck.scdocheck(cb, (varr, parr))

    # This should return false (MPI rank != 0 always returns true)
    varr[1][3] = varr[1][3] * 10.0
    test_stat_02 = ClimateMachine.StateCheck.scdocheck(cb, (varr, parr))
    if MPI.Comm_rank(MPI.COMM_WORLD) != 0
        test_stat_02 = false
    end

    test_stat = test_stat_01 && !test_stat_02

    @test test_stat

end


text:
println("=====================test_pricing_problem=====================")
# P = Problem(
#     [0:1, 0:1],             # Capacity of edges
#     [3, 5],                 # Selling period end of edges
#     5,                      # Number of timesteps (Start at 1)
#     [10, 20],               # Actions (prices)
#     [(1,), (2,), (1, 2)],   # Products (seqeuences of edge indeces)
#     Dict((1,) => 0.5,       # λ: Dictionary of demand intensities for products
#         (2,) => 0.5,
#         (1, 2) => 0.4),
# )

@testset "PricingProblem.jl" begin
using Mdp.PricingProblem

um = price -> prob_user_accept(price, slope_start=5., slope_end=30.)
e1, e2 = Edge(1, 3), Edge(1, 5)
prd1, prd2, prd12 = Product([e1,], 0.5, um),
                    Product([e2,], 0.5, um),
                    Product([e1,e2], 0.4, um)
P = Problem(
    [e1,e2],
    [prd1, prd2, prd12],
    [10,20],
    5
)

@testset "prob_user_accept" begin
    @test prob_user_accept(5.0; slope_start = 5.0, slope_end = 10.0) == 1.0
    @test prob_user_accept(10.0; slope_start = 5.0, slope_end = 10.0) == 0.0
    @test prob_user_accept(7.5; slope_start = 5.0, slope_end = 10.0) == 0.5
    @test prob_user_accept(3.0; slope_start = 5.0, slope_end = 10.0) == 1.0
    @test prob_user_accept(15.0; slope_start = 5.0, slope_end = 10.0) == 0.0
end

@testset "len_product_selling_period" begin
    @test len_product_selling_period(prd1) == 3
    @test len_product_selling_period(prd2) == 5
    @test len_product_selling_period(prd12) == 3
end

@testset "prob_prod_req" begin
    k = 3
    @test prob_prod_req(prd1, k, P) == 0.16666666666666666
    @test prob_prod_req(prd2, k, P) == 0.1
    @test prob_prod_req(prd12, k, P) == 0.13333333333333333
    k = 4
    @test prob_prod_req(prd1, k, P) == 0.
    @test prob_prod_req(prd2, k, P) == 0.1
    @test prob_prod_req(prd12, k, P) == 0.
end

@testset "product_generator" begin
    all_possible_products([e1, e2], 0.5, um)==[Product([e1], 0.5,um),
                                                Product([e2], 0.5,um),
                                                Product([e1,e2], 0.5,um)]
end

end


text:
using BenchmarkTools, FinancialMonteCarlo, ReverseDiff

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 100000;
Nstep = 30;
sigma = 0.2
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

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], Underlying(x[2], d)), ZeroRate(r), mc, EuropeanOption(x[3], K));
#f(x) = pricer(BlackScholesProcess(x[1]),ZeroRate(x[2],r,d),mc,EuropeanOption(T,K),Underlying(S0,d));
#x=Float64[sigma,S0,r,d,T]
x = Float64[sigma, S0, T]
g = x -> ReverseDiff.gradient(f, x);
y0 = g(x);
@btime f(x);
@btime g(x);

#f1_(x) = pricer(BlackScholesProcess(x[1]),ZeroRate(x[2],x[3],x[4]),mc,AmericanOption(x[5],K),Underlying(S0,d));
f1_(x) = pricer(BlackScholesProcess(x[1], Underlying(x[2], d)), ZeroRate(r), mc, AmericanOption(x[3], K));
g_ = x -> ReverseDiff.gradient(f1_, x);
#y0=g(x);
@btime f1_(x);
#@btime g_(x);


text:
#####
# Data-Driven Portfolio Allocation Example
#####
#Solves the w.c. portfolio allocation problem
##  max t
#   s.t. r^T x >= t  for all  r in U
#        1^T x == 1
#        x >= 0
# where U is constructed from data
using JuMP, JuMPeR, DDUS, Distributions

#generates some silly market data for example
function genReturnData(numAssets, numObs)
	#Single factor CAPM model
	z = .1 * randn(numObs) + .2
	betas = linspace(0, 1, numAssets)
	z * betas' + .05*randn(numObs, numAssets)
end

##########
srand(8675309)
mkt_data = genReturnData(5, 500)

##
# Other sets can be constructed by suitably changing this line
oracle = UCSOracle(mkt_data, .1, .2)

#build model
m = RobustModel()
setDefaultOracle!(m, oracle)
@defVar(m, xs[1:5] >= 0)
@defVar(m, t <= 100)
@defUnc(m, us[1:5])
@addConstraint(m, sum(xs) == 1.)
@addConstraint(m, sum{us[i] * xs[i], i=1:5} >= t)
@setObjective(m, Max, t)

#for now, must use cuts to solve
#Notice how set construction takes advantage of the correlations tructure
#So that allocation favors higher indices
println(solve(m, prefer_cuts=true))
println("Obj Value:\t", getObjectiveValue(m))
println("Portfolio:\t", getValue(xs))


text: iscoprime(x::Integer, y::Integer) = gcd(x,y) == 1

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Object-based wrapper for level-2 BLIS routines.
#

macro blis_object_api_level2_mv(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           a, xObj,
                           x, xObj,
                           β, xObj,
                           y, xObj)
    end
end

@blis_object_api_level2_mv gemv
@blis_object_api_level2_mv hemv
@blis_object_api_level2_mv symv

macro blis_object_api_level2_r2(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           x, xObj,
                           y, xObj,
                           a, xObj)
    end
end

@blis_object_api_level2_r2 ger
@blis_object_api_level2_r2 her2
@blis_object_api_level2_r2 syr2

macro blis_object_api_level2_r(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           x, xObj,
                           a, xObj)
    end
end

@blis_object_api_level2_r her
@blis_object_api_level2_r syr

macro blis_object_api_level2_tr(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           α, xObj,
                           a, xObj,
                           x, xObj)
    end
end

@blis_object_api_level2_tr trmv
@blis_object_api_level2_tr trsv


text:
Elements = (
    Hydrogen = Element(1, "H", "Hydrogen", 0.99212, 19.2eV, 8.375E-05g / cm^3),
    Helium = Element(2, "He", "Helium", 0.49968, 41.8eV, 1.663E-04g / cm^3),
    Lithium = Element(3, "Li", "Lithium", 0.43221, 40.0eV, 5.340E-01g / cm^3),
    Beryllium =
        Element(4, "Be", "Beryllium", 0.44384, 63.7eV, 1.848E+00g / cm^3),
    Boron = Element(5, "B", "Boron", 0.46245, 76.0eV, 2.370E+00g / cm^3),
    Carbon = Element(6, "C", "Carbon", 0.49954, 78.0eV, 1.700E+00g / cm^3),
    Nitrogen = Element(7, "N", "Nitrogen", 0.49976, 82.0eV, 1.165E-03g / cm^3),
    Oxygen = Element(8, "O", "Oxygen", 0.50002, 95.0eV, 1.332E-03g / cm^3),
    Fluorine = Element(9, "F", "Fluorine", 0.47372, 115.0eV, 1.580E-03g / cm^3),
    Neon = Element(10, "Ne", "Neon", 0.49555, 137.0eV, 8.385E-04g / cm^3),
    Sodium = Element(11, "Na", "Sodium", 0.47847, 149.0eV, 9.710E-01g / cm^3),
    Magnesium =
        Element(12, "Mg", "Magnesium", 0.49373, 156.0eV, 1.740E+00g / cm^3),
    Aluminum =
        Element(13, "Al", "Aluminum", 0.48181, 166.0eV, 2.699E+00g / cm^3),
    Silicon = Element(14, "Si", "Silicon", 0.49848, 173.0eV, 2.330E+00g / cm^3),
    Phosphorus =
        Element(15, "P", "Phosphorus", 0.48428, 173.0eV, 2.200E+00g / cm^3),
    Sulfur = Element(16, "S", "Sulfur", 0.49897, 180.0eV, 2.000E+00g / cm^3),
    Chlorine =
        Element(17, "Cl", "Chlorine", 0.47951, 174.0eV, 2.995E-03g / cm^3),
    Argon = Element(18, "Ar", "Argon", 0.45059, 188.0eV, 1.662E-03g / cm^3),
    Potassium =
        Element(19, "K", "Potassium", 0.48595, 190.0eV, 8.620E-01g / cm^3),
    Calcium = Element(20, "Ca", "Calcium", 0.49903, 191.0eV, 1.550E+00g / cm^3),
    Scandium =
        Element(21, "Sc", "Scandium", 0.46712, 216.0eV, 2.989E+00g / cm^3),
    Titanium =
        Element(22, "Ti", "Titanium", 0.45948, 233.0eV, 4.540E+00g / cm^3),
    Vanadium =
        Element(23, "V", "Vanadium", 0.45150, 245.0eV, 6.110E+00g / cm^3),
    Chromium =
        Element(24, "Cr", "Chromium", 0.46157, 257.0eV, 7.180E+00g / cm^3),
    Manganese =
        Element(25, "Mn", "Manganese", 0.45506, 272.0eV, 7.440E+00g / cm^3),
    Iron = Element(26, "Fe", "Iron", 0.46556, 286.0eV, 7.874E+00g / cm^3),
    Cobalt = Element(27, "Co", "Cobalt", 0.45815, 297.0eV, 8.900E+00g / cm^3),
    Nickel = Element(28, "Ni", "Nickel", 0.47708, 311.0eV, 8.902E+00g / cm^3),
    Copper = Element(29, "Cu", "Copper", 0.45636, 322.0eV, 8.960E+00g / cm^3),
    Zinc = Element(30, "Zn", "Zinc", 0.45879, 330.0eV, 7.133E+00g / cm^3),
    Gallium = Element(31, "Ga", "Gallium", 0.44462, 334.0eV, 5.904E+00g / cm^3),
    Germanium =
        Element(32, "Ge", "Germanium", 0.44071, 350.0eV, 5.323E+00g / cm^3),
    Arsenic = Element(33, "As", "Arsenic", 0.44046, 347.0eV, 5.730E+00g / cm^3),
    Selenium =
        Element(34, "Se", "Selenium", 0.43060, 348.0eV, 4.500E+00g / cm^3),
    Bromine = Element(35, "Br", "Bromine", 0.43803, 343.0eV, 7.072E-03g / cm^3),
    Krypton = Element(36, "Kr", "Krypton", 0.42959, 352.0eV, 3.478E-03g / cm^3),
    Rubidium =
        Element(37, "Rb", "Rubidium", 0.43291, 363.0eV, 1.532E+00g / cm^3),
    Strontium =
        Element(38, "Sr", "Strontium", 0.43369, 366.0eV, 2.540E+00g / cm^3),
    Yttrium = Element(39, "Y", "Yttrium", 0.43867, 379.0eV, 4.469E+00g / cm^3),
    Zirconium =
        Element(40, "Zr", "Zirconium", 0.43848, 393.0eV, 6.506E+00g / cm^3),
    Niobium = Element(41, "Nb", "Niobium", 0.44130, 417.0eV, 8.570E+00g / cm^3),
    Molybdenum = Element(
        42,
        "Mo",
        "Molybdenum",
        0.43777,
        424.0eV,
        1.022E+01g / cm^3,
    ),
    Technetium = Element(
        43,
        "Tc",
        "Technetium",
        0.43919,
        428.0eV,
        1.150E+01g / cm^3,
    ),
    Ruthenium =
        Element(44, "Ru", "Ruthenium", 0.43534, 441.0eV, 1.241E+01g / cm^3),
    Rhodium = Element(45, "Rh", "Rhodium", 0.43729, 449.0eV, 1.241E+01g / cm^3),
    Palladium =
        Element(46, "Pd", "Palladium", 0.43225, 470.0eV, 1.202E+01g / cm^3),
    Silver = Element(47, "Ag", "Silver", 0.43572, 470.0eV, 1.050E+01g / cm^3),
    Cadmium = Element(48, "Cd", "Cadmium", 0.42700, 469.0eV, 8.650E+00g / cm^3),
    Indium = Element(49, "In", "Indium", 0.42676, 488.0eV, 7.310E+00g / cm^3),
    Tin = Element(50, "Sn", "Tin", 0.42120, 488.0eV, 7.310E+00g / cm^3),
    Antimony =
        Element(51, "Sb", "Antimony", 0.41889, 487.0eV, 6.691E+00g / cm^3),
    Tellurium =
        Element(52, "Te", "Tellurium", 0.40752, 485.0eV, 6.240E+00g / cm^3),
    Iodine = Element(53, "I", "Iodine", 0.41764, 491.0eV, 4.930E+00g / cm^3),
    Xenon = Element(54, "Xe", "Xenon", 0.41130, 482.0eV, 5.485E-03g / cm^3),
    Cesium = Element(55, "Cs", "Cesium", 0.41383, 488.0eV, 1.873E+00g / cm^3),
    Barium = Element(56, "Ba", "Barium", 0.40779, 491.0eV, 3.500E+00g / cm^3),
    Lanthanum =
        Element(57, "La", "Lanthanum", 0.41035, 501.0eV, 6.154E+00g / cm^3),
    Cerium = Element(58, "Ce", "Cerium", 0.41395, 523.0eV, 6.657E+00g / cm^3),
    Praseodymium = Element(
        59,
        "Pr",
        "Praseodymium",
        0.41871,
        535.0eV,
        6.710E+00g / cm^3,
    ),
    Neodymium =
        Element(60, "Nd", "Neodymium", 0.41597, 546.0eV, 6.900E+00g / cm^3),
    Promethium = Element(
        61,
        "Pm",
        "Promethium",
        0.42094,
        560.0eV,
        7.220E+00g / cm^3,
    ),
    Samarium =
        Element(62, "Sm", "Samarium", 0.41234, 574.0eV, 7.460E+00g / cm^3),
    Europium =
        Element(63, "Eu", "Europium", 0.41457, 580.0eV, 5.243E+00g / cm^3),
    Gadolinium = Element(
        64,
        "Gd",
        "Gadolinium",
        0.40699,
        591.0eV,
        7.900E+00g / cm^3,
    ),
    Terbium = Element(65, "Tb", "Terbium", 0.40900, 614.0eV, 8.229E+00g / cm^3),
    Dysprosium = Element(
        66,
        "Dy",
        "Dysprosium",
        0.40615,
        628.0eV,
        8.550E+00g / cm^3,
    ),
    Holmium = Element(67, "Ho", "Holmium", 0.40623, 650.0eV, 8.795E+00g / cm^3),
    Erbium = Element(68, "Er", "Erbium", 0.40655, 658.0eV, 9.066E+00g / cm^3),
    Thulium = Element(69, "Tm", "Thulium", 0.40844, 674.0eV, 9.321E+00g / cm^3),
    Ytterbium =
        Element(70, "Yb", "Ytterbium", 0.40453, 684.0eV, 6.730E+00g / cm^3),
    Lutetium =
        Element(71, "Lu", "Lutetium", 0.40579, 694.0eV, 9.840E+00g / cm^3),
    Hafnium = Element(72, "Hf", "Hafnium", 0.40338, 705.0eV, 1.331E+01g / cm^3),
    Tantalum =
        Element(73, "Ta", "Tantalum", 0.40343, 718.0eV, 1.665E+01g / cm^3),
    Tungsten =
        Element(74, "W", "Tungsten", 0.40250, 727.0eV, 1.930E+01g / cm^3),
    Rhenium = Element(75, "Re", "Rhenium", 0.40278, 736.0eV, 2.102E+01g / cm^3),
    Osmium = Element(76, "Os", "Osmium", 0.39958, 746.0eV, 2.257E+01g / cm^3),
    Iridium = Element(77, "Ir", "Iridium", 0.40058, 757.0eV, 2.242E+01g / cm^3),
    Platinum =
        Element(78, "Pt", "Platinum", 0.39984, 790.0eV, 2.145E+01g / cm^3),
    Gold = Element(79, "Au", "Gold", 0.40108, 790.0eV, 1.932E+01g / cm^3),
    Mercury = Element(80, "Hg", "Mercury", 0.39882, 800.0eV, 1.355E+01g / cm^3),
    Thallium =
        Element(81, "Tl", "Thallium", 0.39631, 810.0eV, 1.172E+01g / cm^3),
    Lead = Element(82, "Pb", "Lead", 0.39575, 823.0eV, 1.135E+01g / cm^3),
    Bismuth = Element(83, "Bi", "Bismuth", 0.39717, 823.0eV, 9.747E+00g / cm^3),
    Polonium =
        Element(84, "Po", "Polonium", 0.40195, 830.0eV, 9.320E+00g / cm^3),
    Astatine =
        Element(85, "At", "Astatine", 0.40479, 825.0eV, 1.000E+01g / cm^3),
    Radon = Element(86, "Rn", "Radon", 0.38736, 794.0eV, 9.066E-03g / cm^3),
    Francium =
        Element(87, "Fr", "Francium", 0.39010, 827.0eV, 1.000E+01g / cm^3),
    Radium = Element(88, "Ra", "Radium", 0.38934, 826.0eV, 5.000E+00g / cm^3),
    Actinium =
        Element(89, "Ac", "Actinium", 0.39202, 841.0eV, 1.007E+01g / cm^3),
    Thorium = Element(90, "Th", "Thorium", 0.38787, 847.0eV, 1.172E+01g / cm^3),
    Protactinium = Element(
        91,
        "Pa",
        "Protactinium",
        0.39388,
        878.0eV,
        1.537E+01g / cm^3,
    ),
    Uranium = Element(92, "U", "Uranium", 0.38651, 890.0eV, 1.895E+01g / cm^3),
)


text:
module MultivariateOrthogonalPolynomials
using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
      QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
      LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
      HarmonicOrthogonalPolynomials

import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size, oneto
import DomainSets: boundary

import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
import ContinuumArrays: @simplify, Weight, grid, plotgrid, TransformFactorization, Expansion

import ArrayLayouts: MemoryLayout
import BlockArrays: block, blockindex, BlockSlice, viewblock
import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
import LinearAlgebra: factorize
import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout
import LazyBandedMatrices: LazyBandedBlockBandedLayout
import InfiniteArrays: InfiniteCardinal

import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted
import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                          PartialDerivative, BlockOneTo, BlockRange1, interlace

export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian,
      MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
      zerniker, zernikez, Weighted, Block, ZernikeWeight, AbsLaplacianPower

include("disk.jl")
include("triangle.jl")

end # module


text:
"""
    OptSummary

Summary of an `NLopt` optimization

# Fields

* `initial`: a copy of the initial parameter values in the optimization
* `lowerbd`: lower bounds on the parameter values
* `ftol_rel`: as in NLopt
* `ftol_abs`: as in NLopt
* `xtol_rel`: as in NLopt
* `xtol_abs`: as in NLopt
* `initial_step`: as in NLopt
* `maxfeval`: as in NLopt
* `final`: a copy of the final parameter values from the optimization
* `fmin`: the final value of the objective
* `feval`: the number of function evaluations
* `optimizer`: the name of the optimizer used, as a `Symbol`
* `returnvalue`: the return value, as a `Symbol`
* `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs
* `REML`: use the REML criterion for LMM fits

The latter field doesn't really belong here but it has to be in a mutable struct in case it is changed.
"""
mutable struct OptSummary{T <: AbstractFloat}
    initial::Vector{T}
    lowerbd::Vector{T}
    finitial::T
    ftol_rel::T
    ftol_abs::T
    xtol_rel::T
    xtol_abs::Vector{T}
    initial_step::Vector{T}
    maxfeval::Int
    final::Vector{T}
    fmin::T
    feval::Int
    optimizer::Symbol
    returnvalue::Symbol
    nAGQ::Integer           # doesn't really belong here but I needed some place to store it
    REML::Bool
end
function OptSummary(initial::Vector{T}, lowerbd::Vector{T},
    optimizer::Symbol; ftol_rel::T=zero(T), ftol_abs::T=zero(T), xtol_rel::T=zero(T),
    initial_step::Vector{T}=T[]) where T <: AbstractFloat
    OptSummary(initial, lowerbd, T(Inf), ftol_rel, ftol_abs, xtol_rel, zero(initial),
        initial_step, -1, copy(initial), T(Inf), -1, optimizer, :FAILURE, 1, false)
end

function Base.show(io::IO, s::OptSummary)
    println(io, "Initial parameter vector: ", s.initial)
    println(io, "Initial objective value:  ", s.finitial)
    println(io)
    println(io, "Optimizer (from NLopt):   ", s.optimizer)
    println(io, "Lower bounds:             ", s.lowerbd)
    println(io, "ftol_rel:                 ", s.ftol_rel)
    println(io, "ftol_abs:                 ", s.ftol_abs)
    println(io, "xtol_rel:                 ", s.xtol_rel)
    println(io, "xtol_abs:                 ", s.xtol_abs)
    println(io, "initial_step:             ", s.initial_step)
    println(io, "maxfeval:                 ", s.maxfeval)
    println(io)
    println(io, "Function evaluations:     ", s.feval)
    println(io, "Final parameter vector:   ", s.final)
    println(io, "Final objective value:    ", s.fmin)
    println(io, "Return code:              ", s.returnvalue)
end

function NLopt.Opt(optsum::OptSummary)
    lb = optsum.lowerbd

    opt = NLopt.Opt(optsum.optimizer, length(lb))
    NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective
    NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective
    NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values
    if length(optsum.xtol_abs) == length(lb)  # not true for the second optimization in GLMM
        NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values
    end
    NLopt.lower_bounds!(opt, lb)
    NLopt.maxeval!(opt, optsum.maxfeval)
    if isempty(optsum.initial_step)
        optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb))
    else
        NLopt.initial_step!(opt, optsum.initial_step)
    end
    opt
end


text:
using Cbc
using QuadraticToBinary

CBC_OPTIMIZER = Cbc.Optimizer()
# MOI.set(CBC_OPTIMIZER, MOI.Silent(), false)
MOI.set(CBC_OPTIMIZER, MOI.TimeLimitSec(), 5)

CBC = CBC_OPTIMIZER

push!(solvers, (opt = CBC, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_sos, (opt = CBC, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_indicator, (opt = CBC, mode = BilevelJuMP.IndicatorMode()))

push!(solvers_fa, (opt = CBC, mode = BilevelJuMP.FortunyAmatMcCarlMode(with_slack = false)))
push!(solvers_fa, (opt = CBC, mode = BilevelJuMP.FortunyAmatMcCarlMode(with_slack = true)))
push!(solvers_fa2, (opt = CBC, mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100)))

CBC_CACHED = MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), CBC_OPTIMIZER)
push!(solvers_cached, (opt = CBC_CACHED, mode = BilevelJuMP.SOS1Mode()))

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(CBC_CACHED),
    mode = BilevelJuMP.StrongDualityMode(1e-9, inequality = true)))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(CBC_CACHED),
    mode = BilevelJuMP.StrongDualityMode(inequality = false)))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

QTB_CBC_BRIDGED = MOI.Bridges.full_bridge_optimizer(
    QuadraticToBinary.Optimizer{Float64}(CBC_CACHED, lb=-100,ub=100),
    Float64)
MOI.Bridges.add_bridge(QTB_CBC_BRIDGED, MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64})
MOI.set(QTB_CBC_BRIDGED, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

push!(solvers_sos_quad_bin, (
    opt = QTB_CBC_BRIDGED,
    mode = BilevelJuMP.SOS1Mode()))

push!(solvers_fa_quad_bin, (
    opt = QTB_CBC_BRIDGED,
    mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100)))

QTB_CBC = QuadraticToBinary.Optimizer{Float64}(CBC_CACHED, lb=-5, ub=5)
QTB_CBC_B = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(QTB_CBC)
push!(solvers_fa_quad_bin_mixed, (
    opt = QTB_CBC_B,
    mode = BilevelJuMP.MixedMode(default=
        BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 10, dual_big_M = 10))))
MOI.set(QTB_CBC_B, QuadraticToBinary.GlobalVariablePrecision(), 1e-3)


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

mutable struct ServiceEndpointPolicyPropertiesFormat <: SwaggerModel
    serviceEndpointPolicyDefinitions::Any # spec type: Union{ Nothing, Vector{ServiceEndpointPolicyDefinition} } # spec name: serviceEndpointPolicyDefinitions
    subnets::Any # spec type: Union{ Nothing, Vector{Subnet} } # spec name: subnets
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ServiceEndpointPolicyPropertiesFormat(;serviceEndpointPolicyDefinitions=nothing, subnets=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(ServiceEndpointPolicyPropertiesFormat, Symbol("serviceEndpointPolicyDefinitions"), serviceEndpointPolicyDefinitions)
        setfield!(o, Symbol("serviceEndpointPolicyDefinitions"), serviceEndpointPolicyDefinitions)
        validate_property(ServiceEndpointPolicyPropertiesFormat, Symbol("subnets"), subnets)
        setfield!(o, Symbol("subnets"), subnets)
        validate_property(ServiceEndpointPolicyPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(ServiceEndpointPolicyPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ServiceEndpointPolicyPropertiesFormat

const _property_map_ServiceEndpointPolicyPropertiesFormat = Dict{Symbol,Symbol}(Symbol("serviceEndpointPolicyDefinitions")=>Symbol("serviceEndpointPolicyDefinitions"), Symbol("subnets")=>Symbol("subnets"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ServiceEndpointPolicyPropertiesFormat = Dict{Symbol,String}(Symbol("serviceEndpointPolicyDefinitions")=>"Vector{ServiceEndpointPolicyDefinition}", Symbol("subnets")=>"Vector{Subnet}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ServiceEndpointPolicyPropertiesFormat }) = collect(keys(_property_map_ServiceEndpointPolicyPropertiesFormat))
Swagger.property_type(::Type{ ServiceEndpointPolicyPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceEndpointPolicyPropertiesFormat[name]))}
Swagger.field_name(::Type{ ServiceEndpointPolicyPropertiesFormat }, property_name::Symbol) =  _property_map_ServiceEndpointPolicyPropertiesFormat[property_name]

function check_required(o::ServiceEndpointPolicyPropertiesFormat)
    true
end

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


text:
using Test
using StatsFuns
using Base.Iterators: take
using Random
using LinearAlgebra
using DynamicIterators: trace, TimeLift
using TransformVariables: transform, as𝕀, inverse

using MeasureTheory
using MeasureTheory: Const

using Aqua
Aqua.test_all(MeasureTheory; ambiguities=false, unbound_args=false)

function draw2(μ)
    x = rand(μ)
    y = rand(μ)
    while x == y
        y = rand(μ)
    end
    return (x,y)
end

@testset "Parameterized Measures" begin
    @testset "Binomial" begin
        D = Binomial{(:n, :p)}
        par = merge((n=20,),transform(asparams(D, (n=20,)), randn(1)))
        d = D(par)
        (n,p) = (par.n, par.p)
        logitp = logit(p)
        probitp = norminvcdf(p)
        y = rand(d)

        ℓ = logdensity(Binomial(;n, p), y)
        @test ℓ ≈ logdensity(Binomial(;n, logitp), y)
        @test ℓ ≈ logdensity(Binomial(;n, probitp), y)

        @test_broken logdensity(Binomial(n,p), CountingMeasure(ℤ[0:n]), x) ≈ binomlogpdf(n,p,x)
    end

    @testset "NegativeBinomial" begin
        D = NegativeBinomial{(:r, :p)}
        par = transform(asparams(D), randn(2))
        d = D(par)
        (r,p) = (par.r, par.p)
        logitp = logit(p)
        λ = p * r / (1 - p)
        logλ = log(λ)
        y = rand(d)

        ℓ = logdensity(NegativeBinomial(;r, p), y)
        @test ℓ ≈ logdensity(NegativeBinomial(;r, logitp), y)
        @test ℓ ≈ logdensity(NegativeBinomial(;r, λ), y)
        @test ℓ ≈ logdensity(NegativeBinomial(;r, logλ), y)

        sample1 = rand(MersenneTwister(123), NegativeBinomial(;r, λ))
        sample2 = rand(MersenneTwister(123), NegativeBinomial(;r, logλ))
        @test sample1 == sample2

        @test_broken logdensity(Binomial(n,p), CountingMeasure(ℤ[0:n]), x) ≈ binomlogpdf(n,p,x)
    end

    @testset "Poisson" begin
        sample1 = rand(MersenneTwister(123), Poisson(;logλ = log(100)))
        sample2 = rand(MersenneTwister(123), Poisson(;λ = 100))
        @test sample1 == sample2
    end

    @testset "Normal" begin
        D = Normal{(:μ,:σ)}
        par = transform(asparams(D), randn(2))
        d = D(par)
        @test params(d) == par

        μ = par.μ
        σ = par.σ
        σ² = σ^2
        τ = 1/σ²
        logσ = log(σ)
        y = rand(d)

        ℓ = logdensity(Normal(;μ,σ), y)
        @test ℓ ≈ logdensity(Normal(;μ,σ²), y)
        @test ℓ ≈ logdensity(Normal(;μ,τ), y)
        @test ℓ ≈ logdensity(Normal(;μ,logσ), y)
    end

    @testset "LKJCholesky" begin
        D = LKJCholesky{(:k,:η)}
        par = transform(asparams(D, (k=4,)), randn(1))
        d = D(merge((k=4,),par))
        # @test params(d) == par

        η  = par.η
        logη = log(η)

        y = rand(d)
        η = par.η
        ℓ = logdensity(LKJCholesky(4,η), y)
        @test ℓ ≈ logdensity(LKJCholesky(k=4,logη=logη), y)
    end
end

@testset "Kernel" begin
    κ = MeasureTheory.kernel(MeasureTheory.Dirac, identity)
    @test rand(κ(1.1)) == 1.1
end

@testset "SpikeMixture" begin
    @test rand(SpikeMixture(Dirac(0), 0.5)) == 0
    @test rand(SpikeMixture(Dirac(1), 1.0)) == 1
    w = 1/3
    m = SpikeMixture(Normal(), w)
    bm = basemeasure(m)
    @test (bm.s*bm.w)*bm.m == 1.0*basemeasure(Normal())
    @test density(m, 1.0)*(bm.s*bm.w) == w*density(Normal(),1.0)
    @test density(m, 0)*(bm.s*(1-bm.w)) ≈ (1-w)
end

@testset "Dirac" begin
    @test rand(Dirac(0.2)) == 0.2
    @test logdensity(Dirac(0.3), 0.3) == 0.0
    @test logdensity(Dirac(0.3), 0.4) == -Inf
end

@testset "For" begin
    FORDISTS = [
        For(1:10) do j Normal(μ=j) end
        For(4,3) do μ,σ Normal(μ,σ) end
        For(1:4, 1:4) do μ,σ Normal(μ,σ) end
        For(eachrow(rand(4,2))) do x Normal(x[1], x[2]) end
        For(rand(4), rand(4)) do μ,σ Normal(μ,σ) end
    ]

    for d in FORDISTS
        @test logdensity(d, rand(d)) isa Float64
    end
end

import MeasureTheory.:⋅
function ⋅(μ::Normal, kernel)
    m = kernel(μ)
    Normal(μ = m.μ.μ, σ = sqrt(m.μ.σ^2 + m.σ^2))
end
struct AffineMap{S,T}
    B::S
    β::T
end
(a::AffineMap)(x) = a.B*x + a.β
(a::AffineMap)(p::Normal) = Normal(μ = a.B*mean(p) + a.β, σ = sqrt(a.B*p.σ^2*a.B'))

@testset "DynamicFor" begin
    mc = Chain(Normal(μ=0.0)) do x Normal(μ=x) end
    r = rand(mc)

    # Check that `r` is now deterministic
    @test logdensity(mc, take(r, 100)) == logdensity(mc, take(r, 100))

    d2 = For(r) do x Normal(μ=x) end

    @test_broken let r2 = rand(d2)
        logdensity(d2, take(r2, 100)) == logdensity(d2, take(r2, 100))
    end
end

@testset "Univariate chain" begin
    ξ0 = 1.
    x = 1.2
    P0 = 1.0

    Φ = 0.8
    β = 0.1
    Q = 0.2

    μ = Normal(μ=ξ0, σ=sqrt(P0))
    kernel = MeasureTheory.kernel(Normal; μ=AffineMap(Φ, β), σ=Const(Q))

    @test (μ ⋅ kernel).μ == Normal(μ = 0.9, σ = 0.824621).μ

    chain = Chain(kernel, μ)

    dyniterate(iter::TimeLift, ::Nothing) = dyniterate(iter, 0=>nothing)
    tr1 = trace(TimeLift(chain), nothing, u -> u[1] > 15)
    tr2 = trace(TimeLift(rand(Random.GLOBAL_RNG, chain)), nothing, u -> u[1] > 15)
    collect(Iterators.take(chain, 10))
    collect(Iterators.take(rand(Random.GLOBAL_RNG, chain), 10))
end

@testset "Transforms" begin
    t = as𝕀
    @testset "Pushforward" begin
        μ = Normal()
        ν = Pushforward(t, μ)
        x = rand(μ)
        @test logdensity(μ, x) ≈ logdensity(Pushforward(inverse(t), ν), x)
    end

    @testset "Pullback" begin
        ν = Uniform()
        μ = Pullback(t,ν)
        y = rand(ν)
        @test logdensity(ν, y) ≈ logdensity(Pullback(inverse(t), μ), y)
    end
end

using TransformVariables

@testset "Likelihood" begin
    dps = [
        (Normal()                             ,    2.0  )
        # (Pushforward(as((μ=asℝ,)), Normal()^1), (μ=2.0,))
    ]

    ℓs = [
        Likelihood(Normal{(:μ,)},              3.0)
        Likelihood(kernel(Normal, x -> (μ=x, σ=2.0)), 3.0)
    ]

    for (d,p) in dps
        for ℓ in ℓs
            @test logdensity(d ⊙ ℓ, p) == logdensity(d, p) + logdensity(ℓ, p)
        end
    end
end

@testset "Reproducibility" begin

    function repro(D, args, nt=NamedTuple())
        t = asparams(D{args}, nt)
        d = D(transform(t, randn(t.dimension)))
        r(d) = rand(Random.MersenneTwister(1), d)
        logdensity(d, r(d)) == logdensity(d, r(d))
    end

    @testset "Bernoulli" begin
        @test repro(Bernoulli, (:p,))
    end
    @testset "Binomial" begin
        @test repro(Binomial, (:n,:p), (n=10,))
    end

    @testset "Beta" begin
        @test repro(Beta, (:α,:β))
    end

    @testset "Cauchy" begin
        @test repro(Cauchy, (:μ,:σ))
    end

    @testset "Dirichlet" begin
        @test_broken repro(Dirichlet, (:p,))
    end

    @testset "Exponential" begin
        @test repro(Exponential, (:λ,))
    end

    @testset "Gumbel" begin
        @test repro(Gumbel, (:μ,:σ))
    end

    @testset "InverseGamma" begin
        @test_broken repro(InverseGamma, (:p,))
    end

    @testset "Laplace" begin
        @test repro(Laplace, (:μ,:σ))
    end

    @testset "LKJCholesky" begin
        @test repro(LKJCholesky, (:k,:η,), (k=3,))
    end

    @testset "Multinomial" begin
        @test_broken repro(Multinomial, (:n,:p,))
    end

    @testset "MvNormal" begin
        @test_broken repro(MvNormal, (:μ,))
    end

    @testset "NegativeBinomial" begin
        @test repro(NegativeBinomial, (:r, :p))
    end

    @testset "Normal" begin
        @test repro(Normal, (:μ,:σ))
    end

    @testset "Poisson" begin
        @test repro(Poisson, (:λ,))
    end

    @testset "StudentT" begin
        @test repro(StudentT, (:ν, :μ))
    end

    @testset "Uniform" begin
        @test repro(Uniform, ())
    end

end


text:
using PackageCompiler
create_sysimage(:Plots;
                precompile_execution_file = "warmup.jl",
                replace_default = true,
                cpu_target = PackageCompiler.default_app_cpu_target())


text:
# ------------------------------- number utils ------------------------------- #

"""
    toint(x)

Compact Float64 -> Int64 conversion
"""
toint(x::Float64) = Int(round(x))
toint(x::AbstractVector) = Int.(round.(x))

"""
Linear interpolation.
"""
lerp(x1::AbstractVector, x2::AbstractVector, p::Float64) = (1 - p) .* x1 .+ p .* x2

"""
Linear interpolation between two colors.
"""
function lerp(c1::AbstractColor, c2::AbstractColor, p::Float64)
    c1, c2  = Vector(RGB(c1)), Vector(RGB(c2))
    RGB(toint(lerp(c1, c2, p)))
end

# ------------------------- number/string conversion ------------------------- #

"""
    int2hex

Converts a single integer to hex string, but with a minimal hex
string length of 2 for hex color codes.
"""
int2hex(i::Int) = i > 15 ? string(i, base=16) : string("0", string(i, base=16))

"""
    hex2ing("ff")

Converts a string with a hex number to an integer
"""
hex2int(hex::AbstractString) = parse(Int, hex, base=16)

# -------------------------- color types conversion -------------------------- #

"""
    rgb2hex(r, g, b)

Converta a RGB color to a Hex color.
"""
rgb2hex(rgb::RGB) = Hex(string("#", int2hex.(rgb)...))

"""
    rgb2hex(r, g, b)

Converta a `rgb` tuple of integers ∈ [0, 255] to a hex color code.
"""
rgb2hex(r::Number, g::Number, b::Number) = rgb2hex(RGB(r, g, b))

"""
    hex(Hex)

Converts a Hex string to a RGB color.
"""
hex2rgb(hex::Hex) = RGB([hex2int(hex.code[i:i+1]) for i in (1, 3, 5)])

"""
    hex("#ffffff")

Converts a hex string color code to a RGB color.
"""
hex2rgb(hex::AbstractString) = hex2rgb(Hex(hex))

# more constructors utils
RGB(hex::Hex) = hex2rgb(hex)
RGB(rgb::RGB) = rgb

Hex(c::RGB) = rgb2hex(c::RGB)
Hex(h::Hex) = h


text:
using LazIO
using LasIO
using FileIO
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using DataFrames
using PointCloudRasterizers

"Store coordinates"
function coordinates(points,h)
    coords = SArray{Tuple{9},Float64,1,9}[]
    for (ind,p) in enumerate(points)
        x = xcoord(p,h)
        y = ycoord(p,h)
        z = zcoord(p,h)
        inten = intensity(p)
        ret = return_number(p)
        nret = number_of_returns(p)
        classificc = classification(p)
        gps = gps_time(p)
        push!(coords,(ind, x, y, z, inten, ret, nret, classificc, gps))
    end
    return coords
end

"Convert Cartesian to Spherical coordinates"
function cart2pol(x::Float64, y::Float64, z::Float64)
    cart = Tuple{Float64,Float64,Float64}
    ρ = sqrt(x^2 + y^2 + z^2) #distance
    θ = atand(y,x) #horizontal angle
    ϕ = acosd(z/ρ) #incidence angle
    cart = (ρ,θ,ϕ)
    return cart
end

"Convert Spherical to Cartesian coordinates"
function pol2cart(ρ::Float64, θ::Float64, ϕ::Float64)
    pol = Tuple{Float64,Float64,Float64}
    x = (ρ * sind(ϕ) * cosd(θ))
    y = (ρ * sind(ϕ) * sind(θ))
    z = (ρ * cosd(ϕ))
    return pol = (x,y,z)
end

const c_air, c_water = 2.99*10^8, 2.25*10^8
const n_air_water = c_water/c_air

"Correct an underwater point according to a watersurface point from the refraction factor by using Snell's Law."
function refraction_correction(watersurface::SArray{Tuple{4},Float64,1,4},underwater::SArray{Tuple{4},Float64,1,4})
    x0, y0, z0 = watersurface[2:4]
    x1, y1, z1 = underwater[2:4]
    Δx, Δy, Δz = (x1-x0), (y1-y0), (z1-z0)
    ind = underwater[1]
    #Calculate the distance of the two points (ρ), the vertical-incidence angle (ϕ) and horizontal angle (θ) between the two points
    ρ,θ,ϕ = cart2pol(Δx, Δy, Δz)

    #Calculate the underwater angle due to the refraction effect.
    #Using Snell's Law: n(air)*sin(ϕ_air) = n(water)*sin(ϕ_water),
    #where ϕ_air: the angle of incidence (I know that the angle of the scan is 20)
    #      ϕ_water: the angle of refraction
    #      n(air): the refraction idex of medium containing the incident ray, value=1.000293
    #      n(water): the refraction idex of medium containing the transmitted ray, value=1.333

    #check the height values of the water-surface and underwater point
    if z0 == z1
        error("the water-surface and underwater point are at the same height level")
    elseif z0 > z1
        ϕ_water = asind(n_air_water*sind(ϕ)) #ϕ_new is the new vertical
    else
        error("the underwater point can not be higher than the watersurface point")
    end

    #Distance under the water has been influenced by the different speed of light in the water.
    #In a specific moment, ρ = c_air * t and ρ_water = c_water * t.
    #So, by diving these two equations, the ratio (ρ1/ρ2) = n_air_water
    ρ_water = (ρ*n_air_water)

    #Calculate the correction in the 3d space using the pol2cart function
    xcor, ycor, zcor = pol2cart(ρ_water, θ, ϕ_water)

    #Define the new corrected coordinates of the points
    xnew = x0 + xcor
    ynew = y0 + ycor
    znew = z0 - zcor

    #Return the new corrected point
    corrections = @SVector[xcor, ycor, zcor]
    correct_point = @SVector[ind, xnew, ynew ,znew]

    return corrections, correct_point
end

"Group the laser pulses based on their GPS time"
function laserpulses(points)
    println("Execution Time of laser pulses " * string(length(points)) * " points")
    @time begin
        laser_pulses = SArray{Tuple{}}[] #specify the type of array
        sorted_gps = sort!(points, by = x -> x[9], rev = false)
        pp = [tuple(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9]) for p in sorted_gps]
        df = DataFrame(pp[:]) #create a dataframe
        groups = groupby(df,:9) #group by gps time
        for gr in groups
            if size(gr,1) > 1
                push!(laser_pulses,[gr])
            end
        end
    end
    return laser_pulses
end

"Calculate the corrected/refracted points for every pulse"
function corrected_points(laser_pulses,coords)
    corrected_laser_pulses = [] #specify the type of array

    INDEXES_corr = [] #indexes of corrected_points
    INDEXES_un_corr = [] #indexes of non-corrected points

    for every_pulse in laser_pulses
        indexes = getindex(every_pulse)[1]
        ind_surf = Int64(indexes[1]) #water surface index
        p_surf = coords[ind_surf] #water surface point
        p_surf = @SVector[ind_surf,p_surf[2],p_surf[3],p_surf[4]] #create Arrays
        corrected_points = SArray{Tuple{4},Float64,1,4}[] #store the corrected points for each pulse!
        not_corrected = SArray{Tuple{4},Float64,1,4}[] #store the un-corrected points due to the Exception ERRORsfor each pulse!

        push!(corrected_points,(p_surf)) # add the water surface point
        push!(not_corrected,(p_surf)) # add the water surface point

        for ind_under = 2:length(indexes) #iterate throught the inder water points
            n = Int64(indexes[ind_under])
            p_under = coords[n]
            p_under = @SVector[n,p_under[2],p_under[3],p_under[4]]
            try
                correction, corre_pt = refraction_correction(p_surf,p_under)
                push!(INDEXES_corr,p_under[1])
                push!(corrected_points,(corre_pt)) # I will classify as un-refracted with classif.code = 19
            catch
                push!(INDEXES_un_corr,p_under[1])
                push!(not_corrected,(p_under)) # I will classify as un-refracted with classif.code = 20
            end
        end

        if length(corrected_points) > 1 #check if there is NOT only water surface point
            push!(corrected_laser_pulses,[corrected_points])
        end
    end
    return corrected_laser_pulses, INDEXES_corr, INDEXES_un_corr
end

#Write a new classified Pointcloud
function write_pointcloud(ds,header,pp,indexes_corr,index_un_corr)
    n = length(indexes_corr) + length(index_un_corr)
    println("Execution Time of writing points " * string(n) * " points")
    @time  begin
        laz_out = joinpath(path, "myoutput_" * string(length(pp)) * "_points.laz")
        LazIO.write(laz_out, ds.header) do io
            for (ind,p) in enumerate(ds)
                if ind in indexes_corr
                    p.classification = UInt8(19)
                elseif ind in indexes_un_corr #this elseif doesn't work!!why???
                    p.classification = UInt8(20)
                else
                    p.classification = UInt8(9)
                end
                LazIO.writepoint(io, p)
            end
        end
    end
end

#MAIN#
const path = "C:/Users/alexandr/OneDrive - Stichting Deltares/Desktop/NL3_subset/"
const filename = "NL3_clipped_water_cropheight_pointformat3_classification.laz" #classification = 9 (water)
const lazinput = path * filename

#Open LAZ pointcloud
open_file = File{format"LAZ_"}(lazinput)
h, points = load(open_file)
n = length(points)

#dataset
ds = LazIO.open(lazinput)

#subset
pp = points[1500000:1520000]
#Array with coordinates
coords = coordinates(pp,h)
#Separate laser pulses
laser_pulses = laserpulses(coords)
#Return corrected points
corrected_laser_pulses, indexes_corr, indexes_un_corr  = corrected_points(laser_pulses,coords)
#write pointcloud
write_pointcloud(ds,h,pp,indexes_corr,indexes_un_corr)


text:
using InstantiateFromURL
activate_github("QuantEcon/QuantEconLectureAllPackages", tag = "v0.9.1") # activate the QuantEcon environment

using LinearAlgebra, Statistics, Compat # load common packages

max_stars_repo_name: yywingliang/lecture-source-jl

max_issues_repo_name: yywingliang/lecture-source-jl

max_forks_repo_name: yywingliang/lecture-source-jl


text:
CODE_EQ(s::State) = if length(s.code) > 1
  push!(s.boolean, pop!(s.code) === pop!(s.code))
end

CODE_APPEND(s::State) = if length(s.code) > 1
  first = pop!(s.code)
  second = peek(s.code)
  s.code[end] = vcat(isa(first, Vector) ? first : [first],
    isa(second, Vector) ? second : [second])
end

CODE_ATOM(s::State) = if !isempty(s.code)
  push!(s.boolean, !isa(pop!(s.code), Vector))
end

CODE_CAR(s::State) = if !isempty(s.code) && isa(s.code[end], Vector)
  s.code[end] = isempty(s.code[end]) ? [] : s.code[end][1]
end

CODE_CDR(s::State) = if !isempty(s.code)
  top = peek(s.code)
  if !isa(top, Vector) || length(top) < 2
    s.code[end] = []
  else
    s.code[end] = top[2:end]
  end
end

CODE_CONS(s::State) = if length(s.code) > 1
  top = pop!(s.code)
  s.code[end] = vcat([peek(s.code)], isa(top, Vector) ? top : [top])
end

CODE_CONTAINER(s::State) = if length(s.code) > 1
  s.code[end] = container(pop!(s.code), peek(s.code))
end

CODE_CONTAINS(s::State) = if length(s.code) > 1
  first = pop!(s.code)
  second = pop!(s.code)
  push!(s.boolean, !isempty(container(second, first)))
end

CODE_DISCREPANCY(s::State) = if length(s.code) > 1
  a = isa(peek(s.code), Vector) ? pop!(s.code) : [pop!(s.code)]
  b = isa(peek(s.code), Vector) ? pop!(s.code) : [pop!(s.code)]
  dis = zero(Int32)

  for v in setdiff(a, b)
    dis += count(x -> x == v, a)
  end
  for v in setdiff(b, a)
    dis += count(x -> x == v, b)
  end

  push!(s.integer, dis)
end

CODE_DO(s::State) = if !isempty(s.code)
  push!(s.exec, Symbol("CODE.POP"), peek(s.code))
end

CODE_DO_STAR(s::State) = if !isempty(s.code)
  push!(s.exec, pop!(s.code))
end

CODE_DO_STAR_COUNT(s::State) = if !isempty(s.code) && !isempty(s.integer) && s.integer[end] > 0
  push!(s.exec, pop!(s.code), Symbol("EXEC.DO*COUNT"))
end

CODE_DO_STAR_TIMES(s::State) = if !isempty(s.code) && !isempty(s.integer) && s.integer[end] > 0
  push!(s.exec, pop!(s.code), Symbol("EXEC.DO*TIMES"))
end

CODE_DO_STAR_RANGE(s::State) = if !isempty(s.code) && length(s.integer) > 1
  push!(s.exec, pop!(s.code), Symbol("EXEC.DO*RANGE"))
end

CODE_DUP(s::State) = if !isempty(s.code)
  push!(s.code, peek(s.code))
end

CODE_EXTRACT(s::State) = if !isempty(s.code) && !isempty(s.integer)
  index = abs(pop!(s.integer))
  code = peek(s.code)

  # Do nothing (other than pop from the INTEGER stack) if the top of the CODE
  # stack is an atom.
  if !isa(code, Vector)
    return
  end

  # Calculate the index modulo the number of points in the given expression.
  index %= num_points(code)

  # Depth-first search for the given index.
  q = Any[code]
  j = 0
  while !isempty(q)
    expr = pop!(q)

    # Put the extracted sub-expression on top of the CODE stack.
    if j == index
      s.code[end] = expr
      break
    end

    # Add sub-expressions to the front of the processing queue, from
    # left to right.
    if isa(expr, Vector)
      for k in length(expr):-1:1
        push!(q, expr[k])
      end
    end

    j += 1
  end
end

CODE_FLUSH(s::State) = empty!(s.code)

CODE_FROM_BOOLEAN(s::State) = if !isempty(s.boolean)
  push!(s.code, pop!(s.boolean))
end

CODE_FROM_FLOAT(s::State) = if !isempty(s.float)
  push!(s.code, pop!(s.float))
end

CODE_FROM_INTEGER(s::State) = if !isempty(s.integer)
  push!(s.code, pop!(s.integer))
end

CODE_FROM_NAME(s::State) = if !isempty(s.name)
  push!(s.code, pop!(s.name))
end

CODE_IF(s::State) = if !isempty(s.boolean) && length(s.code) > 1
  a = pop!(s.code)
  b = pop!(s.code)
  push!(s.exec, pop!(s.boolean) ? b : a)
end

CODE_INSERT(s::State) = if !isempty(s.integer) && length(s.code) > 1
  cntr = pop!(s.code)
  item = pop!(s.code)

  # Calculate the insertion point.
  i = abs(pop!(s.integer)) % num_points(cntr)

  # Push the modified container back onto the stack.
  push!(s.code, insert_at_point!(cntr, item, i))
end

#
# TODO
#
CODE_INSTRUCTIONS(s::State) = return

CODE_LENGTH(s::State) = if !isempty(s.code)
  top = pop!(s.code)
  push!(s.integer, isa(top, Vector) ? length(top) : 1)
end

CODE_LIST(s::State) = if length(s.code) > 1
  s.code[end] = [pop!(s.code), peek(s.code)]
end

CODE_MEMBER(s::State) = if length(s.code) > 1
  first = isa(peek(s.code), Vector) ? pop!(s.code) : [pop!(s.code)]
  push!(s.boolean, in(pop!(s.code), first))
end

CODE_NOOP(s::State) = return

CODE_NTH(s::State) = if !isempty(s.integer) && !isempty(s.code) && isa(s.code[end], Vector)
  if isempty(peek(s.code))
    s.code[end] = []
  else
    s.code[end] = s.code[end][1 + (abs(pop!(s.integer)) % length(peek(s.code)))]
  end
end

# Zeroth CDR is the first CDR in this implementation.
# Seems to be implied by the Push 3.0 language reference.
CODE_NTH_CDR(s::State) = if !isempty(s.code) && !isempty(s.integer)
  i = pop!(s.integer)
  if !isa(peek(s.code), Vector) || isempty(peek(s.code))
    s.code[end] = []
  else
    i = abs(i) % length(peek(s.code))
    deleteat!(s.code[end], 1:i+1)
  end
end

CODE_NULL(s::State) = if !isempty(s.code)
  top = pop!(s.code)
  push!(s.boolean, isa(top, Vector) && isempty(top))
end

CODE_POP(s::State) = !isempty(s.code) && pop!(s.code)

#
# TODO
#
CODE_POSITION(s::State) = return

CODE_QUOTE(s::State) = s.flag_quote_code = true

CODE_ROT(s::State) = if length(s.code) > 2
  s.code[end], s.code[end-1], s.code[end-2] =
    s.code[end-1], s.code[end-2], s.code[end]
end

CODE_SHOVE(s::State) = if !isempty(s.integer) && !isempty(s.code)
  shove!(s.code, pop!(s.integer))
end

CODE_SIZE(s::State) = if !isempty(s.code)
  push!(s.integer, num_points(pop!(s.code)))
end

CODE_STACK_DEPTH(s::State) = push!(s.integer, length(s.code))

#
# TODO
#
CODE_SUBST(s::State) = return

CODE_SWAP(s::State) = if length(s.code) >= 2
  s.code[end], s.code[end-1] = s.code[end-1], s.code[end]
end

CODE_YANK(s::State) = if !isempty(s.integer) && !isempty(s.code)
  yank!(s.code, pop!(s.integer))
end

CODE_YANK_DUP(s::State) = if !isempty(s.integer) && !isempty(s.code)
  yankdup!(s.code, pop!(s.integer))
end

CODE_RAND(s::State) = push!(s.code, random_code(s))

CODE_DEFINE(s::State) = if !isempty(s.code) && !isempty(s.name)
  s.definitions[pop!(s.name)] = pop!(s.code)
end

CODE_DEFINITION(s::State) = if !isempty(s.name) && haskey(s.definitions, pop!(s.name))
  push!(s.code, deepcopy(s.definitions[pop!(s.name)]))
end

Push.register("CODE.=",                CODE_EQ)
Push.register("CODE.APPEND",           CODE_APPEND)
Push.register("CODE.ATOM",             CODE_ATOM)
Push.register("CODE.CAR",              CODE_CAR)
Push.register("CODE.CDR",              CODE_CDR)
Push.register("CODE.CONS",             CODE_CONS)
Push.register("CODE.CONTAINER",        CODE_CONTAINER)
Push.register("CODE.CONTAINS",         CODE_CONTAINS)
Push.register("CODE.DEFINE",           CODE_DEFINE)
Push.register("CODE.DEFINITION",       CODE_DEFINITION)
Push.register("CODE.DISCREPANCY",      CODE_DISCREPANCY)
Push.register("CODE.DO",               CODE_DO)
Push.register("CODE.DO*",              CODE_DO_STAR)
Push.register("CODE.DO*COUNT",         CODE_DO_STAR_COUNT)
Push.register("CODE.DO*RANGE",         CODE_DO_STAR_RANGE)
Push.register("CODE.DO*TIMES",         CODE_DO_STAR_TIMES)
Push.register("CODE.DUP",              CODE_DUP)
Push.register("CODE.EXTRACT",          CODE_EXTRACT)
Push.register("CODE.FLUSH",            CODE_FLUSH)
Push.register("CODE.FROMBOOLEAN",      CODE_FROM_BOOLEAN)
Push.register("CODE.FROMINTEGER",      CODE_FROM_INTEGER)
Push.register("CODE.FROMNAME",         CODE_FROM_NAME)
Push.register("CODE.IF",               CODE_IF)
Push.register("CODE.INSERT",           CODE_INSERT)
Push.register("CODE.INSTRUCTIONS",     CODE_INSTRUCTIONS)
Push.register("CODE.LENGTH",           CODE_LENGTH)
Push.register("CODE.LIST",             CODE_LIST)
Push.register("CODE.MEMBER",           CODE_MEMBER)
Push.register("CODE.NOOP",             CODE_NOOP)
Push.register("CODE.NTH",              CODE_NTH)
Push.register("CODE.NTHCDR",           CODE_NTH_CDR)
Push.register("CODE.NULL",             CODE_NULL)
Push.register("CODE.POP",              CODE_POP)
Push.register("CODE.POSITION",         CODE_POSITION)
Push.register("CODE.QUOTE",            CODE_QUOTE)
Push.register("CODE.RAND",             CODE_RAND)
Push.register("CODE.ROT",              CODE_ROT)
Push.register("CODE.SHOVE",            CODE_SHOVE)
Push.register("CODE.SIZE",             CODE_SIZE)
Push.register("CODE.STACKDEPTH",       CODE_STACK_DEPTH)
Push.register("CODE.SUBST",            CODE_SUBST)
Push.register("CODE.SWAP",             CODE_SWAP)
Push.register("CODE.YANK",             CODE_YANK)
Push.register("CODE.YANKDUP",          CODE_YANK_DUP)


text:
# name = "BARON"
# add BARON

using BARON

BARON_OPT = BARON.Optimizer()
# MOI.set(BARON_OPT, MOI.Silent(), true)
# MOI.set(BARON_OPT, MOI.TimeLimitSec(), 5)
# BARON_OPT = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(BARON_OPT)

push!(solvers_nlp, (opt = BARON_OPT, mode = BilevelJuMP.ProductMode(1e-5)))


text:
#=
Catch all for all the parameters that makes this easy to use;
=#

struct SystemParameters
    mode_tag::String
    unknown_dynamics_fcn
    known_dynamics_fcn
    measurement_noise_dist
    process_noise_dist
    lipschitz_bound::Float64
    dependency_dims::Dict
    n_dims_in::Int
    n_dims_out::Int
end

mutable struct DataParameters
    data_num::Int
    data_frac_optimize::Float64
    bound_type::String
    noise_sigma::Float64
    epsilon::Float64
    eta::Float64
    safety_dims
    local_radius::Float64
    num_neighbors::Int
end

mutable struct ExperimentParameters
    experiment_directory::String
    system_tag::String
    specification_file::String
    domain::Dict
    discretization_step::Dict
    verification_mode::String
    random_seed::Int
    system_params::SystemParameters
    data_params
end

function save_experiment_params(params::ExperimentParameters)
	# TODO: Finish this
end

function load_experiment_params(filename::String)
	# TODO: Finish this
end


text:
"""
    compute_delta!(g::AbstractGraph)

Determines `t₀` and `t₁` for the given evolved `g`

Arguments
* `g` : An evolved instance of type AbstractGraph
Extra Information
* `t₀`: The step in the evolution process where the heterogeneity peaks
* `t₁`: The step in the evolution process where the size of the largest cluster increases the most
Returns
* `nothing`, updates `g` in-place
"""
function compute_delta!(g::AbstractGraph)
    t₀ = argmax(g.observables.heterogeneity) - 1 # because t+1 observations when t edges active
    t₁ = argmax([
        g.observables.largest_cluster_size[i+1] - g.observables.largest_cluster_size[i]
        for i = 1:length(g.observables.largest_cluster_size)-1
    ])

    C₀ = g.observables.largest_cluster_size[t₀+1]
    C₁ = g.observables.largest_cluster_size[t₁+1]

    g.observables.delta = (t₀, t₁, C₀, C₁) ./ g.N
    return nothing
end

"""
    finalize_observables!(g::AbstractGraph)

Finalize the observables for the given evolved `g`

Arguments
* `g` : An evolved instance of type AbstractGraph
Returns
* `nothing`, updates `g` in-place
"""
function finalize_observables!(g::AbstractGraph)
    compute_delta!(g)
    return nothing
end


text:
function sample_interpolation(x,y,xmin,xmax,step)

  nx = length(x)
  n = round(Int64,(xmax-xmin)/step+1)
  xnew = Vector{Float64}(undef,n)
  ynew = Vector{Float64}(undef,n)

  for i in 1:n
   xnew[i] = xmin + (i-1)*step
   ix = findfirst( x -> x > xnew[i], x )
   if ix == 1
     ynew[i] = y[1]
   elseif ix == nothing
     ynew[i] = y[nx]
   else
     ynew[i] = y[ix-1] + (xnew[i]-x[ix-1])*((y[ix]-y[ix-1])/(x[ix]-x[ix-1]))
   end
  end

  return xnew, ynew

end


text:
#####################################################################
# IMPLEMENTACIÓN DEL ALGORITMOS DE INICIALIZACIÓN DE PESOS
# Basado en capítulo 7, algoritmo 6
#   CASO ENTRADA REAL de dimensión d > 1 SALIDA REAL de dimensión s>1
#####################################################################
export nn_from_data
"""
    nn_from_data(X_train,Y_train, n, M=10)
    Devuelve una red neuronal con los pesos ya inicializados
    de acorte a los conjuntos de entrenamiento.
    `n` es el número de neuronas en la capa oculta.
    El tamaño de entrada  y salida de la red neuronal vienen determinados por
    por los propios datos de entranamiento.
    El tamaño de entrada se corresponde con el número de atributos de X_train (su número de columnas)
    y a salida con el número de columnas de `Y_train`.

    M es una constante que depende de la función de activación,
    por lo visto en teoría M=10 funciona para todas las

"""
function nn_from_data(X_train::Matrix,Y_train::Matrix, n::Int, M=10)::AbstractOneLayerNeuralNetwork
    (_ , entry_dimension) = size(X_train)
    (_ , output_dimension) = size(Y_train)
    # inicializamos p
    p = rand(Float32, entry_dimension)

    index :: Int = 1
    tam :: Int = 0
    nodes = [Vector{Float64}(undef, output_dimension) for _ in 1:n]
    y_values = [Vector{Float64}(undef, output_dimension) for _ in 1:n]
    my_keys = zeros(Float64, n)

    while tam < n && index <= n
        new_point = X_train[index, :]
        if notOrtonormal(nodes, p, new_point, tam)
            tam += 1
            ordered_vector = sum(p.*new_point)
            my_keys[tam] =  ordered_vector
            nodes[tam] = new_point
            y_values[tam] = Y_train[index,:]
        end
        index += 1
    end
    ordered_values_index = sortperm(my_keys)
    # Matrices de la red neuronal
    # A = n x d
    # S = n x 1
    # B = s x n
    A = zeros(Float64, n, entry_dimension)
    S = zeros(Float64, n)
    B = zeros(Float64, output_dimension, n)

    # Cálculo del valor de las neuronas
    key =  ordered_values_index[1]
    x_a = nodes[key]
    y_a = y_values[key]
    # valores iniciales
    S[1]=M
    B[:, 1] = y_a

    for index in 2:n
        key =  ordered_values_index[index]
        x_s = nodes[key]
        y_s = y_values[key]

        coeff_aux = 2M / sum(p.* (x_s - x_a))
        S[index] = M -  coeff_aux*sum(p .* x_s)
        A[index,:] = coeff_aux * p
        B[:,index] = y_s - y_a

        x_a = x_s
        y_a = y_s
    end
    return FromMatrixNN(S,A,B)
end


text:
export IterationOptions

"""
Options for iterative solvers.
"""
@with_kw immutable IterationOptions
    max_iter::Int = 100                    # maximum number of iterations
    tol1::Float64 = √eps()                 # 1 norm for change
    tol∞::Float64 = √eps()                 # ∞ norm for change
end

"""
Check convergence.
"""
function _converged(Δ, options::IterationOptions)
    norm(Δ, 1) ≤ options.tol1 && norm(Δ, ∞) ≤ options.tol∞
end


text:
module Utils

export pwf

function pwf(s)#print with fieldnames
    for f in fieldnames(typeof(s))
        println("$f: $(getfield(s,f))")
    end
end

end


text:
```Write a program that prints all prime numbers. (Note: if your programming language
does not support arbitrary size numbers,
 printing all primes up to the largest number you can easily represent is fine too.)```

 for i in 1:10000
     if i % 2 == 1
        println(i)
     end
  end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# Memoisation for Approximate Computation
#
# #######################################

function ApproximateHashingMemoise(fn, hashmap, hashfn, inputs...)
    hashid = hashfn(fn, inputs...)
    if( haskey(hashmap, hashid))
       return hashmap[hashid]
    else
        res = fn(inputs...)
        hashmap[hashid] = res
        return res
    end
end

function TrendingMemoisation(fn, hashfn, hasharray, limit, inputs...)
    index = hashfn(fn, inputs...)

    # Training is complete
    if(hasharray[index][1] > limit)
        return hasharray[index][3]
    end

    # Continue Training
    res = fn(inputs...)
    hasharray[index][1] += 1
    hasharray[index][2] += res

    if(hasharray[index][1] > limit)
        hasharray[index][3] = hasharray[index][2]/Float64.(hasharray[index][1])
        return hasharray[index][3]
    end

    return res
end


text:
module TraitSimulation
using GLM # this already defines some useful distribution and link types
using DataFrames # so we can test it
using Tables
using LinearAlgebra
using Random
using SnpArrays
using OrdinalMultinomialModels
using VarianceComponentModels
using Distributions
using DelimitedFiles
import Base: show
using LinearAlgebra: BlasReal, copytri!

include("simulatematrixnormal.jl")

include("modelframework.jl")

include("simulatesnparray.jl")

include("modelparameterparsers.jl")

include("simulatepower.jl")

include("writetofilepheno.jl")

"Simulate the trait under the given model."
simulate(trait::AbstractTraitModel) = __default_behavior(trait)

"Simulate a trait `n` times independently."
simulate(trait::AbstractTraitModel, n::Integer) = __default_behavior(trait)

"Simulate a trait and store the result in y."
simulate!(y, trait::AbstractTraitModel) = __default_behavior(trait)

# default behavior for UnivariateDistribution from GLM.jl package
function __get_distribution(dist::Type{D}, μ) where D <: UnivariateDistribution
    return dist(μ)
end

# specific to Gamma
function __get_distribution(dist::Type{Gamma}, μ)
    β = 1 / μ # here β is the rate parameter for gamma distribution
    return dist(1, β) # α = 1
end

# specific to NegativeBinomial
function __get_distribution(dist::Type{NegativeBinomial}, μ)
    p = inv(1 + μ)
    return dist(1, p) # r = 1
end

"""
```
simulate(trait::GLMTrait, n::Integer)
```
This simulates a GLM trait once under the desired generalized linear model, specified using the GLMTrait type.
"""
  function simulate(trait::GLMTrait)
      # pre-allocate output
      y = Vector{eltype(trait.dist)}(undef, nsamplesize(trait))
      # do the simulation
      simulate!(y, trait)
      return y
  end

  """
  ```
  simulate(trait::GLMTrait, n::Integer)
  ```
  This simulates a GLM trait once times under the desired generalized linear model, specified using the GLMTrait type and overwrites the field y.
  """
  function simulate!(y, trait::GLMTrait)
      dist = trait.dist
      for i in eachindex(y)
          y[i] = rand(__get_distribution(dist, trait.μ[i]))
      end
      return y
  end

  """
  ```
  simulate(trait::GLMTrait, n::Integer)
  ```
  This simulates a GLM trait n times under the desired generalized linear model, specified using the GLMTrait type.
  """
  function simulate(trait::GLMTrait, n::Integer)
      # pre-allocate output
      Y_n = [zeros(size(trait.μ)) for _ in 1:n] # replicates
      # do the simulation n times, storing each result in a column
      for k in 1:n
          simulate!(Y_n[k], trait)
      end
      return Y_n
  end

  """
  ```
  simulate(OrderedMultinomialModel, n_reps; Logistic = false, threshold == empty)
  ```
    This simulates a OrderedMultinomialTrait trait n times under the desired model, specified using the OrderedMultinomialTrait type.
  By default we simulate the multinomial ordered outcome, but with the specification of the Logistic and threshold arguments, we can do the transformation to ordinal logistic.
  """
  function simulate(trait::OrderedMultinomialTrait; Logistic::Bool = false, threshold::Union{T, Nothing} = nothing) where T <: Real
     y = Vector{Int64}(undef, nsamplesize(trait)) # preallocate
     simulate!(y, trait; Logistic = Logistic, threshold = threshold) # do the simulation
     return y
  end

  """
  ```
  simulate!(y, OrderedMultinomialModel; Logistic = false, threshold == empty)
  ```
    This simulates a OrderedMultinomialTrait trait one time under the desired model and overwrites y, specified using the OrderedMultinomialTrait type.
  By default we simulate the multinomial ordered outcome, but with the specification of the Logistic and threshold arguments, we can do the transformation to ordinal logistic.
  """
  function simulate!(y, trait::OrderedMultinomialTrait; Logistic::Bool = false, threshold::Union{T, Nothing} = nothing) where T <: Real
      y .= rpolr(trait.X, trait.β, trait.θ, trait.link)
      if Logistic
          threshold == nothing && error("I need the cutoff for case/control")
          y .= Int64.(y .> threshold) # need threshold/ cutoff for case/control
      end
      return y
  end

  """
  ```
  simulate(OrderedMultinomialModel, n::Integer; Logistic = false, threshold == empty)
  ```
    This simulates a OrderedMultinomialTrait trait n times under the desired model, specified using the OrderedMultinomialTrait type.
  By default we simulate the multinomial ordered outcome, but with the specification of the Logistic and threshold arguments, we can do the transformation to ordinal logistic.
  """
  function simulate(trait::OrderedMultinomialTrait, n::Integer; Logistic::Bool = false, threshold::Union{T, Nothing} = nothing) where T <: Real
      # pre-allocate output
      Y_n = [zeros(nsamplesize(trait)) for _ in 1:n] # replicates
      # do the simulation n times, storing each result in a column
      for k in 1:n
          simulate!(Y_n[k], trait; Logistic = Logistic, threshold = threshold)
      end
      return Y_n
  end

  """
  ```
  simulate!(Y, trait::VCMTrait)
  ```
  This simulates a VCM trait one time under the desired variance component model, specified using the VCMTrait type.
  """
  function simulate(trait::VCMTrait)
     Y = zeros(size(trait.μ)) # preallocate
     simulate!(Y, trait) # do the simulation
     return Y
  end

  """
  ```
  simulate!(Y, trait::VCMTrait)
  ```
  This simulates a VCM trait one time under the desired variance component model, specified using the VCMTrait type and writes over the field Y.
  """
  function simulate!(Y, trait::VCMTrait)
      fill!(Y, 0.0)
      TraitSimulation.VCM_trait_simulation(Y, trait.Z, trait.μ, trait.vc)
      return Y
  end

  """
  ```
  simulate!(Y, Z, trait::VCMTrait)
  ```
  This simulates a VCM trait one time under the desired variance component model, specified using the VCMTrait type and writes over the fields Y and Z to simulate using threads.
  """
  function simulate!(Y, Z, trait::VCMTrait)
      fill!(Y, 0.0)
      #LinearAlgebra.BLAS.set_num_threads(1) # make simulations more consistent
      TraitSimulation.VCM_trait_simulation(Y, Z, trait.μ, trait.vc)
      return Y
  end

  """
  ```
  simulate(trait::VCMTrait, n::Integer)
  ```
  This simulates a VCM trait n times under the desired variance component model, specified using the VCMTrait type.
  """
  function simulate(trait::VCMTrait, n::Integer)
      # pre-allocate output
      Y_n = [zeros(size(trait.μ)) for _ in 1:n] # replicates
      # do the simulation n times, each thread needs its own Z to collect randome effects
      Z = [zeros(size(trait.μ))  for _ in 1:Threads.nthreads()]
      Threads.@threads for k in 1:n
          simulate!(Y_n[k], Z[Threads.threadid()], trait)
      end
      return Y_n
  end

  """
  ```
  simulate(trait::GLMMTrait)
  ```
  This simulates a GLMM trait one time under the desired generalized linear mixed model, specified using the GLMMTrait type.
  """
  function simulate(trait::GLMMTrait)
      # pre-allocate output
      Y = Matrix{eltype(trait.dist)}(undef, size(trait.μ))
      # do the simulation
      simulate!(Y, trait)
      return Y
  end

  """
  ```
  simulate!(Y, trait::GLMMTrait)
  ```
  This simulates a GLMM trait one time under the desired generalized linear mixed model, specified using the GLMMTrait type and overwrite the field Y.
  """
  function simulate!(Y, trait::GLMMTrait)
      #  simulate random effects
      fill!(trait.Z, 0.0)
      fill!(trait.Y_vcm, 0.0)
      TraitSimulation.VCM_trait_simulation(trait.Y_vcm, trait.Z, trait.η, trait.vc)
      # simulate from the glm with the mean μ and vector of ones for coefficients
      trait.μ .= GLM.linkinv.(trait.link, trait.Y_vcm)
      copyto!(Y, rand.(__get_distribution.(trait.dist, trait.μ)))
      Y
  end

  """
  ```
  simulate!(Y, trait::GLMMTrait)
  ```
  This simulates a GLMM trait one time under the desired generalized linear mixed model, specified using the GLMMTrait type and overwrite the field Y and Z using threads.
  """
  function simulate!(Y, Z, trait::GLMMTrait)
      #  simulate random effects
      fill!(Z, 0.0)
      fill!(trait.Y_vcm, 0.0)
      TraitSimulation.VCM_trait_simulation(trait.Y_vcm, Z, trait.η, trait.vc)
      # simulate from the glm with the mean μ and vector of ones for coefficients
      trait.μ .= GLM.linkinv.(trait.link, trait.Y_vcm)
      copyto!(Y, rand.(__get_distribution.(trait.dist, trait.μ)))
      Y
  end

  """
  ```
  simulate(trait::GLMMTrait, n::Integer)
  ```
  This simulates a trait n times under the desired generalized linear mixed model, specified using the GLMMTrait type.
  """
  function simulate(trait::GLMMTrait, n::Integer)
      # pre-allocate output
      Y_n = [zeros(size(trait.μ)) for _ in 1:n]
      Z = [zeros(size(trait.μ))  for _ in 1:Threads.nthreads()]
      # do the simulation n times, storing each result
      Threads.@threads for k in 1:n
         simulate!(Y_n[k], trait)
      end
      return Y_n
  end

  export mean_formula, VarianceComponent, @vc, vcobjtuple
  export GLMTrait, OrderedMultinomialTrait, VCMTrait, GLMMTrait
  export simulate_effect_size, snparray_simulation, genotype_sim
  export nsamplesize, neffects, noutcomecategories, nvc, ntraits
  export simulate!, simulate, writepheno
  export null_and_alternative_vcm_and_rotate, power_simulation, power, VCM_trait_simulation
end #module


text:
"""
    Exo2Binion{T <: Real} = Parabolic{Exo1Binion{T}}

A 2-exo-binion is a parabolic Cayley-Dickson construct with 1-exo-binions.
"""
const Exo2Binion{T <: Real} = Parabolic{Exo1Binion{T}}

function Exo2Binion(z::Exo1Binion{T}) where T <: Real
    Exo2Binion{T}(z, zero(Exo1Binion{T}))
end

function Exo2Binion(z::Binion{T}) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(z, zero(Binion{T})),
        zero(Exo1Binion{T})
    )
end

function Exo2Binion(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Exo1Binion{T}(
            Binion{T}(f, g),
            Binion{T}(h, j)
        )
    )
end

function Exo2Binion(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Exo1Binion{T}(
            Binion{T}(f, g),
            Binion{T}(h, zero(T))
        )
    )
end

function Exo2Binion(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Exo1Binion{T}(
            Binion{T}(f, g),
            zero(Binion{T})
        )
    )
end

function Exo2Binion(a::T, b::T, c::T, d::T, f::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        Exo1Binion{T}(
            Binion{T}(f, zero(T)),
            zero(Binion{T})
        )
    )
end

function Exo2Binion(a::T, b::T, c::T, d::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, d)
        ),
        zero(Exo1Binion{T})
    )
end

function Exo2Binion(a::T, b::T, c::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            Binion{T}(c, zero(T))
        ),
        zero(Exo1Binion{T})
    )
end

function Exo2Binion(a::T, b::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, b),
            zero(Binion{T})
        ),
        zero(Exo1Binion{T})
    )
end

function Exo2Binion(a::T) where T <: Real
    Exo2Binion{T}(
        Exo1Binion{T}(
            Binion{T}(a, zero(T)),
            zero(Binion{T})
        ),
        zero(Exo1Binion{T})
    )
end

function Exo2Binion(a::Real, b::Real)
    Exo2Binion(promote(a, b)...)
end

function Exo2Binion(a::Real, b::Real, c::Real)
    Exo2Binion(promote(a, b, c)...)
end

function Exo2Binion(a::Real, b::Real, c::Real, d::Real)
    Exo2Binion(promote(a, b, c, d)...)
end

function Exo2Binion(a::Real, b::Real, c::Real, d::Real, f::Real)
    Exo2Binion(promote(a, b, c, d, f)...)
end

function Exo2Binion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    Exo2Binion(promote(a, b, c, d, f, g)...)
end

function Exo2Binion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    Exo2Binion(promote(a, b, c, d, f, g, h)...)
end

function Exo2Binion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    Exo2Binion(promote(a, b, c, d, f, g, h, j)...)
end

function Exo2Binion(z::Exo2Real{T}) where T <: Real
    Exo2Binion(z.l.l, zero(T), z.l.r, zero(T), z.r.l, zero(T), z.r.r)
end

"""
    unreal(z::Exo2Binion)

The unreal part of a 2-exo-binion is a 7-vector.
"""
function unreal(z::Exo2Binion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo2Binion)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", i: ")
    print(io, z.l.l.r)
    print(io, ", a: ")
    print(io, z.l.r.l)
    print(io, ", ia: ")
    print(io, z.l.r.r)
    print(io, ", b: ")
    print(io, z.r.l.l)
    print(io, ", ib: ")
    print(io, z.r.l.r)
    print(io, ", ab: ")
    print(io, z.r.r.l)
    print(io, ", (ia)b: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{Exo2Binion{T}}) where T <: Real
    Exo2Binion{T}(random(Exo1Binion{T}), random(Exo1Binion{T}))
end


text:
module Algopt

include("UniDerivative.jl")
include("Bracketing.jl")
include("LocalDescent.jl")
include("FirstOrder.jl")
include("TestFunctions.jl")

end # module