text:
using Base
using Interpolations

# While the StateID struct and its helper functions may seem a little over-engineered, they
# primarily exist to document the iSLP and iLV quantities used by TOPBase

"""
    StateID(iSLP, iLV)

`StateID` represents a given state (electron configuration) of an ion species. The arguments are
taken from the TOPBase.

# Arguments
- `iSLP::UInt16`: An integer representing the quantum atomic numbers SLπ (total spin, orbital
  angular momentum, and parity). The value is given by `100*(2*S+1) + 10*L + π`. Note that TOPBase
  documentation refers to the parity quantum number as P or π, but for consistency with how TOPBase
  presents the data, we prefer P in our implementation code.
- `iLV::UInt8`: Distinguishes between configurations that have the same values of `iSLP`. TOPBase
  seems to index the configuration based on the energy level of the configuration (with `iLV = 1`
  corresponding to the lowest energy).

# Notes
The combination of this struct and an ion species (or neutral element) should be enough to fully
specify a given electron configuration. For example, consider He I:
- `iSLP = 100, iLV = 1` has the electron configuration: 1s²
- `iSLP = 100, iLV = 2` has the electron configuration: 1s 2s
- `iSLP = 100, iLV = 3` has the electron configuration: 1s 3s
- `iSLP = 111, iLV = 1` has the electron configuration: 1s 2p
- `iSLP = 111, iLV = 2` has the electron configuration: 1s 3p
"""
struct StateID
    iSLP::UInt16
    iLV::UInt8
    function StateID(iSLP, iLV)
        @assert 99 < iSLP < 992 # this upper limit is set based on the data source
        @assert (iSLP % 10) < 2 # trailing digit must be 0 or 1
        @assert 0 < iLV # the upper limit is set based on the data source
        new(iSLP, iLV)
    end
end

Base.:(<)(a::StateID, b::StateID) = (a.iSLP < b.iSLP) || ((a.iSLP == b.iSLP) && (a.iLV < b.iLV))
Base.show(io::IO, obj::StateID) = print(io, "StateID{iSLP = ", obj.iSLP, ", iLV = ", obj.iLV, "}")

_parity_vals = ("even", "odd") # predefine these to avoid allocations

function _quantum_prop_type(state_id::StateID)
    S_code = div(state_id.iSLP, 100) # hundreds place of state_id.iSLP
    twice_S = S_code - 1
    L = div(state_id.iSLP - S_code*100, 10) # tens place of state_id.iSLP
    parity = _parity_vals[(state_id.iSLP % 2) + 1] # 0 or 1 in ones place of state_id.iSLP
    (Int32(twice_S), Int32(L), parity)
end
get_S_quantum_num(state_id::StateID) = _quantum_prop_type(state_id)[1] / 2.0f0
get_L_quantum_num(state_id::StateID) = _quantum_prop_type(state_id)[2]
get_parity(state_id::StateID) = _quantum_prop_type(state_id)[3]

function get_statistical_weight(state_id::StateID)
    twice_S, L, _ = _quantum_prop_type(state_id)
    (twice_S+1)*(2*L+1)
end

"""
    StateProp(ion_energy_ryd, excitation_potential_ryd, statistical_weight, electron_config)

`StateProp` holds data associated with different electron configurations. It's used to hold
data read from the TOPBase cross-section queries and to compute the statistical weight.
"""
struct StateProp
    ion_energy_ryd::Float64 # energy with respect to ionization potential in Rydbergs
    excitation_potential_ryd::Float64 # energy above ground state in Rydbergs
    statistical_weight::UInt16
    electron_config::String
end

# define some generally useful functions related to parsing
function _parse_species_stateid(str::AbstractString, start_Z::Integer,
                                start_numElectrons::Integer, start_iSLP::Integer,
                                start_iLV::Integer, last_iLV::Integer)
    # parses the atomic_num, number of electrons, iSLP, and iLV entries of the table
    atomic_num = parse(UInt8, str[start_Z:start_numElectrons-1])
    num_electrons = parse(UInt8, str[start_numElectrons:start_iSLP-1])
    iSLP = parse(UInt16, str[start_iSLP:start_iLV-1])
    iLV = parse(UInt16, str[start_iLV:last_iLV])
    (Species(Formula(atomic_num), atomic_num - num_electrons), StateID(iSLP, iLV))
end

# Functionallity for reading photoionization cross sections. Since these tables can be massive, we
# use an iterator to do this. This is based on the design used in the Julia standard library for
# implementing the ``EachLine`` iterable object. That code used the following license:
#
# Copyright (c) 2009-2021: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
# and other contributors:
#
# https://github.com/JuliaLang/julia/contributors
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

const _TABLE_BREAK = " ================================================"
const _HEADER_LINE_LENGTH = length(_TABLE_BREAK)

struct EachPhotoIonSubtable{IOT <: IO}
    stream::IOT
    ondone::Function
end

"""
    each_photo_ion_subtable(filename::AbstractString)

Create an iterable `EachPhotoIonSubtable` object for parsing photoionization cross-sections.

The input file is organized as a series of subtables. Each subtable provides the photoionization
energies (make sure this isn't the energy of the photoionized electron) in Ryd, and the associated
cross section (in megabarns) for a different state.

Each entry in the `EachPhotoIonSubtable` iterable provides a named tuple with the keys:
- name: holds the species name
- state_id: holds a `StateID` instance
- photon_energy_ryd: the photoionization values (check), in Ryd, stored in an array of Float32 vals
- cross_section_MegaBarn: the associated photoionization cross-sections. This array should have the
  same shape and datatype as photon_energy_ryd

The underlying file is only closed after the `EachPhotoIonSubtable` is garbage collected.

#Notes
This implementation is not robust enough to deal with blank lines at the start or end of the file.

The table should have been produced by a
[TOPbase Photoionisation Cross Section query](http://cdsweb.u-strasbg.fr/topbase/xsections.html)
that uses the following query options:
- Only a single atomic number and a single number of electrons is included in the query (other
  choices may yield too many results, which causes TOPBase to truncate output).
- There isn't any limitation on the number the quantum numbers, the energy or the levels
- The query request levels ordered in the "Level order in each series"
"""
function each_photo_ion_subtable(filename::AbstractString)
    stream = open(filename, "r")
    # confirm that standard file header is present and move the position of the
    # stream past this point
    @assert readline(stream) == _TABLE_BREAK
    @assert readline(stream) == "       I  NZ  NE  ISLP  ILV        E(RYD)      NP"
    @assert readline(stream) == _TABLE_BREAK
    EachPhotoIonSubtable(stream, ()->close(stream))::EachPhotoIonSubtable
end

function Base.iterate(itr::EachPhotoIonSubtable, state=nothing)
    next_line = readline(itr.stream)

    if next_line == _TABLE_BREAK
        return (itr.ondone(); nothing)
    end

    # first parse the subtable's header
    species, state_id = _parse_species_stateid(next_line, 9, 13, 17, 23, 27)
    ion_energy_ryd = parse(Float32, next_line[28:41])

    num_entries = parse(UInt32, next_line[42:49])

    # this could probably be 32-bit.
    photon_energy_ryd = Array{Float32,1}(undef, num_entries)
    cross_section_MegaBarn = Array{Float32,1}(undef, num_entries)
    for i = 1:num_entries # this properly handles the case when num_entries == 0
        cur_line = readline(itr.stream)
        photon_energy_ryd[i] = parse(Float32, cur_line[1:14])
        cross_section_MegaBarn[i] = parse(Float32, cur_line[15:24])
    end

    # parse the species, and the level
    ((species = species, state_id=state_id, ion_energy_ryd = ion_energy_ryd,
      photon_energy_ryd_arr = photon_energy_ryd,
      cross_section_MegaBarn_arr = cross_section_MegaBarn), nothing)
end

# this and the following definition implement Julia's informal iterator interface
Base.eltype(::Type{<:EachPhotoIonSubtable}) =
    Tuple{String, StateID, Float32, Array{Float32,1}, Array{Float32,1}}
Base.IteratorSize(::Type{<:EachPhotoIonSubtable}) = SizeUnknown()

# In a couple tables, there's a place where the photon energy is listed out of order. This seems to
# be okay given that this is in a section of the table where the cross-section is zero.
const _species_with_problematic_subtables = (Species("He_I"), Species("C_I"), Species("C_II"),
                                             Species("Al_II"), Species("O_II"))

function _interpolator_from_subtable(table_photon_energy_ryd::Vector{Float32},
                                     table_cross_section_MBarn::Vector{Float32},
                                     species::Union{Species,Nothing} = nothing,
                                     extrapolation_bc = 0.0)
    if species in _species_with_problematic_subtables

        for j in 1:length(table_photon_energy_ryd)-1
            if table_photon_energy_ryd[j+1] < table_photon_energy_ryd[j]
                @assert j > 1
                @assert (j + 2) <= length(table_photon_energy_ryd)
                @assert table_photon_energy_ryd[j+2] >= table_photon_energy_ryd[j+1]
                @assert table_photon_energy_ryd[j+2] >= table_photon_energy_ryd[j]
                @assert table_photon_energy_ryd[j+1] >= table_photon_energy_ryd[j-1]
                @assert table_photon_energy_ryd[j] >= table_photon_energy_ryd[j-1]

                @assert table_cross_section_MBarn[j+2] == 0.0
                @assert table_cross_section_MBarn[j+1] == 0.0
                @assert table_cross_section_MBarn[j] == 0.0
                @assert table_cross_section_MBarn[j-1] == 0.0

                # swap the order of the 2 values
                temp = table_photon_energy_ryd[j]
                table_photon_energy_ryd[j] = table_photon_energy_ryd[j+1]
                table_photon_energy_ryd[j+1] = temp
            end
        end
    end
    LinearInterpolation(table_photon_energy_ryd, table_cross_section_MBarn,
                        extrapolation_bc=extrapolation_bc)
end

function _try_env_prefix_path(prefix_env_var, fname)
    prefix = get(ENV, prefix_env_var, nothing)
    if isnothing(prefix)
        e = ErrorException(string("The \"", prefix_env_var, "\" environment variable was not set"))
        throw(e);
    end
    path = joinpath(prefix, fname)
    if !isfile(path) # this properly handles the case where path is a symlink to a file
        throw(ErrorException(string("There is no file at ", path)));
    end
    path
end

# return the minimum iterator over cross_section subtables and the ionization energy of the ground
# state (according to the table)
function _get_cross_sec_data(species::Species,
                             cross_sec_file::Union{AbstractString,Nothing} = nothing)

    # get the photo-ionization subtable iterator:
    my_cross_sec_file = if isnothing(cross_sec_file)
        # determine the string representation of Species
        # we don't directly use repr(species) since the representation could change in the future
        # to not use roman numeral (e.g. if support is added for negatively charged species)
        formula_str = repr("text/plain", species.formula)
        roman_numeral = get_roman_numeral(species)
        fname = string(formula_str, "_", roman_numeral, ".txt")

        println("cross_sec_file isn't specified. Searching for file in the directory specified ",
                "by the KORG_OP_DIR environment variable.")
        _try_env_prefix_path("KORG_OP_DIR", fname)
    else
        cross_sec_file
    end

    # determine the ionization energy of the ground state (from the subtables)
    gs_ion_energy_ryd = begin # this could be optimized significantly!
        cur_min = floatmax(Float32)
        for (cur_species, _, ion_energy_ryd, _, _) in each_photo_ion_subtable(my_cross_sec_file)
            @assert (cur_species == species) # this will catch mistakes in the query that generated
                                             # the cross_section data file
            cur_min = min(ion_energy_ryd, cur_min)
        end
        cur_min
    end

    subtable_iter = each_photo_ion_subtable(my_cross_sec_file)
    subtable_iter, gs_ion_energy_ryd
end

