text:
module Menus

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export menu

register_normal_element("q__menu", context = @__MODULE__)

function menu(
              fieldname::Union{Symbol,Nothing} = nothing,
              args...;
              content::Union{String,Vector} = "",
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)

  wrap() do
    q__menu(args...; attributes([:fieldname => fieldname, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...) do
      join(content)
    end
  end
end

function menu(content::Function,
              fieldname::Union{Symbol,Nothing} = nothing,
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  menu(label, fieldname, args...; wrap = wrap, content = content(), kwargs...)
end

end


text:
include("tests/qccdDev.jl")
using qccdSimulator
using Test
using LightGraphs

@testset "Read JSON" begin
    @test readJSONOK("./testFiles/topology.json")
    @test_throws ArgumentError("Input is not a file") readJSON(".")
    @test_throws ArgumentError readJSON("./testFiles/wrongTopology.json")
    @test readTimeJSONOK("./testFiles/times.json")
    @test readTimeJSONfail(["./testFiles/negativeTimes.json",
                                        "testFiles/wrongTimes.json","testFiles/wrongTimes2.json"])
    @test readTimeJSONnoFile()
end

@testset "QCCDevControl initialization" begin
    @test QCCDevCtrlOKTest()
    # @test nv(QCCDevCtrlTest().graph) == 5
    # @test ne(QCCDevCtrlTest().graph) == 6
    @test initJunctionsTest()
    @test initJunctionsTestRepId()
    @test initJunctionsTestIsolated()
    @test initJunctionsTestWrongType()
    @test initAuxGateZonesTestRepId()
    @test initAuxGateZonesTestInvZone()
    @test initAuxGateZonesTestWithNothing()
    @test initAuxZonesTest()
    @test initGateZoneTest()
    @test_throws ArgumentError("Repeated gate zone with ID: 1.") initGateZoneRepeatedIdTest()
    @test initLoadingZoneTest()
    @test_throws ArgumentError("Repeated loading zone with ID: 1") initLoadingZoneRepeatedIdTest()
    @test checkInitErrorsTest()
    @test checkInitErrorsTestEdgeCases()
end

@testset "Linear transport" begin
    @test isallowedLinearTransportTestTime()
    @test isallowedLinearTransportTestNoIon()
    @test isallowedLinearTransportTestNoZone()
    @test isallowedLinearTransportTestNonAdjacent()
    @test isallowedLinearTransportTestAllGood()
    @test isallowedLinearTransportTestFull()
    @test isallowedLinearTransportTestBlockedEnd0()
    @test isallowedLinearTransportTestBlockedEnd1()
    @test isallowedLinearTransportTestNotBlockedEnd0()
    @test isallowedLinearTransportTestNotBlockedEnd1()
    @test linearTransportTestOK()
    @test linearTransportTest1()
    @test linearTransportTest2()
    @test linearTransportTest3()

end

@testset "Utils" begin
    @test giveZoneTest()
end

@testset "Time checks" begin
    @test_throws OperationNotAllowedException("Time must be higher than " *
                        "10") time_check_timeFailsTest()
    @test_throws OperationNotAllowedException("Time model for test not " *
                        "defined.") time_check_modelFailsTest()
    @test time_checkOKTest()
end

@testset "Load ions `load()` & `isallowed_load()` & `initQubit()`" begin
    @test initQubitTest()
    @test isallowedLoad_OK()
    @test_throws OperationNotAllowedException("Loading zone with id test doesn't" *
                                    " exist.") isallowedLoad_zoneNotExistTest()
    @test_throws OperationNotAllowedException("Loading hole is " *
                                              "busy.") isallowedLoad_loadingHoleBusyTest()
    @test loadOKTest()
end

@testset "Swap & isallowed_swap" begin
    @test_throws OperationNotAllowedException("Qubit with id 99" *
                                              " doesn't exist.") isallowedSwap_qubitNotExist()
    @test_throws OperationNotAllowedException("Qubits with ids 1  " *
                                    "and 2 are not in the same zone.") isallowedSwap_qubitNotSameZone()
    @test_throws OperationNotAllowedException("Swap can only be done " *
                                             "in Gate Zones.") isallowedSwap_qubitNotGateZone()
    @test_throws OperationNotAllowedException("Qubits with ids 1 and 3 are not" *
                                             " adjacents.") isallowedSwap_qubitsNotAdjacents1()
    @test_throws OperationNotAllowedException("Qubits with ids 4 and 6 are not" *
                                             " adjacents.") isallowedSwap_qubitsNotAdjacents2()

    @test_throws OperationNotAllowedException("Qubits with ids 3 and 5 are not in the same" *
                                             " chain.") isallowedSwap_qubitsNotSameChain1()
    @test_throws OperationNotAllowedException("Qubits with ids 3 and 4 are not in the same" *
                                             " chain.") isallowedSwap_qubitsNotSameChain2()
    @test isallowedSwap_OK()
    @test swap_OK1()
    @test swap_OK2()
end

ext: jl

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
function RootsOfUnity()
    return RootsOfUnity(())
end

function compute_roots(obj::RootsOfUnity, arg0::jint)
    return jcall(obj, "computeRoots", void, (jint,), arg0)
end

function get_imaginary(obj::RootsOfUnity, arg0::jint)
    return jcall(obj, "getImaginary", jdouble, (jint,), arg0)
end

function get_number_of_roots(obj::RootsOfUnity)
    return jcall(obj, "getNumberOfRoots", jint, ())
end

function get_real(obj::RootsOfUnity, arg0::jint)
    return jcall(obj, "getReal", jdouble, (jint,), arg0)
end

function is_counter_clock_wise(obj::RootsOfUnity)
    return jcall(obj, "isCounterClockWise", jboolean, ())
end


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/invdet_calc.jl
## (c) 2018 Danley C. Hsu
# Script for producing DR25 FGK planet candidate occurrence rate estimates
#    using both the inverse detection efficiency and the simple
#    Bayesian methods

using ExoplanetsSysSim
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","dr25_gaia_fgk", "christiansen_func.jl"))

global sim_param_closure = setup_sim_param_christiansen()
sim_param_closure = set_test_param(sim_param_closure)

df_koi,usable_koi = read_koi_catalog(sim_param_closure)
println("# Finished reading in KOI data")
df_star = setup_star_table_christiansen(sim_param_closure)
println("# Finished reading in stellar data")
cat_obs = setup_actual_planet_candidate_catalog(df_star, df_koi, usable_koi, sim_param_closure)

#@time inv_det_simp_bayes(cat_obs, sim_param_closure)
@time simp_bayes(cat_obs, sim_param_closure)


text:
module SystemStructures

using DataStructures
using SymbolicUtils: istree, operation, arguments, Symbolic
using ..ModelingToolkit
import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten, value
using ..BipartiteGraphs
using UnPack
using Setfield
using SparseArrays

#=
When we don't do subsitution, variable information is split into two different
places, i.e. `states` and the right-hand-side of `observed`.

eqs = [0 ~ z + x; 0 ~ y + z^2]
states = [y, z]
observed = [x ~ sin(y) + z]
struct Reduced
    var
    expr
    idxs
end
fullvars = [Reduced(x, sin(y) + z, [2, 3]), y, z]
active_𝑣vertices = [false, true, true]
      x   y   z
eq1:  1       1
eq2:      1   1

      x   y   z
eq1:      1   1
eq2:      1   1

for v in 𝑣vertices(graph); active_𝑣vertices[v] || continue

end
=#

export SystemStructure, initialize_system_structure, find_linear_equations
export diffvars_range, dervars_range, algvars_range
export isdiffvar, isdervar, isalgvar, isdiffeq, isalgeq
export DIFFERENTIAL_VARIABLE, ALGEBRAIC_VARIABLE, DERIVATIVE_VARIABLE
export DIFFERENTIAL_EQUATION, ALGEBRAIC_EQUATION
export vartype, eqtype

struct SystemStructure
    dxvar_offset::Int
    fullvars::Vector # [diffvars; dervars; algvars]
    varassoc::Vector{Int}
    algeqs::BitVector
    graph::BipartiteGraph{Int,Nothing}
    solvable_graph::BipartiteGraph{Int,Nothing}
    assign::Vector{Int}
    inv_assign::Vector{Int}
    scc::Vector{Vector{Int}}
    partitions::Vector{NTuple{4, Vector{Int}}}
end

diffvars_range(s::SystemStructure) = 1:s.dxvar_offset
# TODO: maybe dervars should be in the end.
dervars_range(s::SystemStructure) = s.dxvar_offset+1:2s.dxvar_offset
algvars_range(s::SystemStructure) = 2s.dxvar_offset+1:length(s.fullvars)

isdiffvar(s::SystemStructure, var::Integer) = var in diffvars_range(s)
isdervar(s::SystemStructure, var::Integer) = var in dervars_range(s)
isalgvar(s::SystemStructure, var::Integer) = var in algvars_range(s)

@enum VariableType DIFFERENTIAL_VARIABLE ALGEBRAIC_VARIABLE DERIVATIVE_VARIABLE

function vartype(s::SystemStructure, var::Integer)::VariableType
    isdiffvar(s, var) ? DIFFERENTIAL_VARIABLE :
    isdervar(s, var)  ? DERIVATIVE_VARIABLE :
    isalgvar(s, var)  ? ALGEBRAIC_VARIABLE : error("Variable $var out of bounds")
end

@enum EquationType DIFFERENTIAL_EQUATION ALGEBRAIC_EQUATION

isalgeq(s::SystemStructure, eq::Integer) = s.algeqs[eq]
isdiffeq(s::SystemStructure, eq::Integer) = !isalgeq(s, eq)
eqtype(s::SystemStructure, eq::Integer)::EquationType = isalgeq(s, eq) ? ALGEBRAIC_EQUATION : DIFFERENTIAL_EQUATION

function initialize_system_structure(sys)
    sys, dxvar_offset, fullvars, varassoc, algeqs, graph = init_graph(flatten(sys))

    solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph))

    @set sys.structure = SystemStructure(
                                         dxvar_offset,
                                         fullvars,
                                         varassoc,
                                         algeqs,
                                         graph,
                                         solvable_graph,
                                         Int[],
                                         Int[],
                                         Vector{Int}[],
                                         NTuple{4, Vector{Int}}[]
                                        )
end

function Base.show(io::IO, s::SystemStructure)
    @unpack fullvars, dxvar_offset, solvable_graph, graph = s
    algvar_offset = 2dxvar_offset
    print(io, "xvars: ")
    print(io, fullvars[1:dxvar_offset])
    print(io, "\ndxvars: ")
    print(io, fullvars[dxvar_offset+1:algvar_offset])
    print(io, "\nalgvars: ")
    print(io, fullvars[algvar_offset+1:end], '\n')

    S = incidence_matrix(graph, Num(Sym{Real}(:×)))
    print(io, "Incidence matrix:")
    show(io, S)
end

function init_graph(sys)
    iv = independent_variable(sys)
    eqs = equations(sys)
    neqs = length(eqs)
    algeqs = trues(neqs)
    varsadj = Vector{Any}(undef, neqs)
    dervars = OrderedSet()
    diffvars = OrderedSet()

    for (i, eq) in enumerate(eqs)
        vars = OrderedSet()
        vars!(vars, eq)
        isalgeq = true
        for var in vars
            if istree(var) && operation(var) isa Differential
                isalgeq = false
                diffvar = arguments(var)[1]
                @assert !(diffvar isa Differential) "The equation [ $eq ] is not first order"
                push!(dervars, var)
                push!(diffvars, diffvar)
            end
        end
        algeqs[i] = isalgeq
        varsadj[i] = vars
    end

    algvars  = setdiff(states(sys), diffvars)
    fullvars = [collect(diffvars); collect(dervars); algvars]

    dxvar_offset = length(diffvars)
    algvar_offset = 2dxvar_offset

    nvars = length(fullvars)
    idxmap = Dict(fullvars .=> 1:nvars)
    graph = BipartiteGraph(neqs, nvars)

    vs = Set()
    for (i, eq) in enumerate(eqs)
        # TODO: custom vars that handles D(x)
        # TODO: add checks here
        lhs = eq.lhs
        if isdiffeq(eq)
            v = lhs
            haskey(idxmap, v) && add_edge!(graph, i, idxmap[v])
        else
            vars!(vs, lhs)
        end
        vars!(vs, eq.rhs)
        for v in vs
            haskey(idxmap, v) && add_edge!(graph, i, idxmap[v])
        end
        empty!(vs)
    end

    varassoc = Int[(1:dxvar_offset) .+ dxvar_offset; zeros(Int, length(fullvars) - dxvar_offset)] # variable association list
    sys, dxvar_offset, fullvars, varassoc, algeqs, graph
end

function find_linear_equations(sys)
    s = structure(sys)
    @unpack fullvars, graph = s
    is_linear_equations = falses(ndsts(graph))
    eqs = equations(sys)
    eadj = Vector{Int}[]
    cadj = Vector{Int}[]
    coeffs = Int[]
    for (i, eq) in enumerate(eqs); isdiffeq(eq) && continue
        empty!(coeffs)
        linear_term = 0
        all_int_algvars = true

        term = value(eq.rhs - eq.lhs)
        for j in 𝑠neighbors(graph, i)
            if !isalgvar(s, j)
                all_int_algvars = false
                continue
            end
            var = fullvars[j]
            c = expand_derivatives(Differential(var)(term), false)
            # test if `var` is linear in `eq`.
            if !(c isa Symbolic) && c isa Number
                if isinteger(c)
                    c = convert(Integer, c)
                else
                    all_int_algvars = false
                end
                linear_term += c * var
                push!(coeffs, c)
            end
        end

        # Check if there are only algebraic variables and the equation is both
        # linear and homogeneous, i.e. it is in the form of
        #
        #       ``∑ c_i * a_i = 0``,
        #
        # where ``c_i`` ∈ ℤ and ``a_i`` denotes algebraic variables.
        if all_int_algvars && isequal(linear_term, term)
            is_linear_equations[i] = true
            push!(eadj, copy(𝑠neighbors(graph, i)))
            push!(cadj, copy(coeffs))
        else
            is_linear_equations[i] = false
        end
    end
    return is_linear_equations, eadj, cadj
end

end # module


text:
using SafeTestsets
@safetestset "my_func tests" begin include("my_func_tests.jl") end


text:
using FastDMTransform, NPZ, Test, Statistics

function standardize(A; dims=1)
    μ = mean(A; dims=dims)
    σ = std(A; mean=μ, dims=dims)
    return @. (A - μ) / σ
end

@testset "FastDMTransform.jl" begin
    pulse = npzread("../pulse.npz")
    dedispersed = standardize(npzread("../dedispersed.npz"))
    dedisped = standardize(fdmt(pulse, 1500, 1200, 1e-3, 0, 2000))
    @test findmax(dedisped) == findmax(dedispersed)
end


text:
module Base64
import Base: read, write, close
export Base64Pipe, base64

# Base64Pipe is a pipe-like IO object, which converts writes (and
# someday reads?) into base64 encoded (decoded) data send to a stream.
# (You must close the pipe to complete the encode, separate from
# closing the target stream).  We also have a function base64(f,
# args...) which works like sprint except that it produces
# base64-encoded data, along with base64(args...)  which is equivalent
# to base64(write, args...), to return base64 strings.

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

type Base64Pipe <: IO
    io::IO
    # writing works in groups of 3, so we need to cache last two bytes written
    b0::Uint8
    b1::Uint8
    nb::Uint8 # number of bytes in cache: 0, 1, or 2

    function Base64Pipe(io::IO)
        b = new(io,0,0,0)
        finalizer(b, close)
        return b
    end
end

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

# Based on code by Stefan Karpinski from https://github.com/hackerschool/WebSockets.jl (distributed under the same MIT license as Julia)

const b64chars = ['A':'Z','a':'z','0':'9','+','/']

function b64(x::Uint8, y::Uint8, z::Uint8)
  n = int(x)<<16 | int(y)<<8 | int(z)
  b64chars[(n >> 18) + 1],
  b64chars[(n >> 12) & 0b111111 + 1],
  b64chars[(n >> 6) & 0b111111 + 1],
  b64chars[(n ) & 0b111111 + 1]
end

function b64(x::Uint8, y::Uint8)
  a, b, c = b64(x, y, 0x0)
  a, b, c, '='
end

function b64(x::Uint8)
  a, b = b64(x, 0x0, 0x0)
  a, b, '=', '='
end

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

function write(b::Base64Pipe, x::AbstractVector{Uint8})
    n = length(x)
    s = 1 # starting index
    # finish any cached data to write:
    if b.nb == 1
        if n >= 2
            write(b.io, b64(b.b0, x[1], x[2])...)
            s = 3
        elseif n == 1
            b.b1 = x[1]
            b.nb = 2
            return
        else
            return
        end
    elseif b.nb == 2
        if n >= 1
            write(b.io, b64(b.b0, b.b1, x[1])...)
            s = 2
        else
            return
        end
    end
    # write all groups of three bytes:
    while s + 2 <= n
        write(b.io, b64(x[s], x[s+1], x[s+2])...)
        s += 3
    end
    # cache any leftover bytes:
    if s + 1 == n
        b.b0 = x[s]
        b.b1 = x[s+1]
        b.nb = 2
    elseif s == n
        b.b0 = x[s]
        b.nb = 1
    else
        b.nb = 0
    end
end

function write(b::Base64Pipe, x::Uint8)
    if b.nb == 0
        b.b0 = x
        b.nb = 1
    elseif b.nb == 1
        b.b1 = x
        b.nb = 2
    else
        write(b.io, b64(b.b0,b.b1,x)...)
        b.nb = 0
    end
end

function close(b::Base64Pipe)
    if b.nb > 0
        # write leftover bytes + padding
        if b.nb == 1
            write(b.io, b64(b.b0)...)
        else # b.nb == 2
            write(b.io, b64(b.b0, b.b1)...)
        end
        b.nb = 0
    end
end

# like sprint, but returns base64 string
function base64(f::Function, args...)
    s = IOBuffer()
    b = Base64Pipe(s)
    f(b, args...)
    close(b)
    takebuf_string(s)
end
base64(x...) = base64(write, x...)

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

# read(b::Base64Pipe, ::Type{Uint8}) = # TODO: decode base64

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

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function mtdmin(
    latt_vars::LatticeVars, atm_vars::AtomicVars, rmt
)
    v1 = zeros(Float64,3)
    v2 = zeros(Float64,3)
    is = 1
    js = 1
    dmin = 1.0e6

    avec = latt_vars.avec
    epslat = latt_vars.epslat
    atposc = atm_vars.atposc
    nspecies = atm_vars.nspecies
    natoms = atm_vars.natoms

    for i1 in -1:1, i2 in -1:1, i3 in -1:1
        #
        @views v1[:] = i1*avec[:,1] + i2*avec[:,2] + i3*avec[:,3]
        #
        for ks in 1:nspecies, ka in 1:natoms[ks]
            #
            @views v2[:] = v1[:] + atposc[:,ka,ks]
            #
            for ls in 1:nspecies
                #
                t1 = rmt[ks] + rmt[ls]
                #
                for la in 1:natoms[ls]
                    if ( (i1 != 0) || (i2 != 0) || (i3 != 0) || (ks != ls) ||  (ka != la) )
                        t2 = sqrt( (v2[1] - atposc[1,la,ls])^2 +
                                   (v2[2] - atposc[2,la,ls])^2 +
                                   (v2[3] - atposc[3,la,ls])^2 )
                        t3 = t2 - t1
                        if t3 < (dmin - epslat)
                            is = ks
                            js = ls
                            dmin = t3
                        end
                    end
                end
            end
        end
    end

    return dmin, is, js

end


text:
####
#### placeholder for constants
####

export Constant

"""
`Contant(value)`

Placeholder for inserting a constant. Inverse checks equality with `==`.
"""
struct Constant{T} <: VectorTransform
    value::T
end

dimension(::Constant) = 0

function transform_with(logjac_flag::LogJacFlag, t::Constant, x::AbstractVector, index)
    t.value, logjac_zero(logjac_flag, eltype(x)), index
end

inverse_eltype(t::Constant, _) = Union{}

function inverse_at!(x::AbstractVector, index, t::Constant, y)
    @argcheck t.value == y
    index
end


text:
# Automatically generated using Clang.jl

const LSL_IRREGULAR_RATE = 0.0
const LSL_DEDUCED_TIMESTAMP = -1.0
const LSL_FOREVER = 3.2e7
const LSL_NO_PREFERENCE = 0

@cenum lsl_channel_format_t::UInt32 begin
    cft_float32 = 1
    cft_double64 = 2
    cft_string = 3
    cft_int32 = 4
    cft_int16 = 5
    cft_int8 = 6
    cft_int64 = 7
    cft_undefined = 0
end

@cenum lsl_processing_options_t::UInt32 begin
    proc_none = 0
    proc_clocksync = 1
    proc_dejitter = 2
    proc_monotonize = 4
    proc_threadsafe = 8
    proc_ALL = 15
end

@cenum lsl_error_code_t::Int32 begin
    lsl_no_error = 0
    lsl_timeout_error = -1
    lsl_lost_error = -2
    lsl_argument_error = -3
    lsl_internal_error = -4
end

# Skipping MacroDefinition: LIBLSL_C_API __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: LIBLSL_COMPILE_HEADER_VERSION = 113 ;

const lsl_streaminfo_struct_ = Cvoid
const lsl_streaminfo = Ptr{lsl_streaminfo_struct_}
const lsl_outlet_struct_ = Cvoid
const lsl_outlet = Ptr{lsl_outlet_struct_}
const lsl_inlet_struct_ = Cvoid
const lsl_inlet = Ptr{lsl_inlet_struct_}
const lsl_xml_ptr_struct_ = Cvoid
const lsl_xml_ptr = Ptr{lsl_xml_ptr_struct_}
const lsl_continuous_resolver_ = Cvoid
const lsl_continuous_resolver = Ptr{lsl_continuous_resolver_}


text:
module MaxHelpingHandFlagLogicGate

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/FlagLogicGate" FlagLogicGate(x::Integer, y::Integer, inputFlags::String="flag1,!flag2,flag3", outputFlag::String="flag4", func::String="AND", not::Bool=false)

const placements = Ahorn.PlacementDict(
    "Flag Logic Gate (max480's Helping Hand)" => Ahorn.EntityPlacement(
        FlagLogicGate
    )
)

function Ahorn.selection(entity::FlagLogicGate)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function Ahorn.editingOptions(effect::FlagLogicGate)
    return Dict{String, Any}(
        "func" => String["AND", "OR", "XOR"]
    )
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagLogicGate, room::Maple.Room)
    text = entity.func
    if entity.not
        if entity.func == "XOR"
            # this one is known as XNOR in electronics, not NXOR, so...
            text = "XNOR"
        else
            text = "N" * entity.func
        end
    end

    Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/flag_logic_gate", 0, 0)
    Ahorn.drawCenteredText(ctx, text, -12, 1, 24, 8, tint=(0.0, 0.0, 0.0, 1.0))
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
export logdensityof

