text:
# ---------------------------------------------------------------------------
# density

@recipe function f(::Type{Val{:density}}, x, y, z; trim=false)
    newx, newy = violin_coords(y, trim=trim)
    if Plots.isvertical(d)
        newx, newy = newy, newx
    end
    x := newx
    y := newy
    seriestype := :path
    ()
end
Plots.@deps density path

# ---------------------------------------------------------------------------
# cumulative density

@recipe function f(::Type{Val{:cdensity}}, x, y, z; trim=false,
                   npoints = 200)
    newx, newy = violin_coords(y, trim=trim)

    if Plots.isvertical(d)
        newx, newy = newy, newx
    end

    newy = [sum(newy[1:i]) for i = 1:length(newy)] / sum(newy)

    x := newx
    y := newy
    seriestype := :path
    ()
end
Plots.@deps cdensity path


text:
function make_gui_bootloader(cur_module::Module, cur_id::UUID)
  comm_id = "interact-$( string(cur_id) )"

  display(MIME("text/javascript"), """

    var _guiBootup = function (curCallback) {
      Jupyter.notebook.kernel.comm_manager.unregister_target("$(comm_id)")
      Jupyter.notebook.kernel.comm_manager.register_target("$(comm_id)", function (comm, msg) {
        globalComms["$(comm_id)"] = comm;

        comm.on_msg(function(msg) {
          if ( "json" in msg.content.data ) {
            var workPlot = \$("#js-interact__$( string(cur_id) ) .js-plotly-plot")[0];
            if ( typeof workPlot === "undefined" ) { alert("Could not find plotly plot!"); }

            var workJson = msg.content.data.json.plot;
            customPlotlyReact(workPlot, workJson.data, workJson.layout, workJson.config);

            \$("#js-interact__$( string(cur_id) ) .js-new-plot-script").replaceWith(msg.content.data.json.script);
          } else {
            if ( "html" in msg.content.data ) {
              var tmpDisplay = \$("#js-interact__$( string(cur_id) ) .js-display");
              tmpDisplay.html(msg.content.data.html);
            } else {
              if ( "text" in msg.content.data ) {
                console.log(msg.content.data.text);
              } else {
                alert("Unrecognized interact content type!");
              }
            }
          }
        });

        comm.on_close(function(msg) { console.log("Julia close message: " + msg); });
        \$("#js-interact__$( string(cur_id) )").trigger("interact");
      })

      var curCommand = '$(string(cur_module)).SimplePlots.interact_comms["$(comm_id)"][] += 10';
      Jupyter.notebook.kernel.execute(curCommand, {"iopub": {"output": function(tmpMessage) {
        console.log(tmpMessage)
      }}});

      if ( typeof curCallback !== "undefined" ) { curCallback(); }
    }

    if ( typeof Jupyter === "undefined" ) {
      customPlotLoader(function() {
        \$("#js-interact__$( string(cur_id) )").trigger("interact");
      });
    } else {
      if (Jupyter.notebook.kernel) {
        _guiBootup()
      } else {
        Jupyter.notebook.events.one('kernel_ready.Kernel', (e) => {
          _guiBootup(function() {
            customPlotLoader(function() {
              \$("#js-interact__$( string(cur_id) )").trigger("interact");
            });
          });
        });
      }
    }
  """)
end

ext: jl


text:
using ApplicationBuilder

using Compat
using Compat.Pkg

examples_blink = joinpath(@__DIR__, "..", "..", "examples", "blink.jl")

# Allow this file to be called either as a standalone file to build the above
# example, or from runtests.jl using a provided builddir.
Compat.isdefined(:builddir) || (builddir="builddir")

blinkPkg = Pkg.dir("Blink")
httpParserPkg = Pkg.dir("HttpParser")
mbedTLSPkg = Pkg.dir("MbedTLS")

@assert blinkPkg != nothing "Blink is not installed!"

using Blink

ApplicationBuilder.build_app_bundle(examples_blink;
    verbose = true,
    resources = [joinpath(blinkPkg, "deps","Julia.app"),
                 joinpath(blinkPkg, "src","AtomShell","main.js"),
                 joinpath(blinkPkg, "src","content","main.html"),
                 joinpath(blinkPkg, "res")],
    # Get the current library names directly from the packages that use them,
    # which keeps this build script robust against lib version changes.
    libraries = [HttpParser.lib,
                 MbedTLS.libmbedcrypto],
    appname="HelloBlink", builddir=builddir)


text:
#=
Récupère le coût d'un circuit
=#
function getCircuitCost(circuit, costs)

	cost = 0
	for (from, to) in zip(circuit[1:end-1], circuit[2:end])
		cost += costs[(from, to)]
	end

	cost += costs[circuit[1], circuit[end]]

	return cost

end

#=
Récupère le coût d'une liste de circuits
=#
function getCircuitsCost(circuits, costs)

	cost = 0

	for circuit in circuits
		cost += getCircuitCost(circuit, costs)
	end

	return cost

end

#=
Convertie le résultat d'un VRP_MTZ en liste de circuits
=#
function vrpToCircuit(model, params, isPDI=false, t=1)

	circuits = Vector{Int64}[]

	endString = "]"

	if isPDI
		endString = ",$t]"
	end

	for i in 1:params["n"]
		#On cherche tout les chemin commençant en 0 (Début de circuit)
		if round(value(variable_by_name(model, "x[0,$i" * endString))) > 0

			circuit = [0]
			#from deviant le premier noeud après 0 du circuit
			from = i

			#Puis tant qu'on a pas fais une boucle (Un circuit complet)...
			while from != 0

				#On ajout from ici pour ne pas avoir de 0 à la fin
				push!(circuit, from)

				for i in 0:params["n"]
					if i != from

						#...On cherche le noeud suivant from...
						if round(value(variable_by_name(model, "x[$from,$i" * endString))) > 0
							#...Et ce noeud devient from
							from = i
							break
						end

					end
				end

			end

			push!(circuits, circuit)

		end
	end

	return circuits

end

#=
Retourne les VRP mais adapter a la vrai demande après un passage de LSP à chaque temps
=#
function getTrueVRP(params, nodes, demands, costs)

	model = createLSP(params, nodes, demands, costs)

	resolvePlne(model, 0)
	println()

	demandsAtT = Array{Int, 2}(undef, params["n"], params["l"]) #demandsAtT = qté à livrer u pas de temps t
	notEmptyIndexAtT = [[] for i in 1:params["l"]]

	for i in 1:params["n"]
		for t in 1:params["l"]
			demandsAtT[i, t] = value(variable_by_name(model, "q[$i,$t]"))

			if demandsAtT[i, t] != 0
				push!(notEmptyIndexAtT[t], i)
			end

		end
	end

	tToVrp = Any[]

	for t in 1:params["l"]

		#NOTE: CHOISIR LE t A TESTER ICI, PLUS TARD LE RETIRER ET METTRE LA BOUCLE
		#t=2

		#On copie les paramètres mais n prend la valeur du nombre de noeud avec une demande supérieur à 0
		copyParams = copy(params)
		copyParams["n"] = size(notEmptyIndexAtT[t], 1)

		#On copie les informations des noeuds mais dans sans les noeud avec une demande 0
		copyNodes = Array{Dict, 1}(undef, size(notEmptyIndexAtT[t], 1) + 1)

		#On copie les demande à chaque temps mais sans les noeud avec une demande à 0
		copyDemandsAtT = Array{Int64, 2}(undef, size(notEmptyIndexAtT[t], 1), params["l"])

		#Initialise le noeud 0 à part
		copyNodes[1] = nodes[0]

		for i in eachindex(notEmptyIndexAtT[t])
			copyNodes[i+1] = nodes[notEmptyIndexAtT[t][i]]
			copyDemandsAtT[i, :] = demandsAtT[notEmptyIndexAtT[t][i], :]
		end

		#On remet les index en commençant à 0
		copyNodes = OffsetVector(copyNodes, 0:(size(copyNodes, 1) - 1))

		#Un dictionnaire faisant le lien entre les anciens index et les nouveaux (Du au décalage d'index dans les array)
		isPresent = Dict{Int, Int}()

		for (index, elem) in enumerate(notEmptyIndexAtT[t])
			isPresent[elem] = index
		end

		isPresent[0] = 0

		#On copie les coûts de transports mais sans les arêtes passant par un noeud avec une demande à 0
		copyCost = Dict{Tuple{Int, Int}, Float64}()
		for (edge, edgeCost) in costs

			if in(edge[1], keys(isPresent)) && in(edge[2], keys(isPresent))
				copyCost[(isPresent[edge[1]], isPresent[edge[2]])] = edgeCost
			end

		end

		push!(tToVrp, (copyParams, copyNodes, copyDemandsAtT, copyCost))

	end

	return tToVrp

end


text:
# Autogenerated wrapper script for htslib_jll for x86_64-unknown-freebsd
export bgzip, htsfile, libhts, tabix

using Zlib_jll
using Bzip2_jll
using XZ_jll
using LibCURL_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("htslib")
JLLWrappers.@declare_executable_product(bgzip)
JLLWrappers.@declare_executable_product(htsfile)
JLLWrappers.@declare_library_product(libhts, "libhts.so.3")
JLLWrappers.@declare_executable_product(tabix)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Bzip2_jll, XZ_jll, LibCURL_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        bgzip,
        "bin/bgzip",
    )

    JLLWrappers.@init_executable_product(
        htsfile,
        "bin/htsfile",
    )

    JLLWrappers.@init_library_product(
        libhts,
        "lib/libhts.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        tabix,
        "bin/tabix",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@inline function change_t_via_interpolation!{T}(integrator,t,modify_save_endpoint::Type{Val{T}}=Val{false})
  # Can get rid of an allocation here with a function
  # get_tmp_arr(integrator.cache) which gives a pointer to some
  # cache array which can be modified.
  if t < integrator.tprev
    error("Current interpolant only works between tprev and t")
  elseif t != integrator.t

    if typeof(integrator.u) <: AbstractArray
      integrator(integrator.u,t)
    else
      integrator.u = integrator(t)
    end
    integrator.t = t
    integrator.dt = integrator.t - integrator.tprev
    reeval_internals_due_to_modification!(integrator)
    if T
      solution_endpoint_match_cur_integrator!(integrator)
    end
  end
end

@inline function reeval_internals_due_to_modification!(integrator)
  if integrator.opts.calck
    resize!(integrator.k,integrator.kshortsize) # Reset k for next step!
    ode_addsteps!(integrator,integrator.f,Val{true},Val{false})
  end
  integrator.u_modified = false
end

@inline function u_modified!(integrator::ODEIntegrator,bool::Bool)
  integrator.u_modified = bool
end

@inline get_proposed_dt(integrator::ODEIntegrator) = integrator.dtpropose
@inline set_proposed_dt!(integrator::ODEIntegrator,dt) = (integrator.dtpropose = dt)

user_cache(integrator::ODEIntegrator) = user_cache(integrator.cache)
u_cache(integrator::ODEIntegrator) = u_cache(integrator.cache)
du_cache(integrator::ODEIntegrator)= du_cache(integrator.cache)
full_cache(integrator::ODEIntegrator) = chain(user_cache(integrator),u_cache(integrator),du_cache(integrator.cache))
default_non_user_cache(integrator::ODEIntegrator) = chain(u_cache(integrator),du_cache(integrator.cache))
@inline add_tstop!(integrator::ODEIntegrator,t) = push!(integrator.opts.tstops,t)
user_cache(cache::OrdinaryDiffEqCache) = (cache.u,cache.uprev,cache.tmp)

resize!(integrator::ODEIntegrator,i::Int) = resize!(integrator,integrator.cache,i)
function resize!(integrator::ODEIntegrator,cache,i)
  for c in user_cache(integrator)
    resize!(c,i)
  end
  resize_non_user_cache!(integrator,cache,i)
end

resize_non_user_cache!(integrator::ODEIntegrator,i::Int) = resize_non_user_cache!(integrator,integrator.cache,i)
deleteat_non_user_cache!(integrator::ODEIntegrator,i) = deleteat_non_user_cache!(integrator,integrator.cache,i)
addat_non_user_cache!(integrator::ODEIntegrator,i) = addat_non_user_cache!(integrator,integrator.cache,i)

function resize_non_user_cache!(integrator::ODEIntegrator,cache,i)
  for c in default_non_user_cache(integrator)
    resize!(c,i)
  end
end

function resize_non_user_cache!(integrator::ODEIntegrator,cache::Union{Rosenbrock23Cache,Rosenbrock32Cache},i)
  for c in default_non_user_cache(integrator)
    resize!(c,i)
  end
  for c in vecu_cache(integrator.cache)
    resize!(c,i)
  end
  Jvec = vec(cache.J)
  cache.J = reshape(resize!(Jvec,i*i),i,i)
  Wvec = vec(cache.W)
  cache.W = reshape(resize!(Wvec,i*i),i,i)
  resize!(cache.jac_config.duals[1],i)
end
user_cache(cache::Union{Rosenbrock23Cache,Rosenbrock32Cache}) = (cache.u,cache.uprev,cache.jac_config.duals[2])

function resize_non_user_cache!(integrator::ODEIntegrator,cache::Union{ImplicitEulerCache,TrapezoidCache},i)
  for c in default_non_user_cache(integrator)
    resize!(c,i)
  end
  for c in vecu_cache(integrator.cache)
    resize!(c,i)
  end
  for c in dual_cache(integrator.cache)
    resize!(c.du,i)
    resize!(c.dual_du,i)
  end
  if alg_autodiff(integrator.alg)
    cache.adf = autodiff_setup(cache.rhs,cache.uhold,integrator.alg)
  end
end

function deleteat_non_user_cache!(integrator::ODEIntegrator,cache,idxs)
  # ordering doesn't matter in deterministic cache, so just resize
  # to match the size of u
  i = length(integrator.u)
  resize_non_user_cache!(integrator,cache,i)
end

function addat_non_user_cache!(integrator::ODEIntegrator,cache,idxs)
  # ordering doesn't matter in deterministic cache, so just resize
  # to match the size of u
  i = length(integrator.u)
  resize_non_user_cache!(integrator,cache,i)
end

function deleteat!(integrator::ODEIntegrator,idxs)
  for c in user_cache(integrator)
    deleteat!(c,idxs)
  end
  deleteat_non_user_cache!(integrator,cache,idxs)
end

function addat!(integrator::ODEIntegrator,idxs)
  for c in user_cache(integrator)
    addat!(c,idxs)
  end
  addat_non_user_cache!(integrator,cache,idxs)
end

function terminate!(integrator::ODEIntegrator)
  integrator.opts.tstops.valtree = typeof(integrator.opts.tstops.valtree)()
end


text:
@testset "Robots: box" begin
    # TODO: add tests
    q0 = nominal_configuration(box)
    # visualizer
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, box, [q0], Δt=0.1);
end


