text:
struct ROCKETOptions
  λ::Float64
  zeroThreshold::Float64
  refit::Bool
  optionsCD::CDOptions
  optionsScaledLasso::IterLassoOptions
end

ROCKETOptions(;
  λ::Float64=0.3,
  zeroThreshold::Float64=1e-4,
  refit::Bool=true,
  optionsCD::CDOptions=CDOptions(),
  optionsScaledLasso::IterLassoOptions=IterLassoOptions()) =
     ROCKETOptions(λ, zeroThreshold, refit, optionsCD, optionsScaledLasso)

function _stdColumn!{T<:AbstractFloat}(out::Vector{T}, X::AbstractMatrix{T})

  n, p = size(X)
  p == length(out) || throw(DimensionMismatch())

  for j=1:p
    t = zero(T)
    for i=1:n
      @inbounds t += X[i, j]^2.
    end
    @inbounds out[j] = sqrt(t / n)
  end

  out
end

# hOmega = inv(ɛ'ɛ/n)
# eP = hOmega[1,2]
# eVar = (hOmega[1,1]*hOmega[2,2] + hOmega[1,2]*hOmega[1,2]) / n
function _TEinv22{T<:AbstractFloat}(ɛ::Matrix{T})
  n, p = size(ɛ)
  p == 2 || throw(ArgumentError())

  a = zero(T)
  b = zero(T)
  c = zero(T)
  for i=1:n
    a += ɛ[i, 1] * ɛ[i, 1]
    c += ɛ[i, 2] * ɛ[i, 2]
    b += ɛ[i, 1] * ɛ[i, 2]
  end

  dt = a*c - b*b
  eP = -b / dt * n
  eVar = (a * c + b * b) / dt / dt * n
  eP, eVar
end

function _teInferenceGaussian{T<:AbstractFloat}(
  Y::Matrix{T}, a::Int, b::Int, S::Vector{Int})

  n, p = size(Y)
  x = view(Y, :, S)
  y = view(Y, :, [a,b])
  β = x \ y
  ɛ = x * β
  @. ɛ = y - ɛ
  _TEinv22(ɛ)
end

function _teInferenceGaussian{T<:AbstractFloat}(
  Y::Matrix{T}, a::Int, b::Int,
  γa::Union{SparseIterate{T}, Vector{T}},
  γb::Union{SparseIterate{T}, Vector{T}})

  n, p = size(Y)
  I = setdiff(1:p, [a,b])
  x = view(Y, :, I)
  y = view(Y, :, [a,b])

  ɛ = zeros(T, (n, 2))
  @inbounds for i=1:n
    ɛ[i, 1] = HD._row_A_mul_b(x, γa, i)
    ɛ[i, 2] = HD._row_A_mul_b(x, γb, i)
  end
  @. ɛ = y - ɛ
  _TEinv22(ɛ)
end

# Method Type:
#               1 (default)     uses lasso to estimate the neighborhood
#                               of node a and node b
#               2               uses sqrt-lasso to estimate the neighborhood
#                               of node a and node b
#               3               uses scalled-lasso to estimate the neighborhood
#                               of node a and node b
#               4               OLS without model selection
function _teInferenceGaussian{T<:AbstractFloat}(
  Y::Matrix{T}, a::Int, b::Int,
  methodType::Int,
  options::ROCKETOptions=ROCKETOptions())

  n, p = size(Y)
  I = setdiff(1:p, [a,b])

  λ = options.λ
  τ0 = options.zeroThreshold

  if methodType in [1, 2, 3]
    # lasso
    x = view(Y, :, I)
    ya = view(Y, :, a)
    yb = view(Y, :, b)

    scaleX = zeros(p-2)
    _stdColumn!(scaleX, x)

    g = AProxL1(λ, scaleX)
    γa = SparseIterate(p-2)
    γb = SparseIterate(p-2)
    if methodType == 1
      fa = CDLeastSquaresLoss(ya, x)
      fb = CDLeastSquaresLoss(yb, x)
      coordinateDescent!(γa, fa, g, options.optionsCD)
      coordinateDescent!(γb, fb, g, options.optionsCD)
    elseif methodType == 2
      fa = CDSqrtLassoLoss(ya, x)
      fb = CDSqrtLassoLoss(yb, x)
      coordinateDescent!(γa, fa, g, options.optionsCD)
      coordinateDescent!(γb, fb, g, options.optionsCD)
    else # if methodType == 3
      scaledLasso!(γa, x, ya, λ, scaleX, options.optionsScaledLasso)
      scaledLasso!(γb, x, yb, λ, scaleX, options.optionsScaledLasso)
    end

    if options.refit
      Sa = find(x -> abs(x) > τ0, γa)
      Sb = find(x -> abs(x) > τ0, γb)
      S = union(I[Sa], I[Sb])

      return _teInferenceGaussian(Y, a, b, S)
    else
      return _teInferenceGaussian(Y, a, b, γa, γb)
    end
  elseif methodType == 4
    return _teInferenceGaussian(Y, a, b, I)
  else
    throw(ArgumentError("Unknown methodType"))
  end
end

"""
teInference - Inference for an edge in an elliptical coopula model.

Arguments:
              Y               (n x p) data matrix
              a, b            index of an edge for which inference is made
              methodType      number from 1-3 indicating the procedure

Method Type:
              1 (default)     uses lasso to estimate the neighborhood
                              of node a and node b
              2               uses sqrt-lasso to estimate the neighborhood
                              of node a and node b
              3               uses scalled-lasso to estimate the neighborhood
                              of node a and node b
              4               OLS without model selection

Covariance Type:
              1 (default)     rank correlation
              2               pearson correlation
              3               pearson correlation after marginal
                              transformation
              4               empirical covariance X'X/n

Returns:
  (eP, eVar, ISa, ISb)

  eP   - estimated edge parameter
  eVar - estimated variance of the parameter
  ISa  - neighborhood of a
  ISb  - neighborhood of b
"""
function teInference{T<:AbstractFloat}(
  Y::Matrix{T}, a::Int, b::Int,
  methodType::Int=1,
  covarianceType::Int=1,
  options::ROCKETOptions=ROCKETOptions())

  a == b && throw(ArgumentError("a should not be equal to b"))
  n, p = size(Y)
  1 <= a <= p || throw(ArgumentError("1 <= a <= p"))
  1 <= b <= p || throw(ArgumentError("1 <= b <= p"))

  if covarianceType != 4
    throw(ArgumentError("unknown covarianceType"))
  end
  return _teInferenceGaussian(Y, a, b, methodType, options)

  # # compute covariance / correlation
  # if covarianceType == 1
  #   covM = kendallsTau(Y)
  #   throw(ArgumentError("unknown covarianceType"))
  # elseif covarianceType == 2
  #   covM = cor(Y)
  #   throw(ArgumentError("unknown covarianceType"))
  # elseif covarianceType == 3
  #   covM = nonparanormalCorrelation(Y)
  #   throw(ArgumentError("unknown covarianceType"))
  # elseif covarianceType == 4
  #   return _teInferenceGaussian(Y, a, b, methodType, options)
  # else
  #   throw(ArgumentError("unknown covarianceType"))
  # end
end


text:
# ps3000_flash_led.jl

export ps3000_flash_led

function ps3000_flash_led(handle::Int16)
  ps_status = ccall((:ps3000_flash_led, ps3000driver),
                  UInt16,
                  (Int16,),
                  handle)
  if ps_status == 0
    error("ps3000 Error: Handle Not Valid")
  end
  return nothing
end

ext: jl


text:
# The content of this file is copied from https://github.com/JuliaStats/StatsModels.jl/blob/master/test/extension.jl (in December 2021)
# The license hereafter is associated with this content

# Copyright (c) 2016: Dave Kleinschmidt.
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

poly(x, n) = x^n

abstract type PolyModel end
struct PolyTerm <: AbstractTerm
    term::Symbol
    deg::Int
end
PolyTerm(t::Term, deg::ConstantTerm) = PolyTerm(t.sym, deg.n)

StatsModels.apply_schema(t::FunctionTerm{typeof(poly)}, sch, ::Type{<:PolyModel}) =
    PolyTerm(t.args_parsed...)

StatsModels.modelcols(p::PolyTerm, d::NamedTuple) =
    reduce(hcat, [d[p.term].^n for n in 1:p.deg])


text:
#!/usr/bin/env ~/julia-1.6.1/bin/julia
println("hello world")


text:
import Flux.testmode!

abstract type AbstractEnv end

include("vis/utils.jl")

IntOrNothing  = Union{Int,  Nothing}
RealOrNothing = Union{Real, Nothing}

mutable struct EnvWrapper
    done::Bool
    total_reward::RealOrNothing
    steps::Int
    train::Bool
	reward_threshold::RealOrNothing
	max_episode_steps::IntOrNothing
    _env::AbstractEnv
    _ctx::AbstractCtx
end

EnvWrapper(env::AbstractEnv, ctx::AbstractCtx, train::Bool=true;
		   reward_threshold=nothing, max_episode_steps=nothing) =
EnvWrapper(false, 0, 0, train, reward_threshold, max_episode_steps, env, ctx)

function step!(env::EnvWrapper, a)
    s′, r, done, dict = step!(env._env, a)
    env.total_reward = env.total_reward .+ r
    env.steps += 1
    env.done = done
	if !isnothing(env.max_episode_steps)
		env.done |= env.steps ≥ env.max_episode_steps
	end
    return s′, r, done, dict
end

function reset!(env::EnvWrapper)
    env.done = false
    env.total_reward = 0
    env.steps = 0
    reset!(env._env)
end

render!(env::EnvWrapper, ctx::AbstractCtx) = render!(env._env, ctx)
render!(env::EnvWrapper) = render!(env, env._ctx)

Ctx(env::EnvWrapper, mode::Symbol = :human_window) = Ctx(env._env, mode)

"""
Returns the observational state of the environment. The original state can
be accessed by `\`env._env.state\``.
"""
function state(env::EnvWrapper)
	try
		return _get_obs(env._env)
	catch y
		if isa(y, UndefVarError) || isa(y, MethodError)
			return env._env.state
		end
	end
end

function testmode!(env::EnvWrapper, val::Bool=true)
    env.train = !val
end

trainable(env::EnvWrapper) = env.train
game_over(env::EnvWrapper) = env.done


text:
function is_binary_operator(sym::Symbol)
    return  is_plus_operator(sym)    ||
            is_minus_operator(sym)   ||
            is_divide_operator(sym)  ||
            is_product_operator(sym) ||
            is_power_operator(sym)
end
function is_plus_operator(sym::Symbol)
    return sym == :+
end
function is_minus_operator(sym::Symbol)
    return sym == :-
end
function is_divide_operator(sym::Symbol)
    return sym == :/
end
function is_product_operator(sym::Symbol)
    return sym == :*
end
function is_power_operator(sym::Symbol)
    return sym == :^
end
function is_unary_operator(sym::Symbol)
    return  is_sqrt_operator(sym) ||
            is_cbrt_operator(sym) ||
            is_exp_operator(sym) ||
            is_log_operator(sym) ||
            is_log2_operator(sym) ||
            is_log10_operator(sym) ||
            is_sin_operator(sym) ||
            is_cos_operator(sym) ||
            is_tan_operator(sym) ||
            is_abs_operator(sym) ||
            is_sign_operator(sym)
end
function is_sqrt_operator(sym::Symbol)
    return sym == :sqrt
end
function is_cbrt_operator(sym::Symbol)
    return sym == :cbrt
end
function is_exp_operator(sym::Symbol)
    return sym == :exp
end
function is_log_operator(sym::Symbol)
    return sym == :log
end
function is_log2_operator(sym::Symbol)
    return sym == :log2
end
function is_log10_operator(sym::Symbol)
    return sym == :log10
end
function is_sin_operator(sym::Symbol)
    return sym == :sin
end
function is_cos_operator(sym::Symbol)
    return sym == :cos
end
function is_tan_operator(sym::Symbol)
    return sym == :tan
end
function is_abs_operator(sym::Symbol)
    return sym == :abs
end
function is_sign_operator(sym::Symbol)
    return sym == :sign
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"A kernel is a function of two variables."
abstract type Kernel end

"A boundary kernel is associated with a parameterized boundary."
abstract type BoundaryKernel <: Kernel end

# A concrete kernel k should implement the routine k(x,y,param).
# The boundary kernel may not know its own parameterization, but
# the parameterization is passed as an extra argument.

# By default, we add a parameterization by asking the kernel
(kernel::BoundaryKernel)(x, y) = kernel(x, y, parameterization(kernel))
# If two numbers are given, we assume they are in the parameter domain
(kernel::BoundaryKernel)(t::Number, tau::Number, param) =
    kernel(t, tau, param, applymap(param, t), applymap(param, tau))
# the first argument is a field point
(kernel::BoundaryKernel)(x::SVector, tau::Number, param) =
    kernel(x, tau, param, applymap(param, tau))


text:
module EMLtank
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("tank/VesselVolume.jl")
	include("tank/TankVL.jl")
	include("tank/TankL.jl")
	include("tank/SumpTank.jl")
	include("tank/tank_simplified.jl")
end


text:
function evaluate_decision(structure::VerticalStructure, decision::AbstractVector)
    # Evalaute decision stage-wise
    cᵀx = _eval_first_stage(structure, decision)
    𝔼Q = _eval_second_stages(structure, decision, objective_sense(structure.first_stage))
    # Return evaluation result
    return cᵀx + 𝔼Q