using NestedTuples: lazymerge
import MeasureTheory

using Accessors

@inline function logdensityof(cm::AbstractConditionalModel, pars::NamedTuple)
    # TODO: Check insupport
    dynamic(unsafe_logdensityof(cm, pars))
end

@inline function MeasureBase.unsafe_logdensityof(cm::AbstractConditionalModel, pars::NamedTuple; cfg=NamedTuple(), ctx=NamedTuple())
    cfg = merge(cfg, (pars=pars,))
    ctx = merge(ctx, (ℓ = partialstatic(0.0),))
    gg_call(cm, unsafe_logdensityof, cfg, ctx, DropReturn())
end

@inline function tilde(::typeof(unsafe_logdensityof), lens, xname, x, d, cfg, ctx::NamedTuple)
    lm = lazymerge(cfg.obs, cfg.pars)
    x = NestedTuples._get(lm, xname)
    # x =get(lm, xname)
    @reset ctx.ℓ += MeasureBase.unsafe_logdensityof(d, lens(x))
    (x, ctx, ctx.ℓ)
end


text:
module Cribbage

export scorehand, discardtocrib

import IterTools
import Random: shuffle!
import UUIDs
using  PlayingCards

# Define Custom Exceptions
struct UnexpectedPlayerException <: Exception
    message :: String
end

# Define Custom Types
abstract type Player end
abstract type CribbageGame end
abstract type GameState end

function game(gamestate::GameState)
    return gamestate.game
end

function dealer(gamestate::GameState)::Player
    return gamestate.dealer
end

function name(player::Player)::String
    return player.name
end

function nextplayer(players::Array{Player,1}, player::Player)::Player
    i = findfirst(p -> p === player, players)
    if i == nothing
        throw(UnexpectedPlayerException("$(player) is not one of the players in this game"))
    end
    if i == length(players)
        return players[1]
    else
        return players[i+1]
    end
end

include("Rules.jl")
include("Scoring.jl")
include("TwoPlayer.jl")
include("ThreePlayer.jl")
include("FourPlayer.jl")
include("RandomPlay.jl")
include("SmartPlay.jl")
include("InteractivePlay.jl")

using  .Scoring
using  .Rules
import .TwoPlayer: incrementscore, numplayers, players, reportscore, cardsdealt, cardstocrib, isgameover, rotatedealer
using  .RandomPlay

macro returnifgameover(statesarray)
    return quote
        if isgameover($(esc(statesarray))[end])
            return $(esc(statesarray))
        end
    end