@doc raw"""
    weighted_bf_cross_section_TOPBase(λ_vals, T_vals, species; extrapolation_bc = 0.0
                                      cross_sec_file = nothing,
                                      partition_func = partition_funcs[species])

Calculate the combined bound-free cross section (in cm²) of an ion species using data from the
Opacity Project. The result is the mean of the the bound-free cross section of each electron state,
weighted by their Boltzmann probabilities. It includes the LTE correction for stimulated emission.
The result is stored in a 2D array with a size given by `(length(λ_vals), length(T_vals))`.

# Arguments
- `λ_vals::AbstractVector`: wavelengths (in cm) to compute the cross sections at
- `T_vals::AbstractVector`: temperatures (in K) to compute the cross sections at
- `species::Species`: the starting atom or atomic ion of the bound-free reaction.

# Keyword Arguments
- `extrapolation_bc`: Specifies handling of extrapolation during linear interpolation (this is
  passed to `Interpolations.LinearInterpolation`).
- `cross_sec_file`: Optional keyword that specifies the path to the text file that holds the
  cross-section subtables. When `nothing` (the default), this function tries to load the file
  from `$KORG_OP_DIR/<atomic-symbol>_<ion-state>.txt`. `<atomic-symbol>` is replaced with the
  atomic symbol of the species (e.g. `H`, `He`, `Al`) and `<ion-state>` is replaced with the
  capital Roman numeral(s) specifying the ionization state.
- `partition_func`: Specifies the partition function for the current species. This should be a
  callable that accepts temperature as an argument. When this is `nothing`, it falls back to the
  default partition functions.

# Explanation
In more mathematical rigor, this function basically evaluates the following equation:

``\overline{\sigma_{\rm bf}}(\lambda,T) = \sum_{{\bf n}\in{\rm states}} \frac{ g_{\bf n}}{U(T)} e^{-E_{\bf n}/(k_B T)}\left(1 - e^{-h c/(\lambda k_B T)}\right) \sigma_{{\rm bf}, {\bf n}}(\lambda)``

In this equation:
- ``g_{\bf n}`` is the statistical weight of state ``{\bf n}``.
- ``U`` is the partition function.
- ``E_{\bf n}`` is the excitation potential of state ``{\bf n}``. In other words, its the energy
  difference between state ``{\bf n}`` and the ground state.
- ``\sigma_{\lambda,{\rm bf}, {\bf n}}`` is the bound-free cross-section (a.k.a. photo-ionization
  cross section) of state ``{\bf n}``.

The linear aborption coefficient for bound-free absorption is simply:
``\alpha_{\lambda,{\rm bf}} = n_{\rm tot} \overline{\sigma_{\rm bf}}(\lambda,T)``.

# Notes
In the future, it needs to be possible to pass an argument that adjusts the energy levels of
different electron configurations.
"""
function weighted_bf_cross_section_TOPBase(λ_vals::AbstractVector, T_vals::AbstractVector,
                                           species::Species; extrapolation_bc=0.0,
                                           cross_sec_file::Union{AbstractString,Nothing} = nothing,
                                           partition_func = partition_funcs[species])
    # Note: it would not take a lot of work to make this work with 32-bit floats (since the tables
    # only store 32-bit floats). If doing that, it might be wise to accept λ in units of ångstroms
    # (for precision purposes)

    # TODO: make it possible to pass in a dict mapping dict names to the various energy levels
    #       (to allow for overwriting the values in the table with empirical values)

    @assert !ismolecule(species) # species must be an atom or atomic ion

    # precompute Temperature-dependent constant
    inv_partition_func_val = 1.0./partition_func.(T_vals)

    β_Ryd = RydbergH_eV./(kboltz_eV .* T_vals)

    # convert λ_vals to photon energies
    photon_energies = (hplanck_eV * c_cgs / RydbergH_eV) ./ λ_vals

    # prepare the output array where results will be accumulated
    weighted_average = zeros(eltype(photon_energies), (length(photon_energies),length(T_vals)))

    cross_section_subtables, gs_ion_energy_ryd = _get_cross_sec_data(species, cross_sec_file)

    for subtable in cross_section_subtables
        if length(subtable.photon_energy_ryd_arr) == 0
            continue
        end

        cross_section_interp = _interpolator_from_subtable(subtable.photon_energy_ryd_arr,
                                                           subtable.cross_section_MegaBarn_arr,
                                                           species, extrapolation_bc)
        statistical_weight = get_statistical_weight(subtable.state_id)
        excitation_potential_ryd = subtable.ion_energy_ryd - gs_ion_energy_ryd

        # now iterate over Temperatures
        for j in 1:length(T_vals)

            # now compute the weighting of the current state under LTE
            weight = inv_partition_func_val[j] * statistical_weight * exp(-excitation_potential_ryd
                                                                          * β_Ryd[j])

            current_cross_section = ( cross_section_interp.(photon_energies) .*
                                      (1.0 .- exp.(-photon_energies.*β_Ryd[j])) )
            view(weighted_average, :, j) .+= (current_cross_section .* weight)
        end
    end

    # convert weighted_average units from megabarn to cm²
    weighted_average .*= 1e-18
    weighted_average
end

"""
    absorption_coef_bf_TOPBase(λ, T, ndens_species, species; kwargs...)

Computes the bound-free linear absorption coefficient, α, (in units of cm⁻¹) using data from the
Opacity Project (downloaded from TOPBase) and assuming LTE. The result is stored in a 2D array with
a size given by `(length(λ), length(T))`.

# Arguments
- `λ::AbstractVector`: wavelengths (in cm) to compute α at.
- `T::AbstractVector`: temperatures (in K) to compute α at.
- `ndens_species::AbstractVector`: number density (in cm⁻³) of the species at each temperature
  (this should have the same length as T).
- `species::Species`: the starting atom or atomic ion of the bound-free reaction.

# Keyword Arguments
This function accepts all the same keyword arguments as `weighted_bf_cross_section_TOPBase`.
"""
function absorption_coef_bf_TOPBase(λ::AbstractVector, T::AbstractVector,
                                    ndens_species::AbstractVector, species::Species; kwargs...)
    @assert size(T) == size(ndens_species)

    arr = weighted_bf_cross_section_TOPBase(λ, T, species; kwargs...)

    # to save time and memory, let's modify the values of arr in-place
    for (i,n) in enumerate(ndens_species)
        view(arr,:,i) .*= n
    end

    arr
end


text:
Base.@kwdef struct PlanarRocket{T} <: AbstractModel
    g::T = 9.81   # gravity (m/s²)
    m::T = 100.0  # mass (kg)
    ℓ::T = 2.0    # length (m)
    max_angle::T = 10.0  # max thrust angle (deg)
    max_roll::T = 10.0   # max roll angle (deg)
    Tlo::T = 0.75        # min thrust, per unit weight (1/kg⋅m/s²)
    Thi::T = 1.50        # max thrust, per unit weight (1/kg⋅m/s²)
end
state_dim(::PlanarRocket) = 8
control_dim(::PlanarRocket) = 2

umin(model::PlanarRocket) = SA[model.Tlo*model.m*model.g, -deg2rad(model.max_angle)]
umax(model::PlanarRocket) = SA[model.Thi*model.m*model.g, +deg2rad(model.max_angle)]

function dynamics(model::PlanarRocket, x, u)
    mass,g,ℓ = model.m, model.g, model.ℓ

    # clamp the controls
    # u = clamp.(u, umin(model), umax(model))

    ulo = umin(model)
    uhi = umax(model)
    T = clamp(x[7],ulo[1],uhi[1])   # thrust (N)
    ϕ = clamp(x[8],ulo[2],uhi[2])  # gimbal angle (rad)

    J = mass*ℓ^2 / 12  # moment of inertia
    θ = x[3]

    # acceleration
    s,c = sincos(θ + ϕ)
    a = SA[
        (T/mass) * s,
        (T/mass) * c - g,
        (T/J) * 0.5 * ℓ * sin(ϕ)
    ]
    return SA[x[4], x[5], x[6], a[1], a[2], a[3], u[1], u[2]]
end

ext: jl


text:
export data

abstract type TensorStorage end

data(S::TensorStorage) = S.data


text:
push!(LOAD_PATH,"../src/")

using NESSie
using Test

# All tests. Will be used if no command line arguments are given.
const tests = [
    "base/constants.jl",
    "base/model.jl",
    "base/potentials.jl",
    "base/quadrature.jl",
    "base/util.jl",
    "bem/local.jl",
    "bem/nonlocal.jl",
    "format/hmo.jl",
    "format/mcsf.jl",
    "format/msms.jl",
    "format/off.jl",
    "format/pqr.jl",
    "format/skel.jl",
    "format/stl.jl",
    "format/vtk.jl",
    "format/xml3d.jl",
    "testmodel/born.jl",
    "testmodel/xie.jl",
    "radon.jl",
    "rjasanow.jl"
]
const testtypes = (Float64, Float32)

function printusage()
    println("\n\e[1mNESSie.jl test suite\e[0m")
    println("====================\n")
    println("\e[1mUsage:\e[0m\n")
    println("\truntests.jl [--help] [<test>...]\n")
    println("If called without arguments, this program will run all available tests.")
    println("\n\e[1mOptions:\e[0m\n")
    println("\t\e[1m--help\e[0m\tPrint this message and exit")
    println("\n\e[1mAvailable tests:\e[0m\n")
    for test in tests
        println("\t$(test[1:end-3])")
    end
    println()
end

#=
    Validate args and return test files
=#
function checkargs(args)
    ret = String[]
    err = String[]
    for arg in ARGS
        if arg == "--help"
            printusage()
            exit(0)
        end
        if arg[1] == '-'
            push!(err, "\e[1;31mInvalid option \e[39m$arg\e[0m")
            continue
        end
        if arg in tests || "$(arg).jl" in tests
            push!(ret, endswith(arg, ".jl") ? arg : "$(arg).jl")
        else
            push!(err, "\e[1;31mInvalid test file \e[39m$arg\e[0m")
        end
    end
    if length(err) > 0
        println("\e[1;31mERROR:\e[0m\n")
        for e in err
            println(" * $e")
        end
        print("\n\e[1;31mRun\e[39m runtests.jl --help")
        println(" \e[31mfor a list of available options and tests.\e[0m")
        exit(1)
    end
    ret
end

# Run tests
include("testutils.jl")
if length(ARGS) > 0
    @testset "$t" for t in checkargs(ARGS)
        include(t)
        println()
    end
else
    @testset "NESSie.jl" begin
        println()
        @testset "$t" for t in tests
            include(t)
        end
    end
end


text:
immutable ViolinGeometry <: Gadfly.GeometryElement
    order::Int

    function ViolinGeometry(; order=1)
        new(order)
    end
end

const violin = ViolinGeometry

element_aesthetics(::ViolinGeometry) = [:x, :y, :color]

default_statistic(::ViolinGeometry) = Gadfly.Stat.violin()