end

function statistically_evaluate_decision(structure::VerticalStructure, decision::AbstractVector)
    # Evalaute decision stage-wise
    cᵀx = _eval_first_stage(structure, decision)
    𝔼Q, σ² = _stat_eval_second_stages(structure, decision, objective_sense(structure.first_stage))
    return cᵀx + 𝔼Q, sqrt(σ²)
end

function _eval_first_stage(structure::VerticalStructure, decision::AbstractVector)
    # Update decisions (checks handled by first-stage model)
    take_decisions!(structure.first_stage,
                    all_decision_variables(structure.first_stage, 1),
                    decision)
    # Optimize first_stage model
    optimize!(structure.first_stage)
    # Switch on return status
    status = termination_status(structure.first_stage)
    result = if status in AcceptableTermination
        result = objective_value(structure.first_stage)
    else
        result = if status == MOI.INFEASIBLE
            result = objective_sense(structure.first_stage) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            result = objective_sense(structure.first_stage) == MOI.MAX_SENSE ? Inf : -Inf
        else
            error("First-stage model could not be solved, returned status: $status")
        end
    end
    # Revert back to untaken decisions
    untake_decisions!(structure.first_stage, all_decision_variables(structure.first_stage, 1))
    # Return evaluation result
    return result
end

function _eval_second_stages(structure::VerticalStructure{2,1,Tuple{SP}},
                             decision::AbstractVector,
                             sense::MOI.OptimizationSense) where SP <: ScenarioProblems
    update_known_decisions!(structure.decisions[2], decision)
    map(subprob -> update_known_decisions!(subprob), subproblems(structure))
    return outcome_mean(subproblems(structure), probability.(scenarios(structure)), sense)
end
function _eval_second_stages(structure::VerticalStructure{2,1,Tuple{SP}},
                             decision::AbstractVector,
                             sense::MOI.OptimizationSense) where SP <: DistributedScenarioProblems
    Qs = Vector{Float64}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async Qs[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                decision,
                sense) do sp, d, x, sense
                    scenarioproblems = fetch(sp)
                    decisions = fetch(d)
                    num_scenarios(scenarioproblems) == 0 && return 0.0
                    update_known_decisions!(decisions, x)
                    map(subprob -> update_known_decisions!(subprob), subproblems(scenarioproblems))
                    return outcome_mean(subproblems(scenarioproblems),
                                        probability.(scenarios(scenarioproblems)),
                                        sense)
                end
        end
    end
    return sum(Qs)
end

function _stat_eval_second_stages(structure::VerticalStructure{2,1,Tuple{SP}},
                                  decision::AbstractVector,
                                  sense::MOI.OptimizationSense) where SP <: ScenarioProblems
    update_known_decisions!(structure.decisions[2], decision)
    map(subprob -> update_known_decisions!(subprob), subproblems(structure))
    return welford(subproblems(structure), probability.(scenarios(structure)), sense)
end
function _stat_eval_second_stages(structure::VerticalStructure{2,1,Tuple{SP}},
                                  decision::AbstractVector,
                                  sense::MOI.OptimizationSense) where SP <: DistributedScenarioProblems
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                decision,
                sense) do sp, d, x, sense
                    scenarioproblems = fetch(sp)
                    decisions = fetch(d)
                    num_scenarios(scenarioproblems) == 0 && return zero(eltype(x)), zero(eltype(x)), zero(eltype(x)), zero(Int)
                    update_known_decisions!(decisions, x)
                    map(subprob -> update_known_decisions!(subprob), subproblems(scenarioproblems))
                    return welford(subproblems(scenarioproblems), probability.(scenarios(scenarioproblems)), sense)
                end
        end
    end
    𝔼Q, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼Q, σ²
end


text:
@testset "Aqua.jl" begin
    Aqua.test_all(Octavian)
end


text:
using PyCall
using Conda
easy_get = x->try; strip(read(pipeline(`which $x`), String)); catch; "";end

tf = pyimport("tensorflow")
lib = readdir(splitdir(tf.__file__)[1])
TF_LIB_FILE = joinpath(splitdir(tf.__file__)[1],lib[findall(occursin.("libtensorflow_framework", lib))[1]])
TF_INC = tf.sysconfig.get_compile_flags()[1][3:end]
TF_ABI = tf.sysconfig.get_compile_flags()[2][end:end]

CC = easy_get("gcc")
CXX = easy_get("g++")
CMAKE = easy_get("cmake")
MAKE = easy_get("make")
GIT = "LibGit2"
PYTHON = easy_get("python")

s = ""
t = []
function adding(k, v)
    global s
    s *= "$k = \"$v\"\n"
    push!(t, "$k")
end
adding("BINDIR", Conda.BINDIR)
adding("LIBDIR", Conda.LIBDIR)
adding("TF_INC", TF_INC)
adding("TF_ABI", TF_ABI)
adding("EIGEN_INC", joinpath(Conda.LIBDIR,"Libraries"))
adding("CC", CC)
adding("CXX", CXX)
adding("CMAKE", CMAKE)
adding("MAKE", MAKE)
adding("GIT", GIT)
adding("PYTHON", PYTHON)
adding("TF_LIB_FILE", TF_LIB_FILE)

t = "join(["*join(t, ",")*"], \";\")"
s *= "__STR__ = $t"
open("deps.jl", "w") do io
    write(io, s)
end


text:
# This file was generated, do not modify it. # hide
ŷ = ẑ |> inverse_transform(z)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
module OneDimensionalCosineAdvectionDiffusion

using Printf
using Statistics

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Advection
using Oceananigans.Fields: interior

using ConvergenceTests: compute_error

# Advection and diffusion of a cosine.
c(x, y, z, t, U, κ) = exp(-κ * t) * cos(x - U * t)

function run_test(; Nx, Δt, stop_iteration, U = 1, κ = 1e-4,
                  architecture = CPU(), topo = (Periodic, Periodic, Periodic), advection = CenteredSecondOrder())

    #####
    ##### Test advection-diffusion in the x-direction
    #####

    domain = (x=(0, 2π), y=(0, 1), z=(0, 1))
    grid = RectilinearGrid(architecture, topology=topo, size=(Nx, 1, 1), halo=(3, 3, 3); domain...)

    model = NonhydrostaticModel(
                                 timestepper = :RungeKutta3,
                                        grid = grid,
                                   advection = advection,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     closure = IsotropicDiffusivity(ν=κ, κ=κ))

    set!(model, u = U,
                v = (x, y, z) -> c(x, y, z, 0, U, κ),
                w = (x, y, z) -> c(x, y, z, 0, U, κ),
                c = (x, y, z) -> c(x, y, z, 0, U, κ))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, iteration_interval=stop_iteration)

    @info "Running cosine advection diffusion test for vx, wx, and cx with Nz = $Nx and Δt = $Δt ($(typeof(advection)))..."
    run!(simulation)

    x = xnodes(model.tracers.c)
    c_analytical = c.(x, 0, 0, model.clock.time, U, κ)

    # Calculate errors
    vx_simulation = interior(model.velocities.v)[:, 1, 1] |> Array
    vx_errors = compute_error(vx_simulation, c_analytical)

    wx_simulation = interior(model.velocities.w)[:, 1, 1] |> Array
    wx_errors = compute_error(wx_simulation, c_analytical)

    cx_simulation = interior(model.tracers.c)[:, 1, 1] |> Array
    cx_errors = compute_error(cx_simulation, c_analytical)

    #####
    ##### Test advection-diffusion in the y-direction
    #####

    ydomain = (x=(0, 1), y=(0, 2π), z=(0, 1))
    ygrid = RectilinearGrid(topology=topo, size=(1, Nx, 1), halo=(3, 3, 3); ydomain...)

    model = NonhydrostaticModel(architecture = architecture,
                                 timestepper = :RungeKutta3,
                                        grid = ygrid,
                                   advection = advection,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     closure = IsotropicDiffusivity(ν=κ, κ=κ))

    set!(model, v = U,
                u = (x, y, z) -> c(y, x, z, 0, U, κ),
                w = (x, y, z) -> c(y, x, z, 0, U, κ),
                c = (x, y, z) -> c(y, x, z, 0, U, κ))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, iteration_interval=stop_iteration)

    @info "Running cosine advection diffusion test for uy, wy, and cy with Ny = $Nx and Δt = $Δt ($(typeof(advection)))..."
    run!(simulation)

    # Calculate errors
    uy_simulation = interior(model.velocities.u)[1, :, 1] |> Array
    uy_errors = compute_error(uy_simulation, c_analytical)

    wy_simulation = interior(model.velocities.w)[1, :, 1] |> Array
    wy_errors = compute_error(wy_simulation, c_analytical)

    cy_simulation = interior(model.tracers.c)[1, :, 1] |> Array
    cy_errors = compute_error(cy_simulation, c_analytical)

    #####
    ##### Test advection-diffusion in the z-direction
    #####

    zdomain = (x=(0, 1), y=(0, 1), z=(0, 2π))
    zgrid = RectilinearGrid(topology=topo, size=(1, 1, Nx), halo=(3, 3, 3); zdomain...)

    model = NonhydrostaticModel(architecture = architecture,
                                 timestepper = :RungeKutta3,
                                        grid = zgrid,
                                   advection = advection,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     closure = IsotropicDiffusivity(ν=κ, κ=κ))

    set!(model, w = U,
                u = (x, y, z) -> c(z, y, x, 0, U, κ),
                v = (x, y, z) -> c(z, y, x, 0, U, κ),
                c = (x, y, z) -> c(z, y, x, 0, U, κ))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, iteration_interval=stop_iteration)

    @info "Running cosine advection diffusion test for uz, vz, and cz with Nz = $Nx and Δt = $Δt ($(typeof(advection)))..."
    run!(simulation)

    # Calculate errors
    uz_simulation = interior(model.velocities.u)[1, 1, :] |> Array
    uz_errors = compute_error(uz_simulation, c_analytical)

    vz_simulation = interior(model.velocities.v)[1, 1, :] |> Array
    vz_errors = compute_error(vz_simulation, c_analytical)

    cz_simulation = interior(model.tracers.c)[1, 1, :] |> Array
    cz_errors = compute_error(cz_simulation, c_analytical)

    return (

            cx = (simulation = cx_simulation,
                  analytical = c_analytical,
                          L₁ = cx_errors.L₁,
                          L∞ = cx_errors.L∞),

            cy = (simulation = cy_simulation,
                  analytical = c_analytical,
                          L₁ = cy_errors.L₁,
                          L∞ = cy_errors.L∞),

            cz = (simulation = cz_simulation,
                  analytical = c_analytical,
                          L₁ = cz_errors.L₁,
                          L∞ = cz_errors.L∞),

            uy = (simulation = uy_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = uy_errors.L₁,
                          L∞ = uy_errors.L∞),

            uz = (simulation = uz_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = uz_errors.L₁,
                          L∞ = uz_errors.L∞),

            vx = (simulation = vx_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = vx_errors.L₁,
                          L∞ = vx_errors.L∞),

            vz = (simulation = vz_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = vz_errors.L₁,
                          L∞ = vz_errors.L∞),

            wx = (simulation = wx_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = wx_errors.L₁,
                          L∞ = wx_errors.L∞),

            wy = (simulation = wy_simulation,
                  analytical = c_analytical, # same solution as c.
                          L₁ = wy_errors.L₁,
                          L∞ = wy_errors.L∞),

            grid = grid

            )
end

end # module


text:
using HttpServer
using WebSockets

function coroutine(ws)
    while isopen(ws)
        data, = readguarded(ws)
        s = String(data)
        s == "" && break
        println("Received: ", s)
        if s[1] == "P"
            writeguarded(ws, "No, I'm not!")
        else
            writeguarded(ws, "Why?")
        end
    end
end

function gatekeeper(req, ws)
    println("\nOrigin:", origin(req), "    Target:", target(req), "    subprotocol:", subprotocol(req))
    # Non-browser clients don't send Origin. We liberally accept in this case.
    if origin(req) == "" || origin(req) == "http://127.0.0.1:8080" || origin(req) == "http://localhost:8080"
        coroutine(ws)
    else
        println("Inacceptable request")
    end
end

handle(req, res) = Response(200)

server = Server(HttpHandler(handle),
                WebSocketHandler(gatekeeper))

@async run(server, 8080)


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for sampling, sketching or selecting from linear
## systems of equations.
##
## Contents
## - Abstract Types
## - `sample` function documentation
## - Vector Row Sampler/Sketch/Selector
## - Vector Column Sampler/Sketch/Selector
## - Block Row Sampler/Sketch/Selector
## - Block Column Sampler/Sketch/Selector
## - Export Statements
##
############################################################################################

# Dependencies: Random, Distributions, LinearAlgebra

#############################################
# Abstract Types
#############################################

"""
    LinSysSampler

Abstract supertype for sampling, sketching or deterministically selecting components
    of a linear system.

# Aliases
- `LinSysSketch`
- `LinSysSelect`
"""
abstract type LinSysSampler end
LinSysSketch = LinSysSampler
LinSysSelect = LinSysSampler

"""
    LinSysVecRowSampler <: LinSysSampler

Abstract supertype for sampling, sketching or deterministically selecting a row space
    element from a linear system.

# Aliases
- `LinSysVecRowSketch`
- `LinSysVecRowSelect`
"""
abstract type LinSysVecRowSampler <: LinSysSampler end
LinSysVecRowSketch = LinSysVecRowSampler
LinSysVecRowSelect = LinSysVecRowSampler