end

function dotheplay(gamestate::GameState, hands::Dict{Player,Array{Card,1}})
    gamestates = Array{typeof(gamestate),1}()
    # loop through the cards
    gameplayers = players(game(gamestate))
    gomap = Dict{Player,Bool}()

    player = dealer(gamestate)
    while any(length(h) > 0 for h in values(hands))
        for p in gameplayers
            gomap[p] = false
        end
        prior::Array{Card,1} = []
        priortotal = 0
        lastcard = nothing
        while !all(values(gomap))
            player = nextplayer(gameplayers, player)
            playerhand = hands[player]
            card = pickcardtoplay(player, tuple(prior...), tuple(playerhand...))
            if ismissing(card)
                gomap[player] = true
                continue
            end
            @assert card in playerhand
            @assert cardvalue(card) + priortotal ≤ 31
            lastcard = card
            gamestate = incrementscore(gamestate, player, scoreplay(prior, card))
            push!(gamestates, gamestate)
            push!(prior, card)
            if isgameover(gamestates[end])
                return gamestates
            end
            priortotal += cardvalue(card)
            hands[player] = setdiff(playerhand, [card])
        end
        total = sum(cardvalue(c) for c in prior)
        if total ≠ 31
            @debug "$(string(lastcard)): last card"
            gamestate = incrementscore(gamestate, player, 1)
            push!(gamestates, gamestate)
            if isgameover(gamestates[end])
                return gamestates
            end
        end
    end
    return gamestates
end

function dotheshow(gamestate::GameState, hands::Dict{Player,Array{Card,1}}, starter::Card, crib::Array{Card,1})
    gamestates = Array{typeof(gamestate),1}()
    # score each hand in order
    gameplayers = players(game(gamestate))
    player = dealer(gamestate)
    while true
        player = nextplayer(gameplayers, player)
        gamestate = incrementscore(gamestate, player, scorehand(hands[player], starter))
        push!(gamestates, gamestate)
        @returnifgameover gamestates
        if player == dealer(gamestate)
            gamestate = incrementscore(gamestate, player, scorehand(crib, starter))
            push!(gamestates, gamestate)
            @returnifgameover gamestates
            break
        end
    end
    return gamestates
end

function playhand(gamestate::GameState)
    gamestates = Array{typeof(gamestate),1}()

    @debug "Building a deck & shuffling"
    deck = createdeck()
    shuffle!(deck)
    crib::Array{Card,1} = []
    n = numplayers(game(gamestate))

    @debug "$(name(dealer(gamestate))) deals a hand"
    hands = deal!(
        deck, cardsdealt(game(gamestate)), numplayers(game(gamestate))
    )
    for i in 1:cardstocrib(game(gamestate))
        push!(crib, pop!(deck))
    end

    @debug "Forming the crib"
    playhands = Dict{Player,Array{Card,1}}()
    showhands = Dict{Player,Array{Card,1}}()
    gameplayers = players(game(gamestate))
    for (player, hand) in zip(gameplayers, hands)
        handcards, cribcards = discardtocrib(player, tuple(hand...))
        @assert handcards ⊊ hand
        @assert cribcards ⊊ hand
        for c in cribcards
            @assert c ∉ handcards
        end
        append!(crib, cribcards)
        playhands[player] = collect(handcards)
        showhands[player] = collect(handcards)
    end

    starter = pop!(deck)
    @debug "The starter card is $(string(starter))"
    gamestate = incrementscore(gamestate, dealer(gamestate), scoreheels(starter))
    push!(gamestates, gamestate)
    @returnifgameover gamestates

    # The play
    append!(gamestates, dotheplay(gamestates[end], playhands))
    @returnifgameover gamestates

    # The Show
    append!(gamestates, dotheshow(gamestates[end], showhands, starter, crib))
    @returnifgameover gamestates

    @debug "Updating the dealer to the next player"
    push!(gamestates, rotatedealer(gamestates[end]))
    return gamestates
end

function playgame(gamestate::T) where T <: GameState
    @info "Playing a game with $(name(dealer(gamestate))) as the current dealer"
    gamestates = Array{T,1}()
    push!(gamestates, gamestate)
    while !isgameover(gamestates[end])
        append!(gamestates, playhand(gamestates[end]))
    end
    @info reportscore(gamestates[end])
    @info "End of game!"
    return gamestates
end

end  # module Cribbage


text:
function add_deeplearning_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(load_weights_cb,b["dl_load_weights"],"clicked",Void,(),false,(handles,))
    signal_connect(load_labels_cb,b["dl_load_labels"],"clicked",Void,(),false,(handles,))
    signal_connect(dl_save_weights_cb,b["dl_save_weights"],"clicked",Void,(),false,(handles,))
    signal_connect(dl_save_labels_cb,b["dl_save_labels"],"clicked",Void,(),false,(handles,))

    signal_connect(create_training_cb,b["dl_create_model_button"],"clicked",Void,(),false,(handles,))

    signal_connect(training_button_cb,b["dl_train_button"],"clicked",Void,(),false,(handles,))
    signal_connect(epochs_sb_cb,b["dl_epoch_button"],"value-changed",Void,(),false,(handles,))
    signal_connect(confidence_sb_cb,b["dl_confidence_sb"],"value-changed",Void,(),false,(handles,))

    signal_connect(predict_frames_cb,b["dl_predict_button"],"clicked",Void,(),false,(handles,))

    signal_connect(create_config_cb,b["dl_export_training"],"clicked",Void,(),false,(handles,true))

    nothing
end

function load_weights_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    config_path = open_dialog("Load Previous Weights",han.b["win"])

    if config_path != ""

        try
            load_hourglass_to_nn(han.nn,config_path)
            setproperty!(han.b["dl_weights_label"],:label,config_path)
        catch
            println("Could not load weights")
        end
    end

    nothing
end

function load_labels_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    config_path = open_dialog("Load Previous Labels",han.b["win"])

    if config_path != ""

        try
            load_training(han,config_path)
            setproperty!(han.b["dl_labels_label"],:label,config_path)
            han.nn.use_existing_labels = true
        catch
            println("Could not load labeled data")
        end

        println("Previous Session Loaded")
    end

    nothing
end

function dl_save_labels_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    save_path = save_dialog("Save Labels",han.b["win"])

    if save_path != ""
        save_training(han,save_path)
    end

end

function create_training_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    try
        if (!han.nn.use_existing_labels)|(han.nn.norm.mean_img == zeros(Float32,0,0,0))
            set_up_training(han) #heatmaps, labels, normalize, augment
            save_training(han)
            han.nn.use_existing_labels=true
        else
            set_up_training(han,false)
            save_training(han)
        end

        if !han.nn.use_existing_weights
            create_new_weights(han.nn)
        end

        set_gtk_property!(han.b["create_model_label"],:label,string("model created at ", Dates.Time(Dates.now())))
    catch
        println("Could not create new training model")
    end

    nothing
end

function training_button_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    check_hg_features(han.nn)

    try
        dtrn=make_training_batch(han.nn.imgs,han.nn.labels);

        myadam=Adam(lr=1e-3)
        @async begin
            run_training(han.nn.hg,dtrn,myadam,han.b["dl_prog"],han.nn.epochs,han.nn.losses)
        end

    catch

        println("Failed to start training. Will use half the batch size and try again")
        dtrn=make_training_batch(han.nn.imgs,han.nn.labels,4);

        myadam=Adam(lr=1e-3)
        @async begin
            run_training(han.nn.hg,dtrn,myadam,han.b["dl_prog"],han.nn.epochs,han.nn.losses)
        end
    end
    #save_hourglass(string(han.paths.backup,"weights.jld"),han.nn.hg)

    nothing
end

function epochs_sb_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.nn.epochs=getproperty(han.b["dl_epochs_adjustment"],:value,Int)

    nothing
end

function confidence_sb_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.nn.confidence_thres=getproperty(han.b["dl_confidence_adjustment"],:value,Float64)

    nothing
end

function dl_save_weights_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    save_path = save_dialog("Save Weights",han.b["win"])

    if save_path != ""
        save_hourglass(save_path,han.nn.hg)
    end

    nothing
end

function predict_frames_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    @async begin
        han.nn.predicted=calculate_whiskers(han)
    end

    nothing
end

function create_config_cb(w::Ptr,user_data::Tuple{Tracker_Handles,Bool})

    han, training = user_data

    create_config(han,training)

    nothing
end


text:
using LinearAlgebra

"""
	face_mapping( V, FV, sigma [, err=10e-8] )

Compute the `mapping` matrix to transform `sigma` and faces in `sp_idx[sigma]` to ``z=0``.
"""
function face_mapping( V, FV, sigma, err=10e-8 )
	vs = FV[sigma]; i = 1
	# compute affinely independent triple
	n = length(vs)
	succ(i) = i % n + 1
	a = Lar.normalize(V[:,vs[succ(i)]] - V[:,vs[i]])
	b = Lar.normalize(V[:,vs[succ(succ(i))]] - V[:,vs[i]])
	c = cross(a,b)
	while (-err < det([a b c]) < err)
		i += 1
		a = Lar.normalize(V[:,vs[succ(i)]] - V[:,vs[i]])
		b = Lar.normalize(V[:,vs[succ(succ(i))]] - V[:,vs[i]])
		c = LinearAlgebra.cross(a,b)
	end
	# sigma translation
	T = Lar.t(-V[:,vs[succ(i)]]...)
	W = T * [ V[:,vs] ; ones(1,length(vs)) ]
	t_vs = W[1:3,:]
	# translated sigma rotation
	r1,r2,r3 = t_vs[:,i], t_vs[:,succ(succ(i))], cross(t_vs[:,i], t_vs[:,succ(succ(i))])
	R = Matrix{Float64}(I, 4, 4)
	R[1:3,1:3] = inv([r1 r2 r3])
	mapping = R * T
	return mapping
end