text:
#
# OIFITS.jl --
#
# Support for OI-FITS (optical interferometry data format) in Julia.
#
#------------------------------------------------------------------------------

module OIFITS

export
    OIArray,
    OICorrelation,
    OIDataBlock,
    OIFlux,
    OIMaster,
    OIPolarization,
    OIT3,
    OITarget,
    OIVis,
    OIVis2,
    OIWavelength

using Compat
import Base: getindex, setindex!, haskey, keys, show

using FITSIO
using FITSIO.Libcfitsio

using FITSIO: libcfitsio, fits_assert_ok, fits_assert_open
using FITSIO.Libcfitsio: fits_get_errstatus

include("types.jl")
include("utils.jl")
include("objects.jl")
include("builder.jl")
include("files.jl")
include("accessors.jl") # must be *after* oifile.jl and all oiformat*.jl

end # module


text:
## This file implements 2-3 trees for sorted containers.
## The types and functions in this file are not exported; they
## are meant to be used by SortedDict, MultiMap and SortedSet.
## A 2-3 tree is a rooted trees in which all leaves are at the same
## depth (i.e., it is "balanced")
## and in which each internal node of the tree has either 2 or 3
## children.
## All the internal tree
## nodes are stored in an array of TreeNodes.  The bottom layer
## of internal tree nodes,
## called the "leaves" in this file, sit above one more layer
## of data nodes, which are stored in a different array.

## KDRec is one data node:
##  k: the key of the node
##  d: the data of the node
##  parent: the tree leaf that is the parent of this
##    node.  Parent pointers are needed in order
##    to implement indices.
##  There are two constructors, the standard one (first)
##  and the incomplete one (second).  The incomplete constructor
##  is needed because when the data structure is first created,
##  there are no valid K or D values to store in the initial
##  data nodes.

struct KDRec{K,D}
    parent::Int
    k::K
    d::D
    KDRec{K,D}(p::Int, k1::K, d1::D) where {K,D} = new{K,D}(p,k1,d1)
    KDRec{K,D}(p::Int) where {K,D} = new{K,D}(p)
end

## TreeNode is an internal node of the tree.
## child1,child2,child3:
##     These are the three children node numbers.
##     If the node is a 2-node (rather than 3), then child3 == 0.
##     If this is a leaf then child1,child2,child3 are subscripts
##       of data nodes, else they are subscripts  of other tree nodes.
## splitkey1:
##    the minimum key of the subtree at child2.  If this is a leaf
##    then it is the key of child2.
## splitkey2:
##    if child3 > 0 then splitkey2 is the minimum key of the subtree at child3.
##    If this is a leaf, then it is the key of child3.
## Again, there are two constructors for the same reason mentioned above.

struct TreeNode{K}
    child1::Int
    child2::Int
    child3::Int
    parent::Int
    splitkey1::K
    splitkey2::K
    TreeNode{K}(::Type{K}, c1::Int, c2::Int, c3::Int, p::Int) where {K} = new{K}(c1, c2, c3, p)
    TreeNode{K}(c1::Int, c2::Int, c3::Int, p::Int, sk1::K, sk2::K) where {K} =
        new{K}(c1, c2, c3, p, sk1, sk2)
end

## The next two functions are called to initialize the tree
## by inserting a dummy tree node with two children, the before-start
## marker whose index is 1 and the after-end marker whose index is 2.
## These two markers live in dummy data nodes.

function initializeTree!(tree::Array{TreeNode{K},1}) where K
    resize!(tree,1)
    tree[1] = TreeNode{K}(K, 1, 2, 0, 0)
    nothing
end

function initializeData!(data::Array{KDRec{K,D},1}) where {K,D}
    resize!(data, 2)
    data[1] = KDRec{K,D}(1)
    data[2] = KDRec{K,D}(1)
    nothing
end

## Type BalancedTree23{K,D,Ord} is 'base class' for
## SortedDict.
## K = key type, D = data type
## Key type must support an ordering operation defined by Ordering
## object Ord.
## The default is Forward which implies that the ordering function
## is isless (see ordering.jl)
## The fields are as follows.
## ord:: The ordering object.  Often the ordering type
##   is a singleton type, so this field is empty, but it
##   is still necessary to direct the multiple dispatch.
## data: the (key,data) pairs of the tree.
##   The first and second entries of the data array are dummy placeholders
##   for the beginning and end of the sorted order of the keys
## tree: the nodes of a 2-3 tree that sits above the data.
## rootloc: the index of the entry of tree (i.e., a subscript to
##   treenodes) that is the tree's root
## depth: the depth of the tree, (number
##    of tree levels, not counting the level of data at the bottom)
##    depth==1 means that there is a single root node
##      whose children are data nodes.
## freetreeinds: Array of indices of free locations in the
##    tree array (locations are freed due to deletion)
## freedatainds: Array of indices of free locations in the
##    data array (locations are freed due to deletion)
## useddatacells: IntSet (i.e., bit vector) showing which
##    data cells are taken.  The complementary positions are
##    exactly those stored in freedatainds.  This array is
##    used only for error checking (only present at debug level 1 and 2)
## deletionchild and deletionleftkey are two work-arrays
## for the delete function.

mutable struct BalancedTree23{K, D, Ord <: Ordering}
    ord::Ord
    data::Array{KDRec{K,D}, 1}
    tree::Array{TreeNode{K}, 1}
    rootloc::Int
    depth::Int
    freetreeinds::Array{Int,1}
    freedatainds::Array{Int,1}
    useddatacells::IntSet
    # The next two arrays are used as a workspace by the delete!
    # function.
    deletionchild::Array{Int,1}
    deletionleftkey::Array{K,1}
    function BalancedTree23{K,D,Ord}(ord1::Ord) where {K,D,Ord<:Ordering}
        tree1 = Vector{TreeNode{K}}(uninitialized, 1)
        initializeTree!(tree1)
        data1 = Vector{KDRec{K,D}}(uninitialized, 2)
        initializeData!(data1)
        u1 = IntSet()
        push!(u1, 1, 2)
        new{K,D,Ord}(ord1, data1, tree1, 1, 1, Vector{Int}(), Vector{Int}(),
                     u1,
                     Vector{Int}(uninitialized, 3), Vector{K}(uninitialized, 3))
    end
end

## Function cmp2 checks a tree node with two children
## against a given key, and returns 1 if the given key is
## less than the node's splitkey or 2 else.  Special case
## if the node is a leaf and its right child is the end
## of the sorted order.

@inline function cmp2_nonleaf(o::Ordering,
                              treenode::TreeNode,
                              k)
    lt(o, k, treenode.splitkey1) ? 1 : 2
end

@inline function cmp2_leaf(o::Ordering,
                           treenode::TreeNode,
                           k)
    (treenode.child2 == 2) ||
    lt(o, k, treenode.splitkey1) ? 1 : 2
end

## Function cmp3 checks a tree node with three children
## against a given key, and returns 1 if the given key is
## less than the node's splitkey1, 2 if the key is greater than or
## equal to splitkey1 but less than splitkey2, or 3 else.  Special case
## if the node is a leaf and its right child is the end
## of the sorted order.

@inline function cmp3_nonleaf(o::Ordering,
                              treenode::TreeNode,
                              k)
    lt(o, k, treenode.splitkey1) ? 1 :
    lt(o, k, treenode.splitkey2) ? 2 : 3
end

@inline function cmp3_leaf(o::Ordering,
                           treenode::TreeNode,
                           k)
    lt(o, k, treenode.splitkey1) ?                           1 :
    (treenode.child3 == 2 || lt(o, k, treenode.splitkey2)) ? 2 : 3
end

## Function cmp2le checks a tree node with two children
## against a given key, and returns 1 if the given key is
## less than or equal to the node's splitkey or 2 else.  Special case
## if the node is a leaf and its right child is the end
## of the sorted order.

@inline function cmp2le_nonleaf(o::Ordering,
                                treenode::TreeNode,
                                k)
    !lt(o,treenode.splitkey1,k) ? 1 : 2
end

@inline function cmp2le_leaf(o::Ordering,
                             treenode::TreeNode,
                             k)
    treenode.child2 == 2 || !lt(o,treenode.splitkey1,k) ? 1 : 2
end

## Function cmp3le checks a tree node with three children
## against a given key, and returns 1 if the given key is
## less than or equal to the node's splitkey1, 2 if less than or equal
## to splitkey2, or
## 3 else. Special case
## if the node is a leaf and its right child is the end
## of the sorted order.

@inline function cmp3le_nonleaf(o::Ordering,
                                treenode::TreeNode,
                                k)
    !lt(o,treenode.splitkey1, k) ? 1 :
    !lt(o,treenode.splitkey2, k) ? 2 : 3
end

@inline function cmp3le_leaf(o::Ordering,
                             treenode::TreeNode,
                             k)
    !lt(o,treenode.splitkey1,k) ?                            1 :
    (treenode.child3 == 2 || !lt(o,treenode.splitkey2, k)) ? 2 : 3
end

## The empty! function deletes all data in the balanced tree.
## Therefore, it invalidates all indices.

function empty!(t::BalancedTree23)
    resize!(t.data,2)
    initializeData!(t.data)
    resize!(t.tree,1)
    initializeTree!(t.tree)
    t.depth = 1
    t.rootloc = 1
    t.freetreeinds = Vector{Int}()
    t.freedatainds = Vector{Int}()
    empty!(t.useddatacells)
    push!(t.useddatacells, 1, 2)
    nothing
end

## Default implementations of eq for Forward, Reverse
## and default ordering

eq(::ForwardOrdering, a, b) = isequal(a,b)
eq(::ReverseOrdering{ForwardOrdering}, a, b) = isequal(a,b)
eq(o::Ordering, a, b) = !lt(o, a, b) && !lt(o, b, a)

## The findkey function finds the index of a (key,data) pair in the tree
## where the given key lives (if it is present), or
## if the key is not present, to the lower bound for the key,
## i.e., the data item that comes immediately before it.
## If there are multiple equal keys, then it finds the last one.
## It returns the index of the key found and a boolean indicating
## whether the exact key was found or not.