"""
    LinSysVecColSampler <: LinSysSampler

Abstract supertype for sampling, sketching or deterministically selecting a column space
    element from a linear system.

# Aliases
- `LinSysVecColSketch`
- `LinSysVecColSelect`
"""
abstract type LinSysVecColSampler <: LinSysSampler end
LinSysVecColSketch = LinSysVecColSampler
LinSysVecColSelect = LinSysVecColSampler

"""
    LinSysBlkRowSampler <: LinSysSampler

Abstract supertype for sampling, sketching or deterministically selecting a collection of
    row space elements from a linear system.

# Aliases
- `LinSysBlkRowSketch`
- `LinSysBlkRowSelect`
"""
abstract type LinSysBlkRowSampler <: LinSysSampler end
LinSysBlkRowSketch = LinSysBlkRowSampler
LinSysBlkRowSelect = LinSysBlkRowSampler

"""
    LinSysBlkColSampler <: LinSysSampler

Abstract supertype for sampling, sketching or deterministically selecting a collection of
    column space elements from a linear system.

# Aliases
- `LinSysBlkColSketch`
- `LinSysBlkColSelect`
"""
abstract type LinSysBlkColSampler <: LinSysSampler end
LinSysBlkColSketch = LinSysBlkColSampler
LinSysBlkColSelect = LinSysBlkColSampler

#############################################
# `sample` Function Documentation
#############################################
"""
    sample(type::T where T<:LinSysSampler,
        A::AbstractArray,
        b::AbstractVector,
        x::AbstractVector,
        iter::Int64)

A common interface for specifying different strategies for sampling, selecting or sketching
    a linear system specified by `A` and `b`. The `type` argument is used to select the an
    appropriately defined strategy. The argument `x` is the current iterate value for the
    solution. The arguent `iter` is the iteration counter.

The value(s) returned by sample depend on the subtype of `LinSysSampler` being used.
    Specifically,
- For `T<:LinSysVecRowSampler`, a vector in the row space of `A` and constant are returned
- For `T<:LinSysVecColSampler`, a vector of `length(x)`, the matrix `A`,
    and a scalar-valued residual are returned.
"""
function sample(
    type::Nothing,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)
    return nothing
end

#############################################
# Vector Row Sampler/Sketch/Selector
#############################################
# Non-adaptive Cyclic Methods
include("linear_samplers/vec_row_deterministic_cyclic.jl")
include("linear_samplers/vec_row_one_rand_cyclic.jl")
include("linear_samplers/vec_row_hop_rand_cyclic.jl")
include("linear_samplers/vec_row_rand_cyclic.jl")

# Non-adaptive Sampling (with replacement)
include("linear_samplers/vec_row_unid_replace.jl")
include("linear_samplers/vec_row_prop_to_norm_replace.jl")

# Non-adaptive Sketching
#include("linear_samplers/vec_row_uniform.jl")
#include("linear_samplers/vec_row_uniform_sym.jl")
#include("linear_samplers/vec_row_gauss.jl")
#include("linear_samplers/vec_row_uniform_sparse.jl")
#include("linear_samplers/vec_row_uniform_sym_sparse.jl")
#include("linear_samplers/vec_row_gauss_sparse.jl")

# Adaptive Deterministic Selection
#include("linear_samplers/vec_row_max_residual.jl")
#include("linear_samplers/vec_row_max_distance.jl")

# Adaptive Cyclic Methods
#include("linear_samplers/vec_row_residual_cycle.jl")
#include("linear_samplers/vec_row_distance_cycle.jl")

#############################################
# Vector Column Sampler/Sketch/Selector
#############################################
# Non-adaptive Cyclic Methods
include("linear_samplers/vec_col_deterministic_cyclic.jl")

# Non-adaptive Sampling (with replacement)
#Leventhal-Lewis (non-symmetric)

#############################################
# Block Row Sampler/Sketch/Selector
#############################################

#############################################
# Block Column Sampler/Sketch/Selector
#############################################

#############################################
# Compositional Sampler/Sketch/Selector
#############################################
# Wishlist

#############################################
# Export Statements
#############################################
# See RLinearAlgebra.jl
#export LinSysSampler, LinSysSketch, LinSysSelect
#export LinSysVecRowSampler, LinSysVecRowSketch, LinSysVecRowSelect
#export LinSysVecColSampler, LinSysVecColSketch, LinSysVecColSelect
#export LinSysBlkRowSampler, LinSysBlkRowSketch, LinSysBlkRowSelect
#export LinSysBlkColSampler, LinSysBlkColSketch, LinSysBlkColSelect


text:
### Load the required packages ###

using DimensionalData, GeoData, ArchGDAL, NCDatasets, HDF5
using CSV, DataFrames, DataFramesMeta, TableView
using DynamicGrids, Dispersal
using Statistics, Dates, Unitful, Setfield, Flatten, FieldMetadata
using Plots, ColorSchemes, Colors
using DimensionalData: setdims, rebuild, Between

const DG = DynamicGrids

### Override default parameter range bounds ###

# For both interfaces and optimiser

import FieldMetadata: @bounds, bounds, @flattenable, flattenable
using Flatten: flatten

@bounds HumanDispersal begin
    human_exponent  | (0.0, 3.0)
    dist_exponent   | (0.0, 3.0)
    dispersalperpop | (0.0, 1e-8)
    max_dispersers  | (1e1, 1e5)
end
# Constant growth
@bounds ExactLogisticGrowth begin
    intrinsicrate | (0.0, 10.0)
end
# Alee
@bounds AlleeExtinction begin
    minfounders | (2e0, 1e3)
end

# Dont parametrise carrycap
@flattenable ExactLogisticGrowthMap begin
    carrycap | false
end
@flattenable ExactLogisticGrowth begin
    carrycap | false
end
# Don't parametrise λ for main models, it will just be the top bound
# so we need to estimate it. This is reverted below for local-dispersal
# only models.
@flattenable ExponentialKernel begin
    λ | false
end

### Define simulation settings ###

timestep = Month(1)

us_tspan = DateTime(2008, 5):timestep:DateTime(2013, 12)
eu_tspan = DateTime(2008, 5):timestep:DateTime(2016, 12)

aus = Lon(Between(113.3402, 153.9523)), Lat(Between(-43.62234, -10.65125))
us = Lon(Between(-130.332, -66.9)), Lat(Between(18.977, 71.59))
eu = Lon(Between(-13.7, 40.15)), Lat(Between(34.9, 70.30))

aus_incursionpoints = (
    Melbourne=(-37.805896, 144.959527),
    Mildura=(-34.219504, 142.130864),
    Coffs_Harbour=(-30.287245, 153.092991),
    Sydney=(-33.839943, 151.006101),
    Adelaide=(-34.901608, 138.601547),
    Port_Augusta=(-32.466201, 137.813850),
    Devonport=(-41.180545, 146.314887),
    Hobart=(-42.881742, 147.323879),
    Brisbane=(-27.436190, 152.990588),
    Cairns=(-16.937281, 145.747709),
    Perth=(-31.9505, 115.8605),
    Geraldton=(-28.778138, 114.615632)
)

us_incursionpoints = (
    San_Francisco=(36.9453, -122.0695),
)

obs = CSV.File(joinpath(basedir, "data/Oersted_occurrence.csv")) |> DataFrame
eu_obs = @where(obs, :Year .!== missing, :Year .== 2008, :Country .!= "USA")
eu_incursionpoints = collect(zip(eu_obs.Latitude, eu_obs.Longitude))

incursionpointstatekeys = (
    Melbourne = (:Vic),
    Mildura	 = (:Vic),
    Coffs_Harbour =	(:NSW),
    Sydney = (:NSW),
    Adelaide = (:SA),
    Port_Augusta = (:SA),
    Devonport  = (:Tas),
    Hobart = (:Tas),
    Brisbane = (:QLD),
    Cairns = (:QLD),
    Perth =	(:WA),
    Geraldton =	(:WA)
)

### Load numbered states rasters ###

loadstates(path, selectors) = GDALarray(joinpath(basedir, path);
    name=:States, usercrs=EPSG(4326))[Band(1), selectors...] |>
    x->permutedims(x, (Lat, Lon)) |> x -> replace_missing(x, 0) |>
    x->trunc.(Int, x)
aus_states = loadstates(joinpath(basedir, "data/aus_states.tif"), aus)
us_states = loadstates(joinpath(basedir, "data/US_states.tif"), us)
eu_states = loadstates(joinpath(basedir, "data/EU_states.tif"), eu)
us_states = loadstates(joinpath(basedir, "data/US_states.tif"), ())
# aus_states |> plot
# eu_states |> plot
# us_states |> plot

### Define the RuleSets ###

growthratesfilepath = joinpath(basedir, "output/sensitivity/growthrates_standard.ncd")
isfile(growthratesfilepath)
growthrates = NCDarray(growthratesfilepath; crs=crs(aus_states), dimcrs=EPSG(4326)) |> GeoArray |>
    x -> setdims(x, (@set dims(x, Ti).mode.span = Regular(Month(1)))) |>
    x -> setdims(x, convertmode(Projected, dims(x, Lon))) |>
    x -> setdims(x, convertmode(Projected, dims(x, Lat))) |>
    x -> permutedims(x, (Lat, Lon, Ti)) |>
    x -> reverse(x; dims=Lat)

# Make growth rates aux grids, replacing the `missing` values with zero.
# Dispersal.jl `init` can't contain `missing` - it will spread everywhere.

grzm = floattype.(replace_missing(growthrates, 0.0))
aus_growthrates_zeromissing = grzm[aus...]
us_growthrates_zeromissing = grzm[us...]
eu_growthrates_zeromissing = grzm[eu...]
aus_growthrates_zeromissing[Ti(1)] |> plot
us_growthrates_zeromissing[Ti(1)] |> plot
eu_growthrates_zeromissing[Ti(1)] |> plot
aus_aux = (growthrates=aus_growthrates_zeromissing,)
us_aux = (growthrates=us_growthrates_zeromissing,)
eu_aux = (growthrates=eu_growthrates_zeromissing,)

### Define masking layers ###

# These define regions that are not run, specifically, the sea.

aus_boolmask = GeoData.boolmask(aus_states)
aus_missingmask = GeoData.missingmask(aus_states)
aus_missingmask = GeoData.missingmask(aus_states)
us_boolmask = GeoData.boolmask(growthrates[us..., Ti(1)])
us_missingmask = GeoData.missingmask(growthrates[us..., Ti(1)])
eu_boolmask = GeoData.boolmask(eu_states)
eu_missingmask = GeoData.missingmask(eu_states)
# plot(aus_boolmask)
# plot(aus_missingmask)
# plot(us_boolmask)
# plot(us_missingmask)
# plot(eu_boolmask)
# plot(eu_missingmask)

# Growth masks. These are the cells where the mean growth rate is above zero.

aus_growthmask = rebuild(aus_boolmask, (mean(aus_growthrates_zeromissing; dims=Ti)[Ti(1)] .* aus_boolmask) .> 0)
us_growthmask = rebuild(us_boolmask, (mean(us_growthrates_zeromissing; dims=Ti)[Ti(1)] .* us_boolmask) .> 0)
eu_growthmask = rebuild(eu_boolmask, (mean(eu_growthrates_zeromissing; dims=Ti)[Ti(1)] .* eu_boolmask) .> 0)
# aus_growthmask |> plot
# us_growthmask |> plot
# eu_growthmask |> plot

### Load human population

human_pop_path = joinpath(basedir, "data/population_density.tif")
human_pop = GDALarray(human_pop_path; usercrs=EPSG(4326))[Band(1)] |>
            permutedims |> replace_missing
# human_pop |> plot
aus_humanpop = human_pop[aus...]
us_humanpop = human_pop[us...]
eu_humanpop = human_pop[eu...]
# aus_humanpop |> plot
# us_humanpop |> plot
# eu_humanpop |> plot

gdp_path = joinpath(basedir, "data/GDP_USD_per_m2.tif")
gdp = GDALarray(gdp_path; usercrs=EPSG(4326))[Band(1)] .* 1e6 |> permutedims
gdp = (@set gdp.missingval = minimum(gdp)) |> replace_missing
# gdp |> plot
aus_gdp = gdp[aus...]
us_gdp = gdp[us...]
eu_gdp = gdp[eu...]
# aus_gdp |> plot
# us_gdp |> plot
# eu_gdp |> plot

@assert size(eu_humanpop) == size(eu_growthmask[Ti(1)]) == size(eu_states)
@assert size(us_humanpop) == size(us_growthmask[Ti(1)]) == size(us_states)
@assert size(aus_humanpop) == size(aus_growthmask[Ti(1)]) == size(aus_states)

### Define Rules ###

carrycap = floattype(1e9)
growth = ExactLogisticGrowthMap{:population,:population}(
    layerkey=Val(:growthrates),
    carrycap=carrycap,
    timestep=Day(1),
    nsteps=floattype(1.0),
);
growth

### Constant growth ###

constant_growth = ExactLogisticGrowth{:population,:population}(
    intrinsicrate=floattype(0.1),
    timestep=Day(1),
    carrycap=carrycap
)

### Local dispersal ###

λ = floattype(0.0125)
radius = 1