function render(geom::ViolinGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    # TODO: What should we do with the color aesthetic?

    Gadfly.assert_aesthetics_defined("Geom.violin", aes, :y, :width)
    Gadfly.assert_aesthetics_equal_length("Geom.violin", aes, :y, :width)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(RGB{Float32}[theme.default_color])
    default_aes.x = Float64[0.5]
    aes = inherit(aes, default_aes)

    n = length(aes.y)

    # Group y and width by x
    grouped_y     = DefaultDict(eltype(aes.x), typeof(aes.y), () -> similar(aes.y, 0))
    grouped_width = DefaultDict(eltype(aes.x), typeof(aes.width), () -> similar(aes.width, 0))
    for (x, y, w) in zip(cycle(aes.x), aes.y, aes.width)
        push!(grouped_y[x], y)
        push!(grouped_width[x], w)
    end

    ctx = context(order=geom.order)
    compose!(ctx,
             Compose.polygon([vcat([(x - w/2, y) for (y, w) in zip(grouped_y[x], grouped_width[x])],
                                   reverse!([(x + w/2, y) for (y, w) in zip(grouped_y[x], grouped_width[x])]))
                              for x in keys(grouped_y)]))

    return compose!(ctx,
                    fill(theme.default_color),
                    svgclass("geometry"))
end


text:
@testset "reddit" begin
    graph = graphdata(Reddit())
    @test typeof(graph) == SimpleGraph{UInt32}
    @test nv(graph) == 232965
    @test ne(graph) == 57307946

    all_X, all_y = alldata(Reddit())
    @test typeof(all_X) == Matrix{Float32}
    @test size(all_X) == (602, 232965)
    @test typeof(all_y) == Matrix{UInt16}
    @test size(all_y) == (1, 232965)

    raw = rawdata(Reddit())
    @test typeof(raw[:graph]) == SparseMatrixCSC{Int32,Int64}
    @test typeof(raw[:X]) == Matrix{Float32}
    @test typeof(raw[:y]) == Vector{Int32}
    @test typeof(raw[:ids]) == Vector{Int32}
    @test typeof(raw[:types]) == Vector{UInt8}

    meta = metadata(Reddit())
    @test meta.graph.num_V == nv(graph)
    @test meta.graph.num_E == ne(graph)
end


text:
const GDALColorTable      = Ptr{GDAL.GDALColorTableH}
const GDALCoordTransform  = Ptr{GDAL.OGRCoordinateTransformationH}
const GDALDataset         = Ptr{GDAL.GDALDatasetH}
const GDALDriver          = Ptr{GDAL.GDALDriverH}
const GDALFeature         = Ptr{GDAL.OGRFeatureH}
const GDALFeatureDefn     = Ptr{GDAL.OGRFeatureDefnH}
const GDALFeatureLayer    = Ptr{GDAL.OGRLayerH}
const GDALField           = Ptr{GDAL.OGRField}
const GDALFieldDefn       = Ptr{GDAL.OGRFieldDefnH}
const GDALGeometry        = Ptr{GDAL.OGRGeometryH}
const GDALGeomFieldDefn   = Ptr{GDAL.OGRGeomFieldDefnH}
const GDALProgressFunc    = Ptr{GDAL.GDALProgressFunc}
const GDALRasterAttrTable = Ptr{GDAL.GDALRasterAttributeTableH}
const GDALRasterBand      = Ptr{GDAL.GDALRasterBandH}
const GDALSpatialRef      = Ptr{GDAL.OGRSpatialReferenceH}
const GDALStyleManager    = Ptr{GDAL.OGRStyleMgrH}
const GDALStyleTable      = Ptr{GDAL.OGRStyleTableH}
const GDALStyleTool       = Ptr{GDAL.OGRStyleToolH}

const StringList          = Ptr{Ptr{UInt8}}

abstract type AbstractGeometry <: GeoInterface.AbstractGeometry end
    # needs to have a `ptr::GDALGeometry` attribute

abstract type AbstractSpatialRef end
    # needs to have a `ptr::GDALSpatialRef` attribute

mutable struct ColorTable;                    ptr::GDALColorTable         end
mutable struct CoordTransform;                ptr::GDALCoordTransform     end
mutable struct Dataset;                       ptr::GDALDataset            end
mutable struct Driver;                        ptr::GDALDriver             end
mutable struct Feature;                       ptr::GDALFeature            end
mutable struct FeatureDefn;                   ptr::GDALFeatureDefn        end
mutable struct FeatureLayer;                  ptr::GDALFeatureLayer       end
mutable struct Field;                         ptr::GDALField              end
mutable struct FieldDefn;                     ptr::GDALFieldDefn          end
mutable struct Geometry <: AbstractGeometry;  ptr::GDALGeometry           end
mutable struct IGeometry <: AbstractGeometry
    ptr::GDALGeometry

    function IGeometry(ptr::GDALGeometry)
        geom = new(GDAL.clone(ptr))
        finalizer(geom, destroy)
        geom
    end
end
mutable struct GeomFieldDefn;                    ptr::GDALGeomFieldDefn   end
mutable struct RasterAttrTable;                  ptr::GDALRasterAttrTable end
mutable struct RasterBand;                       ptr::GDALRasterBand      end
mutable struct SpatialRef <: AbstractSpatialRef; ptr::GDALSpatialRef      end
mutable struct ISpatialRef <: AbstractSpatialRef
    ptr::GDALSpatialRef

    function ISpatialRef(ptr::GDALSpatialRef)
        spref = new(GDAL.clone(ptr))
        finalizer(spref, destroy)
        spref
    end
end
mutable struct StyleManager;                     ptr::GDALStyleManager    end
mutable struct StyleTable;                       ptr::GDALStyleTable      end
mutable struct StyleTool;                        ptr::GDALStyleTool       end

CPLErr = GDAL.CPLErr
CPLXMLNodeType = GDAL.CPLXMLNodeType
GDALDataType = GDAL.GDALDataType
GDALAsyncStatusType = GDAL.GDALAsyncStatusType
GDALAccess = GDAL.GDALAccess
GDALRWFlag = GDAL.GDALRWFlag
GDALRIOResampleAlg = GDAL.GDALRIOResampleAlg
GDALColorInterp = GDAL.GDALColorInterp
GDALPaletteInterp = GDAL.GDALPaletteInterp
GDALRATFieldType = GDAL.GDALRATFieldType
GDALRATFieldUsage = GDAL.GDALRATFieldUsage
GDALTileOrganization = GDAL.GDALTileOrganization
GDALGridAlgorithm = GDAL.GDALGridAlgorithm
OGRwkbGeometryType = GDAL.OGRwkbGeometryType
OGRwkbVariant = GDAL.OGRwkbVariant
OGRwkbByteOrder = GDAL.OGRwkbByteOrder
OGRFieldType = GDAL.OGRFieldType
OGRFieldSubType = GDAL.OGRFieldSubType
OGRJustification = GDAL.OGRJustification
OGRSTClassId = GDAL.OGRSTClassId
OGRSTUnitId = GDAL.OGRSTUnitId
OGRSTPenParam = GDAL.OGRSTPenParam
OGRSTBrushParam = GDAL.OGRSTBrushParam
OGRSTSymbolParam = GDAL.OGRSTSymbolParam
OGRSTLabelParam = GDAL.OGRSTLabelParam
GDALResampleAlg = GDAL.GDALResampleAlg
GWKAverageOrModeAlg = GDAL.GWKAverageOrModeAlg
OGRAxisOrientation = GDAL.OGRAxisOrientation
OGRDatumType = GDAL.OGRDatumType

"return the corresponding `DataType` in julia"
const _JLTYPE = Dict{GDAL.GDALDataType, DataType}(
    GDAL.GDT_Unknown    => Any,
    GDAL.GDT_Byte       => UInt8,
    GDAL.GDT_UInt16     => UInt16,
    GDAL.GDT_Int16      => Int16,
    GDAL.GDT_UInt32     => UInt32,
    GDAL.GDT_Int32      => Int32,
    GDAL.GDT_Float32    => Float32,
    GDAL.GDT_Float64    => Float64)

const _GDALTYPE = Dict{DataType,GDAL.GDALDataType}(
    Any         => GDAL.GDT_Unknown,
    UInt8       => GDAL.GDT_Byte,
    UInt16      => GDAL.GDT_UInt16,
    Int16       => GDAL.GDT_Int16,
    UInt32      => GDAL.GDT_UInt32,
    Int32       => GDAL.GDT_Int32,
    Float32     => GDAL.GDT_Float32,
    Float64     => GDAL.GDT_Float64)

"return the corresponding `DataType` in julia"
const _FIELDTYPE = Dict{OGRFieldType, DataType}(
    GDAL.OFTInteger         => Int32,
    GDAL.OFTIntegerList     => Void,
    GDAL.OFTReal            => Float64,
    GDAL.OFTRealList        => Void,
    GDAL.OFTString          => String,
    GDAL.OFTStringList      => Void,
    GDAL.OFTWideString      => Void, # deprecated
    GDAL.OFTWideStringList  => Void, # deprecated
    GDAL.OFTBinary          => Void,
    GDAL.OFTDate            => Date,
    GDAL.OFTTime            => Void,
    GDAL.OFTDateTime        => DateTime,
    GDAL.OFTInteger64       => Int64,
    GDAL.OFTInteger64List   => Void)

@enum(GDALOpenFlag,
      OF_ReadOnly = GDAL.GDAL_OF_READONLY, # 0x00
      OF_Update   = GDAL.GDAL_OF_UPDATE,   # 0x01
#     OF_All      = GDAL.GDAL_OF_ALL,      # 0x00
      OF_Raster   = GDAL.GDAL_OF_RASTER,   # 0x02
      OF_Vector   = GDAL.GDAL_OF_VECTOR,   # 0x04
      OF_GNM      = GDAL.GDAL_OF_GNM)      # 0x08
                 # const GDAL_OF_KIND_MASK = 0x1e
                 # const GDAL_OF_SHARED = 0x20
                 # const GDAL_OF_VERBOSE_ERROR = 0x40
                 # const GDAL_OF_INTERNAL = 0x80
                 # const GDAL_OF_DEFAULT_BLOCK_ACCESS = 0
                 # const GDAL_OF_ARRAY_BLOCK_ACCESS = 0x0100
                 # const GDAL_OF_HASHSET_BLOCK_ACCESS = 0x0200
                 # const GDAL_OF_RESERVED_1 = 0x0300
                 # const GDAL_OF_BLOCK_ACCESS_MASK = 0x0300

import Base.|

|(x::GDALOpenFlag,y::UInt8) = UInt8(x) | y
|(x::UInt8,y::GDALOpenFlag) = x | UInt8(y)
|(x::GDALOpenFlag,y::GDALOpenFlag) = UInt8(x) | UInt8(y)

"Get data type size in bits."
typesize(dt::GDALDataType) =
    @gdal(GDALGetDataTypeSize::Cint,
        dt::GDAL.GDALDataType
    )

"name (string) corresponding to GDAL data type"
typename(dt::GDALDataType) =
    unsafe_string(@gdal(GDALGetDataTypeName::Cstring,
        dt::GDAL.GDALDataType
    ))

"Returns GDAL data type by symbolic name."
gettype(name::AbstractString) = GDALDataType(GDAL.getdatatypebyname(name))

"Return the smallest data type that can fully express both input data types."
typeunion(dt1::GDALDataType,dt2::GDALDataType) =
    GDALDataType(@gdal(GDALDataTypeUnion::GDAL.GDALDataType,
        dt1::GDAL.GDALDataType,
        dt2::GDAL.GDALDataType
    ))

"""
`true` if `dtype` is one of `GDT_{CInt16|CInt32|CFloat32|CFloat64}`
"""
iscomplex(dtype::GDALDataType) =
    Bool(@gdal(GDALDataTypeIsComplex::Cint,
        dtype::GDAL.GDALDataType
    ))

"Get name of AsyncStatus data type."
getname(dtype::GDALAsyncStatusType) =
    unsafe_string(@gdal(GDALGetAsyncStatusTypeName::Cstring,
        dtype::GDAL.GDALAsyncStatusType
    ))

"Get AsyncStatusType by symbolic name."
asyncstatustype(name::AbstractString) =
    GDALAsyncStatusType(GDAL.getasyncstatustypebyname(name))

"Return name (string) corresponding to color interpretation"
getname(obj::GDALColorInterp) =
    unsafe_string(@gdal(GDALGetColorInterpretationName::Cstring,
        obj::GDAL.GDALColorInterp
    ))

"Get color interpretation corresponding to the given symbolic name."
colorinterp(name::AbstractString) =
    GDALColorInterp(GDAL.getcolorinterpretationbyname(name))

"Get name of palette interpretation."
getname(obj::GDALPaletteInterp) =
    unsafe_string(@gdal(GDALGetPaletteInterpretationName::Cstring,
        obj::GDAL.GDALPaletteInterp
    ))

"Fetch human readable name for a field type."
getname(obj::OGRFieldType) =
    unsafe_string(@gdal(OGR_GetFieldTypeName::Cstring,
        obj::GDAL.OGRFieldType
    ))

"Fetch human readable name for a field subtype."
getname(obj::OGRFieldSubType) =
    unsafe_string(@gdal(OGR_GetFieldSubTypeName::Cstring,
        obj::GDAL.OGRFieldSubType
    ))

"Return if type and subtype are compatible."
arecompatible(dtype::OGRFieldType, subtype::OGRFieldSubType) =
    Bool(@gdal(OGR_AreTypeSubTypeCompatible::Cint,
        dtype::GDAL.OGRFieldType,
        subtype::GDAL.OGRFieldSubType
    ))


text:
###
# Format and filters

# Reader

support_filter_all(archive::Reader) =
    @_la_call(archive_read_support_filter_all, (Ptr{Void},), archive)
support_filter_bzip2(archive::Reader) =
    @_la_call(archive_read_support_filter_bzip2, (Ptr{Void},), archive)
support_filter_compress(archive::Reader) =
    @_la_call(archive_read_support_filter_compress, (Ptr{Void},), archive)
support_filter_gzip(archive::Reader) =
    @_la_call(archive_read_support_filter_gzip, (Ptr{Void},), archive)
support_filter_grzip(archive::Reader) =
    @_la_call(archive_read_support_filter_grzip, (Ptr{Void},), archive)
support_filter_lrzip(archive::Reader) =
    @_la_call(archive_read_support_filter_lrzip, (Ptr{Void},), archive)
support_filter_lzip(archive::Reader) =
    @_la_call(archive_read_support_filter_lzip, (Ptr{Void},), archive)
support_filter_lzma(archive::Reader) =
    @_la_call(archive_read_support_filter_lzma, (Ptr{Void},), archive)
support_filter_lzop(archive::Reader) =
    @_la_call(archive_read_support_filter_lzop, (Ptr{Void},), archive)
support_filter_none(archive::Reader) =
    @_la_call(archive_read_support_filter_none, (Ptr{Void},), archive)
"""
Data is fed through the specified external program before being dearchived.
Note that this disables automatic detection of the compression format, so it
makes no sense to specify this in conjunction with any other decompression
option.
"""
support_filter_program(archive::Reader, cmd) =
    @_la_call(archive_read_support_filter_program,
              (Ptr{Void}, Cstring), archive, cmd)
support_filter_program_signature(archive::Reader, cmd, sig) =
    @_la_call(archive_read_support_filter_program_signature,
              (Ptr{Void}, Cstring, Ptr{Void}, Csize_t),
              archive, cmd, sig, sizeof(sig))
support_filter_rpm(archive::Reader) =
    @_la_call(archive_read_support_filter_rpm, (Ptr{Void},), archive)
support_filter_uu(archive::Reader) =
    @_la_call(archive_read_support_filter_uu, (Ptr{Void},), archive)
support_filter_xz(archive::Reader) =
    @_la_call(archive_read_support_filter_xz, (Ptr{Void},), archive)

support_format_7zip(archive::Reader) =
    @_la_call(archive_read_support_format_7zip, (Ptr{Void},), archive)
support_format_all(archive::Reader) =
    @_la_call(archive_read_support_format_all, (Ptr{Void},), archive)
support_format_ar(archive::Reader) =
    @_la_call(archive_read_support_format_ar, (Ptr{Void},), archive)
support_format_by_code(archive::Reader, code) =
    @_la_call(archive_read_support_format_by_code,
              (Ptr{Void}, Cint), archive, code)
support_format_cab(archive::Reader) =
    @_la_call(archive_read_support_format_cab, (Ptr{Void},), archive)
support_format_cpio(archive::Reader) =
    @_la_call(archive_read_support_format_cpio, (Ptr{Void},), archive)
support_format_empty(archive::Reader) =
    @_la_call(archive_read_support_format_empty, (Ptr{Void},), archive)
support_format_gnutar(archive::Reader) =
    @_la_call(archive_read_support_format_gnutar, (Ptr{Void},), archive)
support_format_iso9660(archive::Reader) =
    @_la_call(archive_read_support_format_iso9660, (Ptr{Void},), archive)
support_format_lha(archive::Reader) =
    @_la_call(archive_read_support_format_lha, (Ptr{Void},), archive)
support_format_mtree(archive::Reader) =
    @_la_call(archive_read_support_format_mtree, (Ptr{Void},), archive)
support_format_rar(archive::Reader) =
    @_la_call(archive_read_support_format_rar, (Ptr{Void},), archive)
support_format_raw(archive::Reader) =
    @_la_call(archive_read_support_format_raw, (Ptr{Void},), archive)
support_format_tar(archive::Reader) =
    @_la_call(archive_read_support_format_tar, (Ptr{Void},), archive)
support_format_xar(archive::Reader) =
    @_la_call(archive_read_support_format_xar, (Ptr{Void},), archive)
support_format_zip(archive::Reader) =
    @_la_call(archive_read_support_format_zip, (Ptr{Void},), archive)

# Functions to manually set the format and filters to be used. This is
# useful to bypass the bidding process when the format and filters to use
# is known in advance.
set_format(archive::Reader, fmt) =
    @_la_call(archive_read_set_format, (Ptr{Void}, Cint), archive, fmt)
append_filter(archive::Reader, filter) =
    @_la_call(archive_read_append_filter, (Ptr{Void}, Cint), archive, filter)
append_filter_program(archive::Reader, cmd) =
    @_la_call(archive_read_append_filter_program, (Ptr{Void}, Cstring),
              archive, cmd)
append_filter_program_signature(archive::Reader, cmd, sig) =
    @_la_call(archive_read_append_filter_program_signature,
              (Ptr{Void}, Cstring, Ptr{Void}, Csize_t),
              archive, cmd, sig, sizeof(sig))

# Writer

"A convenience function to set the filter based on the code."
add_filter(archive::Writer, filter_code::Integer) =
    @_la_call(archive_write_add_filter, (Ptr{Void}, Cint),
              archive, filter_code)
add_filter(archive::Writer, name::AbstractString) =
    @_la_call(archive_write_add_filter_by_name, (Ptr{Void}, Cstring),
              archive, name)
add_filter_b64encode(archive::Writer) =
    @_la_call(archive_write_add_filter_b64encode, (Ptr{Void},), archive)
add_filter_bzip2(archive::Writer) =
    @_la_call(archive_write_add_filter_bzip2, (Ptr{Void},), archive)
add_filter_compress(archive::Writer) =
    @_la_call(archive_write_add_filter_compress, (Ptr{Void},), archive)
add_filter_grzip(archive::Writer) =
    @_la_call(archive_write_add_filter_grzip, (Ptr{Void},), archive)
add_filter_gzip(archive::Writer) =
    @_la_call(archive_write_add_filter_gzip, (Ptr{Void},), archive)
add_filter_lrzip(archive::Writer) =
    @_la_call(archive_write_add_filter_lrzip, (Ptr{Void},), archive)
add_filter_lzip(archive::Writer) =
    @_la_call(archive_write_add_filter_lzip, (Ptr{Void},), archive)
add_filter_lzma(archive::Writer) =
    @_la_call(archive_write_add_filter_lzma, (Ptr{Void},), archive)
add_filter_lzop(archive::Writer) =
    @_la_call(archive_write_add_filter_lzop, (Ptr{Void},), archive)
add_filter_none(archive::Writer) =
    @_la_call(archive_write_add_filter_none, (Ptr{Void},), archive)
add_filter_program(archive::Writer, cmd::AbstractString) =
    @_la_call(archive_write_add_filter_program,
              (Ptr{Void}, Cstring), archive, cmd)
add_filter_uuencode(archive::Writer) =
    @_la_call(archive_write_add_filter_uuencode, (Ptr{Void},), archive)
add_filter_xz(archive::Writer) =
    @_la_call(archive_write_add_filter_xz, (Ptr{Void},), archive)

"A convenience function to set the format based on the code or name."
set_format(archive::Writer, format_code::Integer) =
    @_la_call(archive_write_set_format, (Ptr{Void}, Cint),
              archive, format_code)
set_format(archive::Writer, name::AbstractString) =
    @_la_call(archive_write_set_format_by_name, (Ptr{Void}, Cstring),
              archive, name)
set_format_7zip(archive::Writer) =
    @_la_call(archive_write_set_format_7zip, (Ptr{Void},), archive)
set_format_ar_bsd(archive::Writer) =
    @_la_call(archive_write_set_format_ar_bsd, (Ptr{Void},), archive)
set_format_ar_svr4(archive::Writer) =
    @_la_call(archive_write_set_format_ar_svr4, (Ptr{Void},), archive)
set_format_cpio(archive::Writer) =
    @_la_call(archive_write_set_format_cpio, (Ptr{Void},), archive)
set_format_cpio_newc(archive::Writer) =
    @_la_call(archive_write_set_format_cpio_newc, (Ptr{Void},), archive)
set_format_gnutar(archive::Writer) =
    @_la_call(archive_write_set_format_gnutar, (Ptr{Void},), archive)
set_format_iso9660(archive::Writer) =
    @_la_call(archive_write_set_format_iso9660, (Ptr{Void},), archive)
set_format_mtree(archive::Writer) =
    @_la_call(archive_write_set_format_mtree, (Ptr{Void},), archive)
set_format_mtree_classic(archive::Writer) =
    @_la_call(archive_write_set_format_mtree_classic, (Ptr{Void},), archive)
set_format_pax(archive::Writer) =
    @_la_call(archive_write_set_format_pax, (Ptr{Void},), archive)
set_format_pax_restricted(archive::Writer) =
    @_la_call(archive_write_set_format_pax_restricted, (Ptr{Void},), archive)
set_format_shar(archive::Writer) =
    @_la_call(archive_write_set_format_shar, (Ptr{Void},), archive)
set_format_shar_dump(archive::Writer) =
    @_la_call(archive_write_set_format_shar_dump, (Ptr{Void},), archive)
set_format_ustar(archive::Writer) =
    @_la_call(archive_write_set_format_ustar, (Ptr{Void},), archive)
set_format_v7tar(archive::Writer) =
    @_la_call(archive_write_set_format_v7tar, (Ptr{Void},), archive)
set_format_xar(archive::Writer) =
    @_la_call(archive_write_set_format_xar, (Ptr{Void},), archive)
set_format_zip(archive::Writer) =
    @_la_call(archive_write_set_format_zip, (Ptr{Void},), archive)
zip_set_compression_deflate(archive::Writer) =
    @_la_call(archive_write_zip_set_compression_deflate, (Ptr{Void},), archive)
zip_set_compression_store(archive::Writer) =
    @_la_call(archive_write_zip_set_compression_store, (Ptr{Void},), archive)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -9.734460520574642 <= q <= 2.116911729869532)

                     add_NL_constraint(m, :(gelu(-0.8046949419530325 + -0.2278494590955038*gelu(0.1770285029406451 + -0.5067338360802731*gelu(-0.7998422425711689 + 0.5883237869246916*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*gelu(-0.7622628423532047 + 0.6518912481623378*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.18111143847071443*gelu(-0.12446319977672715 + 0.7018706953017659*gelu(-0.7998422425711689 + 0.5883237869246916*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*gelu(-0.7622628423532047 + 0.6518912481623378*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + gelu(-0.4025778812055263 + 0.08853333939955377*gelu(0.1770285029406451 + -0.5067338360802731*gelu(-0.7998422425711689 + 0.5883237869246916*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*gelu(-0.7622628423532047 + 0.6518912481623378*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.4269192713141794*gelu(-0.12446319977672715 + 0.7018706953017659*gelu(-0.7998422425711689 + 0.5883237869246916*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*gelu(-0.7622628423532047 + 0.6518912481623378*gelu(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*gelu(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("practic_8.jl")

"""

"""
function  shuttle!(condition::Function, robot, side)
    n=0 # число шагов от начального положения
    while condition() == false # condition(robot) - как вариант
        n += 1
        movements!(() -> !condition(), robot, side, n)
        side = inverse(side)
    end
end

"""
movements!(condition::Function, robot::AbstractBorderRobot2, side, num_steps)

-- делает не более чем num_steps шагов в заданном направлении, пока выпоняется условие condition()==true",
и возвращает число сделанных шагов

-- condition - функция без аргументов, возвращающая логическое значение
"""
function movements!(condition::Function, robot::AbstractBorderRobot2, side, max_num_steps)
    n = 0
    while n < max_num_steps && condition()
        try_move!(robot, side) # try_move! - позволяет обходить прямоугольные перегородки (если они есть), как бы по прямой
        n += 1
    end
    return n
end

struct BorderRobot <: AbstractBorderRobot2{Robot}
    robot::Robot
end

get(r::BorderRobot) = r.robot
action!(r::BorderRobot) = nothing

#= Примеры вызова вункции shuttle:

julia> include("practic_9.jl")
movements!

julia> shuttle!(() -> ismarker(r), r, Ost) # - робот ходит влево-вправо пока не наткнется на маркер

julia> shuttle!(() -> isborder(r,Nord), r, Ost) # - робот ходит влево-вправо пока не наткнется на перегородку сверху

julia> shuttle!(() -> isborder(r,Ost) || isborder(r,West), r, Ost) # - робот ходит влево-вправо пока не наткнется на перегородку по ходу движения

julia>

julia> border_robot = BorderRobot(r);

julia> shuttle!(() -> ismarker(border_robot), border_robot, Ost) # - робот ходит влево-вправо, обходя внутренние перегородки, пока не наткнется на маркер

julia>
=#

function spiral!(condition::Function, robot, side=Nord)
    n=1
    while condition() == false
        movements!(() -> !condition(), robot, side, n)
        side = left(side)
        movements!(() -> !condition(), robot, side, n)
        side = left(side)
        n += 1
    end
end


text:
using DSGE, ModelConstructors
using FileIO, HDF5, Test
using DataFrames

path = dirname(@__FILE__)
writing_output = false

# Set up model for testing
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Read in the data, mode, and hessian
mode = load("$path/../reference/hessian.jld2", "paramsmode")
data = Matrix{Float64}(load("$path/../reference/hessian.jld2", "data")')

# Read in the covariance matrix for Metropolis-Hastings and reference parameter draws
# TODO: check that new MH agrees with old MH. (read in metropolis_hastings.h5)
hessian_inv, ref_draws, ref_cov =
    h5open("$path/../reference/metropolis_hastings_test.h5", "r") do file
        read(file, "hessian_inv"),
        read(file, "mhparams"),
        read(file, "ref_cov")
    end

# Set up and run metropolis-hastings
DSGE.update!(m, mode)
prop_cov = DegenerateMvNormal(mode, hessian_inv)

metropolis_hastings(prop_cov, m, data, .01, .09; verbose=:none)
compute_parameter_covariance(m)

# Read in the parameter draws and covariance just generated from estimate.
test_draws = h5open(rawpath(m, "estimate", "mhsave.h5"), "r") do file
    read(file, "mhparams")
end
test_cov = h5open(workpath(m, "estimate", "parameter_covariance.h5"), "r") do file
    read(file, "mhcov")
end

if writing_output
    h5open("$path/../reference/metropolis_hastings_test.h5", "w") do file
        write(file, "hessian_inv", hessian_inv),
        write(file, "mhparams",    test_draws),
        write(file, "ref_cov",     test_cov)
    end
end

# Test that the parameter draws and covariance matrices are equal
@testset "Check equality of parameter draws and cov matrices in MH (1 block)" begin
    @test @test_matrix_approx_eq ref_draws test_draws
    @test @test_matrix_approx_eq ref_cov test_cov
end

# Set up and run metropolis-hastings with three blocks!
hessian_inv, ref_draws, ref_cov =
    h5open("$path/../reference/metropolis_hastings_test_3_blocks.h5", "r") do file
        read(file, "hessian_inv"),
        read(file, "mhparams"),
        read(file, "ref_cov")
    end

DSGE.update!(m, mode)
prop_cov = DegenerateMvNormal(mode, hessian_inv)
m <= Setting(:n_mh_blocks, 3)

metropolis_hastings(prop_cov, m, data, .01, .09; verbose=:none)
compute_parameter_covariance(m)

# Read in the parameter draws and covariance just generated from estimate.
test_draws = h5open(rawpath(m, "estimate", "mhsave.h5"), "r") do file
    read(file, "mhparams")
end
test_cov = h5open(workpath(m, "estimate", "parameter_covariance.h5"), "r") do file
    read(file, "mhcov")
end

if writing_output
    h5open("$path/../reference/metropolis_hastings_test_3_blocks.h5", "w") do file
        write(file, "hessian_inv", hessian_inv),
        write(file, "mhparams",    test_draws),
        write(file, "ref_cov",     test_cov)
    end
end

# Test that the parameter draws and covariance matrices are equal
@testset "Check equality of parameter draws and cov matrices in MH (3 blocks)" begin
    @test @test_matrix_approx_eq ref_draws test_draws
    @test @test_matrix_approx_eq ref_cov test_cov
end

nothing


text:
@testset "multilevel_PARSDMM" begin

# if nworkers()==1
#   addprocs(3)
# end

#test serial multilevel PARSDMM
    options=PARSDMM_options()
    default_PARSDMM_options(options,Float64)
    options.parallel = false
    options.evol_rel_tol = 10*eps()
    options.maxit=5000

    n_levels=2
    coarsening_factor=1.87

    x=randn(100,201)

    comp_grid=compgrid((1.0, 1.0),(100, 201))
    x         = vec(x)
    m         = deepcopy(x)

    constraint=Dict()

    #total variation
    (TV_OP, AtA_diag, dense, TD_n, banded) = get_TD_operator(comp_grid,"TV",options.FL)
    constraint["use_TD_l1_1"]      = true
    constraint["TD_l1_operator_1"] = "TV"
    constraint["TD_l1_sigma_1"]    = 0.4*norm(TV_OP*x,1)

    (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,constraint_level)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

    (x1,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
    result=Vector{typeof(x[1])}(length(x1))
    for i=1:length(TD_OP_levels[1])-1
      copy!(result,x1)
      @test norm(P_sub_levels[1][i](TD_OP_levels[1][i]*result)-(TD_OP_levels[1][i]*result))/norm((TD_OP_levels[1][i]*result)) <= options.feas_tol
    end

#test parallel multilevel PARSDMM
        options=PARSDMM_options()
        default_PARSDMM_options(options,Float64)
        options.parallel = true

            (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,constraint_level)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

            (x2,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
            result=Vector{typeof(x[1])}(length(x2))
            for i=1:length(TD_OP_levels[1])-1
              copy!(result,x2)
              @test norm(P_sub_levels[1][i](TD_OP_levels[1][i]*result)-(TD_OP_levels[1][i]*result))/norm((TD_OP_levels[1][i]*result)) <= options.feas_tol
            end
@test isapprox(x1,x2,rtol=options.feas_tol)

#run serial PARSDMM and compare results
default_PARSDMM_options(options,Float64)
options.parallel    = false

(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL);
(TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

(x3,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
result=Vector{typeof(x[1])}(length(x3))
for i=1:length(TD_OP)-1
  copy!(result,x3)
  @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= options.feas_tol
end

@test isapprox(x1,x3,rtol=options.feas_tol)

end


text:
"""
    get_classes(mtg)

Compute the mtg classes based on its content. Usefull after having mutating the mtg nodes.
"""
function get_classes(mtg)
    attributes = traverse(mtg, node -> (SYMBOL = node.MTG.symbol, SCALE = node.MTG.scale))
    attributes = unique(attributes)
    df = DataFrame(attributes)

    # Make everything to default values:
    df[!, :DECOMPOSITION] .= "FREE"
    df[!, :INDEXATION] .= "FREE"
    df[!, :DEFINITION] .= "IMPLICIT"
    return df
end

"""
    get_description(mtg)

Returns `nothing`, because we can't really predict the description section from an mtg.
"""
function get_description(mtg)
    return nothing
end

"""
    get_features(mtg)

Compute the mtg features section based on its attributes. Usefull after having computed new attributes
in the mtg.
"""
function get_features(mtg)
    attributes = traverse(mtg, node -> (collect(keys(node.attributes)), [typeof(i) for i in values(node.attributes)]))

    df = DataFrame(
        :NAME => vcat([i[1] for i in attributes]...),
        :TYPE => vcat([i[2] for i in attributes]...)
    )

    # filter-out the attributes that have more than one value inside them (not compatible with
    # the mtg format yet):
    filter!(x -> x.TYPE <: Number || x.TYPE <: AbstractString, df)

    # Remove repeated rows:
    unique!(df)

    new_type = fill("", size(df)[1])
    for (index, value) in enumerate(df.TYPE)
        if value <: AbstractFloat
            new_type[index] = "REAL"
        elseif value <: Int
            new_type[index] = "INT"
            # elseif df.NAME[i] in () # Put reserved keywords here if needed in the future
            # new_type[index] = "ALPHA"
        else
            # All others are parsed as string
            new_type[index] = "STRING"
        end
    end

    df.TYPE = new_type

    return df
end

"""
    scales(mtg)

Get all the scales of an MTG.
"""
function scales(mtg)
    unique(traverse(mtg, node -> node.MTG.scale))
end

function symbols(mtg)
    unique(traverse(mtg, node -> node.MTG.symbol))
end

components = symbols

"""
    symbols(mtg)
    components(mtg)

Get all the symbols names, a.k.a. components of an MTG.
"""
components, symbols

"""
    get_attributes(mtg)

Get all attributes names available on the mtg and its children.
"""
get_attributes(mtg) = unique(vcat(traverse(mtg, node -> collect(keys(node.attributes)))...))

"""
    names(mtg)

Get all attributes names available on the mtg and its children. This is an alias for
[`get_attributes`](@ref).
"""
Base.names(mtg::T) where {T<:MultiScaleTreeGraph.Node} = get_attributes(mtg)

"""
    list_nodes(mtg)

List all nodes IDs in the subtree of `mtg`.
"""
list_nodes(mtg) = traverse(mtg, node -> node.id)

"""
    max_id(mtg)

Returns the maximum id of the mtg
"""
function max_id(mtg)
    maxid = [0]

    function update_maxname(id, maxid)
        id > maxid[1] ? maxid[1] = id : nothing
    end

    traverse!(get_root(mtg), x -> update_maxname(x.id, maxid))

    return maxid[1]
end


text:
using VibronicToolkit: @cat

let A = collect(reshape(1:9, 3, 3))

    @test diag(A, 1, 2) == diag(A)
end

let A = collect(reshape(1:8, 2, 2, 2))

    @test diag(A, 1, 2) == [1 5; 4 8]
    @test diag(A, 1, 3) == [1 3; 6 8]
    @test diag(A, 2, 3) == [1 7; 2 8]
end

let A = collect(reshape(1:72, 2, 3, 2, 3, 2)),
    B = diag(A, 2, 4)

    for i in 1:2, j in 1:3, k in 1:2, l in 1:2
        @test B[i, j, k, l] == A[i, j, k, j, l]
    end
    C = diag(diag(A, 1, 3), 1, 4)
    for i in 1:2, j in 1:3, k in 1:3
        @test C[i, j, k] == A[i, j, i, k, i]
    end
    @test_throws DomainError diag(A, 1, 2)
end

let A = collect(reshape(1:72, 2, 3, 2, 3, 2)),
    B = collect(reshape(73:180, 2, 3, 2, 3, 3)),
    C = collect(reshape(181:216, 2, 3, 2, 3))

    @test A == reshape(1:72, 2, 3, 2, 3, 2)
    @cat(A, B)
    @test A == reshape(1:180, 2, 3, 2, 3, 5)
    @cat(A, C)
    @test A == reshape(1:216, 2, 3, 2, 3, 6)
end


text:
"""
    RemainderModel(main::BalanceLaw, subcomponents::Tuple)

Compute the "remainder" contribution of the `main` model, after subtracting `subcomponents`.
"""
struct RemainderModel{M,S} <: BalanceLaw
  main::M
  subs::S
end

vars_state(rem::RemainderModel, FT) = vars_state(rem.main,FT)
vars_gradient(rem::RemainderModel, FT) = vars_gradient(rem.main,FT)
vars_diffusive(rem::RemainderModel, FT) = vars_diffusive(rem.main,FT)
vars_aux(rem::RemainderModel, FT) = vars_aux(rem.main,FT)
vars_integrals(rem::RemainderModel,FT) = vars_integrals(rem.main,FT)

update_aux!(dg::DGModel, rem::RemainderModel, Q::MPIStateArray, t::Real) =
  update_aux!(dg, rem.main, Q, t)

integrate_aux!(rem::RemainderModel, integ::Vars, state::Vars, aux::Vars) =
  integrate_aux!(rem.main, integ, state, aux)

flux_diffusive!(rem::RemainderModel, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real) =
  flux_diffusive!(rem.main, flux, state, diffusive, aux, t)

gradvariables!(rem::RemainderModel, transform::Vars, state::Vars, aux::Vars, t::Real) =
  gradvariables!(rem.main, transform, state, aux, t)

diffusive!(rem::RemainderModel, diffusive::Vars, ∇transform::Grad, state::Vars, aux::Vars, t::Real) =
  diffusive!(rem.main, diffusive, ∇transform, state, aux, t)

function wavespeed(rem::RemainderModel, nM, state::Vars, aux::Vars, t::Real)
  ref = aux.ref_state
  return wavespeed(rem.main, nM, state, aux, t) - sum(sub -> wavespeed(sub, nM, state, aux, t), rem.subs)
end

boundary_state!(nf, rem::RemainderModel, x...) = boundary_state!(nf, rem.main, x...)

init_aux!(rem::RemainderModel, aux::Vars, geom::LocalGeometry) = nothing
init_state!(rem::RemainderModel, state::Vars, aux::Vars, coords, t) = nothing

function flux_nondiffusive!(rem::RemainderModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  m = getfield(flux, :array)
  flux_nondiffusive!(rem.main, flux, state, aux, t)

  flux_s = similar(flux)
  m_s = getfield(flux_s, :array)

  for sub in rem.subs
    fill!(m_s, 0)
    flux_nondiffusive!(sub, flux_s, state, aux, t)
    m .-= m_s
  end
  nothing
end

function source!(rem::RemainderModel, source::Vars, state::Vars, aux::Vars, t::Real)
  m = getfield(source, :array)
  source!(rem.main, source, state, aux, t)

  source_s = similar(source)
  m_s = getfield(source_s, :array)

  for sub in rem.subs
    fill!(m_s, 0)
    source!(sub, source_s, state, aux, t)
    m .-= m_s
  end
  nothing
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
using Replay

instructions = """
println("julia -q")
"""

replay(instructions, cmd="-q")


text:
function nthroot(n::Integer, r::Real)
    r < 0 || n == 0 && throw(DomainError())
    n < 0 && return 1 / nthroot(-n, r)
    r > 0 || return 0
    x = r / n
    prevdx = r
    while true
        y = x ^ (n - 1)
        dx = (r - y * x) / (n * y)
        abs(dx) ≥ abs(prevdx) && return x
        x += dx
        prevdx = dx
    end
end

@show nthroot.(-5:2:5, 5.0)
@show nthroot.(-5:2:5, 5.0) - 5.0 .^ (1 ./ (-5:2:5))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"A `Token` can only be held by one async process at one time. Use `Base.take!(token)` to claim the token, `Base.put!(token)` to give the token back."
struct Token
    c::Channel{Nothing}
    Token() = let
        c = Channel{Nothing}(1)
        push!(c, nothing)
        new(c)
    end
end

Base.take!(token::Token) = Base.take!(token.c)
Base.put!(token::Token) = Base.put!(token.c, nothing)
Base.isready(token::Token) = Base.isready(token.c)
Base.wait(token::Token) = Base.put!(token.c, Base.take!(token.c))

"Track whether some task needs to be done. `request!(requestqueue)` will make sure that the task is done at least once after calling it. Multiple calls might get bundled into one."
mutable struct RequestQueue
    is_processing::Bool
    c::Channel{Nothing}
    RequestQueue() = new(false, Channel{Nothing}(1))
end

"Give a function (with no arguments) that should be called after a request."
function process(f::Function, requestqueue::RequestQueue)
    @assert !requestqueue.is_processing
    requestqueue.is_processing = true
    while true
        take!(requestqueue.c)
        f()
    end
end

function request!(queue::RequestQueue)
    if isready(queue.c)
        push!(queue.c, nothing)
    end
end


text:
#### Variable Mapper

"""
    get_var_mapper(vars)

Get a `NamedTuple` that maps a variable name
and subdomain ID to a unique index. Example:
```
julia> vars = ((:ρ_0, 1), (:w, 3), (:a, 3), (:α_0, 1));

julia> var_mapper = get_var_mapper(vars)
(ρ_0 = (1,), w = (2, 3, 4), a = (5, 6, 7), α_0 = (8,))
```
"""
function get_var_mapper(vars, domain_set::DomainSet)
    var_names = tuple([v for (v, dss) in vars]...)
    n_sd_per_var = [sum(domain_set, dss) for (v, dss) in vars]
    end_index, n_sd_max = cumsum(n_sd_per_var), max(n_sd_per_var...)
    start_index = [1, [1 + x for x in end_index][1:(end - 1)]...]
    vals = Tuple([Tuple(a:b) for (a, b) in zip(start_index, end_index)])
    var_mapper = NamedTuple{var_names}(vals)
    D = Dict([v => dss for (v, dss) in vars]...)
    dss_per_var = NamedTuple{Tuple(keys(D))}(Tuple(values(D)))
    return var_mapper, dss_per_var, var_names
end

function get_sd_mapped(vars, idx::DomainIdx, domain_set::DomainSet)
    sd_mapped = Dict()
    al = alldomains(idx)
    for (v, dss) in vars
        idx_ss = DomainIdx(domain_set, dss)
        if !(haskey(sd_mapped, v))
            sd_mapped[v] = Int[]
        end
        for i in al
            i == gridmean(idx) && push!(sd_mapped[v], gridmean(idx_ss))
            i == environment(idx) && push!(sd_mapped[v], environment(idx_ss))
            if i in updraft(idx) && !(i in sd_mapped[v])
                push!(sd_mapped[v], updraft(idx_ss)...)
            end
        end
    end
    return sd_mapped
end

function get_sv_a_map(idx::DomainIdx, idx_ss::DomainIdx)
    a = zeros(Int, length(alldomains(idx)))
    for i in alldomains(idx)
        if i ≠ 0 # protect against null map
            if i == gridmean(idx) && has_gridmean(idx_ss)
                a[i] = gridmean(idx_ss)
            elseif i == environment(idx) && has_environment(idx_ss)
                a[i] = environment(idx_ss)
            elseif i in updraft(idx) && has_updraft(idx_ss)
                a[i] = updraft(idx_ss)[i]
            else
                a[i] = 0
            end
        end
    end
    return a
end

function get_sd_unmapped(vars, idx::DomainIdx, domain_set::DomainSet)
    sd_unmapped = Dict{Symbol, Vector{Int}}()
    al = alldomains(idx)
    for (v, dss) in vars
        idx_ss = DomainIdx(domain_set, dss)
        if !(haskey(sd_unmapped, v))
            sd_unmapped[v] = Int[]
        end
        for i in al
            i == gridmean(idx_ss) && push!(sd_unmapped[v], gridmean(idx))
            i == environment(idx_ss) && push!(sd_unmapped[v], environment(idx))
            if i in updraft(idx_ss) && !(i in sd_unmapped[v])
                push!(sd_unmapped[v], updraft(idx)...)
            end
        end
    end
    return sd_unmapped
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module REAliases

using ..RE
using ..Empty
using ..SK
using ..EA
using ..SAT
using ..PercLinear
using ..PercStep
using ..CommStep
using ..CommReLU
using ..CommQu

export Graph0RE, GraphSKRE, GraphEARE, GraphSATRE,
       GraphPercLinearRE, GraphPercStepRE, GraphCommStepRE,
       GraphCommReLURE, GraphCommQuRE

const Graph0RE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphEmpty}

# """
#     Graph0RE(...)
#
# TODO
# """
Graph0RE(Nk::Integer, M::Integer, γ::Float64, β::Float64) = GraphRobustEnsemble(Nk, M, γ, β, GraphEmpty, Nk)

const GraphSKRE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphSK}

# """
#     GraphSKRE(N::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
GraphSKRE(Nk::Integer, M::Integer, γ::Float64, β::Float64) = GraphRobustEnsemble(Nk, M, γ, β, GraphSK, SK.gen_J(Nk))

const GraphEARE{M,γ,β,twoD} = GraphRobustEnsemble{M,γ,β,GraphEANormal{twoD}}

# """
#     GraphEARE(L::Integer, D::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphEARE(L::Integer, D::Integer, M::Integer, γ::Float64, β::Float64)
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    A = EA.gen_EA(L, D)
    N = length(A)
    J = EA.gen_J(Float64, N, A) do
        4*rand() - 2
    end

    GraphRobustEnsemble(N, M, γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEARE(fname::AbstractString, M::Integer, γ::Float64, β::Float64)
    L, D, A, J = EA.gen_AJ(fname)
    N = length(A)
    @assert N == L^D
    GraphRobustEnsemble(N, M, γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEARE(X::GraphEANormal{twoD}, M::Integer, γ::Float64, β::Float64) where {twoD}
    @assert iseven(twoD)
    D = twoD ÷ 2
    N = X.N
    L = round(Int, N^(1/D))
    @assert L^D == N
    GraphRobustEnsemble(N, M, γ, β, GraphEANormal{twoD}, L, X.A, X.J)
end

const GraphSATRE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphSAT}

# """
#     GraphSATRE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphSATRE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, β::Float64)
    A, J = SAT.gen_randomKSAT(N, K, α)
    GraphRobustEnsemble(N, M, γ, β, GraphSAT, N, A, J)
end

function GraphSATRE(X::GraphSAT, M::Integer, γ::Float64, β::Float64)
    N = X.N
    GraphRobustEnsemble(N, M, γ, β, GraphSAT, N, X.A, X.J)
end

const GraphPercLinearRE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphPercLinear}

# """
#     GraphPercLinearRE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphPercLinearRE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64)
    ξ, ξv = PercLinear.gen_ξ(N, P)
    GraphRobustEnsemble(N, M, γ, β, GraphPercLinear, ξ, ξv)
end

function GraphPercLinearRE(X::GraphPercLinear, M::Integer, γ::Float64, β::Float64)
    GraphRobustEnsemble(X.N, M, γ, β, GraphPercLinear, X.ξ, X.ξv)
end

const GraphPercStepRE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphPercStep}

# """
#     GraphPercStepRE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphPercStepRE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64)
    ξ, ξv = PercStep.gen_ξ(N, P)
    GraphRobustEnsemble(N, M, γ, β, GraphPercStep, ξ, ξv)
end

function GraphPercStepRE(X::GraphPercStep, M::Integer, γ::Float64, β::Float64)
    GraphRobustEnsemble(X.N, M, γ, β, GraphPercStep, X.ξ, X.ξv)
end

const GraphCommStepRE{M,γ,β} = GraphRobustEnsemble{M,γ,β,GraphCommStep}

# """
#     GraphCommStepRE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommStepRE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv = CommStep.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphRobustEnsemble(N, M, γ, β, GraphCommStep, K2, ξ, ξv)
end

function GraphCommStepRE(X::GraphCommStep, M::Integer, γ::Float64, β::Float64)
    GraphRobustEnsemble(X.N, M, γ, β, GraphCommStep, X.K2, X.ξ, X.ξv)
end

# """
#     GraphCommReLURE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommReLURE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv, y = CommReLU.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphRobustEnsemble(N, M, γ, β, GraphCommReLU, K2, ξ, ξv, y)
end

function GraphCommReLURE(X::GraphCommReLU, M::Integer, γ::Float64, β::Float64)
    GraphRobustEnsemble(X.N, M, γ, β, GraphCommReLU, X.K2, X.ξ, X.ξv, X.y)
end

# """
#     GraphCommQuRE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommQuRE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv, y = CommQu.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphRobustEnsemble(N, M, γ, β, GraphCommQu, K2, ξ, ξv, y)
end

function GraphCommQuRE(X::GraphCommQu, M::Integer, γ::Float64, β::Float64)
    GraphRobustEnsemble(X.N, M, γ, β, GraphCommQu, X.K2, X.ξ, X.ξv, X.y)
end

end # module


text:
using Test, Pkg, CuArrays

ENV["JULIA_LOG_LEVEL"] = "WARN"

@test CuArrays.functional()

for submodule in [
    #"Common/MoistThermodynamics",
    #"Common/SurfaceFluxes",
    "Arrays",
    #"Numerics/Mesh",
    #"Numerics/DGmethods",
    "Numerics/ODESolvers",
]
    println("Starting tests for $submodule")
    t = @elapsed include(joinpath(submodule, "runtests.jl"))
    println("Completed tests for $submodule, $(round(Int, t)) seconds elapsed")
end


text:
#Covariance functions
export SECov

abstract type CovFunc end

#---------------------------------------
#Squared-exponential covariance function
#    separate length-scales per param
#---------------------------------------

type SECov{N} <: CovFunc
    sig_f::Float64
    lambdas::Vector{Float64}
    function SECov{N}(sf, l::Vector) where N
        @assert isa(N, Integer)
        @assert length(l) == N
        new(sf, l)
    end
end
SECov(p::Vector) = SECov{length(p)-1}(p[1], p[2:end])

function Base.show{N}(io::IO, K::SECov{N})
    print(io,"[sig_f: $(@sprintf("%3.3f",K.sig_f)), lambdas: ")
    _printvec(io, K.lambdas)
    print(io, "]")
end

#Calculate covariance for a pair of points
function (cf::SECov{N}){N}(x1::Vector, x2::Vector)
    @assert length(x1) == N
    @assert length(x2) == N
    cf.sig_f^2*exp(-0.5*dot(x1-x2, diagm(1./cf.lambdas)*(x1-x2)))
end

#Generate covariance matrix from sets of points
function (cf::SECov{N}){N}(X1, X2)
    @assert size(X1,1) == N
    @assert size(X2,1) == N
    D = diagm(1./cf.lambdas)
    K = zeros(size(X1,2),size(X2,2))
    for i=1:size(X1,2)
        for j=1:size(X2,2)
            dx = X1[:,i] - X2[:,j]
            K[i,j] = cf.sig_f^2*exp(-0.5*dot(dx, D*dx))
        end
    end
    K
end
(cf::SECov{N}){N}(X::Matrix) = cf(X,X)

#Return the derivative of a covariance matrix for a data set with respect to a
#single hyperparameter.
#Parameter i == 0: sig_f
function grad{N}(cf::SECov{N}, ::Type{Val{0}}, X)
    D = diagm(1./cf.lambdas)
    K = zeros(size(X,2),size(X,2))
    for i=1:size(X,2)
        for j=1:size(X,2)
            dx = X[:,i] - X[:,j]
            K[i,j] = 2*cf.sig_f*exp(-0.5*dot(dx, D*dx))
        end
    end
    K
end

#Parameter 1 <= i <= n: lambda_i
function grad{Nd,N}(cf::SECov{Nd}, ::Type{Val{N}}, X)
    @assert 1 <= N <= Nd
    D = diagm(1./cf.lambdas)
    K = zeros(size(X,2),size(X,2))
    for i=1:size(X,2)
        for j=1:size(X,2)
            dx = X[:,i] - X[:,j]
            K[i,j] = cf.sig_f^2*exp(-0.5*dot(dx, D*dx))*(dx[N]^2/(2*cf.lambdas[N]^2))
        end
    end
    K
end

"Compute the derivatives of the test point covariance vector w.r.t. the
prediction variables"
function covgrad{Nd}(x::Vector, cf::SECov{Nd}, X)
    @assert length(x) == Nd
    @assert size(X,1) == Nd
    Kd = zeros(size(X,2), Nd)
    D = diagm(1./cf.lambdas)
    for j=1:Nd
        ej = zeros(Nd)
        ej[j] = 1.0
        for i=1:size(X,2)
            dx = X[:,i] - x
            Kd[i,j] = cf.sig_f^2*exp(-0.5*dot(dx,D*dx))*dot(dx,D*ej)
        end
    end
    Kd
end

#The squared-exponential symmetric covariance function
type SECovSym{N} <: CovFunc
    sig_f::Float64
    lambda::Float64
    function SECovSym{N}(sf, l) where N
        @assert isa(N, Integer)
        new(sf, l)
    end
end
SECovSym(N, p::Vector) = SECovSym{N}(p[1], p[2])

function Base.show{N}(io::IO, K::SECovSym{N})
    print(io, @sprintf("[sig_f: %3.3f, lambda: %3.3f]", K.sig_f, K.lambda))
end

#Calculate covariance for a pair of points
function (cf::SECovSym{N}){N}(x1::Vector, x2::Vector)
    @assert length(x1) == N
    @assert length(x2) == N
    cf.sig_f^2*exp(dot(-0.5*(x1-x2),diagm(1./(cf.lambda*ones(N)))*(x1-x2)))
end

#Generate covariance matrix from sets of points
function (cf::SECovSym{N}){N}(X1, X2)
    @assert size(X1,1) == N
    @assert size(X2,1) == N
    D = diagm(1./(cf.lambda*ones(N)))
    K = zeros(size(X1,2),size(X2,2))
    for i=1:size(X1,2)
        for j=1:size(X2,2)
            dx = X1[:,i] - X2[:,j]
            K[i,j] = cf.sig_f^2*exp(-0.5*dot(dx, D*dx))
        end
    end
    K
end
(cf::SECovSym{N}){N}(X::Matrix) = cf(X,X)

#Return the derivative of a covariance matrix for a data set with respect to a
#single hyperparameter.
#Parameter i == 0: sig_f
function grad{N}(cf::SECovSym{N}, ::Type{Val{0}}, X)
    D = diagm(1./(cf.lambda*ones(N)))
    K = zeros(size(X,2),size(X,2))
    for i=1:size(X,2)
        for j=1:size(X,2)
            dx = X[:,i] - X[:,j]
            K[i,j] = 2*cf.sig_f*exp(-0.5*dot(dx, D*dx))
        end
    end
    K
end

#Parameter 1 : lambda
#TODO: recalculate
function grad{N}(cf::SECovSym{N}, ::Type{Val{1}}, X)
    D = diagm(1./(cf.lambda*ones(N)))
    K = zeros(size(X,2),size(X,2))
    for i=1:size(X,2)
        for j=1:size(X,2)
            dx = X[:,i] - X[:,j]
            K[i,j] = cf.sig_f^2*exp(-0.5*dot(dx, D*dx))*(dx[N]^2/(2*cf.lambda^2))
        end
    end
    K
end


text:
@testset "isbinaryor.jl" begin

    @testset "isBinaryOr()" begin
        trailer = SeaPearl.Trailer()

        b = SeaPearl.BoolVar("b", trailer)
        x = SeaPearl.BoolVar("x", trailer)
        y = SeaPearl.BoolVar("y", trailer)

        constraint = SeaPearl.isBinaryOr(b, x, y, trailer)

        @test constraint in b.onDomainChange
        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::isBinaryOr)" begin

        @testset "b false - x,y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test !SeaPearl.assignedValue(x)
            @test !SeaPearl.assignedValue(y)
            @test prunedDomains == SeaPearl.CPModification("x" => [true], "y" => [true])
            @test !constraint.active.value
        end

        @testset "b false - x false, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, false)
            SeaPearl.assign!(x, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test !SeaPearl.assignedValue(y)
            @test prunedDomains == SeaPearl.CPModification("y" => [true])
            @test !constraint.active.value

        end

        @testset "b false - x true, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, false)
            SeaPearl.assign!(x, true)

            @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        end

        @testset "b true - x false, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, true)
            SeaPearl.assign!(x, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test SeaPearl.assignedValue(y)
            @test prunedDomains == SeaPearl.CPModification("y" => [false])
            @test !constraint.active.value

        end

        @testset "b true - x false, y false" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, true)
            SeaPearl.assign!(x, false)
            SeaPearl.assign!(y, false)

            @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        end

        @testset "b true - x true, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, true)
            SeaPearl.assign!(x, true)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test prunedDomains == SeaPearl.CPModification()
            @test !constraint.active.value

        end

        @testset "b unassigned - x true, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(x, true)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test !SeaPearl.isbound(y)
            @test SeaPearl.assignedValue(b)
            @test prunedDomains == SeaPearl.CPModification("b" => [false])
            @test !constraint.active.value

        end

        @testset "b unassigned - x false, y false" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(x, false)
            SeaPearl.assign!(y, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test !SeaPearl.assignedValue(b)
            @test prunedDomains == SeaPearl.CPModification("b" => [true])
            @test !constraint.active.value

        end

        @testset "b unassigned - x false, y unassigned" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.BoolVar("x", trailer)
            y = SeaPearl.BoolVar("y", trailer)

            constraint = SeaPearl.isBinaryOr(b, x, y, trailer)
            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(x, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
            @test !SeaPearl.isbound(b)
            @test !SeaPearl.isbound(y)
            @test prunedDomains == SeaPearl.CPModification()
            @test constraint.active.value

        end

    end
end


text:
# The R Graphics book has a nice example of styling circles with the following code
#```R
# > levels <- round(seq(90, 10, length=25))
# > greys <- paste("grey", c(levels, rev(levels)), sep="")
# > grid.circle(x=seq(0.1, 0.9, length=100),
#               y=0.5 + 0.4*sin(seq(0, 2*pi, length=100)),
#               r=abs(0.1*cos(seq(0, 2*pi, length=100))),
#               gp=gpar(col=greys))
#```
#
# I don't love the vectorization of the circle creation, but I should be able to do
# this effect
using OdinSon
using Colors

# not really the same as getting the gray color names like above, but I should be able to
# tweak this to get a nice version
##TODO: so so ugly, figure out how to cycle values in a nice way
# The key for the R code is that is ramps up for 25 then reverse down for 25 and repeats
ncirc = 100
baseseq = linspace(NC"white", NC"black", 25)
grays = cycle(vcat(baseseq, reverse(baseseq))) # this implicitly gets recycled 4 times in R code
xs = linspace(0.1, 0.9, ncirc)
Θ = linspace(0, 2*π, ncirc)
ys = 0.5 + 0.4*sin(Θ)
rs = abs(0.1*cos(Θ))

circles = Array{Circle}(ncirc)
for (i, col) in enumerate(take(grays, ncirc)) # I do this to get the repeats working, think of cleaner way
    circles[i] = Circle((xs[i], ys[i]), rs[i], style=Style(stroke=col, fill=nothing))
end
render(Canvas(circles))

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
"""
Compute orientations of simplices given the `points` forming the vertices of the simplices
(a n_vertices x dim array) and `simplex_inds` (a n_simplices x (dim+1) array), telling how
the simplices are formed from the vertices.
"""
function orientations(points::AbstractArray{Float64, 2}, simplex_inds::AbstractArray{Int, 2})
    n_simplices = size(simplex_inds, 1)
    dim = size(simplex_inds, 2) - 1
    orientations = Vector{Float64}(n_simplices)

    for i = 1:n_simplices
        orientations[i] = det(hcat(view(points, view(simplex_inds, i, :), :), ones(dim + 1)))
    end

    return orientations
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Order

import ..@__MODULE__, ..parentmodule
const Base = parentmodule(@__MODULE__)
import .Base:
    AbstractVector, @propagate_inbounds, isless, identity, getindex,
    +, -, !, &, <, |

## notions of element ordering ##

export # not exported by Base
    Ordering, Forward, Reverse,
    By, Lt, Perm,
    ReverseOrdering, ForwardOrdering,
    DirectOrdering,
    lt, ord, ordtype

abstract type Ordering end

struct ForwardOrdering <: Ordering end
struct ReverseOrdering{Fwd<:Ordering} <: Ordering
    fwd::Fwd
end

ReverseOrdering(rev::ReverseOrdering) = rev.fwd
ReverseOrdering(fwd::Fwd) where {Fwd} = ReverseOrdering{Fwd}(fwd)

const DirectOrdering = Union{ForwardOrdering,ReverseOrdering{ForwardOrdering}}

const Forward = ForwardOrdering()
const Reverse = ReverseOrdering(Forward)

struct By{T} <: Ordering
    by::T
end

struct Lt{T} <: Ordering
    lt::T
end

struct Perm{O<:Ordering,V<:AbstractVector} <: Ordering
    order::O
    data::V
end

lt(o::ForwardOrdering,       a, b) = isless(a,b)
lt(o::ReverseOrdering,       a, b) = lt(o.fwd,b,a)
lt(o::By,                    a, b) = isless(o.by(a),o.by(b))
lt(o::Lt,                    a, b) = o.lt(a,b)

@propagate_inbounds function lt(p::Perm, a::Integer, b::Integer)
    da = p.data[a]
    db = p.data[b]
    lt(p.order, da, db) | (!lt(p.order, db, da) & (a < b))
end

ordtype(o::ReverseOrdering, vs::AbstractArray) = ordtype(o.fwd, vs)
ordtype(o::Perm,            vs::AbstractArray) = ordtype(o.order, o.data)
# TODO: here, we really want the return type of o.by, without calling it
ordtype(o::By,              vs::AbstractArray) = try typeof(o.by(vs[1])) catch; Any end
ordtype(o::Ordering,        vs::AbstractArray) = eltype(vs)

_ord(lt::typeof(isless), by::typeof(identity), order::Ordering) = order
_ord(lt::typeof(isless), by,                   order::Ordering) = By(by)
_ord(lt,                 by::typeof(identity), order::Ordering) = Lt(lt)
_ord(lt,                 by,                   order::Ordering) = Lt((x,y)->lt(by(x),by(y)))

ord(lt, by, rev::Nothing, order::Ordering=Forward) = _ord(lt, by, order)

function ord(lt, by, rev::Bool, order::Ordering=Forward)
    o = _ord(lt, by, order)
    return rev ? ReverseOrdering(o) : o
end

end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# Compare a local instant to a UTC transition instant by using the offset to make them both
# into local time. We could just as easily convert both of them into UTC time.
lt_local(local_dt::DateTime, t::Transition) = isless(local_dt, t.utc_datetime + t.zone.offset)
lt_local(t::Transition, local_dt::DateTime) = isless(t.utc_datetime + t.zone.offset, local_dt)

lt_utc(utc_dt::DateTime, t::Transition) = isless(utc_dt, t.utc_datetime)
lt_utc(t::Transition, utc_dt::DateTime) = isless(t.utc_datetime, utc_dt)

function transition_range(local_dt::DateTime, tz::VariableTimeZone, ::Type{Local})
    # To understand the logic in this function some background on transitions is needed:
    #
    # A transition (`t[i]`) is applicable to a given UTC instant that occurs on or after the
    # transition start (`t[i].utc_datetime`). The transition (`t[i]`) ends at the start of
    # the next transition in the list (`t[i + 1].utc_datetime`).
    #
    # Any UTC instant that occurs prior to the first transition (`t[1].utc_datetime`) has no
    # associated transitions. Any UTC instant that occurs on or after the last transition
    # (`t[end].utc_datetime`) is associated, at a minimum, with the last transition.

    # Determine the latest transition that applies to `local_dt`. If the `local_dt`
    # preceeds all transitions `finish` will be zero and produce the empty range `1:0`.
    finish = searchsortedlast(tz.transitions, local_dt, lt=lt_local)

    # Usually we'll begin by having `start` be larger than `finish` to create an empty
    # range by default. In the scenario where last transition applies to the `local_dt` we
    # can avoid a bounds by setting `start = finish`.
    len = length(tz.transitions)
    start = finish < len ? finish + 1 : len

    # To determine the first transition that applies to the `local_dt` we will work
    # backwards. Typically, this loop will only use single iteration as multiple iterations
    # only occur when local times are ambiguous.
    @inbounds for i in (start - 1):-1:1
        # Compute the end of the transition in local time. Note that this instant is not
        # included in the implicitly defined transition interval (known as right-open in
        # interval parlance).
        transition_end = tz.transitions[i + 1].utc_datetime + tz.transitions[i].zone.offset

        # If the end of the transition occurs after the `local_dt` then this transition
        # applies to the `local_dt`.
        if transition_end > local_dt
            start = i
        else
            break
        end
    end

    return start:finish
end

function transition_range(utc_dt::DateTime, tz::VariableTimeZone, ::Type{UTC})
    finish = searchsortedlast(tz.transitions, utc_dt, lt=lt_utc)
    start = max(finish, 1)
    return start:finish
end

"""
    transition_range(dt::DateTime, tz::VariableTimeZone, context::Type{Union{Local,UTC}}) -> UnitRange

Finds the indexes of the `tz` transitions which may be applicable for the `dt`. The given
DateTime is expected to be local to the time zone or in UTC as specified by `context`. Note
that UTC context will always return a range of length one.
"""
transition_range(::DateTime, ::VariableTimeZone, ::Type{Union{Local,UTC}})

function interpret(local_dt::DateTime, tz::VariableTimeZone, ::Type{Local})
    t = tz.transitions
    r = transition_range(local_dt, tz, Local)

    possible = (ZonedDateTime(local_dt - t[i].zone.offset, tz, t[i].zone) for i in r)
    return IndexableGenerator(possible)
end

function interpret(utc_dt::DateTime, tz::VariableTimeZone, ::Type{UTC})
    t = tz.transitions
    r = transition_range(utc_dt, tz, UTC)
    length(r) == 1 || error("Internal TimeZones error: A UTC DateTime should only have a single interpretation")

    possible = (ZonedDateTime(utc_dt, tz, t[i].zone) for i in r)
    return IndexableGenerator(possible)
end

"""
    interpret(dt::DateTime, tz::VariableTimeZone, context::Type{Union{Local,UTC}}) -> Array{ZonedDateTime}

Produces a list of possible `ZonedDateTime`s given a `DateTime` and `VariableTimeZone`.
The result will be returned in chronological order. Note that `DateTime`s in the local
context typically return 0-2 results while the UTC context will always return 1 result.
"""
interpret(::DateTime, ::VariableTimeZone, ::Type{Union{Local,UTC}})

"""
    shift_gap(local_dt::DateTime, tz::VariableTimeZone) -> Tuple

Given a non-existent local `DateTime` in a `TimeZone` produces a tuple containing two valid
`ZonedDateTime`s that span the gap. Providing a valid local `DateTime` returns an empty
tuple. Note that this function does not support passing in a UTC `DateTime` since there are
no non-existent UTC `DateTime`s.

Aside: the function name refers to a period of invalid local time (gap) caused by daylight
saving time or offset changes (shift).
"""
function shift_gap(local_dt::DateTime, tz::VariableTimeZone)
    r = transition_range(local_dt, tz, Local)
    boundaries = if isempty(r) && last(r) > 0
        t = tz.transitions
        i, j = last(r), first(r)  # Empty range has the indices we want but backwards
        tuple(
            ZonedDateTime(t[i + 1].utc_datetime - eps(local_dt), tz, t[i].zone),
            ZonedDateTime(t[j].utc_datetime, tz, t[j].zone),
        )
    else
        tuple()
    end

    return boundaries
end

"""
    first_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)

Construct a valid `ZonedDateTime` by adjusting the local `DateTime`. If the local `DateTime`
is non-existent then it will be adjusted using the `step` to be *after* the gap. When the
local `DateTime` is ambiguous the *first* ambiguous `DateTime` will be returned.
"""
function first_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)
    possible = interpret(local_dt, tz, Local)

    # Skip all non-existent local datetimes.
    while isempty(possible)
        local_dt += step
        possible = interpret(local_dt, tz, Local)
    end

    return first(possible)
end

"""
    last_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)

Construct a valid `ZonedDateTime` by adjusting the local `DateTime`. If the local `DateTime`
is non-existent then it will be adjusted using the `step` to be *before* the gap. When the
local `DateTime` is ambiguous the *last* ambiguous `DateTime` will be returned.
"""
function last_valid(local_dt::DateTime, tz::VariableTimeZone, step::Period)
    possible = interpret(local_dt, tz, Local)

    # Skip all non-existent local datetimes.
    while isempty(possible)
        local_dt -= step
        possible = interpret(local_dt, tz, Local)
    end

    return last(possible)
end

function first_valid(local_dt::DateTime, tz::VariableTimeZone)
    possible = interpret(local_dt, tz, Local)
    return isempty(possible) ? last(shift_gap(local_dt, tz)) : first(possible)
end

function last_valid(local_dt::DateTime, tz::VariableTimeZone)
    possible = interpret(local_dt, tz, Local)
    return isempty(possible) ? first(shift_gap(local_dt, tz)) : last(possible)
end


text:
function form_lcp(J, v, body, physics)
    μ = body.μ
    M = body.M
    Mi = body.M_inv
    ϵ = body.ϵ
    dt = physics.dt
    g = physics.g

    Jₙ = J[:,2]
    Jₜ = J[:,1]

    V = zeros((4,4))
    V[1,1] = Jₙ'* Mi * (dt*Jₙ)
    V[1,2] = -Jₙ'* Mi * (dt*Jₜ)
    V[1,3] = Jₙ'* Mi * (dt*Jₜ)
    V[2,1] = -Jₜ'* Mi * (dt*Jₙ)
    V[2,2] = Jₜ'* Mi * (dt*Jₜ)
    V[2,3] = -Jₜ'* Mi * (dt*Jₜ)
    V[2,4] = 1.
    V[3,1] = Jₜ'* Mi * (dt*Jₙ)
    V[3,2] = -Jₜ'* Mi * (dt*Jₜ)
    V[3,3] = Jₜ'* Mi * (dt*Jₜ)
    V[3,4] = 1.
    V[4,1] = μ
    V[4,2] = -1.
    V[4,3] = -1.
    fₑ = M*g
    p = zeros(4)
    p[1] = Jₙ' * ((1+ϵ)*v + dt*Mi*fₑ)
    p[2] = -Jₜ' * (v + dt*Mi*fₑ)
    p[3] = Jₜ' * (v + dt*Mi*fₑ)

    return V, p
end

function solve_lcp(V, p)
    lcp = Model(Ipopt.Optimizer)
    N = length(p)
    @variable(lcp, x[1:N])
    @constraint(lcp,  V*x + p .>= 0.)
    @constraint(lcp, x .>= 0.)
    @objective(lcp, Min, x'*(V*x + p))
    optimize!(lcp)

    return value.(x)

    # result = solve!(LCP(M, q))
    # return result.sol
end