function findkey(t::BalancedTree23, k)
    curnode = t.rootloc
    for depthcount = 1 : t.depth - 1
        @inbounds thisnode = t.tree[curnode]
        cmp = thisnode.child3 == 0 ?
                         cmp2_nonleaf(t.ord, thisnode, k) :
                         cmp3_nonleaf(t.ord, thisnode, k)
        curnode = cmp == 1 ? thisnode.child1 :
                  cmp == 2 ? thisnode.child2 : thisnode.child3
    end
    @inbounds thisnode = t.tree[curnode]
    cmp = thisnode.child3 == 0 ?
                cmp2_leaf(t.ord, thisnode, k) :
                cmp3_leaf(t.ord, thisnode, k)
    curnode = cmp == 1 ? thisnode.child1 :
              cmp == 2 ? thisnode.child2 : thisnode.child3
    @inbounds return curnode, (curnode > 2 && eq(t.ord, t.data[curnode].k, k))
end

## The findkeyless function finds the index of a (key,data) pair in the tree that
## with the greatest key that is less than the given key.  If there is no
## key less than the given key, then it returns 1 (the before-start node).

function findkeyless(t::BalancedTree23, k)
    curnode = t.rootloc
    for depthcount = 1 : t.depth - 1
        @inbounds thisnode = t.tree[curnode]
        cmp = thisnode.child3 == 0 ?
               cmp2le_nonleaf(t.ord, thisnode, k) :
               cmp3le_nonleaf(t.ord, thisnode, k)
        curnode = cmp == 1 ? thisnode.child1 :
                  cmp == 2 ? thisnode.child2 : thisnode.child3
    end
    @inbounds thisnode = t.tree[curnode]
    cmp = thisnode.child3 == 0 ?
            cmp2le_leaf(t.ord, thisnode, k) :
            cmp3le_leaf(t.ord, thisnode, k)
    curnode = cmp == 1 ? thisnode.child1 :
              cmp == 2 ? thisnode.child2 : thisnode.child3
    curnode
end

## The following are helper routines for the insert! and delete! functions.
## They replace the 'parent' field of either an internal tree node or
## a data node at the bottom tree level.

function replaceparent!(data::Array{KDRec{K,D},1}, whichind::Int, newparent::Int) where {K,D}
    data[whichind] = KDRec{K,D}(newparent, data[whichind].k, data[whichind].d)
    nothing
end

function replaceparent!(tree::Array{TreeNode{K},1}, whichind::Int, newparent::Int) where K
    tree[whichind] = TreeNode{K}(tree[whichind].child1, tree[whichind].child2,
                                 tree[whichind].child3, newparent,
                                 tree[whichind].splitkey1,
                                 tree[whichind].splitkey2)
    nothing
end

## Helper function for either grabbing a brand new location in a
## array (if there are no free locations) or else grabbing a free
## location and marking it as used.  The return value is the
## index of the data just inserted into the vector.

function push_or_reuse!(a::Vector, freelocs::Array{Int,1}, item)
    if isempty(freelocs)
        push!(a, item)
        return length(a)
    end
    loc = pop!(freelocs)
    a[loc] = item
    return loc
end

## Function insert! inserts a new data item into the tree.
## The arguments are the (K,D) pair to insert.
## The return values are a bool and an index.  The
## bool indicates whether the insertion inserted a new record (true) or
## whether it replaced an existing record (false).
## The index returned is the subscript in t.data where the
## inserted value sits.
## "allowdups" (i.e., "allow duplicate keys") means that no check is
## done whether the iterm
## is already in the tree, so insertion of a new item always succeeds.

function insert!(t::BalancedTree23{K,D,Ord}, k, d, allowdups::Bool) where {K,D,Ord <: Ordering}

    ## First we find the greatest data node that is <= k.
    leafind, exactfound = findkey(t, k)
    parent = t.data[leafind].parent

    ## The following code is necessary because in the case of a
    ## brand new tree, the initial tree and data entries were incompletely
    ## initialized by the constructor.  In this case, the call to insert!
    ## underway carries
    ## valid K and D values, so these valid values may now be
    ## stored in the dummy placeholder nodes so that they no
    ## longer hold undefined references.

    if size(t.data,1) == 2
        # @assert(t.rootloc == 1 && t.depth == 1)
        t.tree[1] = TreeNode{K}(t.tree[1].child1, t.tree[1].child2,
                                t.tree[1].child3, t.tree[1].parent,
                                k, k)
        t.data[1] = KDRec{K,D}(t.data[1].parent, k, d)
        t.data[2] = KDRec{K,D}(t.data[2].parent, k, d)
    end

    ## If we have found exactly k in the tree, then we
    ## replace the data associated with k and return.

    if exactfound && !allowdups
        t.data[leafind] = KDRec{K,D}(parent, k,d)
        return false, leafind
    end

    # We get here if k was not already found in the tree or
    # if duplicates are allowed.
    # In this case we insert a new node.
    depth = t.depth
    ord = t.ord

    ## Store the new data item in the tree's data array.  Later
    ## go back and fix the parent.

    newind = push_or_reuse!(t.data, t.freedatainds, KDRec{K,D}(0,k,d))
    p1 = parent
    oldchild = leafind
    newchild = newind
    minkeynewchild = k
    splitroot = false
    curdepth = depth

    ## This loop ascends the tree (i.e., follows the path from a leaf to the root)
    ## starting from the parent p1 of
    ## where the new key k would go.  For each 3-node we encounter
    ## during the ascent, we add a new child, which requires splitting
    ## the 3-node into two 2-nodes.  Then we keep going until we hit the root.
    ## If we encounter a 2-node, then the ascent can stop; we can
    ## change the 2-node to a 3-node with the new child. Invariants
    ## during this loop are:
    ##     p1: the parent node (a tree node index) where the insertion must occur
    ##     oldchild,newchild: the two children of the parent node; oldchild
    ##          was already in the tree; newchild was just added to it.
    ##     minkeynewchild:  This is the key that is the minimum value in
    ##         the subtree rooted at newchild.

    while t.tree[p1].child3 > 0
        isleaf = (curdepth == depth)
        oldtreenode = t.tree[p1]

        ## Node p1 index a 3-node. There are three cases for how to
        ## insert new child.  All three cases involve splitting the
        ## existing node (oldtreenode, numbered p1) into
        ## two new nodes.  One keeps the index p1; the other has
        ## has a new index called newparentnum.

        cmp = isleaf ? cmp3_leaf(ord, oldtreenode, minkeynewchild) :
                      cmp3_nonleaf(ord, oldtreenode, minkeynewchild)

        if cmp == 1
            lefttreenodenew = TreeNode{K}(oldtreenode.child1, newchild, 0,
                                          oldtreenode.parent,
                                          minkeynewchild, minkeynewchild)
            righttreenodenew = TreeNode{K}(oldtreenode.child2, oldtreenode.child3, 0,
                                           oldtreenode.parent, oldtreenode.splitkey2,
                                           oldtreenode.splitkey2)
            minkeynewchild = oldtreenode.splitkey1
            whichp = 1
        elseif cmp == 2
            lefttreenodenew = TreeNode{K}(oldtreenode.child1, oldtreenode.child2, 0,
                                          oldtreenode.parent,
                                          oldtreenode.splitkey1, oldtreenode.splitkey1)
            righttreenodenew = TreeNode{K}(newchild, oldtreenode.child3, 0,
                                           oldtreenode.parent,
                                           oldtreenode.splitkey2, oldtreenode.splitkey2)
            whichp = 2
        else
            lefttreenodenew = TreeNode{K}(oldtreenode.child1, oldtreenode.child2, 0,
                                          oldtreenode.parent,
                                          oldtreenode.splitkey1, oldtreenode.splitkey1)
            righttreenodenew = TreeNode{K}(oldtreenode.child3, newchild, 0,
                                           oldtreenode.parent,
                                           minkeynewchild, minkeynewchild)
            minkeynewchild = oldtreenode.splitkey2
            whichp = 2
        end
        # Replace p1 with a new 2-node and insert another 2-node at
        # index newparentnum.
        t.tree[p1] = lefttreenodenew
        newparentnum = push_or_reuse!(t.tree, t.freetreeinds, righttreenodenew)
        if isleaf
            par = (whichp == 1) ? p1 : newparentnum
            # fix the parent of the new datanode.
            replaceparent!(t.data, newind, par)
            push!(t.useddatacells, newind)
            replaceparent!(t.data, righttreenodenew.child1, newparentnum)
            replaceparent!(t.data, righttreenodenew.child2, newparentnum)
        else
            # If this is not a leaf, we still have to fix the
            # parent fields of the two nodes that are now children
            ## of the newparent.
            replaceparent!(t.tree, righttreenodenew.child1, newparentnum)
            replaceparent!(t.tree, righttreenodenew.child2, newparentnum)
        end
        oldchild = p1
        newchild = newparentnum
        ## If p1 is the root (i.e., we have encountered only 3-nodes during
        ## our ascent of the tree), then the root must be split.
        if p1 == t.rootloc
            # @assert(curdepth == 1)
            splitroot = true
            break
        end
        p1 = t.tree[oldchild].parent
        curdepth -= 1
    end

    ## big loop terminated either because a 2-node was reached
    ## (splitroot == false) or we went up the whole tree seeing
    ## only 3-nodes (splitroot == true).
    if !splitroot

        ## If our ascent reached a 2-node, then we convert it to
        ## a 3-node by giving it a child3 field that is >0.
        ## Encountering a 2-node halts the ascent up the tree.

        isleaf = curdepth == depth
        oldtreenode = t.tree[p1]
        cmpres = isleaf ? cmp2_leaf(ord, oldtreenode, minkeynewchild) :
                         cmp2_nonleaf(ord, oldtreenode, minkeynewchild)

        t.tree[p1] = cmpres == 1 ?
                         TreeNode{K}(oldtreenode.child1, newchild, oldtreenode.child2,
                                     oldtreenode.parent,
                                     minkeynewchild, oldtreenode.splitkey1) :
                         TreeNode{K}(oldtreenode.child1, oldtreenode.child2, newchild,
                                     oldtreenode.parent,
                                     oldtreenode.splitkey1, minkeynewchild)
        if isleaf
            replaceparent!(t.data, newind, p1)
            push!(t.useddatacells, newind)
        end
    else
        ## Splitroot is set if the ascent of the tree encountered only 3-nodes.
        ## In this case, the root itself was replaced by two nodes, so we need
        ## a new root above those two.

        newroot = TreeNode{K}(oldchild, newchild, 0, 0,
                              minkeynewchild, minkeynewchild)
        newrootloc = push_or_reuse!(t.tree, t.freetreeinds, newroot)
        replaceparent!(t.tree, oldchild, newrootloc)
        replaceparent!(t.tree, newchild, newrootloc)
        t.rootloc = newrootloc
        t.depth += 1
    end
    true, newind
end

## nextloc0: returns the next item in the tree according to the
## sort order, given an index i (subscript of t.data) of a current
## item.
## The routine returns 2 if there is no next item (i.e., we started
## from the last one in the sorted order).

function nextloc0(t, i::Int)
    ii = i
    # @assert(i != 2 && i in t.useddatacells)
    @inbounds p = t.data[i].parent
    nextchild = 0
    depthp = t.depth
    @inbounds while true
        if depthp < t.depth
            p = t.tree[ii].parent
        end
        if t.tree[p].child1 == ii
            nextchild = t.tree[p].child2
            break
        end
        if t.tree[p].child2 == ii && t.tree[p].child3 > 0
            nextchild = t.tree[p].child3
            break
        end
        ii = p
        depthp -= 1
    end
    @inbounds while true
        if depthp == t.depth
            return nextchild
        end
        p = nextchild
        nextchild = t.tree[p].child1
        depthp += 1
    end
end

## prevloc0: returns the previous item in the tree according to the
## sort order, given an index i (subscript of t.data) of a current
## item.
## The routine returns 1 if there is no previous item (i.e., we started
## from the first one in the sorted order).

function prevloc0(t::BalancedTree23, i::Int)
    # @assert(i != 1 && i in t.useddatacells)
    ii = i
    @inbounds p = t.data[i].parent
    prevchild = 0
    depthp = t.depth
    @inbounds while true
        if depthp < t.depth
            p = t.tree[ii].parent
        end
        if t.tree[p].child3 == ii
            prevchild = t.tree[p].child2
            break
        end
        if t.tree[p].child2 == ii
            prevchild = t.tree[p].child1
            break
        end
        ii = p
        depthp -= 1
    end
    @inbounds while true
        if depthp == t.depth
            return prevchild
        end
        p = prevchild
        c3 = t.tree[p].child3
        prevchild = c3 > 0 ? c3 : t.tree[p].child2
        depthp += 1
    end
end

## This function takes two indices into t.data and checks which
## one comes first in the sorted order by chasing them both
## up the tree until a common ancestor is found.
## The return value is -1 if i1 precedes i2, 0 if i1 == i2
##, 1 if i2 precedes i1.