@time hood = DispersalKernel{radius}(
    formulation=ExponentialKernel(λ),
    distancemethod=AreaToArea(30),
    cellsize=floattype(1.0)
)
localdisp = InwardsPopulationDispersal{:population,:population}(hood)
log.(hood.kernel) |> heatmap
savefig(joinpath(basedir, "output/dispersal_kernel.png"))

### Allee effects

allee = AlleeExtinction{:population,:population}(
    minfounders=floattype(22.0)
);

### Human Dispersal ###

scale = 10
hkwargs = (
    read=:population,
    scale=scale,
    human_exponent=floattype(1.075),
    dist_exponent=floattype(1.429),
    dispersalperpop=floattype(8.703e-9),
    max_dispersers=floattype(3.264e4),
    nshortlisted=100,
)

if humanactivity == :pop
    println("human activity from population")
    aus_humandisp = HumanDispersal(; human_pop=parent(aus_humanpop), hkwargs...)
    us_humandisp = HumanDispersal(; human_pop=parent(us_humanpop), hkwargs...)
    eu_humandisp = HumanDispersal(; human_pop=parent(eu_humanpop), hkwargs...)
elseif humanactivity == :gdp
    println("human activity from gdp")
    aus_humandisp = HumanDispersal(; human_pop=parent(aus_gdp), hkwargs...)
    us_humandisp = HumanDispersal(; human_pop=parent(us_gdp), hkwargs...)
    eu_humandisp = HumanDispersal(; human_pop=parent(eu_gdp), hkwargs...)
end

I = (DimensionalData.dims2indices(aus_states, (Lat(Contains(-36.805896)), Lon(Contains(143.959527)))) .-1) .÷ scale .+ 1
melb_dests = Dispersal.populate!(zero(aus_growthrates_zeromissing[Ti(1)]), aus_humandisp, I...)
plot(melb_dests, title="Melbourne destinations")
savefig(joinpath(basedir, "output/human_dispersal_destinations.png"))

### Define initialisation grids

assign_coord!(popgrid, (lat, lon), initsize) = begin
   coord = Lat(Contains(lat)), Lon(Contains(lon))
   popgrid[coord...] = initsize
end
init_popgrid!(popgrid, (lat, lon), initsize) = begin
   popgrid .= floattype(0.0)
    assign_coord!(popgrid, (lat, lon), initsize)
   nothing
end

aus_populationgrid = zero(aus_growthrates_zeromissing[Ti(1)])
init_popgrid!(aus_populationgrid, aus_incursionpoints[:Brisbane], carrycap)
us_populationgrid = zero(us_growthrates_zeromissing[Ti(1)])
init_popgrid!(us_populationgrid, us_incursionpoints[:San_Francisco], carrycap)
eu_populationgrid = zero(eu_growthrates_zeromissing[Ti(1)])
map(ip -> assign_coord!(eu_populationgrid, ip, carrycap), eu_incursionpoints)
# p = plot(eu_boolmask, label="")
# scatter!(p, reverse.(eu_incursionpoints); title="")
# eu_populationgrid |> plot

### Define the combined rulesets ###

# HumanDispersal has location specific precalculated
# data that the optimiser needs to update, making this a little
# messier than it otherwise would be

chain = Chain(localdisp, allee, growth)

# Test rule type stability
output = ArrayOutput((population=aus_populationgrid,);
    mask=collect(aus_boolmask),
    aux=(growthrates=aus_growthrates_zeromissing,),
    tspan=DateTime(2020, 1):timestep:DateTime(2026, 1)
)
DynamicGrids.isinferred(output, localdisp)
DynamicGrids.isinferred(output, allee)
DynamicGrids.isinferred(output, growth)
DynamicGrids.isinferred(output, chain)
DynamicGrids.isinferred(output, aus_humandisp)

ruleset_kwargs = (
    timestep=timestep,
    opt=SparseOpt(),
)
aus_ruleset = Ruleset(
    aus_humandisp, chain;
    ruleset_kwargs...
);
us_ruleset = Ruleset(
    us_humandisp, chain;
    ruleset_kwargs...
);
eu_ruleset = Ruleset(
    eu_humandisp, chain;
    ruleset_kwargs...
);
DynamicGrids.isinferred(output, us_ruleset)

### Output ###

# Define some color processors to use in live simulations.

zerocolor = ARGB32(0.6)
maskcolor = ARGB32(0.3)
textconfig = TextConfig(font="cantarel", bcolor=maskcolor, fcolor=ARGB32(0.8))

greyscale = ColorProcessor(Greyscale(), zerocolor, maskcolor, textconfig)
oranges = ColorProcessor(ColorSchemes.Oranges_3, zerocolor, maskcolor, textconfig)
jet = ColorProcessor(ColorSchemes.jet, zerocolor, maskcolor, textconfig)
viridis = ColorProcessor(ColorSchemes.viridis, zerocolor, maskcolor, textconfig)
inferno = ColorProcessor(ColorSchemes.inferno, zerocolor, maskcolor, textconfig)
magma = ColorProcessor(ColorSchemes.magma, zerocolor, maskcolor, textconfig)

max_stars_repo_name: rafaqz/DSuzukiiInvasionPaper

max_issues_repo_name: rafaqz/DSuzukiiInvasionPaper


text:
using .StaticArrays: SMatrix, SVector, MMatrix, MVector, SA

eval(load_projection_matrix_static())
eval(load_copy_finalize_static())