"""
	sigmamodel(V::Lar.Points, copEV::Lar.ChainOP, FV::Lar.Cells, copFE::Lar.ChainOP,
		sigma::Int, sp_idx::Array)

Transform `sigma` and `sp_idx[sigma]` faces to ``z=0`` space.
"""
function sigmamodel(V, copEV, FV, copFE, sigma, sp_idx)
	bigpi = sp_idx[sigma]
	Q = face_mapping(V', FV, sigma)
	EV = [findnz(copEV[e,:])[1] for e=1:size(copEV,1)]
	FE = [findnz(copFE[e,:])[1] for e=1:size(copFE,1)]
	sigma_edges = EV[FE[sigma]]
	sigma_verts = union(sigma_edges...)
	sigma_vdict =  Dict(zip(sigma_verts, 1:length(sigma_verts)))
	sigma_lines = [[sigma_vdict[v] for v in edge] for edge in sigma_edges]
	# sigma face on z=0 plane
	S = V'[:,sigma_verts] # sigma vertices
	Z = (Q * [S; ones(1,size(S,2))])[1:3,:] # sigma mapped in z=0
	v,ev = Z,sigma_lines
	#Plasm.view(Plasm.numbering(0.5)((v,[[[k] for k=1:size(v,2)],ev])))
	return v,ev,Q
end

"""
	sigma_intersect(V, EV, FE, sigma, Q, bigpi)

Intersect `sigma` plane `z=0` with edges in transformed `bigpi`.
"""
function sigma_intersect(V, EV, FE, sigma, Q, bigpi)
	sigma_edges = EV[FE[sigma]]
	sigma_verts = union(sigma_edges...)
	sigma_vdict =  Dict(zip(sigma_verts, 1:length(sigma_verts)))
	sigma_lines = [[sigma_vdict[v] for v in edge] for edge in sigma_edges]
	# sigma face on z=0 plane
	S = V[:,sigma_verts] # sigma vertices
	Z = (Q * [S; ones(1,size(S,2))])[1:3,:] # sigma mapped in z=0
	#Plasm.view(Z, sigma_lines)
	# initialization of line storage (to be intersected)
	linestore = [[Z[:,v1] Z[:,v2]] for (v1,v2) in sigma_lines]
	linenum = length(sigma_lines)
	# reindexing of ∏(σ) chain
	bigpi_edges = [EV[fe] for fe in FE[bigpi]]
	bigpi_verts = union(union(bigpi_edges...)...)
	bigpi_vdict =  Dict(zip(bigpi_verts, 1:length(bigpi_verts)))
	bigpi_lines = [[sort([bigpi_vdict[v] for v in edge]) for edge in faceedges]
		for faceedges in bigpi_edges]
	# bigpi trasformed in 3D according to Q mapping  ==> BUG ??!! :  NO
	P = V[:,bigpi_verts] # bigpi vertices
	W = (Q * [P; ones(1,size(P,2))])[1:3,:] # bigpi mapped by Q
	#Plasm.view(W, union(bigpi_lines...))
	# filter on bigpi_lines that do not cross z=0
	filtered_edges = [[[v1,v2] for (v1,v2) in face_lines
		if (sign(W[3,v1]) * sign(W[3,v2])) <= 0]
			for face_lines in bigpi_lines]
	filtered_edges = [edge for edge in filtered_edges if edge!=[]]
	#Plasm.view(W, union(filtered_edges...))
	#Plasm.view(Plasm.numbering()((W,[[[k] for k=1:size(W,2)],union(filtered_edges...)])))
	# computation of ordered z=0 points by face
	facepoints = []
	for face in filtered_edges
		points = Array{Float64,1}[]
		for (v1,v2) in face
			#v1[3] + (v2[3]-v1[3])*t = 0 # compute parameter of intersection point with z=0
			if (W[3,v2] - W[3,v1]) != 0
				t = -W[3,v1] / (W[3,v2] - W[3,v1])  # z1 + (z2-z1)*t = 0
				point = W[:,v1] + (W[:,v2] - W[:,v1]) * t
				push!(points,point)
			end
		end
		if length(Set(points)) > 1
			push!(facepoints,Set(points))
		end
		# TODO: linearly order facepoints[face]
	end
	#Plasm.view(Lar.Struct([ (hcat(collect(edge)...), [[1,2]]) for edge in facepoints ]))
	# compute z_lines
	c = collect
	z0_lines = [[[c(ps)[k] c(ps)[k+1]] for k=1:2:(length(ps)-1)] for ps in facepoints]

	# intersecting lines upload in linestore
	linestore = [[Z[:,v1] Z[:,v2]] for (v1,v2) in sigma_lines]
	if z0_lines != []
		linestore = append!(linestore, union(z0_lines...))
	end
	return linestore,linenum
end

"""
	computeparams(linestore,linenum)

Intersect all (non-simmetric) pairs of `line` in `linestore`.
"""
function computeparams(linestore,linenum)
	m = length(linestore)
	params = [[] for i=1:m]
	for h=1:(m-1)
		(x1,y1,x2,y2) = reshape(linestore[h][1:2,:],1,4)
		for k=(h+1):m
			(x3,y3,x4,y4) = reshape(linestore[k][1:2,:],1,4)
			# intersect lines h,k
			# http://www.cs.swan.ac.uk/~cssimon/line_intersection.html
			det = (x4-x3)*(y1-y2)-(x1-x2)*(y4-y3)
			α,β = 99,99
			if det != 0.0
				a = 1/det
				b = [y1-y2 x2-x1; y3-y4 x4-x3]  # x1-x2 => x2-x1 bug in the source link !!
				c = [x1-x3; y1-y3]
				(β,α) = a * b * c
			else
				if (y1==y2) == (y3==y4) # segments collinear
					α = -x1/(x2-x1)
					β = -x3/(x4-x3)
				elseif (x1==x2) == (x3==x4)
					α = -y1/(y2-y1)
					β = -y3/(y4-y3)
				else
					 # segments parallel: no intersection
					 (β,α) = 999
				end
			end
			if 0<α<1 && 0<β<1
				push!(params[h], α)
				push!(params[k], β)
			end
		end
	end
	out = Array{Float64,1}[]
	for line in params
		push!(line, 0.0)
		push!(line, 1.0)
		line = sort(collect(Set(line)))
		push!(out, line)
	end
	return out
end

"""
	fragface(V::Lar.Points, EV::Lar.Cells, FV::Lar.Cells, FE::Lar.Cells,
		sp_idx::Array, sigma::Int)::Lar.LAR

"""
function fragface(V, cop_EV, cop_FE, sp_idx, sigma)
	# format conversion of function parameters
	EV = [findnz(cop_EV[k,:])[1] for k=1:size(cop_EV,1)]
	FE = [findnz(cop_FE[k,:])[1] for k=1:size(cop_FE,1)]
	FV = [collect(Set(cat(EV[e] for e in FE[f]))) for f=1:length(FE)]
	V = convert(Lar.Points, V')

	Q = face_mapping(V, FV, sigma)
	bigpi = sp_idx[sigma]
	linestore, linenum = sigma_intersect(V, EV, FE, sigma, Q, bigpi)
	lineparams = computeparams(linestore,linenum)
	pairs = collect(zip(lineparams,linestore))

	linepoints = Array{Array{Float64,1},1}[]
	number = 0
	for (k,(params,line)) in enumerate(pairs)
		v1,v2 = line[1:2,1],line[1:2,2]
		if params==[0.0, 1.0]
			points = [v1,v2]
			if k<=linenum number += 1 end
		elseif length(params)>2
			points = [v1+α*(v2-v1) for α ∈ params]
			if k<=linenum number += length(params)-1 end
		end
		push!(linepoints,points)
	end
	verts = []
	edges = Array{Array{Int},1}()
	offset = 0
	for (h,pointline) in enumerate(linepoints)
		offset
		append!(verts, pointline)
		append!(edges, [[k+offset,k+1+offset] for k=1:length(pointline)-1])
		offset = length(verts)
	end
	#verts = hcat(verts...)
	#verts = convert(Lar.Points,verts')
	copEV = Lar.coboundary_0(convert(Lar.Cells,edges))
	copFE = Lar.coboundary_1(verts,faces,edges)
	return verts, copEV, copFE
end


text:
module ERA5Reanalysis

## Base Modules Used
using Base64
using DelimitedFiles
using Logging
using Printf

## Modules Used
using HTTP
using JSON
using NCDatasets
using PrettyTables
using Statistics

## Reexporting exported functions within these modules
using Reexport
@reexport using Dates
@reexport using GeoRegions

import Base: show, read

## Exporting the following functions:
export
        ERA5Module, ERA5Hourly, ERA5Monthly, ERA5MonthlyHour,

        ERA5Variable,
        SingleVariable, SingleCustom, PressureVariable, PressureCustom,
        listSingles,   isSingle,   rmSingle,   resetSingles,   tableSingles,
        listPressures, isPressure, rmPressure, resetPressures, tablePressures,
        resetERA5Variables, addERA5Variables, rmERA5Variable,  tableERA5Variables,

        ERA5Region,
        isinERA5Region, coordERA5Region, ERA5RegionGrid,

        era5Pressures

## Abstract SuperTypes
"""
    ERA5Module

Abstract supertype for ERA5 dataset modules.

All `ERA5Module` Types contain the following fields:
- `modID` : The module ID, that also acts as a prefix to filenames
- `eroot` : The specified directory in which to save the data
- `dtbeg` : The date for which downloads/analysis begins
- `dtend` : The date for which downloads/analysis finishes

!!! note
    The `ERA5MonthlyHour` subType also has the field `hours` that specifies the hour(s) of day for which monthly data is downloaded
"""
abstract type ERA5Module end

"""
    ERA5Variable

Abstract supertype for ERA5 variables.

All `ERA5Variable` Types contain the following fields:
- `varID` : The variable ID, that is also the identifier in the NetCDF files
- `lname` : The variable long-name, which is used to specify retrievals from CDS
- `vname` : The full-name of the variable
- `units` : The units of the variable
"""
abstract type ERA5Variable end

## ERA5Reanalysis.jl logging preface

modulelog() = "$(now()) - ERA5Reanalysis.jl"

function __init__()
    jfol = joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis"); mkpath(jfol);
    fvar = ["singlevariable.txt","singlecustom.txt","pressurevariable.txt","pressurecustom.txt"]
    for fname in fvar
        if !isfile(joinpath(jfol,fname))
            copyera5variables(fname)
            @info "$(modulelog()) - $(fname) does not exist in $(jfol), copying ..."
        end
    end
end

## Including other files in the module

include("module/module.jl")
include("module/show.jl")

include("variable/common.jl")
include("variable/single.jl")
include("variable/pressure.jl")
include("variable/show.jl")

include("region/region.jl")
include("region/functions.jl")

include("downloads/cdsapi.jl")
include("downloads/downloads.jl")
include("downloads/cdsretrieve.jl")
include("downloads/pythonprint.jl")

end


text:
include("util.jl")
include("locate.jl")

p = "StatPlots"
p = "Plots"
p = "PlotRecipes"

replace_d_pkg(p)

# pkgs = sort!(unique(mapreduce(Pkg.dependents, vcat, ["Plots", "RecipesBase"])))
#
# const PKGDB = joinpath(Pkg.dir("ReplaceD"), "pkgdb")
# ispath(PKGDB) || mkpath(PKGDB)
# pkgdbpath(p) = joinpath(PKGDB, p)
#
# for p in pkgs
#     target = pkgdbpath(p)
#     if !ispath(target)
#         cmd = `git clone $(url(p)) $target`
#         run(cmd)
#     end
# end
#
#
# for pkg in pkgs
#     root = pkgdbpath(pkg)
#     println(pkg)
#     for path in recipe_files(root)
#         # parse("begin $(readstring(path)) end")
#         println(" "^4, splitdir(path)[2])
#         search_token_in_recipes(tok -> string(tok) == "d", path)
#     end
# end


text:
# ========== Transforms On A Group Of Signals ==========
# ----- ACDWT on a set of signals -----
"""
    acdwtall(x, wt[, L])

Computes the autocorrelation discrete wavelet transform (ACDWT) on each slice of signal.

# Arguments
- `x::AbstractArray{T,2} where T<:Number`: Input signals, where each column corresponds to a
  signal.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T,3}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACDWT on all signals in x
xw = acdwtall(x, wt)
```

**See also:** [`acdwt`](@ref)
"""
function acdwtall(x::AbstractArray{T,2},
                  wt::OrthoFilter,
                  L::Integer = maxtransformlevels(x,1)) where T<:Number
    # Allocate space for transforms
    n,N = size(x)
    xw = Array{T,3}(undef, (n,L+1,N))
    # Dimension to slice
    dim_xw = 3
    dim_x = 2
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acdwt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacdwtall(xw[, wt])

Computes the inverse autocorrelation discrete wavelet transform (IACDWT) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T,3} where T<:Number`: ACDWT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
- `::Array{T,2}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACDWT on all signals in x
xw = acdwtall(x, wt)

# IACDWT on all signals in xw
x̂ = iacdwtall(xw)
```

**See also:** [`iacdwt`](@ref)
"""
function iacdwtall(xw::AbstractArray{T,3}, wt::Union{OrthoFilter,Nothing} = nothing) where
                   T<:Number
    # Allocate space for transforms
    n,_,N = size(xw)
    x = Array{T,2}(undef, (n,N))
    # Dimension to slice
    dim_x = 2
    dim_xw = 3
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacdwt!(xᵢ, xwᵢ)
        end
    end
    return x
end

# ----- ACWPT on a set of signals -----
"""
    acwptall(x, wt[, L])

Computes the autocorrelation wavelet packet transform (ACWPT) on each slice of signal.

# Arguments
- `x""AbstractArray{T,2} where T<:Number`: Input signals, where each column corresponds to a
  signal.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T,3}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPT on all signals in x
xw = acwptall(x, wt)
```

**See also:** [`acwpt`](@ref)
"""
function acwptall(x::AbstractArray{T,2},
                  wt::OrthoFilter,
                  L::Integer = maxtransformlevels(x,1)) where T<:Number
    # Allocate space for transforms
    n,N = size(x)
    xw = Array{T,3}(undef, (n,1<<L,N))
    # Dimension to slice
    dim_xw = 3
    dim_x = 2
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acwpt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacwptall(xw[, wt])

Computes the inverse autocorrelation wavelet packet transform (IACWPT) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T,3} where T<:Number`: ACWPT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.

# Returns
- `::Array{T,2}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPT on all signals in x
xw = acwptall(x, wt)

# IACWPT on all signals in xw
x̂ = iacwptall(xw)
```

**See also:** [`iacwpt`](@ref)
"""
function iacwptall(xw::AbstractArray{T,3}, wt::Union{OrthoFilter,Nothing} = nothing) where
                   T<:Number
    # Allocate space for transforms
    n,_,N = size(xw)
    x = Array{T,2}(undef, (n,N))
    # Dimension to slice
    dim_xw = 3
    dim_x = 2
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacwpt!(xᵢ, xwᵢ)
        end
    end
    return x
end

# ----- ACWPD on a set of signals -----
"""
    acwpdall(x, wt[, L])

Computes the autocorrelation wavelet packet decomposition (ACWPD) on each slice of signal.

# Arguments
- `x""AbstractArray{T,2} where T<:Number`: Input signals, where each column corresponds to a
  signal.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T,3}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPD on all signals in x
xw = acwpdall(x, wt)
```

**See also:** [`acwpd`](@ref)
"""
function acwpdall(x::AbstractArray{T,2},
                  wt::OrthoFilter,
                  L::Integer = maxtransformlevels(x,1)) where T<:Number
    # Allocate space for transforms
    n,N = size(x)
    xw = Array{T,3}(undef, (n,1<<(L+1)-1,N))
    # Dimension to slice
    dim_xw = 3
    dim_x = 2
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            acwpd!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iacwpdall(xw[, wt, L])
    iacwpdall(xw, L)
    iacwpdall(xw, wt, tree)
    iacwpdall(xw, tree)

Computes the inverse autocorrelation wavelet packet decomposition (IACWPD) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T,3} where T<:Number`: ACWPT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.
- `tree::BitVector`: Binary tree for inverse transform to be computed accordingly.

# Returns
- `::Array{T,2}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# ACWPD on all signals in x
xw = acwpdall(x, wt)