function compareInd(t::BalancedTree23, i1::Int, i2::Int)
    @assert(i1 in t.useddatacells && i2 in t.useddatacells)
    if i1 == i2
        return 0
    end
    i1a = i1
    i2a = i2
    p1 = t.data[i1].parent
    p2 = t.data[i2].parent
    curdepth = t.depth
    while true
        @assert(curdepth > 0)
        if p1 == p2
            if i1a == t.tree[p1].child1
                @assert(t.tree[p1].child2 == i2a || t.tree[p1].child3 == i2a)
                return -1
            end
            if i1a == t.tree[p1].child2
                if (t.tree[p1].child1 == i2a)
                    return 1
                end
                @assert(t.tree[p1].child3 == i2a)
                return -1
            end
            @assert(i1a == t.tree[p1].child3)
            @assert(t.tree[p1].child1 == i2a || t.tree[p1].child2 == i2a)
            return 1
        end
        i1a = p1
        i2a = p2
        p1 = t.tree[i1a].parent
        p2 = t.tree[i2a].parent
        curdepth -= 1
    end
end

## beginloc, endloc return the index (into t.data) of the first, last item in the
## sorted order of the tree.  beginloc works by going to the before-start marker
## (data node 1) and executing a next operation on it.  endloc is the opposite.

beginloc(t::BalancedTree23) = nextloc0(t,1)
endloc(t::BalancedTree23) = prevloc0(t,2)

## delete! routine deletes an entry from the balanced tree.

function delete!(t::BalancedTree23{K,D,Ord}, it::Int) where {K,D,Ord<:Ordering}

    ## Put the cell indexed by 'it' into the deletion list.
    ##
    ## Create the following data items maintained in the
    ## upcoming loop.
    ##
    ## p is a tree-node ancestor of the deleted node
    ## The children of p are stored in
    ## t.deletionchild[..]
    ## The number of these children is newchildcount, which is 1, 2 or 3.
    ## The keys that lower bound the children
    ## are stored in t.deletionleftkey[..]
    ## There is a special case for t.deletionleftkey[1]; the
    ## flag deletionleftkey1_valid indicates that the left key
    ## for the immediate right neighbor of the
    ## deleted node has not yet been been stored in the tree.
    ## Once it is stored, t.deletionleftkey[1] is no longer needed
    ## or used.
    ## The flag mustdeleteroot means that the tree has contracted
    ## enough that it loses a level.

    p = t.data[it].parent
    newchildcount = 0
    c1 = t.tree[p].child1
    deletionleftkey1_valid = true
    if c1 != it
        deletionleftkey1_valid = false
        newchildcount += 1
        t.deletionchild[newchildcount] = c1
        t.deletionleftkey[newchildcount] = t.data[c1].k
    end
    c2 = t.tree[p].child2
    if c2 != it
        newchildcount += 1
        t.deletionchild[newchildcount] = c2
        t.deletionleftkey[newchildcount] = t.data[c2].k
    end
    c3 = t.tree[p].child3
    if c3 != it && c3 > 0
        newchildcount += 1
        t.deletionchild[newchildcount] = c3
        t.deletionleftkey[newchildcount] = t.data[c3].k
    end
    # @assert(newchildcount == 1 || newchildcount == 2)
    push!(t.freedatainds, it)
    pop!(t.useddatacells,it)
    defaultKey = t.tree[1].splitkey1
    curdepth = t.depth
    mustdeleteroot = false
    pparent = -1

    ## The following loop ascends the tree and contracts nodes (reduces their
    ## number of children) as
    ## needed.  If newchildcount == 2 or 3, then the ascent is terminated
    ## and a node is created with 2 or 3 children.
    ## If newchildcount == 1, then the ascent must continue since a tree
    ## node cannot have one child.

    while true
        pparent = t.tree[p].parent
        ## Simple cases when the new child count is 2 or 3
        if newchildcount == 2
            t.tree[p] = TreeNode{K}(t.deletionchild[1],
                                    t.deletionchild[2], 0, pparent,
                                    t.deletionleftkey[2], defaultKey)

            break
        end
        if newchildcount == 3
            t.tree[p] = TreeNode{K}(t.deletionchild[1], t.deletionchild[2],
                                    t.deletionchild[3], pparent,
                                    t.deletionleftkey[2], t.deletionleftkey[3])
            break
        end
        # @assert(newchildcount == 1)
        ## For the rest of this loop, we cover the case
        ## that p has one child.

        ## If newchildcount == 1 and curdepth==1, this means that
        ## the root of the tree has only one child.  In this case, we can
        ## delete the root and make its one child the new root (see below).

        if curdepth == 1
            mustdeleteroot = true
            break
        end

        ## We now branch on three cases depending on whether p is child1,
        ## child2 or child3 of its parent.

        if t.tree[pparent].child1 == p
            rightsib = t.tree[pparent].child2

            ## Here p is child1 and rightsib is child2.
            ## If rightsib has 2 children, then p and
            ## rightsib are merged into a single node
            ## that has three children.
            ## If rightsib has 3 children, then p and
            ## rightsib are reformed so that each has
            ## two children.

            if t.tree[rightsib].child3 == 0
                rc1 = t.tree[rightsib].child1
                rc2 = t.tree[rightsib].child2
                t.tree[p] = TreeNode{K}(t.deletionchild[1],
                                        rc1, rc2,
                                        pparent,
                                        t.tree[pparent].splitkey1,
                                        t.tree[rightsib].splitkey1)
                if curdepth == t.depth
                    replaceparent!(t.data, rc1, p)
                    replaceparent!(t.data, rc2, p)
                else
                    replaceparent!(t.tree, rc1, p)
                    replaceparent!(t.tree, rc2, p)
                end
                push!(t.freetreeinds, rightsib)
                newchildcount = 1
                t.deletionchild[1] = p
            else
                rc1 = t.tree[rightsib].child1
                t.tree[p] = TreeNode{K}(t.deletionchild[1], rc1, 0,
                                        pparent,
                                        t.tree[pparent].splitkey1,
                                        defaultKey)
                sk1 = t.tree[rightsib].splitkey1
                t.tree[rightsib] = TreeNode{K}(t.tree[rightsib].child2,
                                               t.tree[rightsib].child3,
                                               0,
                                               pparent,
                                               t.tree[rightsib].splitkey2,
                                               defaultKey)
                if curdepth == t.depth
                    replaceparent!(t.data, rc1, p)
                else
                    replaceparent!(t.tree, rc1, p)
                end
                newchildcount = 2
                t.deletionchild[1] = p
                t.deletionchild[2] = rightsib
                t.deletionleftkey[2] = sk1
            end

            ## If pparent had a third child (besides p and rightsib)
            ## then we add this to t.deletionchild

            c3 = t.tree[pparent].child3
            if c3 > 0
                newchildcount += 1
                t.deletionchild[newchildcount] = c3
                t.deletionleftkey[newchildcount] = t.tree[pparent].splitkey2
            end
            p = pparent
        elseif t.tree[pparent].child2 == p

            ## Here p is child2 and leftsib is child1.
            ## If leftsib has 2 children, then p and
            ## leftsib are merged into a single node
            ## that has three children.
            ## If leftsib has 3 children, then p and
            ## leftsib are reformed so that each has
            ## two children.

            leftsib = t.tree[pparent].child1
            lk = deletionleftkey1_valid ?
                      t.deletionleftkey[1] :
                      t.tree[pparent].splitkey1
            if t.tree[leftsib].child3 == 0
                lc1 = t.tree[leftsib].child1
                lc2 = t.tree[leftsib].child2
                t.tree[p] = TreeNode{K}(lc1, lc2,
                                        t.deletionchild[1],
                                        pparent,
                                        t.tree[leftsib].splitkey1,
                                        lk)
                if curdepth == t.depth
                    replaceparent!(t.data, lc1, p)
                    replaceparent!(t.data, lc2, p)
                else
                    replaceparent!(t.tree, lc1, p)
                    replaceparent!(t.tree, lc2, p)
                end
                push!(t.freetreeinds, leftsib)
                newchildcount = 1
                t.deletionchild[1] = p
            else
                lc3 = t.tree[leftsib].child3
                t.tree[p] = TreeNode{K}(lc3, t.deletionchild[1], 0,
                                        pparent, lk, defaultKey)
                sk2 = t.tree[leftsib].splitkey2
                t.tree[leftsib] = TreeNode{K}(t.tree[leftsib].child1,
                                              t.tree[leftsib].child2,
                                              0, pparent,
                                              t.tree[leftsib].splitkey1,
                                              defaultKey)
                if curdepth == t.depth
                    replaceparent!(t.data, lc3, p)
                else
                    replaceparent!(t.tree, lc3, p)
                end
                newchildcount = 2
                t.deletionchild[1] = leftsib
                t.deletionchild[2] = p
                t.deletionleftkey[2] = sk2
            end

            ## If pparent had a third child (besides p and leftsib)
            ## then we add this to t.deletionchild

            c3 = t.tree[pparent].child3
            if c3 > 0
                newchildcount += 1
                t.deletionchild[newchildcount] = c3
                t.deletionleftkey[newchildcount] = t.tree[pparent].splitkey2
            end
            p = pparent
            deletionleftkey1_valid = false
        else
            ## Here p is child3 and leftsib is child2.
            ## If leftsib has 2 children, then p and
            ## leftsib are merged into a single node
            ## that has three children.
            ## If leftsib has 3 children, then p and
            ## leftsib are reformed so that each has
            ## two children.

            # @assert(t.tree[pparent].child3 == p)
            leftsib = t.tree[pparent].child2
            lk = deletionleftkey1_valid ?
                       t.deletionleftkey[1] :
                       t.tree[pparent].splitkey2
            if t.tree[leftsib].child3 == 0
                lc1 = t.tree[leftsib].child1
                lc2 = t.tree[leftsib].child2
                t.tree[p] = TreeNode{K}(lc1, lc2,
                                        t.deletionchild[1],
                                        pparent,
                                        t.tree[leftsib].splitkey1,
                                        lk)
                if curdepth == t.depth
                    replaceparent!(t.data, lc1, p)
                    replaceparent!(t.data, lc2, p)
                else
                    replaceparent!(t.tree, lc1, p)
                    replaceparent!(t.tree, lc2, p)
                end
                push!(t.freetreeinds, leftsib)
                newchildcount = 2
                t.deletionchild[1] = t.tree[pparent].child1
                t.deletionleftkey[2] = t.tree[pparent].splitkey1
                t.deletionchild[2] = p
            else
                lc3 = t.tree[leftsib].child3
                t.tree[p] = TreeNode{K}(lc3, t.deletionchild[1], 0,
                                        pparent, lk, defaultKey)
                sk2 = t.tree[leftsib].splitkey2
                t.tree[leftsib] = TreeNode{K}(t.tree[leftsib].child1,
                                              t.tree[leftsib].child2,
                                              0, pparent,
                                              t.tree[leftsib].splitkey1,
                                              defaultKey)
                if curdepth == t.depth
                    replaceparent!(t.data, lc3, p)
                else
                    replaceparent!(t.tree, lc3, p)
                end
                newchildcount = 3
                t.deletionchild[1] = t.tree[pparent].child1
                t.deletionchild[2] = leftsib
                t.deletionchild[3] = p
                t.deletionleftkey[2] = t.tree[pparent].splitkey1
                t.deletionleftkey[3] = sk2
            end
            p = pparent
            deletionleftkey1_valid = false
        end
        curdepth -= 1
    end
    if mustdeleteroot
        # @assert(!deletionleftkey1_valid)
        # @assert(p == t.rootloc)
        t.rootloc = t.deletionchild[1]
        t.depth -= 1
        push!(t.freetreeinds, p)
    end

    ## If deletionleftkey1_valid, this means that the new
    ## min key of the deleted node and its right neighbors
    ## has never been stored in the tree.  It must be stored
    ## as splitkey1 or splitkey2 of some ancestor of the
    ## deleted node, so we continue ascending the tree
    ## until we find a node which has p (and therefore the
    ## deleted node) as its descendent through its second
    ## or third child.
    ## It cannot be the case that the deleted node is
    ## is a descendent of the root always through
    ## first children, since this would mean the deleted
    ## node is the leftmost placeholder, which
    ## cannot be deleted.

    if deletionleftkey1_valid
        while true
            pparentnode = t.tree[pparent]
            if pparentnode.child2 == p
                t.tree[pparent] = TreeNode{K}(pparentnode.child1,
                                              pparentnode.child2,
                                              pparentnode.child3,
                                              pparentnode.parent,
                                              t.deletionleftkey[1],
                                              pparentnode.splitkey2)
                break
            elseif pparentnode.child3 == p
                t.tree[pparent] = TreeNode{K}(pparentnode.child1,
                                              pparentnode.child2,
                                              pparentnode.child3,
                                              pparentnode.parent,
                                              pparentnode.splitkey1,
                                              t.deletionleftkey[1])
                break
            else
                p = pparent
                pparent = pparentnode.parent
                curdepth -= 1
                # @assert(curdepth > 0)
            end
        end
    end
    nothing