text:
function trainfunction2d(func,x,y)
    if func == 1
        a1 = 3 #rand(Uniform(-10,10))
        b1 = -4#rand(Uniform(-10,10))
        u = a1*x + b1*y
    elseif func == 2
        a = [0.9, 0.1, -0.4]#rand(Uniform(-1,1),3)
        b = [0.2, -0.6, -0.3]#rand(Uniform(-1,1),3)
        u = a[2]*sin(1*pi*y)+cos(1*pi*x)
    elseif func == 3
        a = [0.9464328087730216; 0.17162310348664844; 0.07857856828523646]    #rand(Uniform(-1,1),3)
        b = [0.9355991072384002, 0.5353806651162376, -0.0613925812939633]
        u = a[3]*sin(2*pi*y)+b[3]*cos(2*pi*x)
    elseif func == 4
        a = [-0.8752648836596824, -0.2937424545901659, 0.5352037919227479]#rand(Uniform(-1,1),3)
        b = [-0.9137179533411737, -0.4640307310161864, -0.8663071119500332]#rand(Uniform(-1,1),3)
        u1 = a[3]*sin(pi*x)+b[3]*cos(pi*y)
        u = u1 + a[1]*sin(2*pi*x)+b[1]*cos(2*pi*y)
    elseif func == 5
        a = [-0.6867267253726719, 0.21059347965868014, -0.728510882976293]#rand(Uniform(-1,1),3)
        b = [0.676235507814718, 0.8294240477938528, -0.3998500839840293]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[2]*cos(pi*y)
        u = u1 + a[3]*sin(2*pi*x)+b[3]*cos(2*pi*y)
    elseif func == 6
        a = [0.7983982959430316, 0.903381400725598, 0.6022370326216002]#rand(Uniform(-1,1),3)
        b = [-0.7513545425595294, -0.7714624763400932, -0.8408910576188569]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 7
        a = [0.4456995188427575, -0.760694865555275, 0.5341392645364422]#rand(Uniform(-1,1),3)
        b = [0.6038471708245794, -0.9293109017054246, -0.03067895573440227]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 8
        a1 = 8 #rand(Uniform(-10,10))
        b1 = -3 #rand(Uniform(-10,10))
        u = a1*x + b1*y
    elseif func == 9
        a =[0.2653910349699058, 0.4700076088116263, 0.09668422924668407, -0.046234710248600486, 0.20474920422233067, 0.5826924980829493]
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 10
        a =[-0.9017561219736625, -0.03447706434100306, 0.3305030118808854, 0.684511518915635, -0.661785410962842, 0.877713860698909]
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 11
        a =[-0.4, 1.0, 0.2, 0.6, -0.9, 0.8]   # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 13
        a = [0.8977676302413324, 0.5503759548450637, -0.6451820733536873]#rand(Uniform(-1,1),3)
        b = [0.8, -0.6, -0.3]#rand(Uniform(-1,1),3)
        u = a[2]*sin(1*pi*x)+cos(1*pi*y)
    elseif func == 12
        a =[-0.21587992491755914, -0.041287018187041724, -0.8635531658151585, -0.040103752643044555, -0.0632059673237162, -0.09508661768127391]
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 14
        a = [-0.9, 0.1, -0.4]#rand(Uniform(-1,1),3)
        b = [0.8, -0.6, -0.3]#rand(Uniform(-1,1),3)
        u = a[2]*sin(pi*x)+b[2]*cos(pi*y)
    elseif func == 15
        a1 = -1 #rand(Uniform(-10,10))
        b1 = 3 #rand(Uniform(-10,10))
        u = a1*x + b1*y
    elseif func == 16
        a = [0.2, 0.9, -0.6]#rand(Uniform(-1,1),3)
        b = [-0.1, 0.7, 0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(2*pi*x)+b[1]*cos(2*pi*y)
    elseif func == 17
        a = [-0.6, 0.2, -0.9]#rand(Uniform(-1,1),3)
        b = [-1.0, -0.8, 0.3]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 18
        a =[-0.4, 1.0, 0.2, 0.6, -0.9, 0.8]   # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 19 #10
        a = [0.2, 0.1, -0.4]#rand(Uniform(-1,1),3)
        b = [0.1, -0.6, -0.3]#rand(Uniform(-1,1),3)
        u = a[2]*sin(pi*x)+b[2]cos(pi*y)
    elseif func == 20
        a = [-0.6, 0.2, -0.9]#rand(Uniform(-1,1),3)
        b = [-1.0, -0.8, 0.3]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 21
        a = [0.3, 0.4, -1.0]#rand(Uniform(-1,1),3)
        b = [1.0, -0.4, 0.6]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 22
        a = [0.5, -0.9, -0.2]#rand(Uniform(-1,1),3)
        b = [0.1, 0.2, 0.8]#rand(Uniform(-1,1),3)
        u = a[1]*sin(2*pi*x)+b[1]*cos(2*pi*y)
    elseif func == 23
        a = [-0.2, -0.4, 0.9]#rand(Uniform(-1,1),3)
        b = [0.8, 0.6, -0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(pi*x)+b[1]cos(pi*y)
    elseif func == 24
        a = [-0.2, -0.4, 0.9]#rand(Uniform(-1,1),3)
        b = [0.8, 0.6, -0.1]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 25
        a = [-0.3, -0.2, 1.0]#rand(Uniform(-1,1),3)
        b = [1.0, -0.2, -0.6]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 26
        a = [-0.1, -0.4, 1.0]#rand(Uniform(-1,1),3)
        b = [-1.0, 0.1, 0.8]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
    elseif func == 27
        a = [0.8, 0.9, -0.2]#rand(Uniform(-1,1),3)
        b = [0.1, 0.2, -0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(2*pi*x)+b[1]*cos(2*pi*y)
    elseif func == 28
        a = [0.2, -0.4, 0.9]#rand(Uniform(-1,1),3)
        b = [-0.8, 0.6, -0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(pi*x)+b[1]cos(pi*y)
    elseif func == 29  #20
        a = [0.9, 0.4, -0.7]#rand(Uniform(-1,1),3)
        b = [0.9, -0.5, -0.3]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 30
        a = [0.7, 0.0, -1.0]#rand(Uniform(-1,1),3)
        b = [0.0, -0.2, -0.6]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 31
        a = [0.1, 1.0, 1.0]#rand(Uniform(-1,1),3)
        b = [-1.0, -0.3, 0.8]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
    elseif func == 32
        a = [0.1, 0.9, -0.2]#rand(Uniform(-1,1),3)
        b = [0.2, 0.2, -0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(2*pi*x)+b[1]*cos(2*pi*y)
    elseif func == 33
        a = [-0.3, -0.4, 0.9]#rand(Uniform(-1,1),3)
        b = [0.8, 0.6, -0.1]#rand(Uniform(-1,1),3)
        u = a[1]*sin(pi*x)+b[1]cos(pi*y)
    elseif func == 34
        a = [0.9, 0.1, -0.8]#rand(Uniform(-1,1),3)
        b = [0.9, 0.3, -0.3]#rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]*cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]*cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]*cos(3*pi*y)
    elseif func == 35
        a =[1.0, 0.3, -0.7, -0.2, 0.1, -0.9]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 36
        a =[-0.6, 1.0, 0.6, -0.8, -0.1, 0.7]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 37
        a =[0.9, -0.1, 0.1, -0.8, 0.4, -0.9]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 38
        a =[-1.0, 0.5, 0.2, -0.8, 0.0, -0.9]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 39
        a =[-0.4, 1.0, 0.2, 0.6, -0.9, 0.8]   # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 40
        a =[0.9, -1.0, 0.4, -0.2, 0.6, -0.2]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 41
        a =[0.4, -0.6, 1.0, 0.5, -0.5, 0.2]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 42
        a =[-0.2, 0.6, 0.3, 0.9, -0.7, 0.9]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 43
        a =[-1.0, -0.1, 0.1, 1.0, 0.1, 0.9]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 44
        a =[0.4, -0.2, 0.8, 0.6, -0.9, 0.8]   # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 45
        a =[1.0, 0.5, 0.9, -1.0, -0.1, 0.3]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 46
        a =[-0.4, -0.4, 0.9, 0.6, -0.9, -0.8]   # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    end
    return u
end

function troubledcellfunctionabs2d(x, y, a, m, x0, y0)
    u =  a*abs((y-y0)-m*(x-x0))
    return u
end

function troubledcellfunctionstep2d(x, y, ui, m, x0, y0)
    h1 = y0+m*(x-x0)
    h2 = y0-(1/m)*(x-x0)

    if y <= h1 && y < h2
        u=ui[1]
    elseif y >= h1 && y > h2
        u=ui[2]
    elseif y < h1 && y >= h2
        u=ui[3]
    elseif y > h1 && y <= h2
        u=ui[4]
    elseif y == h1 && y == h2
        u=ui[1]
    end
    return u
end

function roundstep2d(x, y, ui, r0, x0, y0)
    inicenter = SVector(x0, y0)
    x_norm = x - inicenter[1]
    y_norm = y - inicenter[2]
    r = sqrt(x_norm^2 + y_norm^2)

    # Calculate primitive variables
    u = r > r0 ? ui[1] : ui[2]
    return u
end

function good_cell2d(node_coord, length, func, m, x0, y0)
    x1 = node_coord[1]
    y1 = node_coord[2]
    x2 = node_coord[1] + length
    y2 = node_coord[2] + length

    if func == 1
        y1_abs = y0 + m*(x1-x0)
        y2_abs = y0 + m*(x2-x0)
        x1_abs = x0 + (1/m)*(y1-y0)
        x2_abs = x0 + (1/m)*(y2-y0)
        if  y1_abs <=y2 && y1_abs >= y1
            false
        elseif y2_abs <=y2 && y2_abs >= y1
            false
        elseif x1_abs <= x2 && x1_abs >= x1
            false
        elseif x2_abs <= x2 && x2_abs >= x1
            false
        else
            true
        end
    elseif func == 2 || func == 3
        y1_step = y0 + m*(x1-x0)
        y2_step = y0 + m*(x2-x0)
        x1_step = x0 + (1/m)*(y1-y0)
        x2_step = x0 + (1/m)*(y2-y0)
        y1_step2 = y0 -(1/m)*(x1-x0)
        y2_step2 = y0 -(1/m)*(x2-x0)
        x1_step2 = x0 - m*(y1-y0)
        x2_step2 = x0 - m*(y2-y0)
        if  (y1_step <=y2 && y1_step >= y1) ||  (y1_step2 <=y2 && y1_step2 >= y1 )
            false
        elseif (y2_step <=y2 && y2_step >= y1 ) || (y2_step2 <=y2 && y2_step2 >= y1  )
            false
        elseif (x1_step <= x2 && x1_step >= x1 ) ||  (x1_step2 <= x2 && x1_step2 >= x1)
            false
        elseif (x2_step <= x2 && x2_step >= x1)  || (x2_step2 <= x2 && x2_step2 >= x1)
            false
        else
            true
        end
    elseif func == 4
        inicenter = SVector(x0, y0)
        x_norm1 = x1 - inicenter[1]
        x_norm2 = x2 - inicenter[1]
        y_norm1 = y1 - inicenter[2]
        y_norm2 = y2 - inicenter[2]

        r1 = sqrt(x_norm1^2 + y_norm1^2)
        r2 = sqrt(x_norm2^2 + y_norm1^2)
        r3 = sqrt(x_norm1^2 + y_norm2^2)
        r4 = sqrt(x_norm2^2 + y_norm2^2)
        if (r1 <= m && r2 >=m) || (r1>= m && r2 <=m)
            false
        elseif (r1 <= m && r3 >=m) || (r1>= m && r3 <=m)
            false
        elseif (r2 <= m && r4 >=m) || (r2>= m && r4 <=m)
            false
        elseif (r3 <= m && r4 >=m) || (r3>= m && r4 <=m)
            false
        else
            true
        end
    end
end

function validfunction2d(func,x,y)
    if func == 1
        a = [-0.2, 0.8, -0.7    ]#rand(Uniform(-1,1),3)
        b = [-1.0, 0.4, -0.1]   #rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]cos(3*pi*y)
    elseif func == 2
        a = [0.4, -0.8, 0.1]    #rand(Uniform(-1,1),3)
        b = [0.1, -0.6, 0.9]    #rand(Uniform(-1,1),3)
        u1 = a[1]*sin(pi*x)+b[1]cos(pi*y)
        u2 = u1 + a[2]*sin(2*pi*x)+b[2]cos(2*pi*y)
        u = u2 + a[3]*sin(3*pi*x)+b[3]cos(3*pi*y)
    elseif func == 3
        a = [0.4]
        b = [0.1]
        u = a[1]*sin(pi*x)+b[1]cos(pi*y)
    elseif func == 4
        a = [-0.4]
        b = [0.9]
        u = a[1]*sin(pi*x)+b[1]sin(pi*y)
    elseif func == 5
        a = [0.4]
        b = [-0.7]
        u = a[1]*sin(2*pi*x)+b[1]sin(2*pi*y)
    elseif func == 6
        a =[0.3, -0.1]# rand(Uniform(-1,1),6)
        u = a[1]*x + a[2]*x^2
    elseif func == 7
        a =[-0.1, 0.6, 0.8]# rand(Uniform(-1,1),6)
        u = a[1]*x + a[2]*x^2 + a[3]*x^3
    elseif func == 8
        #a =[-0.5, -0.1, -0.8, 0.8]# rand(Uniform(-1,1),6)
        #u = a[1]*x + a[2]*x^2 + a[3]*x^3 + a[4]*x^4
        a =[1.0, 0.5, 0.9, -1.0, -0.1, 0.3]    # rand(Uniform(-1,1),6)
        u = exp(a[1]*((x-a[2])^2+(y-a[3])^2)) + exp(a[4]*((x-a[5])^2+(y-a[6])^2))
    elseif func == 9
        a =[0.8, -0.4, 0.1]# rand(Uniform(-1,1),6)
        u = a[1]*x + a[2]*x^2 + a[3]*x^3
    end
    return u
end


text:
println("Test mcmc")

@testset "mcmc.jl" begin
  @testset "subsetnamedtuple" begin
    state = MCMC.make_init_state(MySimpleModel(10))
    sample = MCMC.subsetnamedtuple(state, (:a, :b))
    @test Set(keys(sample)) == Set((:a, :b))
  end

  @testset "mcmc" begin
    state = MCMC.make_init_state(MySimpleModel(10))
    nsamps = 1000
    nburn = 44
    thin = 33
    output = mcmc(MySimpleModel(10), nsamps, init=state, nburn=nburn, thin=thin)
    chain = output.chain

    @test length(chain) == nsamps
    @test chain[end].a == nsamps * thin + nburn
    @test chain[1].a == nburn + thin
    @test all(chain[end].c .== nsamps * thin + nburn)
    @test all(chain[1].c .== nburn + thin)
  end

  @testset "thin, nburn" begin
    model = MySimpleModel(10)
    state = MCMC.make_init_state(model)
    nsamps = 1000
    nburn = 44
    thin = 33

    Random.seed!(1234)
    ref_chain = mcmc(model, nsamps * thin + nburn, init=state, exclude=[:c]).chain

    Random.seed!(1234)
    thinned_chain = mcmc(model, nsamps, init=state, thin=thin,
                         exclude=[:c]).chain

    Random.seed!(1234)
    burned_chain = mcmc(model, nsamps, init=state, nburn=nburn,
                        exclude=[:c]).chain

    Random.seed!(1234)
    thinned_and_burned_chain = mcmc(model, nsamps, init=state, nburn=nburn,
                                    thin=thin, exclude=[:c]).chain

    @test length(ref_chain) == nsamps * thin + nburn
    @test length(thinned_chain) == nsamps
    @test length(burned_chain) == nsamps
    @test length(thinned_and_burned_chain) == nsamps
    @test all(thinned_chain .== ref_chain[thin:thin:nsamps*thin])
    @test all(burned_chain .== ref_chain[nburn+1:nburn+nsamps])
    @test all(thinned_and_burned_chain .== ref_chain[nburn+thin:thin:nburn+thin*nsamps])
  end

  @testset "make_init_state" begin
    model = MySimpleModel(10)
    init = MCMC.make_init_state(model)

    Random.seed!(1234)
    ref_chain = mcmc(model, 20, init=init).chain

    Random.seed!(1234)
    uninitialized_chain = mcmc(model, 20).chain

    @test all(ref_chain .== uninitialized_chain)
  end

  @testset "gibbs" begin
    Random.seed!(0)
    spl = Gibbs(M1(),
                Conditional(:mu, (m, s) -> s.mu + 1),
                Conditional(:sigma, (m, s) -> s.sigma + 2))
    nsamps = 50000
    chain = mcmc(spl, nsamps).chain
    @test chain[end].mu == nsamps
    @test chain[end].sigma == nsamps * 2

    spl = Gibbs(M1(),
                Conditional(:mu, (m, s) -> s.mu + 1),
                Conditional(:sigma, (m, s) -> randn()))
    chain = mcmc(spl, nsamps, init=(mu=0.0, sigma=0.0)).chain
    @test chain[end].mu == nsamps
    @test isapprox(std(getindex.(chain, :sigma)), 1, atol=1e-2)
  end

  @testset "gibbs multi-step" begin
    spl = Gibbs(M2(),
                Conditional((:mu, :sigma),
                            (m, s) -> (mu = s.mu + 1, sigma = s.sigma + 2)),
                Conditional(:eta, (m, s) -> s.eta .- 1))
    nsamps = 1000
    chain = mcmc(spl, nsamps).chain
    @test chain[end].mu == nsamps
    @test chain[end].sigma == nsamps * 2
    @test all(chain[end].eta .== -nsamps)
  end

  # TODO: Add tests for callback.
end


text:
#=
run examples tests from the examples folder
=#

include(joinpath(@__DIR__, "Examples.jl"))
using Test

# uncomment path for writing to results CSV
csv = nothing
# csv = joinpath(mkpath(joinpath(@__DIR__, "..", "benchmarks", "raw")), "tests.csv")

import MathOptInterface
const MOI = MathOptInterface

# import GLPK
# glpk = MOI.OptimizerWithAttributes(
#     GLPK.Optimizer,
#     MOI.Silent() => true,
#     "tol_int" => 1e-10,
#     "tol_bnd" => 1e-10,
#     "mip_gap" => 1e-10,
# )

import Gurobi
gurobi = MOI.OptimizerWithAttributes(
    Gurobi.Optimizer,
    MOI.Silent() => true,
    "IntFeasTol" => 1e-9,
    "FeasibilityTol" => 1e-9,
    "MIPGap" => 1e-9,
    "DualReductions" => 0, # fix infeasible or unbounded statuses
    "InfUnbdInfo" => 1, # get ray of the primal OA continuous relaxation
)

# default MOIPajarito options
options = (;
    iteration_limit = 500,
    time_limit = 120.0,
    # verbose = false,
    oa_solver = gurobi,
    # oa_solver = glpk,
    use_iterative_method = true,
    # use_extended_form = false,
    # solve_relaxation = false,
    # solve_subproblems = false,
    # use_init_fixed_oa = false,
)

# list of names of JuMP examples to run
examples = [
    "completablepsd",
    "polyfacilitylocation",
    "polyregression",
    "twostagestochastic",
    "matrixcompletion",
    "matrixdecomposition",
    "matrixregression",
    "experimentdesign",
    "inversecovariance",
    "knapsack",
    "vectorregression",
    "ballpacking",
    "modulardesign",
]

Examples.load_examples(examples)

function run_examples_tests()
    @testset verbose = true "examples tests" begin
        Examples.run_examples(examples, ["test"], options, csv, false)
    end
    println()
    return
end

run_examples_tests()


text:
# include(raw"D:\programming\src\julia\dev\utilities\utilities.jl")

mutable struct DelaunayMeshTriangle
    alive::Bool

    # indices into the vertex array
    index_1::Int
    index_2::Int
    index_3::Int

    # indices into the triangle array
    parent::Int
    child_1::Int
    child_2::Int
    child_3::Int

    neighbor_12::Int
    neighbor_23::Int
    neighbor_31::Int

    function DelaunayMeshTriangle()
        obj = new()
        obj.alive = false
        obj.index_1 = 0
        obj.index_2 = 0
        obj.index_3 = 0

        obj.parent = 0
        obj.child_1 = 0
        obj.child_2 = 0
        obj.child_3 = 0

        obj.neighbor_12 = 0
        obj.neighbor_23 = 0
        obj.neighbor_31 = 0

        return obj
    end
    function DelaunayMeshTriangle(
        alive,
        index_1,
        index_2,
        index_3,
        parent,
        child_1,
        child_2,
        child_3,
        neighbor_12,
        neighbor_23,
        neighbor_31)

        obj = new()
        obj.alive = alive
        obj.index_1 = index_1
        obj.index_2 = index_2
        obj.index_3 = index_3

        obj.parent = parent
        obj.child_1 = child_1
        obj.child_2 = child_2
        obj.child_3 = child_3

        obj.neighbor_12 = neighbor_12
        obj.neighbor_23 = neighbor_23
        obj.neighbor_31 = neighbor_31

        return obj
    end
end

function Base.show(io::IO, obj::DelaunayMeshTriangle)
    str = util.matlab_display(obj)
    print(io, str)
end

mutable struct DelaunayMesh
    vertex_array::Vector{<:AbstractPoint}
    triangle_array::datastructs.Stack{DelaunayMeshTriangle}
    triangle_index::Matrix{Int}
    # triangle_array::Vector{DelaunayMeshTriangle}

    function DelaunayMesh(va::Vector{<:AbstractPoint})
        obj = new()
        L = length(va)
        obj.vertex_array = Vector{eltype(va)}(undef, L + 3)
        copyto!(obj.vertex_array, 4, va, 1, L)
        # obj.triangle_array = datastructs.Stack{DelaunayMeshTriangle}(div(L,3)+1)
        obj.triangle_array = datastructs.Stack{DelaunayMeshTriangle}(L)
        return obj
    end
end

DelaunayMesh() = DelaunayMesh( Vector{MutablePoint}(undef, 0) )

function create_initial_triangle(dmesh::DelaunayMesh)
    PointType = eltype(dmesh.vertex_array)
    tan_a = 1.0 # alpha = 45 degrees
    scale_factor = 1.12 # 1.12
    scale_L = 0.1
    (min_x, max_x, min_y, max_y) = extrema(dmesh.vertex_array[4:length(dmesh.vertex_array)])
    Lx = max_x.x - min_x.x
    Ly = max_y.y - min_y.y
    if Lx > Ly
        b_y = min_y.y - Ly * scale_L # make the triangle 10% bigger than necessary
        b_x = min_x.x + Lx/2
        B = PointType(b_x, b_y)
        a = Lx/2 + Ly / tan_a
        h = Ly + Lx/2 * tan_a
        a = a * scale_factor
        h = h * scale_factor
        X = PointType(B.x - a, B.y)
        Y = PointType(B.x + a, B.y)
        Z = PointType(B.x, B.y + h)
    else
        b_y = min_y.y + Ly/2 # make the triangle 10% bigger than necessary
        b_x = min_x.x - Lx * scale_L
        B = PointType(b_x, b_y)
        a = Ly/2 + Lx * tan_a
        h = Lx + Ly/2 * 1 / tan_a
        a = a * scale_factor
        h = h * scale_factor
        X = PointType(B.x, B.y + a)
        Y = PointType(B.x, B.y - a)
        Z = PointType(B.x + h, B.y)
    end
    return (X, Y, Z)
end

function insert_first_point!(dmesh::DelaunayMesh)
    p_idx = 4

    tri_pxy = DelaunayMeshTriangle()
    tri_pyz = DelaunayMeshTriangle()
    tri_pzx = DelaunayMeshTriangle()

    push!(dmesh.triangle_array, tri_pxy)
    push!(dmesh.triangle_array, tri_pyz)
    push!(dmesh.triangle_array, tri_pzx)

    tri_parent_idx = 1
    tri_parent = dmesh.triangle_array.data[tri_parent_idx]
    tri_parent.alive = false

    tri_pxy_idx = 2
    tri_pyz_idx = 3
    tri_pzx_idx = 4

    tri_pxy.parent = tri_parent_idx
    tri_pxy.alive = true

    tri_pyz.parent = tri_parent_idx
    tri_pyz.alive = true

    tri_pzx.parent = tri_parent_idx
    tri_pzx.alive = true

    tri_parent.child_1 = tri_pxy_idx
    tri_parent.child_2 = tri_pyz_idx
    tri_parent.child_3 = tri_pzx_idx

    tri_pxy.index_1 = p_idx
    tri_pxy.index_2 = 1
    tri_pxy.index_3 = 2

    tri_pxy.neighbor_12 = tri_pzx_idx
    tri_pxy.neighbor_23 = 0
    tri_pxy.neighbor_31 = tri_pyz_idx

    tri_pyz.index_1 = p_idx
    tri_pyz.index_2 = 2
    tri_pyz.index_3 = 3

    tri_pyz.neighbor_12 = tri_pxy_idx
    tri_pyz.neighbor_23 = 0
    tri_pyz.neighbor_31 = tri_pzx_idx

    tri_pzx.index_1 = p_idx
    tri_pzx.index_2 = 3
    tri_pzx.index_3 = 1

    tri_pzx.neighbor_12 = tri_pyz_idx
    tri_pzx.neighbor_23 = 0
    tri_pzx.neighbor_31 = tri_pxy_idx

    return dmesh
end

function create_triangle(dmesh::DelaunayMesh, tri_idx::Int)
    tri = dmesh.triangle_array.data[tri_idx]

    p1 = dmesh.vertex_array[tri.index_1]
    p2 = dmesh.vertex_array[tri.index_2]
    p3 = dmesh.vertex_array[tri.index_3]

    return Triangle(p1, p2, p3)
end

function is_in_triangle(dmesh::DelaunayMesh, tri_idx::Int, p_idx::Int)::Bool
    # tri = create_triangle(dmesh, tri_idx)
    # return is_in_triangle(tri, dmesh.vertex_array[p_idx])
    query_point = dmesh.vertex_array[p_idx]

    tri = dmesh.triangle_array.data[tri_idx]
    p1 = dmesh.vertex_array[tri.index_1]
    p2 = dmesh.vertex_array[tri.index_2]
    p3 = dmesh.vertex_array[tri.index_3]
    return is_in_triangle(p1, p2, p3, query_point)
end

function find_containing_triange(dmesh::DelaunayMesh, p_idx::Int)
    tri_idx = 1
    tri = dmesh.triangle_array.data[tri_idx]
    succes = true
    while !tri.alive
        if tri.child_1 != 0 && is_in_triangle(dmesh, tri.child_1, p_idx)
            child_tri_idx = tri.child_1
            succes = true
        elseif tri.child_2 != 0 &&  is_in_triangle(dmesh, tri.child_2, p_idx)
            child_tri_idx = tri.child_2
            succes = true
        elseif tri.child_3 != 0 &&  is_in_triangle(dmesh, tri.child_3, p_idx)
            child_tri_idx = tri.child_3
            succes = true
        else
            succes = false
            # try to shuffle the query point a bit
            x_ = dmesh.vertex_array[p_idx].x
            y_ = dmesh.vertex_array[p_idx].y

            dmesh.vertex_array[p_idx].x += sqrt(eps(x_))
            dmesh.vertex_array[p_idx].y += sqrt(eps(y_))

            parent_p1 = dmesh.vertex_array[tri.index_1]
            parent_p2 = dmesh.vertex_array[tri.index_2]
            parent_p3 = dmesh.vertex_array[tri.index_3]

            println("The requested point: [$(dmesh.vertex_array[p_idx].x), $(dmesh.vertex_array[p_idx].y)]")

            println("The parent triangle P1: [$(parent_p1.x), $(parent_p1.y)]")
            println("The parent triangle P2: [$(parent_p2.x), $(parent_p2.y)]")
            println("The parent triangle P3: [$(parent_p3.x), $(parent_p3.y)]")

            println("The triangle child indices: ($(tri.child_1), $(tri.child_2), $(tri.child_3))")

            if tri.child_1 != 0
                tri_ch = dmesh.triangle_array[tri.child_1]
                p1 = dmesh.vertex_array[tri_ch.index_1]
                p2 = dmesh.vertex_array[tri_ch.index_2]
                p3 = dmesh.vertex_array[tri_ch.index_3]

                println("The child 1 triangle P1: [$(p1.x), $(p1.y)]")
                println("The child 1 triangle P1: [$(p2.x), $(p2.y)]")
                println("The child 1 triangle P1: [$(p3.x), $(p3.y)]")
            end

            if tri.child_2 != 0
                tri_ch = dmesh.triangle_array[tri.child_2]
                p1 = dmesh.vertex_array[tri_ch.index_1]
                p2 = dmesh.vertex_array[tri_ch.index_2]
                p3 = dmesh.vertex_array[tri_ch.index_3]

                println("The child 2 triangle P1: [$(p1.x), $(p1.y)]")
                println("The child 2 triangle P1: [$(p2.x), $(p2.y)]")
                println("The child 2 triangle P1: [$(p3.x), $(p3.y)]")
            end

            if tri.child_3 != 0
                tri_ch = dmesh.triangle_array[tri.child_3]
                p1 = dmesh.vertex_array[tri_ch.index_1]
                p2 = dmesh.vertex_array[tri_ch.index_2]
                p3 = dmesh.vertex_array[tri_ch.index_3]

                println("The child 3 triangle P1: [$(p1.x), $(p1.y)]")
                println("The child 3 triangle P1: [$(p2.x), $(p2.y)]")
                println("The child 3 triangle P1: [$(p3.x), $(p3.y)]")
            end

            # error("Something went wrong with the point location")
        end
        if succes
            tri_idx = child_tri_idx
            tri = dmesh.triangle_array.data[tri_idx]
        end
    end
    return (tri, tri_idx)
end

function point_vector_in_triangle(dmesh, tri)
    return [
            dmesh.vertex_array[tri.index_1],
            dmesh.vertex_array[tri.index_2],
            dmesh.vertex_array[tri.index_3] ]
end

function __create_triangles_in_containing!(dmesh::DelaunayMesh, tri::DelaunayMeshTriangle, tri_idx::Int, p_idx::Int)
    # create 3 triangles in the containtment triangle
    tri_pxy_idx = dmesh.triangle_array.ptr + 1
    tri_pyz_idx = tri_pxy_idx + 1
    tri_pzx_idx = tri_pyz_idx + 1

    tri.alive = false
    tri.child_1 = tri_pxy_idx
    tri.child_2 = tri_pyz_idx
    tri.child_3 = tri_pzx_idx

    tri_pxy = DelaunayMeshTriangle()
    tri_pyz = DelaunayMeshTriangle()
    tri_pzx = DelaunayMeshTriangle()

    push!(dmesh.triangle_array, tri_pxy)
    push!(dmesh.triangle_array, tri_pyz)
    push!(dmesh.triangle_array, tri_pzx)

    tri_pxy.parent = tri_idx
    tri_pyz.parent = tri_idx
    tri_pzx.parent = tri_idx

    tri_pxy.alive = true
    tri_pyz.alive = true
    tri_pzx.alive = true

    tri_pxy.index_1 = p_idx
    tri_pxy.index_2 = tri.index_1
    tri_pxy.index_3 = tri.index_2

    tri_pyz.index_1 = p_idx
    tri_pyz.index_2 = tri.index_2
    tri_pyz.index_3 = tri.index_3

    tri_pzx.index_1 = p_idx
    tri_pzx.index_2 = tri.index_3
    tri_pzx.index_3 = tri.index_1

    tri_pxy.neighbor_12 = tri_pzx_idx
    tri_pxy.neighbor_23 = tri.neighbor_12
    tri_pxy.neighbor_31 = tri_pyz_idx

    tri_pyz.neighbor_12 = tri_pxy_idx
    tri_pyz.neighbor_23 = tri.neighbor_23
    tri_pyz.neighbor_31 = tri_pzx_idx

    tri_pzx.neighbor_12 = tri_pyz_idx
    tri_pzx.neighbor_23 = tri.neighbor_31
    tri_pzx.neighbor_31 = tri_pxy_idx

    # update the neighborhood relationships

    # X = dmesh.triangle_array.data[tri.index_1]
    # Y = dmesh.triangle_array.data[tri.index_2]
    # Z = dmesh.triangle_array.data[tri.index_3]

    X = dmesh.vertex_array[tri.index_1]
    Y = dmesh.vertex_array[tri.index_2]
    Z = dmesh.vertex_array[tri.index_3]

    if tri.neighbor_12 != 0
        tri_ayx = dmesh.triangle_array.data[tri.neighbor_12]
        x_local_idx = find_point_local_index(dmesh, tri_ayx, X)
        if x_local_idx == 1
            tri_ayx.neighbor_31 = tri_pxy_idx
        elseif x_local_idx == 2
            tri_ayx.neighbor_12 = tri_pxy_idx
        elseif x_local_idx == 3
            tri_ayx.neighbor_23 = tri_pxy_idx
        else
            error("Something went wrong")
        end
    end

    if tri.neighbor_23 != 0
        tri_bzy = dmesh.triangle_array.data[tri.neighbor_23]
        y_local_idx = find_point_local_index(dmesh, tri_bzy, Y)
        if y_local_idx == 1
            tri_bzy.neighbor_31 = tri_pyz_idx
        elseif y_local_idx == 2
            tri_bzy.neighbor_12 = tri_pyz_idx
        elseif y_local_idx == 3
            tri_bzy.neighbor_23 = tri_pyz_idx
        else
            error("Something went wrong")
        end
    end

    if tri.neighbor_31 != 0
        tri_cxz = dmesh.triangle_array.data[tri.neighbor_31]
        z_local_idx = find_point_local_index(dmesh, tri_cxz, Z)
        if z_local_idx == 1
            tri_cxz.neighbor_31 = tri_pzx_idx
        elseif z_local_idx == 2
            tri_cxz.neighbor_12 = tri_pzx_idx
        elseif z_local_idx == 3
            tri_cxz.neighbor_23 = tri_pzx_idx
        else
            error("Something went wrong")
        end
    end

    return dmesh
end

function neighbor_triangle_index(dmesh::DelaunayMesh, tri::DelaunayMeshTriangle, not_edge_point)::Int
    A = dmesh.vertex_array[tri.index_1]
    B = dmesh.vertex_array[tri.index_2]
    C = dmesh.vertex_array[tri.index_3]

    neighbor_index = 0
    if A == not_edge_point
        # edge BC
        neighbor_index = tri.neighbor_23
    elseif B == not_edge_point
        # edge CA
        neighbor_index = tri.neighbor_31
    elseif C == not_edge_point
        # edge AB
        neighbor_index = tri.neighbor_12
    else
        println("Point is not one of the triangle vertices")
    end
    return neighbor_index
end

function find_point_local_index(dmesh, tri, point)
    A = dmesh.vertex_array[tri.index_1]
    B = dmesh.vertex_array[tri.index_2]
    C = dmesh.vertex_array[tri.index_3]

    local_index = 0
    if A == point
        # edge BC
        local_index = 1
    elseif B == point
        # edge CA
        local_index = 2
    elseif C == point
        # edge AB
        local_index = 3
    else
        println("Point is not one of the triangle vertices")
    end

    return local_index
end

function point_index_from_local_index(tri, local_index::Int)
    idx = 0
    if local_index == 1
        idx = tri.index_1
    elseif local_index == 2
        idx = tri.index_2
    elseif local_index == 3
        idx = tri.index_3
    else
        println("Index out of bound")
    end
    return idx
end

function insert_existing_point!(dmesh::DelaunayMesh, p_idx::Int)
    L = length(dmesh.vertex_array)
    if p_idx > L
        println("Requested index greater than the length")
        return
    end
    P = dmesh.vertex_array[p_idx]
    (tri_parent, tri_parent_idx) = find_containing_triange(dmesh, p_idx)
    __create_triangles_in_containing!(dmesh, tri_parent, tri_parent_idx, p_idx)

    # loop over the neighboring triangles for checking the delaunay criterion
    # tri_pab_idx = dmesh.triangle_array.ptr-2
    # tri_pbc_idx = tri_pab_idx + 1
    # tri_pca_idx = tri_pbc_idx + 1

    tri_pab_idx = tri_parent.child_1
    tri_pbc_idx = tri_parent.child_2
    tri_pca_idx = tri_parent.child_3

    tri_pab = dmesh.triangle_array.data[tri_pab_idx]
    tri_pbc = dmesh.triangle_array.data[tri_pbc_idx]
    tri_pca = dmesh.triangle_array.data[tri_pca_idx]

    tri_stack = datastructs.Stack{DelaunayMeshTriangle}(3)
    push!(tri_stack, tri_pca)
    push!(tri_stack, tri_pbc)
    push!(tri_stack, tri_pab)

    while !isempty(tri_stack)
        act_tri = pop!(tri_stack)
        neighbor_index = neighbor_triangle_index(dmesh, act_tri, P)
        if neighbor_index == 0
            continue
        end
        neighbor_tri = dmesh.triangle_array.data[neighbor_index]

        x_global_idx = act_tri.index_2
        y_global_idx = act_tri.index_3
        X = dmesh.vertex_array[x_global_idx]
        Y = dmesh.vertex_array[y_global_idx]

        local_idx_x = find_point_local_index(dmesh, neighbor_tri, X)
        local_idx_y = find_point_local_index(dmesh, neighbor_tri, Y)
        if local_idx_x == 0
            error("Local index X cannot be zero")
        end
        if local_idx_y == 0
            error("Local index Y cannot be zero")
        end

        if local_idx_x == 1 && local_idx_y == 3
            local_idx_z = 2
        elseif local_idx_x == 2 && local_idx_y == 1
            local_idx_z = 3
        elseif local_idx_x == 3 && local_idx_y == 2
            local_idx_z = 1

        elseif local_idx_x == 1 && local_idx_y == 2
            local_idx_z = 3
        elseif local_idx_x == 2 && local_idx_y == 3
            local_idx_z = 1
        elseif local_idx_x == 3 && local_idx_y == 1
            local_idx_z = 2

        else
            error("I dont think the triangle is in ccw order")
        end
        z_global_idx = point_index_from_local_index(neighbor_tri, local_idx_z)
        Z = dmesh.vertex_array[z_global_idx]

        det_mat = in_circle(X, Z, Y, P)
        if abs(det_mat) <= 1e-10
            # fuzz the point
            pm = rand(1)
            s = +1.0
            if pm[1] < 0.5
                s = -1.0
            end
            fuzz_value = s * eps(P.x) # s * sqrt( eps(P.x) )
            fuzz_value = -s * eps(P.y) # -s * sqrt( eps(P.y) )
            P.x = P.x + fuzz_value
            P.y = P.y + fuzz_value
        end
        if det_mat > 0
            # harder part:
            #   - make the current triagngle dead
            #   - create two new triangles
            #   - push them to the stack
            #   - set the act_tri and neighbor_tri children 1,2 to this new triangles

            #   - reproduce neighborhood indices

            neighbor_tri.alive = false
            act_tri.alive = false

            tri_pzy = DelaunayMeshTriangle()
            tri_pxz = DelaunayMeshTriangle()
            push!(tri_stack, tri_pzy)
            push!(tri_stack, tri_pxz)

            new_tri_pzy_idx = dmesh.triangle_array.ptr + 1
            new_tri_pxz_idx = dmesh.triangle_array.ptr + 2
            push!(dmesh.triangle_array, tri_pzy)
            push!(dmesh.triangle_array, tri_pxz)

            tri_pzy.alive = true
            tri_pxz.alive = true

            tri_pzy.parent = neighbor_index
            tri_pxz.parent = neighbor_index

            neighbor_tri.child_1 = new_tri_pzy_idx
            neighbor_tri.child_2 = new_tri_pxz_idx
            act_tri.child_1 = new_tri_pzy_idx
            act_tri.child_2 = new_tri_pxz_idx

            tri_pzy.index_1 = p_idx
            tri_pzy.index_2 = z_global_idx
            tri_pzy.index_3 = y_global_idx

            tri_pxz.index_1 = p_idx
            tri_pxz.index_2 = x_global_idx
            tri_pxz.index_3 = z_global_idx

            # update neighbor indices
            if local_idx_x == 1
                tri_pzy_neighbor_23 = neighbor_tri.neighbor_23
                tri_pxz_neighbor_23 = neighbor_tri.neighbor_12
            elseif local_idx_x == 2
                tri_pzy_neighbor_23 = neighbor_tri.neighbor_31
                tri_pxz_neighbor_23 = neighbor_tri.neighbor_23
            elseif local_idx_x == 3
                tri_pzy_neighbor_23 = neighbor_tri.neighbor_12
                tri_pxz_neighbor_23 = neighbor_tri.neighbor_31
            else
                error("Something went wrong, for debugging only")
            end

            tri_pxz.neighbor_12 = act_tri.neighbor_12 # new_tri_pyz_idx
            tri_pxz.neighbor_23 = tri_pxz_neighbor_23 # ok
            tri_pxz.neighbor_31 = new_tri_pzy_idx # tri_pbc_idx

            tri_pzy.neighbor_12 = new_tri_pxz_idx # tri_pab_idx
            tri_pzy.neighbor_23 = tri_pzy_neighbor_23 # ok
            tri_pzy.neighbor_31 = act_tri.neighbor_31 # new_tri_pzx_idx

            # get the neighbors if exist
            # find the local indices in the neighbors
            # update the common edge neighbor
            if tri_pxz.neighbor_12 != 0
                # neighbor exists
                nei_tri = dmesh.triangle_array.data[tri_pxz.neighbor_12]
                local_idx_p_in_px_neighbor = find_point_local_index(dmesh, nei_tri, P)
                if local_idx_p_in_px_neighbor == 1
                    nei_tri.neighbor_31 = new_tri_pxz_idx
                elseif local_idx_p_in_px_neighbor == 2
                    nei_tri.neighbor_12 = new_tri_pxz_idx
                elseif local_idx_p_in_px_neighbor == 3
                    nei_tri.neighbor_23 = new_tri_pxz_idx
                else
                    error("Something went wrong")
                end
            end

            if tri_pxz.neighbor_23 != 0
                # neighbor exists
                nei_tri = dmesh.triangle_array.data[tri_pxz.neighbor_23]
                local_idx_x_in_xz_neighbor = find_point_local_index(dmesh, nei_tri, X)
                if local_idx_x_in_xz_neighbor == 1
                    nei_tri.neighbor_31 = new_tri_pxz_idx
                elseif local_idx_x_in_xz_neighbor == 2
                    nei_tri.neighbor_12 = new_tri_pxz_idx
                elseif local_idx_x_in_xz_neighbor == 3
                    nei_tri.neighbor_23 = new_tri_pxz_idx
                else
                    error("Something went wrong")
                end
            end

            if tri_pzy.neighbor_23 != 0
                # neighbor exists
                nei_tri = dmesh.triangle_array.data[tri_pzy.neighbor_23]
                local_idx_z_in_zy_neighbor = find_point_local_index(dmesh, nei_tri, Z)
                if local_idx_z_in_zy_neighbor == 1
                    nei_tri.neighbor_31 = new_tri_pzy_idx
                elseif local_idx_z_in_zy_neighbor == 2
                    nei_tri.neighbor_12 = new_tri_pzy_idx
                elseif local_idx_z_in_zy_neighbor == 3
                    nei_tri.neighbor_23 = new_tri_pzy_idx
                else
                    error("Something went wrong")
                end
            end

            if tri_pzy.neighbor_31 != 0
                # neighbor exists
                nei_tri = dmesh.triangle_array.data[tri_pzy.neighbor_31]
                local_idx_y_in_yp_neighbor = find_point_local_index(dmesh, nei_tri, Y)
                if local_idx_y_in_yp_neighbor == 1
                    nei_tri.neighbor_31 = new_tri_pzy_idx
                elseif local_idx_y_in_yp_neighbor == 2
                    nei_tri.neighbor_12 = new_tri_pzy_idx
                elseif local_idx_y_in_yp_neighbor == 3
                    nei_tri.neighbor_23 = new_tri_pzy_idx
                else
                    error("Something went wrong")
                end
            end

        else
            # easy part: do nothing
        end
    end

    return dmesh
end

using LinearAlgebra
function is_in_circle(tri_A::P, tri_B::P, tri_C::P, D::P)::Bool where {P <: AbstractPoint}
    det_mat = in_circle(tri_A, tri_B, tri_C, D)
    return det_mat >= 0
end

function in_circle(tri_A::P, tri_B::P, tri_C::P, D::P) where {P <: AbstractPoint}
    # mat = [
    # tri_A.x     tri_A.y     tri_A.x^2 + tri_A.y^2   1;
    # tri_B.x     tri_B.y     tri_B.x^2 + tri_B.y^2   1;
    # tri_C.x     tri_C.y     tri_C.x^2 + tri_C.y^2   1;
    # D.x         D.y         D.x^2 + D.y^2           1
    # ]

    sa = (tri_A.x - D.x)^2 + (tri_A.y - D.y)^2
    sb = (tri_B.x - D.x)^2 + (tri_B.y - D.y)^2
    sc = (tri_C.x - D.x)^2 + (tri_C.y - D.y)^2

    dx_ad = tri_A.x - D.x
    dy_ad = tri_A.y - D.y
    dx_bd = tri_B.x - D.x
    dy_bd = tri_B.y - D.y
    dx_cd = tri_C.x - D.x
    dy_cd = tri_C.y - D.y

    det_a = dx_bd * dy_cd - dx_cd * dy_bd
    det_b = dx_ad * dy_cd - dx_cd * dy_ad
    det_c = dx_ad * dy_bd - dx_bd * dy_ad

    det_mat = sa * det_a - sb * det_b + sc * det_c

    # det_mat = det(mat)
    return det_mat
end

function build_delaunay_mesh!(point_array)
    L = length(point_array)
    if L <= 3
        dmesh = DelaunayMesh()
        for kk = 1:L
            dmesh.vertex_array[kk] = point_array[kk]
        end
        return dmesh
    end
    # sort the point_array -> do not sort them!
    # sort!(point_array, lt = (p1, p2) -> p1.x < p2.x)

    dmesh = DelaunayMesh(point_array)

    # make an initial reshuffling
    # reshuffling_indices = randperm(L)
    # dmesh.vertex_array = dmesh.vertex_array[reshuffling_indices]
    # shuffle!( view(dmesh.vertex_array, 4:L) )

    (X, Y, Z) = create_initial_triangle(dmesh)

    dmesh.vertex_array[1] = X
    dmesh.vertex_array[2] = Y
    dmesh.vertex_array[3] = Z

    root_triangle = DelaunayMeshTriangle()
    root_triangle.index_1 = 1
    root_triangle.index_2 = 2
    root_triangle.index_3 = 3

    push!(dmesh.triangle_array, root_triangle)
    insert_first_point!(dmesh)

    p_idx_min = 2
    p_idx_max = length(dmesh.vertex_array)-3

    for p_idx = p_idx_min:p_idx_max
        insert_existing_point!(dmesh, 3 + p_idx)
    end

    # shrink to fit
    # dmesh.triangle_array

    return dmesh
end

function triangle_vector_for_plot(dmesh; init_triangle = false)
    tri_stack = datastructs.Stack{Triangle}()
    for kk = 1:dmesh.triangle_array.ptr
        tri_kk = dmesh.triangle_array.data[kk]
        if tri_kk.alive

            if !init_triangle
                if  tri_kk.index_1 == 1 || tri_kk.index_1 == 2 || tri_kk.index_1 == 3 ||
                    tri_kk.index_2 == 1 || tri_kk.index_2 == 2 || tri_kk.index_2 == 3 ||
                    tri_kk.index_3 == 1 || tri_kk.index_3 == 2 || tri_kk.index_3 == 3
                    continue
                end
            end

            push!(tri_stack, Triangle(
                dmesh.vertex_array[tri_kk.index_1],
                dmesh.vertex_array[tri_kk.index_2],
                dmesh.vertex_array[tri_kk.index_3]))
        end
    end
    tri_vec = tri_stack.data[1:tri_stack.ptr]
    return tri_vec
end

function triangle_indices(delmesh::geom.DelaunayMesh; init_triangle = false)

    L = length(delmesh.triangle_array)
    triidx = zeros(Int64, L, 3)

    last_idx = 0
    for ii = 1:L

        if !delmesh.triangle_array.data[ii].alive
            continue
        end

        if init_triangle

            last_idx += 1
            triidx[last_idx, 1] = delmesh.triangle_array.data[ii].index_1
            triidx[last_idx, 2] = delmesh.triangle_array.data[ii].index_2
            triidx[last_idx, 3] = delmesh.triangle_array.data[ii].index_3

        else

            tri_kk = delmesh.triangle_array.data[ii]
            if  tri_kk.index_1 == 1 || tri_kk.index_1 == 2 || tri_kk.index_1 == 3 ||
                tri_kk.index_2 == 1 || tri_kk.index_2 == 2 || tri_kk.index_2 == 3 ||
                tri_kk.index_3 == 1 || tri_kk.index_3 == 2 || tri_kk.index_3 == 3
                continue
            end

            last_idx += 1
            triidx[last_idx, 1] = delmesh.triangle_array.data[ii].index_1 - 3
            triidx[last_idx, 2] = delmesh.triangle_array.data[ii].index_2 - 3
            triidx[last_idx, 3] = delmesh.triangle_array.data[ii].index_3 - 3
        end

    end

    triidx = triidx[1:last_idx, :]

    return triidx

end

function envelope_points(delmesh::DelaunayMesh)
    return delmesh.vertex_array[1:3]
end

function valid_points(delmesh::DelaunayMesh)
    return delmesh.vertex_array[4:end]
end

function plot_triangulation(delmesh::DelaunayMesh; init_triangle = false)

    px = x(delmesh.vertex_array[4:end])
    py = y(delmesh.vertex_array[4:end])

    # px = point_array[:,1]
    # py = point_array[:,2]

    if isdefined(delmesh, :triangle_index)
        triangle_index = delmesh.triangle_index
    else
        triangle_index = triangle_indices(delmesh)
    end

    PyPlot.grid()
    PyPlot.plot( px, py,
        marker = :., linestyle = :none, markersize = 10)

    for tt = 1:size(triangle_index, 1)
        t = triangle_index[tt,:]
        tx = px[ [t[1], t[2], t[3], t[1]]]
        ty = py[ [t[1], t[2], t[3], t[1]]]
        PyPlot.plot(tx, ty)
    end

    #=
    p_x = x(delmesh.vertex_array)
    p_y = y(delmesh.vertex_array)

    tri_vec = triangle_vector_for_plot(delmesh; init_triangle = init_triangle)

    PyPlot.plot(p_x, p_y,
        color = "k",
        alpha = 0.7,
        linestyle = "none",
        marker = ".",
        markersize = 10)

    if init_triangle
        (init_tri_x, init_tri_y) = discretize_for_plot(
            [delmesh.vertex_array[1],
            delmesh.vertex_array[2],
            delmesh.vertex_array[3]], true)

        PyPlot.plot(init_tri_x, init_tri_y,
            marker = "o", linewidth = 2)
    end

    for kk = 1:length(tri_vec)
        (p_x, p_y) = discretize_for_plot(tri_vec[kk])
        PyPlot.plot(p_x, p_y,
        marker = "o",
        linewidth = 1,
        color = "r",
        alpha = 0.7,)
    end

    =#
end


text:
module BlogPlots

include("utils.jl")
# include("cpu.jl")
include("rockets.jl")
include("geerthofstede.jl")
include("petrolstates.jl")

end


text:
module MacroSpec

import Docile: @doc, meta

TESTCASE = "testcase"

## Testing methods. ##

@doc meta("docs_and_meta/1", key = :docs_and_meta) ->
function docs_and_meta(x)
end

@doc "docs_no_meta/1" ->
function docs_no_meta(x)
end

@doc meta(file = "test_case.md", key = :test_case) ->
function extern_docs_meta(x)
end

## Testing single line methods. ##

@doc meta("inline_docs_meta/1", key = :inline_docs_meta) ->
inline_docs_meta(x) = x

@doc "inline_docs_no_meta/1" ->
inline_docs_no_meta(x) = x

@doc meta(file = "test_case.md", key = :test_case) ->
inline_extern_docs_meta(x) = x

## Testing generic functions. ##

@doc meta("inline_docs_meta", key = :inline_docs_meta) ->
inline_docs_meta

@doc "inline_docs_no_meta" ->
inline_docs_no_meta

@doc meta(file = "test_case.md", key = :test_case) ->
inline_extern_docs_meta

## Testing types. ##

@doc meta("DocsAndMeta", key = :DocsAndMeta) ->
type DocsAndMeta
end

@doc "DocsNoMeta" ->
type DocsNoMeta
end

@doc meta(file = "test_case.md", key = :test_case) ->
type ExternDocsMeta
end

## Testing immutable types. ##

@doc meta("DocsAndMetaImm", key = :DocsAndMetaImm) ->
immutable DocsAndMetaImm
end

@doc "DocsNoMetaImm" ->
immutable DocsNoMetaImm
end

@doc meta(file = "test_case.md", key = :test_case) ->
immutable ExternDocsMetaImm
end

## Testing abstract types. ##

@doc meta("DocsAndMetaAbs", key = :DocsAndMetaAbs) ->
abstract DocsAndMetaAbs

@doc "DocsNoMetaAbs" ->
abstract DocsNoMetaAbs

@doc meta(file = "test_case.md", key = :test_case) ->
abstract ExternDocsMetaAbs

## Testing constants. ##

@doc meta("DocsAndMetaConst", key = :DocsAndMetaConst) ->
const DocsAndMetaConst = 1

@doc "DocsNoMetaConst" ->
const DocsNoMetaConst = 1

@doc meta(file = "test_case.md", key = :test_case) ->
const ExternDocsMetaConst = 1

## Testing globals. ##

@doc meta("DocsAndMetaGlobal", key = :DocsAndMetaGlobal) ->
DocsAndMetaGlobal = 1

@doc "DocsNoMetaGlobal" ->
DocsNoMetaGlobal = 1

@doc meta(file = "test_case.md", key = :test_case) ->
ExternDocsMetaGlobal = 1

## Testing macros. ##

@doc meta("macro_docs_and_meta", key = :macro_docs_and_meta) ->
macro docs_and_meta(x)
end

@doc "macro_docs_no_meta" ->
macro docs_no_meta(x)
end

@doc meta(file = "test_case.md", key = :test_case) ->
macro extern_docs_meta(x)
end

## Testing modules. ##

@doc meta("MacroSpec_meta", key = :MacroSpec_meta) ->
MacroSpec

@doc "MacroSpec_no_meta" ->
MacroSpec

@doc meta(file = "test_case.md", key = :test_case) ->
MacroSpec

end


text:
################################################################################################
#
# This comes from Knitro

# Default configuration.
const config = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                               optimal_status=MOI.OPTIMAL,
                               query=false,
                               infeas_certificates=false, # Do not ask for infeasibility certificates.
                               modify_lhs=false)

const config_local = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                                     optimal_status=MOI.LOCALLY_SOLVED,
                                     query=false,
                                     duals=false, # well presolve give the solution, hence no duals ...
                                     infeas_certificates=false, # Do not ask for infeasibility certificates.
                                     modify_lhs=false)

const config_noduals = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                                       optimal_status=MOI.OPTIMAL,
                                       query=false,
                                       duals=false,
                                       infeas_certificates=false,
                                       modify_lhs=false)

configs = Dict(
               "local" => config_local,
               "nodual" => config_noduals
              )

optimizer_kw = Dict(:rtol => 1e-9)

const OPTIMIZER = ReSHOP.Optimizer(; optimizer_kw...)

# MOIT.nlptest does not support :ExprGraph
#@testset "MOI NLP tests" begin
#    @testset "with $nlp_solver" for nlp_solver in nlp_solvers
#        MOIT.nlptest(ReSHOP.Optimizer(solver=nlp_solver; optimizer_kw...), get(config_solver, nlp_solver, config))
#    end
#end

@testset "MOI utils" begin
    @testset "SolverName" begin
        optimizer = ReSHOP.Optimizer()
        @test MOI.get(optimizer, MOI.SolverName()) == "ReSHOP"
    end
    @testset "supports_default_copy_to" begin
        optimizer = ReSHOP.Optimizer()
        @test MOIU.supports_default_copy_to(optimizer, false)
        # Use `@test !...` if names are not supported
        @test MOIU.supports_default_copy_to(optimizer, true)
    end
    @testset "MOI.Silent" begin
        optimizer = ReSHOP.Optimizer()
        @test MOI.supports(optimizer, MOI.Silent())
        MOI.set(optimizer, MOI.Silent(), true)
        @test MOI.get(optimizer, MOI.Silent()) == true
        MOI.set(optimizer, MOI.Silent(), false)
        @test MOI.get(optimizer, MOI.Silent()) == false
    end
#    @testset "MOI.TimeLimitSec" begin
#        optimizer = ReSHOP.Optimizer()
#        @test MOI.supports(optimizer, MOI.TimeLimitSec())
#        # TimeLimitSec is set to 1e8 by default in Knitro.
#        @test MOI.get(optimizer, MOI.TimeLimitSec()) == 1e8
#        my_time_limit = 10.
#        MOI.set(optimizer, MOI.TimeLimitSec(), my_time_limit)
#        @test MOI.get(optimizer, MOI.TimeLimitSec()) == my_time_limit
#    end
end

#@testset "MOI Basic constraint test" begin
#    MOIT.basic_constraint_tests(ReSHOP.Optimizer(), config, get_constraint_function=false, get_constraint_set=false)
#end

@testset "MOI unit test" begin

    model = ReSHOP.Optimizer(solver="cplex")
    exclude = ["number_threads",
			  "solve_with_lowerbound", # c1 is not registered
				"solve_singlevariable_obj", # problem with dual constraint value
			  ]

    MOI.Test.unittest(model, config, exclude)
end

@testset "MOI Linear tests" begin
    @testset "using $lp_solver" for lp_solver in lp_solvers
    optimizer = ReSHOP.Optimizer(solver=lp_solver; optimizer_kw...)
    exclude = [
               "linear1",
               "linear2", # DualObjectivevalue not supported
               "linear4", # ConstraintSet not supported
               "linear6", # ConstraintSet not supported
               "linear10", # No support of SAF in  Interval for now
               "linear10b", # No support of SAF in  Interval for now
               "linear11", # problem accessing constraint function
               "linear12", # Same as above.
               "linear14", # Delete not allowed
               "linear15", # DualObjectivevalue not supported
               ]
    append!(exclude, get(lp_exclude, lp_solver, []))
    MOIT.contlineartest(optimizer, config, exclude)
    # Test linear2 and linear15 without querying the dual solution.
    MOIT.linear2test(optimizer, config_noduals)
    # TODO URG this test now fails a lot
#    MOIT.linear15test(optimizer, config_noduals)
    end
end

@testset "MOI QP/QCQP tests" begin
    @testset "with $quad_solver" for (quad_solver, config_str) in quad_solvers
    # Exclude NCQP
    # Exclude qp2 and qp3, no solver supports it
    exclude = String["qp2", "qp3", "ncqcp1", "ncqcp2"]
    append!(exclude, get(quad_exclude, quad_solver, []))
    solver_config = get(configs, config_str, config)
    MOIT.contquadratictest(ReSHOP.Optimizer(solver=quad_solver; optimizer_kw...), get(config_solver, quad_solver, solver_config), exclude)
    end
end

@testset "MOI non-convex QP/QCQP tests" begin
    @testset "with $global_solver" for (global_solver, config_str) in global_solvers
    solver_config = get(configs, config_str, config)
    moi_solver = ReSHOP.Optimizer(solver=global_solver; optimizer_kw...)
    MOIT.ncqcp1test(moi_solver, get(config_solver, global_solver, solver_config))
    MOIT.ncqcp2test(moi_solver, get(config_solver, global_solver, solver_config))
    end
end

#= @testset "MOI SOCP tests" begin =#
#=     # TODO: DualObjectivevalue not supported =#
#=     # Presolve must be switch off to get proper dual variables. =#
#=     config2 = MOIT.TestConfig(atol=1e-4, rtol=1e-4, infeas_certificates=false, =#
#=                               optimal_status=MOI.LOCALLY_SOLVED, query=false) =#
#=     # Behavior in infeasible case doesn't match test. =#
#=     exclude = ["lin4"] =#
#=     BRIDGED2 = MOIB.full_bridge_optimizer(ReSHOP.Optimizer(outlev=0, presolve=0), Float64) =#
#=     MOIT.lintest(BRIDGED2, config2) =#
#= end =#

# We don't support ConstraintFunction for now
#@testset "MOI MILP test with $mip_solver" for mip_solver in mip_solvers
#    exclude = [
#               "int1",        # ObjectiveBound()
#               "int3",        # Need SAF in Interval
#               "indicator1",  # Sadly not support yet
#               "indicator2",  # Sadly not support yet
#               "indicator3",  # Need SAF in Interval
#              ]
#    MOIT.intlineartest(ReSHOP.Optimizer(solver=mip_solver; optimizer_kw...), config, exclude)
#end