# IACWPD on all signals in xw
x̂ = iacwpdall(xw)
x̂ = iacwpdalll(xw, maketree(x))
x̂ = iacwpdall(xw, 5)
```

**See also:** [`iacwpd`](@ref)
"""
function iacwpdall(xw::AbstractArray{T,3},
                   wt::Union{OrthoFilter,Nothing} = nothing,
                   L::Integer = maxtransformlevels(x,1)) where T<:Number
    return iacwpdall(xw, L)
end

function iacwpdall(xw::AbstractArray{T,3}, L::Integer) where T<:Number
    return iacwpdall(xw, maketree(size(xw,1), L, :full))
end

function iacwpdall(xw::AbstractArray{T,3},
                   wt::Union{OrthoFilter, Nothing},
                   tree::BitVector) where T<:Number
    return iacwpdall(xw, tree)
end

function iacwpdall(xw::AbstractArray{T,3}, tree::BitVector) where T<:Number
    # Allocate space for transforms
    n,_,N = size(xw)
    x = Array{T,2}(undef, (n,N))
    # Dimension to slice
    dim_xw = 3
    dim_x = 2
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iacwpd!(xᵢ, xwᵢ)
        end
    end
    return x
end


text:
function choixInstance()
    println("\n Choix de la bibliothèque :")
    println(" --------------------------")
    println("  1) Agatz")
    println("  2) Poikonen")
    println("  3) Nous")
    println("  4) TSPLIB : all_tsp")
    println("  5) TSPLIB : vlsi")
    println("  6) TSPLIB : tnm")
    println("  7) Sur mesure")
    choix = choixListe(7)
    if choix == "1"
        instance = parserAgatz()
    elseif choix == "2"
        instance = parserPoikonen()
    elseif choix == "3"
        instance = parserNous()
    elseif choix == "4"
        instance = parserTSPLIB("all_tsp/")
    elseif choix == "5"
        instance = parserTSPLIB("vlsi/")
    elseif choix == "6"
        instance = parserTSPLIB("tnm/")
    elseif choix == "7"
        instance = instanceSurMesure()
    end
end

function test()
    # Parsing de l'instance
    f = open("instances/agatz/base/doublecenter-53-n10.txt")
    lignes = readlines(f) # Lecture de l'instance choisie
    vitesseCamion = parse(Float64, lignes[2])
    vitesseDrone = parse(Float64, lignes[4])
    nbPoints = parse(Int, lignes[6])
    listePoints = Vector{Point}(undef, nbPoints)
    xDepot, yDepot = parse.(Float64, split(lignes[8])[1:2])
    idDepot = "1"
    depot = Point(xDepot, yDepot, idDepot)
    listePoints[1] = depot
    compteur = 2
    for i in 10:length(lignes) # Pour tous les autres points
        xPoint, yPoint = parse.(Float64, split(lignes[i])[1:2])
        idPoint = string(compteur)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[compteur] = point
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(1, 0.5, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserAgatz()
    # Choix du type d'instance
    println("\n Choix du type d'instance :")
    println(" --------------------------")
    println("  1) Test ")
    println("  2) Uniforme ")
    println("  3) Simplement centrée ")
    println("  4) Doublement centrée ")
    choix1 = parse(Int64, choixListe(4))
    liste = ["base", "uniform", "singlecenter", "doublecenter"]
    listeInstances = readdir("instances/agatz/"*liste[choix1])

    # Choix de l'instance
    println("\n Liste des instances :")
    println(" ---------------------")
    for i in 1:length(listeInstances)
        println("   ",i,") ", listeInstances[i])
    end
    choix2 = parse(Int64, choixListe(length(listeInstances)))

    # Parsing de l'instance
    f = open("instances/agatz/"*liste[choix1]*"/"*listeInstances[choix2])
    lignes = readlines(f) # Lecture de l'instance choisie
    vitesseCamion = parse(Float64, lignes[2])
    vitesseDrone = parse(Float64, lignes[4])
    nbPoints = parse(Int, lignes[6])
    listePoints = Vector{Point}(undef, nbPoints)
    xDepot, yDepot = parse.(Float64, split(lignes[8])[1:2])
    idDepot = "1"
    depot = Point(xDepot, yDepot, idDepot)
    listePoints[1] = depot
    compteur = 2
    for i in 10:length(lignes) # Pour tous les autres points
        xPoint, yPoint = parse.(Float64, split(lignes[i])[1:2])
        idPoint = string(compteur)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[compteur] = point
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserPoikonen()
    # Choix de la taille de l'instance
    println("\n Choix de la taille de l'instance :")
    println(" ----------------------------------")
    println("  1) 10 noeuds ")
    println("  2) 15 noeuds ")
    println("  3) 20 noeuds ")
    println("  4) 30 noeuds ")
    println("  5) 40 noeuds ")
    println("  6) 50 noeuds ")
    println("  7) 60 noeuds ")
    println("  8) 70 noeuds ")
    println("  9) 80 noeuds ")
    println("  10) 90 noeuds ")
    println("  11) 100 noeuds ")
    println("  12) 200 noeuds ")
    choix = parse(Int64, choixListe(12))
    liste = ["10_points", "15_points", "20_points", "30_points",
             "40_points", "50_points", "60_points", "70_points",
             "80_points", "90_points", "100_points", "200_points",]

    # Choix de l'instance
    lignes = CSV.read("instances/poikonen/instances_"*liste[choix]*".csv", header=false) #openquotechar='[', closequotechar=']')
    println("\n Nombre d'instances : " * string(size(lignes,1)))
    choix = parse(Int64, choixListe(size(lignes,1)))

    # Parsing de l'instance
    temp = lignes[choix,2]
    temp = replace(temp, '['=>"")
    temp = split(temp[1:end-2], "],")
    nbPoints = length(temp)-1
    listePoints = Vector{Point}(undef, nbPoints)
    for i in 1:length(temp)-1
        couple = split(temp[i], ",")
        xPoint, yPoint = parse(Float64, couple[1]), parse(Float64, couple[2])
        idPoint = string(i)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[i] = point
    end
    D = creationDistancier(nbPoints, listePoints)

    # Paramètres à renseigner manuellement
    println("\n========== Paramètres d'instance ==========")
    print("\n --> Vitesse du camion : ")
    vitesseCamion = parse(Float64, readline())
    print(" --> Vitesse du drone : " )
    vitesseDrone = parse(Float64, readline())
    println("\n========== Paramètres d'instance ==========")

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserNous()
    # Choix de la taille de l'instance
    println("\n Choix de la taille de l'instance :")
    println(" ----------------------------------")
    println("  1) 10-100 noeuds")
    println("  2) 100-1000 noeuds")
    println("  3) Autres")
    choix1 = parse(Int64, choixListe(3))
    liste = ["10_100", "100_1000", "autres"]
    listeInstances = readdir("instances/nous/bibliotheque_"*liste[choix1])

    # Choix de l'instance
    println("\n Liste des instances :")
    println(" ---------------------")
    for i in 1:length(listeInstances)
        println("   ",i,") ", listeInstances[i])
    end
    choix2 = parse(Int64, choixListe(length(listeInstances)))

    # Parsing de l'instance
    f = open("instances/nous/bibliotheque_"*liste[choix1]*"/"*listeInstances[choix2])
    lignes = readlines(f) # Lecture de l'instance choisie
    nbPoints = parse(Int64, lignes[1])
    listePoints = Vector{Point}(undef, nbPoints)
    vitesseCamion, vitesseDrone = parse.(Float64, split(lignes[2]))
    compteur = 1
    for i in 3:length(lignes)
        xPoint, yPoint = parse.(Float64, split(lignes[i])[1:2])
        idPoint = string(compteur)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[compteur] = point
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserTSPLIB(librairie)
   # Choix de l'instance
    listeInstances = readdir("instances/tsplib/"*librairie)
    println("\n Liste des instances :")
    println(" ---------------------")
    for i in 1:length(listeInstances)
        println("   ",i,") ", listeInstances[i])
    end
    choix = parse(Int64, choixListe(length(listeInstances)))

    # Parsing de l'instance
    f = open("instances/tsplib/"*librairie*listeInstances[choix])
    lignes = readlines(f)
    i = 1
    while lignes[i] != "NODE_COORD_SECTION"
        # On cherche le début des coordonées
        i += 1
    end
    nbPoints = length(lignes) - (i + 1) # Il y a i lignes au début + une ligne avec le EOF à la fin

    listePoints = Vector{Point}(undef, nbPoints)
    compteur = 1
    for i in (i+1):length(lignes)-1
        xClient, yClient = parse.(Float64, split(lignes[i])[2:3])
        pointClient = Point(xClient, yClient, string(compteur))
        listePoints[compteur] = pointClient
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Paramètres à renseigner manuellement
    println("\n========== Paramètres d'instance ==========")
    print("\n --> Vitesse du camion : ")
    vitesseCamion = parse(Float64, readline())
    print(" --> Vitesse du drone : " )
    vitesseDrone = parse(Float64, readline())
    println("\n========== Paramètres d'instance ==========")

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

########## Parser rapides (non verbeux) pour faire des expérimentations ##########

function parserNousBis(nomInstance) # 0 verbosité
    # Parsing de l'instance
    f = open("experimentations/"*nomInstance)
    lignes = readlines(f) # Lecture de l'instance choisie
    nbPoints = parse(Int64, lignes[1])
    listePoints = Vector{Point}(undef, nbPoints)
    vitesseCamion, vitesseDrone = parse.(Float64, split(lignes[2]))
    compteur = 1
    for i in 3:length(lignes)
        xPoint, yPoint = parse.(Float64, split(lignes[i])[1:2])
        idPoint = string(compteur)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[compteur] = point
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserAgatzBis(nomInstance) # 0 verbosité
    # Parsing de l'instance
    f = open("experimentations/"*nomInstance)
    lignes = readlines(f) # Lecture de l'instance choisie
    vitesseCamion = parse(Float64, lignes[2])
    vitesseDrone = parse(Float64, lignes[4])
    nbPoints = parse(Int, lignes[6])
    listePoints = Vector{Point}(undef, nbPoints)
    xDepot, yDepot = parse.(Float64, split(lignes[8])[1:2])
    idDepot = "1"
    depot = Point(xDepot, yDepot, idDepot)
    listePoints[1] = depot
    compteur = 2
    for i in 10:length(lignes) # Pour tous les autres points
        xPoint, yPoint = parse.(Float64, split(lignes[i])[1:2])
        idPoint = string(compteur)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[compteur] = point
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(vitesseCamion, vitesseDrone, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserTsplibBis(nomInstance)
    # Parsing de l'instance
    f = open("experimentations/"*nomInstance)
    lignes = readlines(f)
    i = 1
    while lignes[i] != "NODE_COORD_SECTION"
        # On cherche le début des coordonées
        i += 1
    end
    nbPoints = length(lignes) - (i + 1) # Il y a i lignes au début + une ligne avec le EOF à la fin

    listePoints = Vector{Point}(undef, nbPoints)
    compteur = 1
    for i in (i+1):length(lignes)-1
        xClient, yClient = parse.(Float64, split(lignes[i])[2:3])
        pointClient = Point(xClient, yClient, string(compteur))
        listePoints[compteur] = pointClient
        compteur += 1
    end
    close(f)
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(1, 3, nbPoints, listePoints, D, solution, 0)
    return instance
end

function parserPoikoBis(lignes, cpt)
    # Parsing de l'instance
    temp = lignes[cpt,2]
    temp = replace(temp, '['=>"")
    temp = split(temp[1:end-2], "],")
    nbPoints = length(temp)-1
    listePoints = Vector{Point}(undef, nbPoints)
    for i in 1:length(temp)-1
        couple = split(temp[i], ",")
        xPoint, yPoint = parse(Float64, couple[1]), parse(Float64, couple[2])
        idPoint = string(i)
        point = Point(xPoint, yPoint, idPoint)
        listePoints[i] = point
    end
    D = creationDistancier(nbPoints, listePoints)

    # Création de l'instance
    solution = Solution([], [], 0, 0) # Solution vide pour le moment
    instance = Instance(1, 3, nbPoints, listePoints, D, solution, 0)
    return instance
end

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.Architectures: arch_array
using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization
using Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_vertically_integrated_volume_flux!
using Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_implicit_free_surface_right_hand_side!
using Oceananigans.Models.HydrostaticFreeSurfaceModels: pressure_correct_velocities!

@testset "Immersed boundaries test divergent flow solve with hydrostatic free surface models" begin
    @info "Testing immersed boundaries divergent flow solve"

    for arch in archs
        Nx = 11
        Ny = 11
        Nz = 1

        underlying_grid = RectilinearGrid(arch,
                                          size = (Nx, Ny, Nz),
                                          extent = (Nx, Ny, 1),
                                          topology = (Periodic, Periodic, Bounded))

        imm1=Int( floor((Nx+1)/2)   )
        imp1=Int( floor((Nx+1)/2)+1 )
        jmm1=Int( floor((Ny+1)/2)   )
        jmp1=Int( floor((Ny+1)/2)+1 )

        bottom = [-1. for i=1:Nx, j=1:Ny ]
        bottom[imm1-1:imp1+1, jmm1-1:jmp1+1] .= 0

        B = arch_array(arch, bottom)
        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(B))

        free_surfaces = [ImplicitFreeSurface(solver_method=:MatrixIterativeSolver, gravitational_acceleration=1.0),
                         ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient, gravitational_acceleration=1.0),
                         ImplicitFreeSurface(gravitational_acceleration=1.0)]

        sol = ()
        f = ()

        for free_surface in free_surfaces

            model = HydrostaticFreeSurfaceModel(grid = grid,
                                                free_surface = free_surface,
                                                tracer_advection = WENO5(),
                                                buoyancy = nothing,
                                                tracers = nothing,
                                                closure = nothing)

            # Now create a divergent flow field and solve for
            # pressure correction
            u, v, w     = model.velocities
            u[imm1,jmm1,1:Nz ] .=  1.
            u[imp1,jmm1,1:Nz ] .= -1.
            v[imm1,jmm1,1:Nz ] .=  1.
            v[imm1,jmp1,1:Nz ] .= -1.

              η = model.free_surface.η
              g = model.free_surface.gravitational_acceleration
            ∫ᶻQ = model.free_surface.barotropic_volume_flux
            rhs = model.free_surface.implicit_step_solver.right_hand_side
            solver = model.free_surface.implicit_step_solver;
                Δt = 1.
            compute_vertically_integrated_volume_flux!(∫ᶻQ, model)
            rhs_event = compute_implicit_free_surface_right_hand_side!(rhs, solver, g, Δt, ∫ᶻQ, η)
            wait(device(arch), rhs_event)

            solve!(η, solver, rhs, g, Δt)

            fill_halo_regions!(η, arch)

            # println(typeof(model))

            println("model.free_surface.gravitational_acceleration = ",model.free_surface.gravitational_acceleration)
            println("∫ᶻQ.u")
            show(stdout,"text/plain",∫ᶻQ.u.data[1:Nx,1:Ny,1])
            println("")

            println("η")
            show(stdout,"text/plain",η.data[1:Nx,1:Ny,1])
            println("")

            pressure_correct_velocities!(model, Δt)
            fill_halo_regions!(u, arch)

            println("u")
            show(stdout,"text/plain",u.data[1:Nx,1:Ny,1])
            println("")

            fs = model.free_surface
            vertically_integrated_lateral_areas = fs.implicit_step_solver.vertically_integrated_lateral_areas
            ∫Axᶠᶜᶜ = vertically_integrated_lateral_areas.xᶠᶜᶜ
            ∫Ayᶜᶠᶜ = vertically_integrated_lateral_areas.yᶜᶠᶜ

            sol = (sol..., model.free_surface.η)
            f  = (f..., model.free_surface)
        end

        @test all(interior(sol[1]) .≈ interior(sol[2]) .≈ interior(sol[3]))
    end
end


text:
# tests transformations from GeoData <=> Cartesian
using Test
using GeophysicalModelGenerator

# Create 3D volume with some fake data
Lon,Lat,Depth   =   LonLatDepthGrid(5:25,20:50,(-1300:100:0)km);
Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Depth*2 + Lon*km,LonData=Lon))

proj            =   ProjectionPoint(Lon=20,Lat=35)

# Convert this 3D dataset to a Cartesian dataset (the grid will not be orthogonal)
Data_set3D_Cart =   Convert2CartData(Data_set3D, proj)
@test sum(abs.(Data_set3D_Cart.x.val)) ≈ 5.293469089428514e6km

# Create Cartesian grid
X,Y,Z           =   XYZGrid(-500:100:500,-900:200:900,(-500:100:0)km);
Data_Cart       =   CartData(X,Y,Z,(Z=Z,))

# Project values of Data_set3D to the cartesian data
Data_Cart       =   ProjectCartData(Data_Cart, Data_set3D, proj)
@test sum(Data_Cart.fields.Depthdata) ≈ -316834.28716859705km
@test sum(Data_Cart.fields.LonData) ≈ 13165.712831402916

# Next, 3D surface (like topography)
Lon,Lat,Depth   =   LonLatDepthGrid(5:25,20:50,0);
Depth           =   cos.(Lon/5).*sin.(Lat)*10
Data_surf       =   GeoData(Lon,Lat,Depth,(Z=Depth,))
Data_surf_Cart  =   Convert2CartData(Data_surf, proj)

# Cartesian surface
X,Y,Z           =   XYZGrid(-500:10:500,-900:20:900,0);
Data_Cart       =   CartData(X,Y,Z,(Z=Z,))

Data_Cart       =   ProjectCartData(Data_Cart, Data_surf, proj)
@test sum(Data_Cart.z.val) ≈ 1858.2487019158766km
@test sum(Data_Cart.fields.Z) ≈ 1858.2487019158766

# Cartesian surface when UTM data is used
WE,SN,depth     =   XYZGrid(420000:1000:430000, 4510000:1000:4520000, 0);

Data_surfUTM    =   UTMData(WE, SN, depth, 33, true, (Depth = WE,));
Data_Cart       =   CartData(X,Y,Z,(Z=Z,))
Data_Cart       =   ProjectCartData(Data_Cart, Data_surfUTM, proj)

@test sum(Data_Cart.z.val) ≈ 0.0km
@test sum(Data_Cart.fields.Depth) ≈ 3.9046959539921126e9


text:
module LogGaussianCoxProcessWithLaplaceApproximation

import LambertW
import Optim
import LinearAlgebra
import SparseArrays
using LinearAlgebra: cholesky, cholesky!, logdet, norm, Symmetric
using LinearAlgebra: diag
using SpecialFunctions: logabsgamma
using PointPatternStatistics: PointPattern, area
#import FillArrays
import Distributions
#using GaussianRandomFields: CovarianceFunction, Matern, GaussianRandomField, CirculantEmbedding
import GaussianRandomFields

# For Poisson correction
using FFTW
using QuadGK: quadgk

struct PoissonProcess{WT}
    intensity::Matrix{Float64}
    window::WT
end

struct LogGaussianCoxProcess{MT} # with laplace approximation and Matern covariance
    # Parameters for matern covariance
    # mean field
    nu::Int
    range::Float64
    sd::Float64
    mean::MT # Matrix
end
LogGaussianCoxProcess(nu, range, sd, mean) = LogGaussianCoxProcess(nu, float(range), float(sd), mean)

# Could check the length of x here
Distributions.insupport(d::LogGaussianCoxProcess, x) = true

include("discretize.jl")
include("matern_precision.jl")
using .MaternQ

include("util.jl")
include("laplace_approximation.jl")
include("log_dens.jl")
include("poissoncorrection.jl")

export LogGaussianCoxProcess,
    Discretization,
    PoissonProcess

end


text:
elementsAbundances = JSON.parsefile(Base.Filesystem.joinpath(dirname(@__FILE__),"elements.json"))

function cart(x::Array,y::Array)
    #return [repmat(x,1,length(y))'[:] repmat(y,length(x),1)[:]]
    leny=length(y)
    lenx=length(x)
    m=leny*lenx
    OUT = zeros(Float64, m, 2)
    c=1
    for i = 1:lenx
        for j = 1:leny
            OUT[c,1] = x[i]
            OUT[c,2] = y[j]
            c+=1
        end
    end
    return OUT
end

function cartSum(x::Array,y::Array)
    return sum(cart(x, y), dims=2)
end

function cartProd(x::Array,y::Array)
    return prod(cart(x, y), dims=2)
end

function parseFormula(formula::AbstractString)
  # support for any kind of brackets
  formula = replace(formula, "[" => "(")
  formula = replace(formula, "]" => ")")
  formula = replace(formula, "{" => "(")
  formula = replace(formula, "}" => ")")
  # check for any existing brackets
  formula_no_brackets = formula
  while occursin(r"(.*)\((.*?)\)(\d*)", formula_no_brackets)
    # find rightmost left parenthesis
    formula_split = match(r"(.*)\((.*?)\)(\d*)", formula_no_brackets);
    leftpart = formula_split.captures[1];
    bracket_content = formula_split.captures[2];
    if (formula_split.captures[3] == "")
      bracket_content_mult = 1
    else
      bracket_content_mult = parse(Int,formula_split.captures[3])
    end
    rightpart = replace(formula_no_brackets, formula_split.match => "")
    formula_no_brackets = leftpart*repeat(bracket_content, bracket_content_mult)*rightpart
  end

  elementsRawIDs = collect((m.match for m = eachmatch(r"([A-Z][a-z]{0,1})(\d*)", string(formula_no_brackets))))
  #elementsRawIDs = matchall(r"([A-Z](?:[a-z])?)(\d*)", string(formula))

  elementsDict = Dict{String, Int}()
  for i in elementsRawIDs
    if occursin(r"[0-9]{1,}", string(i))
      elementID = match(r"[A-Za-z]*", i).match
      elementMult = parse(Int, match(r"\d*$", i).match)
      try
        elementsDict[elementID] = elementsDict[elementID] + elementMult
      catch
        elementsDict[elementID] = elementMult
      end
    else
      try
        elementsDict[i] = elementsDict[i] + 1
      catch
        elementsDict[i] = 1
      end
    end
  end
  return OrderedDict(elementsDict)
end

function parseCharge(charge::AbstractString)
  # support for any kind of brackets
  charge = replace(charge, "[" => "(")
  charge = replace(charge, "]" => ")")
  charge = replace(charge, "{" => "(")
  charge = replace(charge, "}" => ")")

  if occursin(r"(\d)\((\+*\-*)\)", charge)
    # find rightmost left parenthesis
    charge_split = match(r"(\d)\((\+*\-*)\)", charge);
    charge_multiple = charge_split.captures[1];
    charge_polarity = charge_split.captures[2];
  else
    charge_split = match(r"(\d)*(\+*\-*)", charge);
    charge_multiple = charge_split.captures[1];
    charge_polarity = charge_split.captures[2];
  end

  return (charge_multiple, charge_polarity)
end

function convoluteDict(sumForm::OrderedDict)
  cD = Array{String,1}()
  for i in keys(sumForm)
    for j in 1:sumForm[i]
      push!(cD, i)
    end
  end
  return cD
end

function elementsStringify!(stringBase::AbstractString, a::AbstractString, b::Number)
  if b == 1
    stringBase = stringBase*string(a)
  else
    stringBase = stringBase*string(a)*string(b)
  end
  stringBase
end

function sumFormula(sumFormDict::OrderedDict)
  sf = ""
  sf_c_prefix = ""
  sf_h_prefix = ""
  for i in keys(sumFormDict)
    if string(i) == "C"
      sf_c_prefix = elementsStringify!(sf_c_prefix, "C", sumFormDict[i])
    elseif string(i) == "H"
      sf_h_prefix = elementsStringify!(sf_h_prefix, "H", sumFormDict[i])
    else
      sf = elementsStringify!(sf, string(i), sumFormDict[i])
    end
  end
  return sf_c_prefix*sf_h_prefix*sf
end

"""
Usage: isotopicPattern(formula::AbstractString; adduct=\"H\", charge=\"+\", abundanceLimit = 1e-5, niceOutput = true)")
"""
function isotopicPattern(formula::String; kwargs...)
    kwargs = Dict(kwargs)
    if haskey(kwargs, :adduct)
      adduct = kwargs[:adduct]
    else
      adduct = ""
    end

    if haskey(kwargs, :charge)
      charge = kwargs[:charge]
    else
      charge = ""
    end

    if haskey(kwargs, :abundanceLimit)
      abundanceLimit = kwargs[:abundanceLimit]
    else
      abundanceLimit = 1e-5
    end

    if haskey(kwargs, :niceOutput)
      niceOutput = kwargs[:niceOutput]
    else
      niceOutput = true
    end

    if haskey(kwargs, :result)
      exportBool = kwargs[:result]
      if exportBool
        niceOutput = false
      end
    else
      exportBool = false
    end

    start =  Base.time()
    elementsDict = parseFormula(formula)
    composition = convoluteDict(elementsDict)
    #println(composition)

    if length(adduct) > 0
        composition = [composition; convoluteDict(parseFormula(adduct))]
    end

    elem = string(composition[1])
    masses = elementsAbundances[elem]["Relative Atomic Mass"]
    abund = elementsAbundances[elem]["Isotopic Composition"]

    for i in 2:length(composition)
        elemNext = string(composition[i])
        massesNext = elementsAbundances[elemNext]["Relative Atomic Mass"]
        abundNext = elementsAbundances[elemNext]["Isotopic Composition"]

        massesTemp = round.(cartSum(masses, massesNext), digits = 6)
        abundTemp = cartProd(abund, abundNext)

        masses = unique(massesTemp)

        abund = Array{Float64}(undef, size(masses))
        for j in 1:length(masses)
            index = findall((in)(masses[j]), massesTemp)
            abund[j] = sum(abundTemp[index])
        end

        toKeep = findall(abund.>abundanceLimit)
        masses = masses[toKeep]
        abund = abund[toKeep]
    end
    #println(abund)

    SF = sumFormula(elementsDict)
    SF_add = ""
    if length(adduct) > 0
      SF_add = ".$(sumFormula(parseFormula(adduct)))"
    end

    if length(charge) > 0
      (charge_multiple, charge_polarity) = parseCharge(charge)

      if charge_polarity in ["+", "-"]
        try
          parse(Int,charge_multiple)
        catch
          charge_multiple = "1"
        end
        for i in 1:parse(Int,charge_multiple)
          masses = [eval(Meta.parse(charge_polarity))(masses[j], -1*elementsAbundances["e"]["Relative Atomic Mass"][1]) for j in 1:length(masses)]
        end
        if parse(Int,charge_multiple) > 1
          masses = masses/parse(Int,charge_multiple)
        end
      end
    else
      charge_multiple = ""
      charge_polarity = ""
    end

    charge_multiple_label = ""
    if charge_multiple != "1"
      try
        charge_multiple_label = Base.REPLCompletions.latex_symbols["\\^$(charge_multiple)"]
      catch
        charge_multiple_label = charge_multiple
      end
    else
      charge_multiple_label = ""
    end

    charge_polarity_label = ""
    try
      charge_polarity_label = Base.REPLCompletions.latex_symbols["\\^$(charge_polarity)"]
    catch
      charge_polarity_label = charge_polarity
    end

    SFx = SF*SF_add*charge_multiple_label*charge_polarity_label

    if niceOutput
      println("\nFormula: \e[1;32m $SFx")
    end

    IPCresult = Dict( SFx => zeros(Float64, length(masses), 2))
    #println(result)
    IPCresult[SFx][:,2] = abund/maximum(abund)
    IPCresult[SFx][:,1] = masses

    IPCresult[SFx] = IPCresult[SFx][sortperm(IPCresult[SFx][:,2], rev=true),:]

    IPCresult[SFx] = IPCresult[SFx][findall(IPCresult[SFx][:,2].>abundanceLimit),:]
    stop = time()
    #println("\e[1;37mevaluated in ", round(stop-start,3), " seconds." )
    if niceOutput
      println("\e[m-------------------------------")
      println("Mass [amu]\tAbbundance [%]")
      println("-------------------------------")
      for i in 1:size(IPCresult[SFx])[1]
          println(string(round(IPCresult[SFx][i,1], digits=6))*"\t"*string(round(IPCresult[SFx][i,2]*100,digits=6)))
        end
        println("Found ", length(masses), " isotopic masses for $abundanceLimit abundance limit.")
        println()
      end
    # if exportBool
    #   return IPCresult[SFx]
    # end
    return IPCresult[SFx]
end;

isopat(x) = isotopicPattern(x)

isotopicPatternProtonated(x) = isotopicPattern(x, adduct="H", charge="+")
isopatprot(x) = isotopicPatternProtonated(x)

monoisotopicMass(x) = isotopicPattern(x;niceOutput=false)[1,1]
monoisomass(x) = monoisotopicMass(x)


text:
using Base.Test

abstract CurveContainer

type StaticContainer <: CurveContainer
	parameters_id::Vector{Int} # for interpolation methods, id stores days_to_maturity on curve's daycount convention.
	parameters_values::Vector{Float64}

	StaticContainer(parameters_id::Vector{Int}, parameters_values::Vector{Float64}) = begin

		if length(parameters_id) != length(parameters_values)
			error("parameters_id and parameters_values must have the same length")
		end

		if !issorted(parameters_id)
			error("parameters_id should be sorted before creating IRCurve instance")
		end

		new(parameters_id, parameters_values)
	end
end

type DictAnyContainer <: CurveContainer
	parameters::Dict{Any, Any}

	DictAnyContainer(parameters_id::Vector{Int}, parameters_values::Vector{Float64}) = begin

		if length(parameters_id) != length(parameters_values)
			error("parameters_id and parameters_values must have the same length")
		end

		if !issorted(parameters_id)
			error("parameters_id should be sorted before creating IRCurve instance")
		end

		d = Dict{Any, Any}()
		d["days_to_maturity"] = parameters_id
		d["zero_rates"] = parameters_values

		new(d)
	end
end

type DictSymbolContainer <: CurveContainer
	parameters::Dict{Symbol, Any}

	DictSymbolContainer(parameters_id::Vector{Int}, parameters_values::Vector{Float64}) = begin

		if length(parameters_id) != length(parameters_values)
			error("parameters_id and parameters_values must have the same length")
		end

		if !issorted(parameters_id)
			error("parameters_id should be sorted before creating IRCurve instance")
		end

		d = Dict{Symbol, Any}()
		d[:days_to_maturity] = parameters_id
		d[:zero_rates] = parameters_values

		new(d)
	end
end

get_dtm(c::DictSymbolContainer) = c.parameters[:days_to_maturity]::Vector{Int}
get_zero_rates(c::DictSymbolContainer) = c.parameters[:zero_rates]::Vector{Float64}

get_dtm(c::DictAnyContainer) = c.parameters["days_to_maturity"]::Vector{Int}
get_zero_rates(c::DictAnyContainer) = c.parameters["zero_rates"]::Vector{Float64}

get_dtm(c::StaticContainer) = c.parameters_id
get_zero_rates(c::StaticContainer) = c.parameters_values

#Xa, Ya, Xb, Yb = _interpolationpoints(curve.parameters_id, curve.parameters_values, x_out)
function _interpolationpoints(x::Vector{Int}, y::Vector{Float64}, x_out::Int)
	index_a::Int
	index_b::Int

	if x_out <= x[1]
		# Interpolation point is before first vertice
		# Slope will be determined by the 1st and 2nd vertices
		index_a = 1
		index_b = 2
	elseif x_out >= x[end]
		# Interpolation point is after last vertice
		# Slope will be determined by the last and last-1 vertices
		index_b = length(x)
		index_a = index_b - 1
	else
		# Inner point
		index_a =  findlast(a -> a < x_out, x) # last element before x_out on x
		index_b = findfirst(b -> b >=  x_out, x) # first element after x_out on x
	end

	Xa = x[index_a]
	Xb = x[index_b]
	Ya = y[index_a]
	Yb = y[index_b]

	return Xa, Ya, Xb, Yb
end

function run(c::CurveContainer, x_out_vec::Vector{Int})
	for i = 1:100000
		for x_out in x_out_vec
			index_a, index_b = _interpolationpoints(get_dtm(c), get_zero_rates(c), x_out)
		end
	end
end

const vert_x = [11, 15, 19, 23]
const vert_y = [0.10, 0.15, 0.20, 0.19]
const x_out_v = collect(1:100)

curve_static = StaticContainer(vert_x, vert_y)
curve_dict_any = DictAnyContainer(vert_x, vert_y)
curve_dict_sym = DictSymbolContainer(vert_x, vert_y)

run(curve_static, x_out_v)
println("static")
@time run(curve_static, x_out_v)

run(curve_dict_any, x_out_v)
println("dict any")
@time run(curve_dict_any, x_out_v)

run(curve_dict_sym, x_out_v)
println("dict sym")
@time run(curve_dict_sym, x_out_v)


text:
# Embed distance matrix
function embedmatrix1(x, delay, metric="max")
    dist = getmetric(metric)
    fun = Dict(Euclidean()=>+, Chebyshev()=>max)[dist]
    nr, nc = size(x)
    x = abs2(x)
    sqrt(fun(x[1:nr-delay, 1:nc-delay], x[1+delay:nr, 1+delay:nc]))
end

"""
    fnn(x, mbounds, delay, thresholds; metric="max")

Calculate the number of false nearest neighbours (FNN) by Kennel's algorithm.

FNN are calculated for the embedding dimensions between `mbounds[1]` and
`mbounds[2]`, and for the given delay, taking `thresholds[1]` as the parameter
Rtol, and `thresholds[2]` as the parameter Atol. The argument `metric` defines
what kind of distance is calculated between pairs of points, as in
`distancematrix`.
"""
function fnn(x, mbounds, delay, thresholds; metric="max")
    Rtol2 = thresholds[1]^2
    Atol = thresholds[2]
    Ra = std(x)
    m1, m2 = mbounds
    nfnn = zeros(m2-m1+1)
    # DM - with maximum value in the diagonal to rule it out in neighbour search
    # Start with matrix embedded at m=m1
    dm = distancematrix(embed(x, m1, delay), metric)
    n = size(dm)[1]
    dm[1:n+1:end] = maximum(dm)
    for m = 1:length(nfnn)
        # Look for nearest neighbours within the first n-delay
        n -= delay
        nnval, nnpos = findmin(dm[1:n,1:n], 2)
        dm = embedmatrix1(dm, delay, metric)
        nnv1 = nnval[1:n]
        nnv2 = dm[1:n,1:n][nnpos]
        fnn1 = ( ((nnv2.^2)./(nnv1.^2)-1) .> Rtol2 )
        fnn2 = ( (nnv2/Ra) .> Atol )
        nfnn[m] = sum(fnn1 | fnn2)
    end
    nfnn
end

"""
    afnn(x, mbounds, delay; metric="max")

Calculate the ratios E1 and E2 of Cao's method of "averaged false nearest neighbours".

The ratios are calculated for the embedding dimensions between `mbounds[1]` and
`mbounds[2]`, and for the given delay. The argument `metric` defines
what kind of distance is calculated between pairs of points, as in
`distancematrix`.
"""

function afnn(x, mbounds, delay; metric="max")
    m1, m2 = mbounds
    dm = distancematrix(embed(x, m1, delay), metric)
    n = size(dm)[1]
    dm[1:n+1:end] = maximum(dm)
    mean_ratio = zeros(m2-m1+2)
    mean_increment = zeros(m2-m1+2)
    for m = 1:length(mean_ratio)
        n -= delay
        nnval, nnpos = findmin(dm[1:n,1:n], 2)
        dm = embedmatrix1(dm, delay, metric)
        nnv1 = nnval[1:n]
        nnv2 = dm[1:n,1:n][nnpos]
        mean_ratio[m] = mean(nnv2./nnv1)
        # Project nnpos in original series
        nnx = div(nnpos, n) + 1
        d = (m1+m-1)*delay
        mean_increment[m] = mean(abs(x[(1:n)+d]-x[nnx+d]))
    end
    e1 = mean_ratio[2:end]./mean_ratio[1:end-1]
    e2 = mean_increment[2:end]./mean_increment[1:end-1]
    e1, e2
end

"""
    ffnn(x, mbounds, delay; metric="max")

Calculate the ratio of false first nearest neighbours (FFNN) by Krakovská's algorithm.

The ratios are calculated for the embedding dimensions between `mbounds[1]` and
`mbounds[2]`, and for the given delay. The argument `metric` defines
what kind of distance is calculated between pairs of points, as in
`distancematrix`.
"""
function ffnn(x, mbounds, delay; metric="max")
    m1, m2 = mbounds
    dm = distancematrix(embed(x, m1, delay), metric)
    n = size(x)[1]
    dm[1:n+1:end] = maximum(dm)
    fnn_ratio = zeros(m2-m1+1)
    for m = 1:length(fnn_ratio)
        n -= delay
        nnval, nnpos1 = findmin(dm[1:n,1:n], 2)
        dm = embedmatrix1(dm, delay, metric)
        nnval, nnpos2 = findmin(dm, 2)
        fnn_ratio[m] = sum(nnpos1 .!= nnpos2)/n
    end
    fnn_ratio
end


text:
using OMEinsum: StaticEinCode, DynamicEinCode
using OMEinsum.ContractionOrder: IncidenceList, ContractionTree
export merge_vectors, embed_simplifier
export merge_greedy

struct NetworkSimplifier
    operations::Vector{NestedEinsum}
end

function merge_vectors(code::EinCode)
    ET = code isa DynamicEinCode ? typeof(code) : StaticEinCode
    ixs = OMEinsum.getixs(code)
    mask = trues(length(ixs))
    ops = [NestedEinsum{ET}(i) for i=1:length(ixs)]
    for i in 1:length(ixs)
        if length(ixs[i]) == 1
            for j in 1:length(ixs)
                if i!=j && mask[j] && ixs[i][1] ∈ ixs[j]  # merge i to j
                    mask[i] = false
                    ops[j] = NestedEinsum([ops[i], ops[j]],
                        _similar(code, (ixs[i], ixs[j]), ixs[j]))
                    break
                end
            end
        end
    end
    newcode = _similar(code, ixs[mask], OMEinsum.getiy(code))
    return NetworkSimplifier(ops[mask]), newcode
end
_similar(::DynamicEinCode, ixs, iy) = DynamicEinCode(collect(collect.(ixs)), collect(iy))
_similar(::StaticEinCode, ixs, iy) = StaticEinCode{(Tuple.(ixs)...,), (iy...,)}()

function merge_greedy(code::EinCode, size_dict; threshhold=-1e-12)
    ixs, iy, L = getixsv(code), getiyv(code), OMEinsum.labeltype(code)
    ET = code isa DynamicEinCode ? typeof(code) : StaticEinCode
    log2_edge_sizes = Dict{L,Float64}()
    for (k, v) in size_dict
        log2_edge_sizes[k] = log2(v)
    end
    incidence_list = IncidenceList(Dict([i=>ixs[i] for i=1:length(ixs)]); openedges=iy)
    n = ContractionOrder.nv(incidence_list)
    if n == 0
        return nothing
    elseif n == 1
        return collect(ContractionOrder.vertices(incidence_list))[1]
    end
    tree = Dict{Int,NestedEinsum}([v=>NestedEinsum{ET}(v) for v in ContractionOrder.vertices(incidence_list)])
    cost_values = ContractionOrder.evaluate_costs(MinSpaceDiff(), incidence_list, log2_edge_sizes)
    while true
        if length(cost_values) == 0
            return _buildsimplifier(code, tree, incidence_list)
        end
        v, pair = findmin(cost_values)
        if v <= threshhold
            _, _, c = ContractionOrder.contract_pair!(incidence_list, pair..., log2_edge_sizes)
            tree[pair[1]] = NestedEinsum((tree[pair[1]], tree[pair[2]]), _similar(code, c.first, c.second))
            if ContractionOrder.nv(incidence_list) <= 1
                return _buildsimplifier(code, tree, incidence_list)
            end
            ContractionOrder.update_costs!(cost_values, pair..., MinSpaceDiff(), incidence_list, log2_edge_sizes)
        else
            return _buildsimplifier(code, tree, incidence_list)
        end
    end
end
function _buildsimplifier(code, tree, incidence_list)
    vertices = sort!(collect(keys(incidence_list.v2e)))
    ixs = [incidence_list.v2e[v] for v in vertices]
    iy = incidence_list.openedges
    NetworkSimplifier([tree[v] for v in vertices]), _similar(code, ixs, iy)
end

function apply_simplifier(s::NetworkSimplifier, xs)
    map(s.operations) do op
        return op(xs...)
    end
end
(s::NetworkSimplifier)(xs...) = apply_simplifier(s, xs)

function embed_simplifier(code::NestedEinsum, simplifier)
    if OMEinsum.isleaf(code)
        op = simplifier.operations[code.tensorindex]
        return op
    else
        return NestedEinsum(map(code.args) do arg
            embed_simplifier(arg, simplifier)
        end, code.eins)
    end
end

embed_simplifier(code::SlicedEinsum, simplifier) = SlicedEinsum(code.slicing, embed_simplifier(code.eins, simplifier))

export optimize_permute
optimize_permute(se::SlicedEinsum, level=0) = SlicedEinsum(se.slicing, se.eins isa EinCode ? se.eins : optimize_permute(se.eins, level))
function optimize_permute(ne::NestedEinsum, level=0)
    if isleaf(ne)
        return ne
    else
        args = optimize_permute.(ne.args, level+1)
        ixs0 = getixsv(ne.eins)
        ixs = [isleaf(x) ? ixs0[i] : getiyv(x.eins) for (i, x) in enumerate(args)]
        iy = level == 0 ? getiyv(ne.eins) : optimize_output_permute(ixs, getiyv(ne.eins))
        return NestedEinsum(args, _similar(ne.eins, ixs, iy))
    end
end

function optimize_output_permute(ixs::AbstractVector{<:AbstractVector{LT}}, iy::AbstractVector{LT}) where LT
    if length(ixs) != 2
        return iy
    else
        iA, iB = ixs
        batchdim = LT[]
        outerA = LT[]
        outerB = LT[]
        bcastdim = LT[]
        for l in iy
            if l ∈ iA
                if l ∈ iB
                    push!(batchdim, l)
                else
                    push!(outerA, l)
                end
            else
                if l ∈ iB
                    push!(outerB, l)
                else
                    push!(bcastdim, l)
                end
            end
        end
        return vcat(
            sort!(outerA, by=l->findfirst(==(l), iA)),
            sort!(outerB, by=l->findfirst(==(l), iB)),
            sort!(batchdim, by=l->findfirst(==(l), iA)),
            bcastdim)
    end
end


text:
using LoopVectorization, Test
# T = Float32

macro outline(ex)
  quote
    (() -> begin
       Base.@_noinline_meta
       $(esc(ex))
     end)()
  end
end

function test_broadcast(::Type{T}) where {T}
  M, N = 37, 47
  @show T, @__LINE__
  R = T <: Integer ? (T(-100):T(100)) : T
  a = rand(R,99,99,99);
  b = rand(R,99,99,1);
  bl = LowDimArray{(true,true,false)}(b);
  @test size(bl) == size(b)
  @test LoopVectorization.ArrayInterface.size(bl) === (size(b,1),size(b,2),LoopVectorization.StaticInt(1))

  br = reshape(b, (99,99));
  c1 = a .+ b;
  c2 = @outline @turbo a .+ bl;
  @test c1 ≈ c2
  fill!(c2, 99999); @outline @turbo c2 .= a .+ br;
  @test c1 ≈ c2
  fill!(c2, 99999); @outline @turbo c2 .= a .+ b;
  @test c1 ≈ c2
  br = reshape(b, (99,1,99));
  bl = LowDimArray{(true,false,true)}(br);
  @test size(bl) == size(br)
  @test LoopVectorization.ArrayInterface.size(bl) === (size(br,1),LoopVectorization.StaticInt(1),size(br,3))
  @. c1 = a + br;
  fill!(c2, 99999); @outline @turbo @. c2 = a + bl;
  @test c1 ≈ c2
  fill!(c2, 99999); @outline @turbo @. c2 = a + br;
  @test c1 ≈ c2
  br = reshape(b, (1,99,99));
  bl = LowDimArray{(false,)}(br);
  @test size(bl) == size(br)
  @test LoopVectorization.ArrayInterface.size(bl) === (LoopVectorization.StaticInt(1),size(br,2),size(br,3))
  @. c1 = a + br;
  fill!(c2, 99999);
  @test c1 ≈ @outline @turbo @. c2 = a + bl
  # @test c1 ≈ c2
  br = reshape(rand(R,99), (1,99,1));
  bl = LowDimArray{(false,)}(br);
  @test size(bl) == size(br)
  @. c1 = a + br;
  fill!(c2, 99999);
  @outline @turbo @. c2 = a + bl;
  @test c1 ≈ c2

  if T <: Integer
    xs = rand(-T(100):T(100), M);
  else
    xs = rand(T, M);
  end
  max_ = maximum(xs, dims=1);
  @test (@outline @turbo exp.(xs .- LowDimArray{(false,)}(max_))) ≈ exp.(xs .- LowDimArray{(false,)}(max_))
  @test size(LowDimArray{(false,)}(max_)) == size(max_)

  if T === Int32
    a = rand(T(1):T(100), 73, 1);
    @test sqrt.(Float32.(a)) ≈ @outline @turbo sqrt.(a)
  elseif T === Int64
    a = rand(T(1):T(100), 73, 1);
    @test sqrt.(a) ≈ @outline @tturbo sqrt.(a)
  else
    a = rand(T, 73, 1);
    @test sqrt.(a) ≈ @outline @turbo sqrt.(a)
  end

  a = rand(R, M); B = rand(R, M, N); c = rand(R, N); c′ = c';
  d1 =       @. a + B * c′;
  d2 = @outline @tturbo @. a + B * c′;
  @test d1 ≈ d2

  @.      d1 = a + B * c′;
  @outline @turbo @. d2 = a + B * c′;
  @test d1 ≈ d2

  d3 = a .+ B * c;
  d4 = @outline @turbo a .+ B *ˡ c;
  @test d3 ≈ d4

  fill!(d3, -1000.0);
  fill!(d4, 91000.0);

  d3 .= a .+ B * c;
  @outline @turbo d4 .= a .+ B *ˡ c;
  @test d3 ≈ d4

  fill!(d4, 91000.0);
  @outline @turbo @. d4 = a + B *ˡ c;
  @test d3 ≈ d4

  M, K, N = 77, 83, 57;
  A = rand(R,M,K); B = rand(R,K,N); C = rand(R,M,N);
  At = copy(A');
  D1 = C .+ A * B;
  D2 = @outline @tturbo C .+ A .*ˡ B;
  @test D1 ≈ D2
  if RUN_SLOW_TESTS
    fill!(D2, -999999); D2 = @outline @turbo C .+ At' *ˡ B;
    @test D1 ≈ D2
    fill!(D2, -999999); @test A * B ≈ (@outline @turbo @. D2 = A *ˡ B)
    D1 .= view(C, 1, :)' .+ A * B;
    fill!(D2, -999999);
    @outline @turbo D2 .= view(C, 1, :)' .+ A .*ˡ B;
    @test D1 ≈ D2
    C3d = rand(R,3,M,N);
    D1 .= view(C3d, 1, :, :) .+ A * B;
    fill!(D2, -999999);
    @outline @tturbo D2 .= view(C3d, 1, :, :) .+ A .*ˡ B;
    @test D1 ≈ D2
  end
  D1 .= 9999;
  @outline @turbo D2 .= 9999;
  @test D1 == D2
  D1 .= -99999;
  @outline @tturbo D2' .= -99999;
  @test D1 == D2

  b = rand(T,K); x = rand(R,N);
  D1 .= C .+ A * (b .+ x');
  @outline @tturbo @. D2 = C + A *ˡ (b + x');
  @test D1 ≈ D2
  D2 = @outline @turbo @. C + A *ˡ (b + x');
  @test D1 ≈ D2
  if T === Int64
    xd = rand(-1_000_000_000_000:1_000_000_000_000,89);
  elseif T === Int32
    xd = rand(-Int32(10_000_000):Int32(10_000_000),89);
  else
    xd = rand(R,89);
  end;
  yd = rand(R,89);
  yd[yd .== 0] .= 77;
  @test xd .÷ yd == @turbo xd .÷ yd

  if T <: Union{Float32,Float64}
    D3 = cos.(B');
    D4 = @outline @turbo cos.(B');
    @test D3 ≈ D4

    fill!(D3, -1e3); fill!(D4, 9e9);
    Bt = transpose(B);
    @. D3 = exp(Bt);
    @outline @tturbo @. D4 = exp(Bt);
    @test D3 ≈ D4

    D1 = similar(B); D2 = similar(B);
    D1t = transpose(D1);
    D2t = transpose(D2);
    @. D1t = exp(Bt);
    @outline @turbo @. D2t = exp(Bt);
    @test D1t ≈ D2t

    fill!(D1, -1e3);
    fill!(D2, 9e9);
    @. D1' = exp(Bt);
    lset = @outline @tturbo @. D2' = exp(Bt);

    @test D1 ≈ D2

    a = rand(137);
    b1 = @outline @turbo @. 3*a + sin(a) + sqrt(a);
    b2 =      @. 3*a + sin(a) + sqrt(a);
    @test b1 ≈ b2
    three = 3; fill!(b1, -9999);
    @outline @tturbo @. b1 = three*a + sin(a) + sqrt(a);
    @test b1 ≈ b2

    C = rand(100,10,10);
    D1 = C .^ 0.3;
    D2 = @outline @tturbo C .^ 0.3;
    @test D1 ≈ D2
    @. D1 = C ^ 2;
    @outline @turbo @. D2 = C ^ 2;
    @test D1 ≈ D2
    @turbo view(C,1:100,1:10,1:10) .= 0;
    @test all(==(0), C)
    @turbo view(B, axes(B,1), axes(B,2))' .= 2
    @test all(==(2), B)
  end
end

@testset "broadcast_float" begin
  for T ∈ (Float32, Float64, Int32, Int64)
    @time test_broadcast(T)
  end
end