end


text:
@testset "Integration tests" begin
    include("integration/json.jl")
    include("integration/basic-read.jl")
end


text:
using PolyChaos, Test

outliers = [-1.3, 1.3]

@test_throws DomainError PolyChaos._throwError(3)

densities = [w_legendre,
			 w_uniform01,
			 w_uniform_11,
			 build_w_jacobi(1.2,3.4),
			 build_w_beta(1.2,3.4)
]

@testset "Supports of densities" begin
	for ρ in densities, x in outliers
		@test_throws DomainError ρ(x)
	end
end

@test_throws DomainError w_laguerre(-1)
@test_throws DomainError build_w_gamma(1)(-1)


text:
# Convenience functions for working with MPI
import MPI

"""
Number of processors used in MPI. Can be called without ensuring initialization.
"""
mpi_nprocs(comm=MPI.COMM_WORLD) = (MPI.Init(); MPI.Comm_size(comm))
mpi_master(comm=MPI.COMM_WORLD) = (MPI.Init(); MPI.Comm_rank(comm) == 0)

mpi_sum(  arr, comm::MPI.Comm) = MPI.Allreduce( arr,   +, comm)
mpi_sum!( arr, comm::MPI.Comm) = MPI.Allreduce!(arr,   +, comm)
mpi_min(  arr, comm::MPI.Comm) = MPI.Allreduce( arr, min, comm)
mpi_min!( arr, comm::MPI.Comm) = MPI.Allreduce!(arr, min, comm)
mpi_max(  arr, comm::MPI.Comm) = MPI.Allreduce( arr, max, comm)
mpi_max!( arr, comm::MPI.Comm) = MPI.Allreduce!(arr, max, comm)
mpi_mean( arr, comm::MPI.Comm) = mpi_sum(arr, comm) ./ mpi_nprocs(comm)
mpi_mean!(arr, comm::MPI.Comm) = (mpi_sum!(arr, comm); arr ./= mpi_nprocs(comm))


text:
using PyCall
using PyPlot
function matplotlibstyle()
    py"""
    from matplotlib import rcParams
    import numpy as np
    def mplfigsize(hscale, publicationcols=1,columnsinches=3.1):
        bonuscol=-0.12
        W=columnsinches*(publicationcols+bonuscol)
        return [W, 0.6180469715698392*hscale*W]
    def setfonts(fs):
        smaller=fs*9/11
        rcParams.update({
        "font.size": fs,
        'xtick.labelsize': smaller,
        'axes.labelsize': fs,
        'ytick.labelsize': smaller,
        'legend.fontsize': smaller,
        })
    def setstyle():
        rcParams.update({
        #    "pgf.texsystem": "pdflatex",
        #    "font.family": "serif",
            "font.family": "STIXGeneral",
            "mathtext.fontset": "stix",
        #    "text.usetex": True,
        #    "pgf.rcfonts": False,
            "figure.figsize": mplfigsize(1),
            "figure.autolayout": True,
            "figure.dpi": 300,
            'axes.spines.left': True,
            'axes.spines.bottom': True,
            'axes.spines.right': False,
            'axes.spines.top': False,
            'ytick.left': True,
            'xtick.bottom': True,
            'ytick.right': False,
            'xtick.top': False,
            'ytick.direction': 'in',
            'xtick.direction': 'in',
            'xtick.color': 'black',
            'ytick.color': 'black',
            'axes.edgecolor': 'black',
            'legend.frameon': False,
            'legend.columnspacing': 0.0,
            'legend.labelspacing': 0.2,
            'legend.borderaxespad': 0.5,
            'legend.borderpad': 0.0,
            'legend.handletextpad': 0.2,
            'legend.handlelength': 1.5,
            'lines.linewidth': 1.,
        })
        setfonts(9)
    setstyle()
    """
end
matplotlibstyle()
setstyle()=py"setstyle"()
mplfigsize(hscale, pubcols=1,colinches=3.1)=py"mplfigsize"(hscale, publicationcols=pubcols,columnsinches=colinches)
setfonts(x=9)=py"setfonts"(x)


text:
function sqlite3_step(stmtptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_step, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        stmtptr
    )
    @unix_only return ccall(
        (:sqlite3_step, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        stmtptr
    )
end

function sqlite3_reset(stmtptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_reset, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, ),
        stmtptr
    )
    @unix_only return ccall(
        (:sqlite3_reset, sqlite3_lib),
        Cint,
        (Ptr{Void}, ),
        stmtptr
    )
end

# TODO: Create all of the bind functions in a macro loop
function sqlite3_bind_double(
    stmtptr::Ptr{Void},
    index::Int,
    value::Float64
)
    @windows_only return ccall(
        (:sqlite3_bind_double, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, Cint, Float64),
        stmtptr,
        index,
        value
    )
    @unix_only return ccall(
        (:sqlite3_bind_double, sqlite3_lib),
        Cint,
        (Ptr{Void}, Cint, Float64),
        stmtptr,
        index,
        value
    )
end

function sqlite3_bind_int(
    stmtptr::Ptr{Void},
    index::Int,
    value::Int32
)
    @windows_only return ccall(
        (:sqlite3_bind_int, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, Cint, Int32),
        stmtptr,
        index,
        value
    )
    @unix_only return ccall(
        (:sqlite3_bind_int, sqlite3_lib),
        Cint,
        (Ptr{Void}, Cint, Int32),
        stmtptr,
        index,
        value
    )
end

function sqlite3_bind_int64(
    stmtptr::Ptr{Void},
    index::Int,
    value::Int64
)
    @windows_only return ccall(
        (:sqlite3_bind_int64, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, Cint, Int64),
        stmtptr,
        index,
        value
    )
    @unix_only return ccall(
        (:sqlite3_bind_int64, sqlite3_lib),
        Cint,
        (Ptr{Void}, Cint, Int64),
        stmtptr,
        index,
        value
    )
end

function sqlite3_bind_null(
    stmtptr::Ptr{Void},
    index::Int
)
    @windows_only return ccall(
        (:sqlite3_bind_null, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, Cint),
        stmtptr,
        index
    )
    @unix_only return ccall(
        (:sqlite3_bind_null, sqlite3_lib),
        Cint,
        (Ptr{Void}, Cint),
        stmtptr,
        index
    )
end

function sqlite3_bind_text(
    stmtptr::Ptr{Void},
    index::Int,
    value::String,
    len::Int,
    cb::Ptr{Void}
)
    @windows_only return ccall(
        (:sqlite3_bind_text, sqlite3_lib),
        stdcall,
        Cint,
        (Ptr{Void}, Cint, Ptr{Uint8}, Cint, Ptr{Void}),
        stmtptr,
        index,
        value,
        len,
        cb
    )
    @unix_only return ccall(
        (:sqlite3_bind_text, sqlite3_lib),
        Cint,
        (Ptr{Void}, Cint, Ptr{Uint8}, Cint, Ptr{Void}),
        stmtptr,
        index,
        value,
        len,
        cb
    )
end


text:
using FactCheck
using Distributions
using PDMats
using StatsFuns
using MacroTools

using VinDsl
srand(56778)

facts("Utility functions") do
    context("advi_rand: distributions") do
        p = 10
        n = 6
        d = MvNormal(rand(p))
        x = VinDsl.advi_rand(d)
        @fact eltype(x) --> Float64
        @fact size(x) --> (p,)
        y = VinDsl.advi_rand(d, n)
        @fact size(y) --> (p, n)

        d = Normal(1, 2.)
        x = VinDsl.advi_rand(d)
        @fact size(x) --> ()
        y = VinDsl.advi_rand(d, n)
        @fact size(y) --> (n,)
    end

    context("advi_rand: vector") do
        p = 10
        n = 6
        x = Array{BigFloat}(rand(2p))
        r = VinDsl.advi_rand(x)
        @fact eltype(r) --> BigFloat
        @fact size(r) --> (p,)

        r = VinDsl.advi_rand(x, n)
        @fact size(r) --> (p, n)
    end

    context("randn!") do
        x = Array{ForwardDiff.Dual}(rand(10, 7))
        randn!(x)
        @fact eltype(x) <: ForwardDiff.Dual --> true
    end

    context("H: entropy") do
        p = 5
        vv = rand(p * (p + 3) ÷ 2)
        @fact H(vv, true) --> p * (log2π + 1)/2 + sum(vv[p + [1, 6, 10, 13, 15]])
        @fact H([1.1, 2.2]) --> (log2π + 1)/2 + 2.2/2
    end

end


text:
#
# Return the coordinates of the atoms
#

struct MaxMinCoords
  xmin::Vector{Float64}
  xmax::Vector{Float64}
  xlength::Vector{Float64}
end

"""

```
maxmin(atoms::Vector{Atom}; selection)
```

Returns the maximum and minimum coordinates of an atom vector, and the length (maximum minus minimum) in each direction.

### Example

```julia-repl
julia> protein = wget("1LBD");

julia> maxmin(protein)

 Minimum atom coordinates: xmin = [-29.301, 57.178, 45.668]
 Maximum atom coordinates: xmax = [47.147, 99.383, 86.886]
 Length in each direction: xlength = [76.448, 42.205, 41.217999999999996]

```

"""
function maxmin(atoms::AbstractVector{Atom}, selection::String)
  query = parse_query(selection)
  return maxmin(atoms, only = atom -> apply_query(query,atom))
end

function maxmin(atoms::AbstractVector{Atom}; only=all)
  xmin = [ +Inf, +Inf, +Inf ]
  xmax = zeros(3)
  for at in atoms
    if only(at)
      xmin .= (min(at.x,xmin[1]), min(at.y,xmin[2]), min(at.z,xmin[3]))
      xmax .= (max(at.x,xmax[1]), max(at.y,xmax[2]), max(at.z,xmax[3]))
    end
  end
  xlength = @. xmax - xmin
  return MaxMinCoords(xmin,xmax,xlength)
end

function Base.show(io::IO, m::MaxMinCoords)
  println(" ")
  println(" Minimum atom coordinates: xmin = ", m.xmin)
  println(" Maximum atom coordinates: xmax = ", m.xmax)
  println(" Length in each direction: xlength = ", m.xlength)
end


text:
# Autogenerated wrapper script for RLIBM_32_jll for powerpc64le-linux-gnu
export floatMathLib

JLLWrappers.@generate_wrapper_header("RLIBM_32")
JLLWrappers.@declare_library_product(floatMathLib, "floatMathLib.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        floatMathLib,
        "lib/floatMathLib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Day12

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day12PuzzleInfo <: PuzzleInfo end

"""
Day 12 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day12PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 12 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day12PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(cavename)

Whether cave with name `cavename` is big (returns `True`) or small (returns `False`).
"""
isbigcave(cavename::String) = uppercase(cavename) == cavename

"""
    $(FUNCTIONNAME)(input)

Read cave network from input stream `input`.
"""
function readcavenetwork(input::IO)
    connectedcaves = Dict{String, Set{String}}()
    for line in eachline(input)
        cave1, cave2 = split(line, '-')
        connectedcaves[cave1] = push!(get(connectedcaves, cave1, Set{String}()), cave2)
        connectedcaves[cave2] = push!(get(connectedcaves, cave2, Set{String}()), cave1)
    end
    connectedcaves
end

"""
    $(FUNCTIONNAME)(currentpath, connectedcaves, canrevisitsmallcave)

Recursively explore cave network described by `connectedcaves` (mapping from cave name
to names of connected caves), with current path of caves visited `currentpath` and
whether small caves are currently allowed to be revisited if already visited once
specified by `canrevisitsmallcave`.
"""

function explorecavenetwork!(
    currentpath::Vector{String},
    connectedcaves::Dict{String, Set{String}},
    canrevisitsmallcave::Bool,
)
    paths::Vector{Vector{String}} = []
    for cave in connectedcaves[currentpath[end]]
        if cave == "end"
            push!(paths, [currentpath; cave])
        elseif (
            isbigcave(cave)
            || cave ∉ currentpath
            || (canrevisitsmallcave && cave != "start")
        )
            append!(
                paths,
                explorecavenetwork!(
                    [currentpath; cave],
                    connectedcaves,
                    canrevisitsmallcave && (isbigcave(cave) || cave ∉ currentpath)
                )
            )
        end
    end
    return paths
end

function solve(info::Day12PuzzleInfo)
    connectedcaves = open(readcavenetwork, info.datapath)
    paths = explorecavenetwork!(["start"], connectedcaves, isa(info, Part2PuzzleInfo))
    length(paths)
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day12-1.txt"))) == 10
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day12-2.txt"))) == 19
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day12-3.txt"))) == 226
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day12-1.txt"))) == 36
end

end # module


text:
"""
Converts and returns an array of string to an array of Symbol
"""
function strarr_to_symarr!(arr)
    return arr = [ Symbol(str) for str in arr ]
end

"""
Converts and returns an equation string to an array of variables as string
"""
function equation_str_to_strarr!(equation)
    if occursin("~", equation)
        vars = split(replace(equation, r"\s+|\s+$/g" => " "), "~")
        equation = [String(strip(var)) for var in vcat(vars[1], split(vars[2], "+"))]
    else
        equation = [String(strip(var)) for var in split(replace(equation, r"\s+|\s+$/g" => ","), ",")]
    end
    equation = filter(x->length(x) > 0, equation)
    return equation
end

"""
Converts and returns a strarr equation with wildcards to a strarr with wildcards processed
"""
function equation_converts_wildcards!(equation, names)
    new_equation = []
    for e in equation
        e = replace(e, "." => "*")
        if e[end] == '*'
            datanames_arr = vec([String(key)[1:length(e[1:end - 1])] == e[1:end - 1] ? String(key) : nothing for key in names])
            append!(new_equation, filter!(x->x != nothing, datanames_arr))
        else
            append!(new_equation, [e])
        end
    end
    equation = unique(new_equation)
    return equation
end

"""
Gets datanames from data DataFrame or data Tuple
"""
function get_datanames_from_data(data, datanames)
    if isa(data, DataFrames.DataFrame)
        datanames = names(data)
    elseif isa(data, Tuple)
        datanames = data[2]
        if !isa(datanames, Vector)
            datanames = vec(datanames)
        end
    end
    return strarr_to_symarr!(datanames)
end

"""
Gets datanames from data DataFrame or data Tuple
"""
function get_data_from_data(data)
    if isa(data, DataFrames.DataFrame)
        data = convert(Matrix{Union{Missing, Float64}}, data)
    elseif isa(data, Tuple)
        data = data[1]
    end
    return data
end

"""
Sorts data
"""
function sort_data(data, datanames; time=nothing, panel=nothing)
    time_pos = ModelSelection.get_column_index(time, datanames)
    panel_pos = ModelSelection.get_column_index(panel, datanames)

    if time_pos != nothing && panel_pos != nothing
        data = sortslices(data, by=x->(x[panel_pos], x[time_pos]), dims=1)
    elseif panel_pos != nothing
        data = sortslices(data, by=x->(x[panel_pos]), dims=1)
    elseif time_pos != nothing
        data = sortslices(data, by=x->(x[time_pos]), dims=1)
    end
    return data
end

"""
Filter data by selected columns
"""
function filter_data_by_selected_columns(data, equation, datanames)
    columns = []
    for i = 1:length(equation)
        append!(columns, ModelSelection.get_column_index(equation[i], datanames))
    end
    data = data[:,columns]
    datanames = datanames[columns]

    return (data, datanames)
end

"""
Validates if there are panel gaps
"""
function validate_panel(data, datanames; panel=nothing)
    if panel != nothing
        return !any(ismissing, data[:, ModelSelection.get_column_index(panel, datanames)])
    end
    return true
end

"""
Validates if there are time gaps
"""
function validate_time(data, datanames; panel=nothing, time=nothing)
    nobs = size(data, 1)

    # TODO: Merge solutions
    if panel == nothing
        previous_value = data[1, ModelSelection.get_column_index(time, datanames)]
        for value in data[2:end, ModelSelection.get_column_index(time, datanames)]
            if previous_value + 1 != value
                return false
            end
            previous_value = value
        end
    else
        panel_index = ModelSelection.get_column_index(panel, datanames)
        csis = unique(data[:, panel_index])
        time_index = ModelSelection.get_column_index(time, datanames)
        for csi in csis
            rows = findall(x->x == csi, data[:,panel_index])
            previous_value = data[rows[1], time_index]
            for row in rows[2:end]
                value = data[row, time_index]
                if previous_value + 1 != value
                    return false
                end
                previous_value = value
            end
        end
    end
    return true
end

"""
Remove outliers from data
"""
function remove_outliers(data)
    for column in 1:size(data, 2)
        remove_outlier(data, column)
    end
end

"""
Remove outliers from a data column
"""
function remove_outlier(data, column)
    threshold = 3

    col = @view(data[:,column])
    aux_col = Array{Union{Int64, Float64, Missing}}(undef, size(col, 1), 2)
    for i in keys(col)
        aux_col[i, 1] = i
        aux_col[i, 2] = col[i]
    end

    valid_data = deleteat!(aux_col[:, 2], findall(ismissing, aux_col[:, 2]))

    mean_ = mean(valid_data)
    std_ = std(valid_data)

    for i in keys(col)
        if !ismissing(col[i])
            z_score = (col[i] - mean_) / std_
            if abs(z_score) > threshold
                col[i] = missing
            end
        end
    end
    return data
end

"""
Remove outliers from data
"""
function seasonal_adjustments(data, factor_dict, datanames)
    for column in factor_dict
        seasonal_adjustment(data, column[1], column[2], datanames)
    end
end

"""
Seasonal adjustment from a data column
"""
function seasonal_adjustment(data, name, factor, datanames)
    column = ModelSelection.get_column_index(name, datanames)
    nobs = size(data, 2)
    col = @view(data[:, column])
    L = Int(round(nobs / 2 / factor)) * factor
    yt, ys = analyze(col, L)
    seasonal_component = sum(ys, dims=2)
    col = col - seasonal_component
    return data
end

"""
Add values to extras
"""
function addextras(data, method, seasonaladjustment, removeoutliers)
    data.extras[ModelSelection.generate_extra_key(PREPROCESSING_EXTRAKEY, data.extras)] = Dict(
        :datanames => vcat(data.depvar, data.expvars),
        :depvar => data.depvar,
        :expvars => data.expvars,
        :data => DEFAULT_DATANAME,
        :method => method,
        :intercept => data.intercept,
        :panel => data.panel,
        :time => data.time,
        :seasonaladjustment => seasonaladjustment,
        :removeoutliers => removeoutliers,
        :removemissings => data.removemissings
    )
    return data
end


text:
# insert some elements: useful when one wants to modify θ by inserting
# new parameters
#

export insert_some_elements

# ****************************************************************

#
# How to do that with issorted is not clear:
#
# https://discourse.julialang.org/t/understanding-issorteds-lt-keyword/60918/42?page=3
#
# ... to avoid any missuses of issorted we define our own test
#
# This function returns true in case of _strict_ ordering, false otherwise
#
function is_strictly_ordered(indices::AbstractVector{Int})

    iter = iterate(indices)
    if iter == nothing return true end # empty sequence case

    iter_next = iterate(indices,last(iter))

    while iter_next != nothing
        first(iter) >= first(iter_next) && return false
        iter=iter_next
        iter_next=iterate(indices,last(iter_next))
    end

    true
end

# A priori, can easily be converted for arbitrary axes:
# - remplace i_dest by iterate(X_dest)
# - remplace i_src by iterate(X)
# then
# - remplace i_dest += 1 by iterate(X_dest,last(i_dest))
# - same idea for i_src += 1
#
@doc raw"""
```julia
insert_some_elements(X::AbstractVector{T},
                     indices::AbstractVector{Int},
                     elements::AbstractVector) -> AbstractVector{T}
```

Insert `elements` in `X` at positions `indices`.

Note: `indices` must be strictly ordered.

Example: this functions does a work inverse to the `deleteat!`

```jldoctest
julia> X = Any[1:10;]
10-element Vector{Any}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> indices = Int[2, 3, 5, 6, 9, 10];

julia> elements = Float64[2, 3, 5, 6, 9, 10];

julia> deleteat!(X, indices)
4-element Vector{Any}:
 1
 4
 7
 8

julia> insert_some_elements(X,indices,elements)
10-element Vector{Any}:
  1
  2.0
  3.0
  4
  5.0
  6.0
  7
  8
  9.0
 10.0

```
"""
function insert_some_elements end

# We first define a helper where elements type is free
# This allows to use the methods with *iterable* elements
#
function _insert_some_elements(X::AbstractVector{T},
                              indices::AbstractVector{Int},
                              elements) where {T} # elements must be iterable

    @assert is_strictly_ordered(indices)
    @assert length(indices)==0 || 1 ≤ minimum(indices) ≤ length(X)+length(indices)

    X_dest = Vector{T}(undef,length(X)+length(indices))

    i_dest = 1
    i_src = 1
    for (i,e) in zip(indices,elements)
        while i_dest<i
            X_dest[i_dest]=X[i_src]
            i_dest += 1
            i_src += 1
        end
        X_dest[i_dest]=e
        i_dest += 1
    end

    n = length(X)
    while i_src<=n
        X_dest[i_dest]=X[i_src]
        i_dest += 1
        i_src += 1
    end

    X_dest
end

# Specialize when elements is an AbstractVector
#
function insert_some_elements(X::AbstractVector{T},
                              indices::AbstractVector{Int},
                              elements::AbstractVector) where {T}

    @assert all(V->axes(V,1) isa Base.OneTo,(X,indices,elements)) # V[i], with i=1..n
    @assert length(indices)==length(elements)

    _insert_some_elements(X,
                          indices,
                          elements)
end

# Specialize when elements is a constant
#
function insert_some_elements(X::AbstractVector{T},
                              indices::AbstractVector{Int},
                              elements::T) where {T}

    _insert_some_elements(X,
                          indices,
                          Base.Iterators.repeated(elements))
end


text:
using EcosystemCore

g = Grass(1,5)
s = Sheep(2,1,1,1,1)
w = Wolf(3,1,1,1,1)
world = World([g,s,w])


text:
using Dates

@testset "Exp2.jl" begin
    filename = "data\\test.exp2"
    df = DDR2import.Exp2.read(filename)
    #EXP
    @test df[2,:ADEP] == "EHAM"
    @test df[2,:ADES] == "UUEE"
    @test df[2,:ACTYPE] == "A321"
    @test df[2,:RFL] == 350
    @test df[2,:ZONE_ORIG] == "EHAM"
    @test df[2,:ZONE_DEST] == "UUEE"
    @test df[2,:FLIGHT_ID] == 213793846
    @test Dates.year(df[2,:ETD_DATETIME]) == 2017
    @test Dates.month(df[2,:ETD_DATETIME]) == 12
    @test Dates.day(df[2,:ETD_DATETIME]) == 15
    @test Dates.hour(df[2,:ETD_DATETIME]) == 21
    @test Dates.minute(df[2,:ETD_DATETIME]) == 49
    @test df[2,:ETA_DATETIME] === missing
    @test Dates.day(df[1,:ETA_DATETIME]) == 15
    @test Dates.minute(df[1,:ETA_DATETIME]) == 0
    @test Dates.day(df[3,:ETA_DATETIME]) == 15
    @test Dates.minute(df[3,:ETA_DATETIME]) == 0
    @test df[2,:CALLSIGN] == "AFL2193"
    @test df[2,:COMPANY] == "AFL"

    #UUID
    @test df[2,:UUID] == "AFL2193-EHAM-UUEE-20171215213500"
    @test df[2,:FIPS_CLONED] == "N"

    #FLF
    #TODO test FLIGHT_SAAM_ID
    #TODO test FLIGHT_SAMAD_ID
    @test df[2,:TACT_ID] == 512114
    #TODO test SSR_CODE
    #TODO test REGISTRATION
    @test Dates.year(df[2,:PTD_DATETIME]) == 2017
    @test Dates.minute(df[2,:PTD_DATETIME]) == 35
    @test df[2,:ATFM_DELAY] == 0
    #TODO test REROUTING_STATE
    @test df[2,:MOST_PEN_REG] == "X"
    #TODO test TYPE_OF_FLIGHT
    #TODO test EQUIPMENT
    #TODO test ICAO_EQUIP
    #TODO test COM_EQUIP
    #TODO test NAV_EQUIP
    #TODO test SSR_EQUIP
    #TODO test SURVIVAL_EQUIP
    #TODO test PERSONS_ON_BOARD
    #TODO test TOP_FL
    @test df[2,:MAX_RFL] == 350
    #TODO test FLT_PLN_SOURCE

    #ALLFT
    @test Dates.minute(df[2,:AOBT]) == 40
    @test df[2, :IFPSID] == "AA70883446"
    #TODO test IOBT
    #TODO test ORIGFLIGHTDATAQUALITY
    #TODO test FLIGHTDATAQUALITY
    #TODO test SOURCE
    @test df[2, :EXEMPTREASON] == "NEXE"
    @test df[2, :EXEMPTREASONDIST] == "NEXE"
    #TODO test LATEFILER
    #TODO test LATEUPDATER
    #TODO test NORTHATLANTIC
    @test df[2, :COBT] === missing
    #TODO test COBT
    @test Dates.minute(df[2,:EOBT]) == 35
    #TODO test FLIGHTSTATE
    #TODO test PREV2ACTIVATIONFLIGHTSTATE
    #TODO test SUSPENSIONSTATUS
    #TODO test SAMCTOT
    #TODO test SAMSENT
    #TODO test SIPCTOT
    #TODO test SIPSENT
    #TODO test SLOTFORCED
    #TODO test MOSTPENALIZINGREGID
    #TODO test REGAFFECTEDBYNROFINST
    #TODO test EXCLFROMNROFINST
    #TODO test LASTRECEIVEDATFMMESSAGETITLE
    #TODO test LASTRECEIVEDMESSAGETITLE
    #TODO test LASTSENTATFMMESSAGETITLE
    #TODO test MANUALEXEMPREASON
    #TODO test SENSITIVEFLIGHT
    #TODO test READYFORIMPROVEMENT
    #TODO test READYFORDEP
    @test df[2, :REVISEDTAXITIME] == 840
    #TODO test TIS
    #TODO test TRS
    #TODO test TOBESENTSLOTMESSAGE
    #TODO test TOBESENTPROPMESSAGETITLE
    #TODO test LASTSENTSLOTMESSAGETITLE
    #TODO test LASTSENTPROPMESSAGETITLE
    #TODO test LASTSENTSLOTMESSAGE
    #TODO test LASTSENTPROPMESSAGE
    #TODO test FLIGHTCOUNTOPTION
    #TODO test NORMALFLIGHTTACT_ID
    #TODO test PROPFLIGHTTACT_ID
    @test df[2, :OPERATINGACOPERICAOID] == "AFL"
    #TODO test REROUTINGWHY
    #TODO test REROUTINGLFIGHTSTATE
    #TODO test RVR
    #TODO test FTFMAIRAC
    #TODO test FTFMENVBASELINENUM
    #TODO test RTFMAIRAC
    #TODO test RTFMENVBASELINENUM
    #TODO test CTFMAIRAC
    #TODO test CTFMENVBASELINENUM
    #TODO test LASTRECEIVEDPROGRESSMESSAGE

    #LAST TEST
    @test df[5,:ADEP] == "UUEE"
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# HitRecord for store the result of an intersection between a Ray and a Shape

"""
    HitRecord

A struct representing the result of an intersection between a [`Ray`](@ref) and a [`Shape`](@ref).

# Fields

- `world_point::Point`: a [`Point`](@ref) representing the world coordinates of the hit point.
- `normal::Normal`: a [`Normal`](@ref) representing the orientation of the normal to the surface where the hit happened.
- `surface_point::Vec2D`: a [`Vec2D`](@ref) representing the position of the hit point on the surface of the object.
- `t::Float32`: distance from the origin of the ray where the hit happened.
- `ray::Ray`: a [`Ray`](@ref) representing the the ray that hit the surface.
- `material::Material`: a [`Material`](@ref) representing the material of the point where the hit happened.
"""
struct HitRecord
    world_point::Point
    normal::Normal
    surface_point::Vec2D
    t::Float32
    ray::Ray
    material::Material
end

"""
    HitOrMiss

Alias for `Union{HitRecord, Nothing}`.
"""
const HitOrMiss = Union{HitRecord, Nothing}

function show(io::IO, ::MIME"text/plain", hr::T) where {T <: HitRecord}
    print(io, T)
    fns = fieldnames(T)
    n = maximum(fns .|> String .|> length)
    for fieldname ∈ fns
        println(io)
        print(io, " ↳ ", rpad(fieldname, n), " = ", getfield(hr, fieldname))
    end
end

"""
    ≈(hr1::HitRecord, hr2::HitRecord)

Check if two [`HitRecord`](@ref) represent the same hit event or not.
"""
function (≈)(hr1::HitRecord, hr2::HitRecord)
    sp1, sp2 = hr1.surface_point, hr2.surface_point
    hr1.world_point     ≈  hr2.world_point     &&
    hr1.normal          ≈  hr2.normal          &&
    findall(isnan, sp1) == findall(isnan, sp2) &&
    filter(!isnan, sp1) ≈  filter(!isnan, sp2) &&
    hr1.t               ≈  hr2.t               &&
    hr1.ray             ≈  hr2.ray
end

function Base.isless(hr1::HitRecord, hr2::HitRecord)
    hr1.t/norm(hr1.ray.dir) < hr2.t/norm(hr2.ray.dir)
end

function Base.isequal(hr1::HitRecord, hr2::HitRecord)
    isequal(hr1.t/norm(hr1.ray.dir), hr2.t/norm(hr2.ray.dir))
end

function Base.:(==)(hr1::HitRecord, hr2::HitRecord)
    hr1.t/norm(hr1.ray.dir) == hr2.t/norm(hr2.ray.dir)
end


text:
function zerotol!(M::AbstractArray; atol=eps(eltype(M)))
    @inbounds for i in eachindex(M)
        if abs(M[i]) < atol
            M[i] = zero(M[i])
        end
    end
    return M
end


text:
using Merlin

type MomentumSGD
    rate::Float64
    mu::Float64 # momentum coefficient
    momentum

    function MomentumSGD(rate, mu)
        new(rate, mu)
    end
end

function update!{T}(opt::MomentumSGD, value::Array{T}, grad::Array{T})
    if isdefined(opt, :momentum)
        opt.momentum = BLAS.scal(length(grad), T(opt.mu), grad, 1)
    else
        BLAS.axpy!(T(opt.mu), grad, opt.momentum)
    end
    BLAS.axpy!(-T(opt.rate), opt.momentum, value)
    fill!(grad, T(0))
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
export Airy

"""
	Airy(param;mesh,label)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the linear ([Airy](https://en.wikipedia.org/wiki/Airy_wave_theory)) water waves equations.

# Arguments
- `param::NamedTuple` must contain
    - the shallowness parameter `μ` (set the infinite-layer case if `μ=Inf`);
    - optionally, `ν` the shallow/deep water scaling factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise;
    - numerical parameters to construct the mesh of collocation points, if `mesh` is not provided.
- `mesh  :: Mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`.
- `label :: String`: a label for future references (default is `"linear (Airy)"`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `Airy.f!` to be called in explicit time-integration solvers;
2. a function `Airy.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `Airy.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct Airy <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info 	:: String

    function Airy(param::NamedTuple; # param is a NamedTuple containing all necessary parameters
		mesh = Mesh(param),
		IL	    = false,
		label = "linear (Airy)"  # using a keyword argument allows the user to supersede the default label.
		 )

		# Set up
		μ = param.μ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used thereafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "$label model.\n"
		if IL == true
			info *= "└─Infinite depth case.\n"
		else
			info *= "└─Shallowness parameter μ=$μ, \
					scaling parameter ν=$nu.\n"
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Collocation points and Fourier modes
		x, k = mesh.x, mesh.k
		# Fourier multipliers
    	∂ₓ	 = 1im * k            # Differentiation
		if IL == true
			∂ₓF₁ = 1im * sign.(k)
		else
			∂ₓF₁ = 1im * tanh.(√μ*k)/(√μ*ν)
		end
		# Pre-allocation
		fftη = zeros(Complex{Float64}, mesh.N)
        fftv = zeros(Complex{Float64}, mesh.N)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

			fftη .= U[:,1]
			fftv .= U[:,2]

		   	U[:,1] .= -∂ₓF₁.*fftv
		   	U[:,2] .= -∂ₓ.*fftη

		end

		# Build raw data from physical data (discrete Fourier transform)
		function mapto(data::InitialData)
			U = [fft(data.η(x)) fft(data.v(x))]
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2])),mesh.x
		end

        new(label, f!, mapto, mapfro, info )

    end
end


text:
"""
    check_return_value(value, location="")

Check the return value `value` of calls to C bindings of the SEAL library and throw an exception if
there were errors. `location` can be used to pass information about were the error occurred and will
be printed as well.

See also: [`@check_return_value`](@ref)
"""
function check_return_value(value, location="")
  if isempty(location)
    loc = ""
  else
    loc = " at $location"
  end

  if value == 0
    return true
  elseif value == 0x80004003
    error("invalid pointer/null pointer error in SEAL function (E_POINTER)", loc)
  elseif value == 0x80070057
    error("invalid argument error (E_INVALIDARG)", loc)
  elseif value == 0x8007000E
    error("out of memory error in SEAL function (E_OUTOFMEMORY)", loc)
  elseif value == 0x8000FFFF
    error("unexpected error in SEAL function (E_UNEXPECTED)", loc)
  elseif value == 0x80131620
    error("I/O error in SEAL function (COR_E_IO)", loc)
  elseif value == 0x80131509
    error("invalid operation error in SEAL function (COR_E_INVALIDOPERATION)", loc)
  else
    error("unknown error", loc)
  end
end

"""
    @check_return_value value

Call `check_return_value(value, location)` with `location` being a string of the format
`<filename>:<line_number>`.

See also: [`check_return_value`](@ref)
"""
macro check_return_value(value)
  return quote
    check_return_value($(esc(value)),
                       join([$(string(__source__.file)), $(string(__source__.line))], ":"))
  end
end


text:
import IJulia

let
    ijulia_kernel_file = joinpath(dirname(pathof(IJulia)), "kernel.jl")
    this_kernel_file = @__FILE__
    if ijulia_kernel_file != this_kernel_file
        @warn "this kernel.jl is different from the one provided by IJulia" this_kernel_file ijulia_kernel_file
    end
end

using InteractiveUtils

# workaround #60:
if Sys.isapple()
    ENV["PATH"] = Sys.BINDIR*":"*ENV["PATH"]
end

# the size of truncated output to show should not depend on the terminal
# where the kernel is launched, since the display is elsewhere
ENV["LINES"] = get(ENV, "LINES", 30)
ENV["COLUMNS"] = get(ENV, "COLUMNS", 80)

IJulia.init(ARGS)

let startupfile = abspath(homedir(), ".julia", "config", "startup_ijulia.jl")
    isfile(startupfile) && Base.JLOptions().startupfile != 2 && Base.include(Main, startupfile)
end

# import things that we want visible in IJulia but not in REPL's using IJulia
import IJulia: ans, In, Out, clear_history

pushdisplay(IJulia.InlineDisplay())

ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

println(IJulia.orig_stdout[], "Starting kernel event loops.")
IJulia.watch_stdio()

# workaround JuliaLang/julia#4259
delete!(task_local_storage(),:SOURCE_PATH)

# workaround JuliaLang/julia#6765
Core.eval(Base, :(is_interactive = true))

# check whether Revise is running and as needed configure it to run before every prompt
if isdefined(Main, :Revise)
    let mode = get(ENV, "JULIA_REVISE", "auto")
        mode == "auto" && IJulia.push_preexecute_hook(Main.Revise.revise)
    end
end

IJulia.waitloop()


text:
function reset!(dynamics_model::Simulator; r_tol=1e-5, κ_tol = 1e-3, undercut = 5.0)
    dynamics_model.ip.opts.r_tol = r_tol
    dynamics_model.ip.opts.κ_tol = κ_tol
    dynamics_model.ip.opts.undercut = undercut
    return nothing
end

function continuation_callback!(solver::Solver, dymamics_model::RoboDojo.Simulator)
    dynamics_model.ip.opts.r_tol = max(dynamics_model.ip.opts.r_tol/3, 1e-7)
    dynamics_model.ip.opts.κ_tol = max(dynamics_model.ip.opts.κ_tol/3, 1e-5)
    println("r_tol $(scn(dynamics_model.ip.opts.r_tol))  " *
        "κ_tol $(scn(dynamics_model.ip.opts.κ_tol))")

    # check that we are not too far from a solution
    solver_kickout = solver.data.max_violation[1] > 1e3
    return solver_kickout
end


text:
int
main() {
	printInt(fac(10));
	printInt(rfac(10));
	printInt(mfac(10));
        printInt(ifac(10));
        double r ; // just to test blocks
	{
	  int n = 10;
	  int r = 1;
	  while (n>0) {
	    r = r * n;
	    n = n - 1;
	  }
	  printInt(r);
	}
	printDouble(dfac(10.0));
	printString ("hello */");
        printString ("/* world") ;
        return 0 ;
}

int
fac(int a) {
	int r;
	int n;

	r = 1;
	n = a;
	while (n > 0) {
		r = r * n;
		n = n - 1;
	}
	return r;
}

int
rfac(int n)
{
	if (n == 0)
		return 1;
	else
		return n * rfac(n-1);
}

int
mfac(int n)
{
	if (n == 0)
		return 1;
	else
		return n * nfac(n-1);
}

int
nfac(int n)
{
	if (n != 0)
		return mfac(n-1) * n;
	else
		return 1;
}

double
dfac(double n)
{
	if (n == 0.0)
		return 1.0;
	else
		return n * dfac(n-1.0);
}

int ifac(int n) { return ifac2f(1,n); }

int ifac2f(int l, int h) {
        if (l == h)
                return l;
        if (l > h)
                return 1;
        int m;
        m = (l+h) / 2;
        return ifac2f(l,m) * ifac2f(m+1,h);
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fd8cab10-5458-11eb-3deb-dfba8e8d5152
begin
	using Pkg
	Pkg.activate()
	using HDF5
	using Plots
	using PlutoUI
	using LaTeXStrings
	using Printf
	import Pda
end

# ╔═╡ 8fa90742-5522-11eb-230f-89826748854e
file = "../examples/U1.000tauR0.2001D.h5"

# ╔═╡ c7b74ea0-5469-11eb-0ee4-73b2f96e3b07
@bind frame PlutoUI.Slider(0:99, show_value=true)

# ╔═╡ b336edc8-5522-11eb-009c-717fef97aff1
begin
	x = Pda.get_h5data(file, "config/$frame/x")
	q = Pda.get_h5data(file, "config/$frame/direction")
end

# ╔═╡ 5a678c92-55ff-11eb-1300-b1148bd3b9c8
begin
	bins=100
	gridx = LinRange(-0.5, 0.5, bins+1)
	bins = length(gridx) -1
    # initialize array to store number
    m = zeros(bins)
    # in a cell [a,b], particles that has a<=x<b is counted.
    # for the last cell, x==b is also counted.
    for i in 1:bins
        for (position, localq) in zip(x,q)
            if gridx[i] <=position <gridx[i+1]
                m[i] += localq
            end
        end
    end
    # last cell: add particles on the right boundary
	#m[end] += sum(x .== gridx[end])
    # use bin centers as representative locations.
    dx = diff(gridx)
    bin_centers = gridx[1:end-1] + dx/2
    # if normalize
    #     n ./= dx *length(x)
    # end
end

# ╔═╡ 305b5bb0-5459-11eb-1d3e-21bbe247b5ef
scatter(bin_centers[2:end-1], m[2:end-1],framestyle=:box,label=nothing)

# ╔═╡ Cell order:
# ╠═fd8cab10-5458-11eb-3deb-dfba8e8d5152
# ╠═8fa90742-5522-11eb-230f-89826748854e
# ╠═c7b74ea0-5469-11eb-0ee4-73b2f96e3b07
# ╠═b336edc8-5522-11eb-009c-717fef97aff1
# ╟─5a678c92-55ff-11eb-1300-b1148bd3b9c8
# ╠═305b5bb0-5459-11eb-1d3e-21bbe247b5ef


text:
using JuMP

"""
    hrep(model::JuMP.Model)

Builds an H-representation from the feasibility set of the JuMP model `model`.
Note that if non-linear constraint are present in the model, they are ignored.
"""
hrep(model::JuMP.Model) = LPHRepresentation(model)

function LPHRepresentation(model::JuMP.Model)
    # Inspired from Joey Huchette's code in ConvexHull.jl
    A = JuMP.prepConstrMatrix(model)
    c = JuMP.prepAffObjective(model)
    lb, ub = JuMP.prepConstrBounds(model)
    l, u = model.colLower, model.colUpper

    m, n = size(A)
    @assert m == length(lb) == length(ub)
    @assert model.nlpdata == nothing
    @assert isempty(model.quadconstr)
    @assert isempty(model.sosconstr)

    LPHRepresentation(A, l, u, lb, ub)
end

function polyhedron(model::JuMP.Model, lib::PolyhedraLibrary=default_library(FullDim{model.numCols}(), Float64))
    polyhedron(LPHRepresentation(model), lib)
end


text:
"""
    MapIn(layers...)
Maps input through multiple layers / functions, so it is broadcasted
over.
# Examples
```jldoctest
julia> m = MapIn(Chain(Dense(1, 4), Dense(4, 1)));
julia> size(m(Float32[1, 2, 3])) == (3,)
true
```
"""
struct MapIn{T<:Tuple}
    layers::T
    MapIn(xs...) = new{typeof(xs)}(xs)
end
applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(Base.tail(fs), first(fs)(x))
(a::MapIn)(x) = broadcast(y -> applychain(a.layers, y), x)
(a::MapIn)(x::Array{Float32,1}) = broadcast(first, a(map(y -> [y], x)))
Flux.trainable(a::MapIn) = a.layers

#Base.getindex(x::MapIn, args...) = Base.getindex(x.layers)
Base.length(x::MapIn, args...) = Base.length(x.layers)
#Base.first(x::MapIn, args...) = Base.first(x.layers)
#Base.last(x::MapIn, args...) = Base.last(x.layers)
#Base.iterate(x::MapIn, args...) = Base.iterate(x.layers)
#Base.lastindex(x::MapIn, args...) = Base.lastindex(x.layers)

function Base.show(io::IO, c::MapIn)
  print(io, "MapIn(")
  join(io, c.layers, ", ")
  print(io, ")")
end

max_stars_repo_name: williamjsdavis/ThermoChroNN

max_issues_repo_name: williamjsdavis/ThermoChroNN

max_forks_repo_name: williamjsdavis/ThermoChroNN


text:
function add!(window::GtkWindow, ::Type{Area})
    function callback_w(widget::Gtk.GtkGLArea, width::Int32, height::Int32)
        rect = window[Area]
        window[Area] = IRect(minimum(rect), width, height)
        return true
    end
    signal_connect(callback_w, window[Window], "Scroll")
    add_events(window[Window], GConstants.GdkEventMask.SCROLL)
    return
end

function add!(window::GtkWindow, ::Type{Mouse.Scroll})
    function callback(widget::Gtk.GtkGLArea, s::Gtk.GdkEventScroll)
        window[Scroll] = (s.x, s.y)
        return true
    end
    signal_connect(callback, window[NativeWindow], "scroll-event")
    add_events(window[NativeWindow], GConstants.GdkEventMask.SCROLL)
    return
end

function add!(window::GtkWindow, ::Type{Mouse.Position})
    function callback(widget::Gtk.GtkGLArea, s::Gtk.GdkEventMotion)
        window[Mouse.Position] = (s.x, s.y)
        return true
    end
    add_events(window[NativeWindow], GConstants.GdkEventMask.POINTER_MOTION)
    signal_connect(callback, window[NativeWindow], "motion-notify-event")
    return true
end

function to_mouse_button(x)
    if x == 1
        Mouse.left
    elseif x == 2
        Mouse.middle
    elseif x == 3
        Mouse.right
    else
        # TODO turn into error
        warn("Button is $x, while $(Gtk.GdkModifierType.BUTTON1)")
        Mouse.left
    end
end
function add!(window::GtkWindow, ::Type{Mouse.Buttons})
    function callback(widget::Gtk.GtkGLArea, event::Gtk.GdkEventButton)
        button = to_mouse_button(event.button)
        action = event.event_type
        set = window[Mouse.Buttons]
        if action in (GdkEventType.BUTTON_PRESS:GdkEventType.TRIPLE_BUTTON_PRESS)
            push!(set, button)
        elseif action == GdkEventType.BUTTON_RELEASE
            delete!(set, button)
        else
            warn("unknown action: $(action)")
        end
        window[Mouse.Buttons] = set # trigger setfield event!
        return true
    end
    add_events(window[NativeWindow],
        GConstants.GdkEventMask.GDK_BUTTON_PRESS_MASK |
        GConstants.GdkEventMask.GDK_BUTTON_RELEASE_MASK
    )
    signal_connect(callback, window[NativeWindow], "button_press_event")
    signal_connect(callback, window[NativeWindow], "button_release_event")
    return
end


text:
# Command: Maximize Window
"""
	maximize!(session::Session; window::AbstractString = "current")::Dict{String,Int}

The Maximize Window command invokes the window manager-specific "maximize" operation, if any, on the window containing the current top-level browsing context.
This typically increases the window to the maximum available size without going full-screen.

!!! note

    This command is implemented using the JSON Wire Protocol.
"""
function maximize!(session::Session; window::AbstractString = "current")
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/window/$window/maximize",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    output = JSON3.read(response.body).value
    (
        width = output.width,
        height = output.height,
        x = output.x,
        y = output.y,
    )::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}
end


text:
using Gadfly

set_default_plot_size(6inch,3inch)

x, y = [1, 2, 3, 4], [3, 7, 5, 1]
xmin1, xmax1 = x.-0.1, x.+0.1
xmin2, xmax2 = x.-0.5, x.-0.3

lyr1 = layer(xmin=xmin1, xmax=xmax1, y=y, Geom.bar, color=[1])
lyr2 = layer(xmin=xmin2, xmax=xmax2, y=y, Geom.bar, color=[2])
cs = repeat([colorant"violet", colorant"hotpink"], outer=2)

p1 = plot(Theme(point_size=3pt), lyr1, lyr2,
    layer(x=x, y=y.+0.9, Geom.point, color=cs),
    Scale.color_discrete
)
p2 = plot(Theme(point_size=3pt), lyr1, lyr2,
    layer(x=x, y=y.+0.9, Geom.point, color=cs),
    Scale.color_continuous
)

hstack(p1, p2)


text:
using HITRAN
using Test

@testset "HITRAN.jl" begin
    @testset "Total internal partition sums" begin
        for i in HITRAN.isotopologues[:, :global_id]
            q_t0 = HITRAN.isotopologues[HITRAN.isotopologues.global_id .== i, :q_t0][1]
            tolerance = 10.0^(ceil(Int, log10(q_t0))-3)
            @test isapprox(tips(i, HITRAN.c_T_ref), q_t0, atol=tolerance)
        end
    end
    @testset "Database and absorption calculation" begin
        @test_nowarn fetch!("StdAtm", iso_id(7, 1), 13000, 13150, [:standard])
        @test isequal(iso_id(["H2O"])[1], 1)
        @test isequal(iso_id(7, 1)[1], 36)
        @test_nowarn α(["StdAtm"])
        @test_nowarn α(["StdAtm"], :voigt)
        @test_nowarn α(["StdAtm"], :sdvoigt)
        @test_nowarn α(["StdAtm"], :lorentz)
        @test_nowarn α(["StdAtm"], :gauss)
        ν, absorption = α(["StdAtm"])
        @test_nowarn apply_instrument_function(ν, absorption, :gaussian, 1.0, 0.1)
        @test isapprox(absorption_spectrum([0.], 1)[1], 0.0)
        @test isapprox(transmittance_spectrum([0.], 1)[1], 1.0)
        @test isapprox(optical_depth([0.], 1)[1], 0.0)
    end

    @testset "Instrument functions" begin
        x = -10:0.001:10
        for (s, fn) in HITRAN.instrument_functions
            @test isapprox(sum(fn(x, 0.1)*0.001), 1, atol=1e-2)
        end
    end

    @testset "Environments" begin
        @test isapprox(HITRAN.p_s_h2o(273.15), 611.2911778902558)
        @test isapprox(HITRAN.kelvin_to_celsius(273.15), 0)
        @test isapprox(HITRAN.celsius_to_kelvin(0), 273.15)
        @test_nowarn moist_air(100)
        @test_nowarn moist_air(100, HITRAN.c_p_ref, 240.0)
    end

    @testset "Utility" begin
        @test isapprox(wavelength_to_wavenumber(1e-6), 10000.0)
        @test isapprox(frequency_to_wavenumber(30e9), 1, atol=1e-3)
    end
end


text:
result = 0

for i = 1:1000
    if ( i % 3 == true ) || ( i % 5 == true )
       global result = result + i
    end
end
println("Result: $result")

max_stars_repo_name: marvins/ProjectEuler

max_issues_repo_name: marvins/ProjectEuler

max_forks_repo_name: marvins/ProjectEuler