text:
module Mesh

using DocStringExtensions

export Mesh1D,
    Mesh2D,
    Mesh3D,
    AbstractMesh,
    equispaced_line_mesh,
    equispaced_rectangular_mesh,
    sphere_mesh,
    cube_panel_mesh,
    AbstractSphereWarp,
    EquiangularSphereWarp,
    EquidistantSphereWarp,
    unfold_cube_panel_to_plane

abstract type AbstractMesh{I,FT} end

abstract type AbstractSphereWarp end
struct EquiangularSphereWarp <: AbstractSphereWarp end
struct EquidistantSphereWarp <: AbstractSphereWarp end
"""
    Mesh1D{I,IA1D,IA2D,FT,FTA2D,SNT,NB} <: AbstractMesh{I,FT}

Conformal mesh for a 1D manifold. The manifold can be embedded in a
higher dimensional space.

For 1D meshes, faces are identical to meshes.
"""
struct Mesh1D{I,IA1D,IA2D,FT,FTA2D,SNT,NB} <: AbstractMesh{I,FT}
    "number of unique vertices in the mesh"
    nverts::I
    "number of elements in the mesh"
    nelems::I
    "number of zones in the mesh"
    nbndry::I
    "x₁, x₂, ... coordinates of vertices `(nverts, dim)`, dim can be greater than 1 for 1D manifolds embedded in higher dimensional space"
    coordinates::FTA2D
    "unique vertices `(n_uniquevertices)`"
    unique_verts::IA1D
    "vertex neighbors"
    vertex_neighbors::IA2D
    "vertex numbers on each boundary"
    vertex_boundary::IA1D
    "boundary tags"
    boundary_tags::IA1D
    "boundary tag names"
    boundary_tag_names::SNT
    "vertex boundary offset"
    vertex_boundary_offset::IA1D
    "vertices numbers for each elem `(nelems, 2)`"
    elem_verts::IA2D
end

Mesh1D(
    nverts,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    vertex_neighbors,
    vertex_boundary,
    boundary_tags,
    boundary_tag_names,
    vertex_boundary_offset,
    elem_verts,
) = Mesh1D{
    eltype(nverts),
    typeof(unique_verts),
    typeof(elem_verts),
    eltype(coordinates),
    typeof(coordinates),
    typeof(boundary_tag_names),
    length(boundary_tag_names),
}(
    nverts,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    vertex_neighbors,
    vertex_boundary,
    boundary_tags,
    boundary_tag_names,
    vertex_boundary_offset,
    elem_verts,
)
"""
    Mesh2D{I,IA2D,FT,FTA2D} <: AbstractMesh{FT}

Conformal mesh for a 2D manifold. The manifold can be
embedded in a higher dimensional space.

                        Quadrilateral

                v4            f3           v3
                  o------------------------o
                  |                        |      face    vertices
                  |       <----------      |
                  |                 |      |        f1 =>  v1 v2
               f4 |                 |      | f2     f2 =>  v2 v3
                  |                 |      |        f3 =>  v3 v4
                  |                 |      |        f4 =>  v4 v1
                  |       -----------      |
                  |                        |
                  o------------------------o
                 v1           f1           v2

Reference:
https://gsjaardema.github.io/seacas-docs/html/element_types.html#ordering

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Mesh2D{I,IA1D,IA2D,FT,FTA2D,SNT,NB} <: AbstractMesh{I,FT}
    "number of unique vertices in the mesh"
    nverts::I
    "number of unique faces in the mesh"
    nfaces::I
    "number of elements in the mesh"
    nelems::I
    "number of zones in the mesh"
    nbndry::I
    "x₁, x₂, ... coordinates of vertices `(nverts, dim)`, dim can be greater than 2 for 2D manifolds embedded in higher dimensional space"
    coordinates::FTA2D
    "unique vertices `(n_uniquevertices)`"
    unique_verts::IA1D
    "connectivity information for unique vertices"
    uverts_conn::IA1D
    "offset information for uverts_conn `(n_unique_verts + 1)`"
    uverts_offset::IA1D
    "face vertices numbers `(nfaces, 2)`"
    face_verts::IA2D
    "boundary elems for each face `(nfaces, 5)` -> [elem1, localface1, elem2, localface2, relative orientation]"
    face_neighbors::IA2D
    "face numbers on each boundary `(nfaces)`"
    face_boundary::IA1D
    "boundary tags"
    boundary_tags::IA1D
    "boundary tag names"
    boundary_tag_names::SNT
    "face boundary offset for each boundary for face_boundary array"
    face_boundary_offset::IA1D
    "vertices numbers for each elem `(nelems, 4)`"
    elem_verts::IA2D
    "face numbers for each elem `(nelems, 4)`"
    elem_faces::IA2D
end

Mesh2D(
    nverts,
    nfaces,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    uverts_conn,
    uverts_offset,
    face_verts,
    face_neighbors,
    face_boundary,
    boundary_tags,
    boundary_tag_names,
    face_boundary_offset,
    elem_verts,
    elem_faces,
) = Mesh2D{
    eltype(nverts),
    typeof(face_boundary),
    typeof(face_verts),
    eltype(coordinates),
    typeof(coordinates),
    typeof(boundary_tag_names),
    length(boundary_tag_names),
}(
    nverts,
    nfaces,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    uverts_conn,
    uverts_offset,
    face_verts,
    face_neighbors,
    face_boundary,
    boundary_tags,
    boundary_tag_names,
    face_boundary_offset,
    elem_verts,
    elem_faces,
)

include("line_mesh.jl")
include("box_mesh.jl")
include("sphere_mesh.jl")
include("mesh_utilities.jl")
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
using Eratosthenes
using Test
using LinearAlgebra

include(joinpath(dirname(pathof(Eratosthenes)), "EratosthenesRotations.jl"))
using .EratosthenesRotations

########################
# Individual Rotations #
########################

# Test that qcomp and qdiff are consistent.
a = normalize([1., 2., 3., 4.])
b = normalize([2., 2., 1., -2.])
Δθ, r = q2aa(qpos(qdiff(qdiff(a, b), qcomp(a, qinv(b)))))
@test Δθ < 1e-12

##############
# Simulation #
##############

# While this isn't much of a unit test, the following exercises a huge amount of
# the system with nonetheless easily-predictable results. Let's see if things
# work.

# Create an actuator that provides a constant torque.
function Constorquer()
    function effects(t, constants, state, draws, implicit, inputs, effects, effects_bus)
        return (BodyTorque(constants),)
    end
    DynamicalModel("constorquer",
                   effects = effects,
                   constants = [2.; 0.; 0.], # Torque (must be about principal axis)
                   timing = ModelTiming(0.01))
end

# Create a default simulation with only a constorquer. We'll be able to predict
# the end result.
scenario = Scenario()
scenario.sim.t_end = 5.
scenario.vehicles[1].components = [Constorquer()]
simulate(scenario)

# Test that we rolled just as much as expected.
τ             = scenario.vehicles[1].components[1].constants
Δt            = scenario.sim.t_end
I_B           = scenario.vehicles[1].body.constants.I_B
q_expected    = aa2q(0.5 * Δt^2 * norm(I_B \ τ), normalize(τ))
ω_expected    = Δt * (I_B \ τ)

# Calculate the errors.
q     = normalize(scenario.vehicles[1].body.state.q_BI)
ω     = scenario.vehicles[1].body.state.ω_BI_B
Δθ, r = q2aa(qpos(qdiff(q, q_expected)))
Δω    = ω - ω_expected

# How'd we do?
@test Δθ < 1e-12
@test norm(Δω) < 1e-12

###########
# Example #
###########

# If the out directory doesn't exist, make it.
if !isdir("out")
    mkdir("out")
end

# min-scenario
success = false;
try
    include("run-min-scenario-ci.jl")
    global success = true;
catch err
    println("There was an error running the basic example.")
    rethrow(err)
end

@test success == true

# coupled reaction wheel scenario
success = false;
try
    local scenario = setup(joinpath(dirname(pathof(Eratosthenes)), "..", "examples", "scenarios", "reaction-wheel.yaml"))
    scenario.sim.t_end = 1.
    simulate(scenario)
    global success = true;
catch err
    println("There was an error running the reaction wheel example.")
    rethrow(err)
end

@test success == true

ext: jl


text:
module GtkApps

include("Viewer.jl")

end # module


text:
module BLP

import LinearAlgebra as LA
using LinearAlgebra: LowerTriangular, ⋅, mul!

include("iteration.jl")
include("softmax.jl")
include("market.jl")

end


text:
using UnitfulColors
using Plots, Unitful, UnitfulRecipes, Colors, Latexify, UnitfulLatexify

function showUnitColor(title, x, varargs...)
    return heatmap(
        x,
        [-1, 1],
        repeat(x, 2, 1);
        background=:black,
        color=colormatch.(x, varargs...),
        legend=false,
        colorbar=false,
        xlabel="\\textrm{$title}",
        yticks=[],
        yshowaxis=false,
        fontfamily="Computer Modern",
        unitformat=latexify,
    )
end

length = 1000;
pl = plot(
    showUnitColor("Wavelength", range(350, 750; length) .* u"nm"),
    showUnitColor("Frequency", range(400, 800; length) .* u"THz"),
    showUnitColor("Energy", range(1.5, 3.5; length) .* u"eV"),
    showUnitColor("Temperature (peak)", range(4000, 8000; length) .* u"K"),
    showUnitColor("Temperature (spectrum)", range(1000, 27000; length) .* u"K", :spectrum);
    layout=(5, 1),
);

savefig(pl, "$(pkgdir(UnitfulColors))/doc/examples.png")


text:
using OrdinaryDiffEq
using Trixi

equations = CompressibleEulerEquations2D(1.4) # set gas gamma = 1.4

# freestream flow state with Ma_inf = 0.3
@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D)

  # set the freestream flow parameters
  rho_freestream = 1.0
  u_freestream = 0.3
  p_freestream = inv(equations.gamma)

  theta = 0.0 # zero angle of attack
  si, co = sincos(theta)
  v1 = u_freestream * co
  v2 = u_freestream * si

  prim = SVector(rho_freestream, v1, v2, p_freestream)
  return prim2cons(prim, equations)
end

# initial condition
initial_condition = uniform_flow_state

# boundary condition types
boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)

# boundary condition dictionary
boundary_conditions = Dict( :Bottom => boundary_condition_uniform_flow,
                            :Top    => boundary_condition_uniform_flow,
                            :Right  => boundary_condition_uniform_flow,
                            :Left   => boundary_condition_uniform_flow,
                            :Ellipse => boundary_condition_slip_wall,
                            :Circle => boundary_condition_slip_wall )

# DGSEM solver.
#    1) polydeg must be >= the polynomial order set in the HOHQMesh control file to guarantee
#       freestream preservation. As a extra task try setting polydeg=3
#    2) VolumeIntegralFluxDifferencing with central volume flux is activated
#       for dealiasing
volume_flux = flux_ranocha
solver = DGSEM(polydeg=4, surface_flux=flux_hll,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# create the unstructured from your mesh file
#mesh_file = joinpath(@__DIR__, "out/box_with_object.mesh")
mesh_file = joinpath(@__DIR__, "out/box_with_two_objects.mesh")
mesh = UnstructuredMesh2D(mesh_file)

# Create semidiscretization with all spatial discretization-related components
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

# Create ODE problem from semidiscretization with time span from 0.0 to 2.0
tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

# Create the callbacks to output solution files and adapt the time step
summary_callback = SummaryCallback()
save_solution = SaveSolutionCallback(interval=10,
                                     save_initial_solution=true,
                                     save_final_solution=true)
stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback)

# Evolve ODE problem in time using `solve` from OrdinaryDiffEq
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
# print the timer summary
summary_callback()


text:
function bracket_minimum(f, x=0.0, s=1e-2, k=2)

end

function minimize(f, a, b)
    # Brent-Dekker method

end

function line_search(f, x, d)
    obj = α -> f(x + α*d);
    a, b = bracket_minimum(obj);
    α = minimize(obj, a, b);
    return x + α * d;
end


text:
using Plots
using LaTeXStrings
pgfplots()

c = 1
a₀ = 1
ω = 1

x(τ) = c*(a₀^2)/4 *(τ - 1/2* ω*sin(2ω*τ))
y(τ) = c*a₀/ω *(1-cos(ω*τ))

p1 = plot(x, y, 0, 15,
  xlabel=L"x", ylabel=L"y",
  xticks=nothing, yticks=nothing,
  legend=:none,
  framestyle=:box,
  tex_output_standalone=true
)

v_drift = c*a₀^2/4

p2 = plot(τ->x(τ)-v_drift*τ, y, 1, 10,
  xlabel=L"x", ylabel=L"y",
  xticks=nothing, yticks=nothing,
  legend=:none,
  framestyle=:box,
  tex_output_standalone=true
)

p = plot(p1, p2,
  layout=grid(1,2,widths=[0.7,0.3]),
  framestyle=:box,
  tex_output_standalone=true
)

savefig(p, "figures/electron-in-plane-wave.tex")

max_stars_repo_name: SebastianM-C/MasterThesis

max_issues_repo_name: SebastianM-C/MasterThesis

max_forks_repo_name: SebastianM-C/MasterThesis


text:
# Evaluating heat capacities

using WaterData
# needed so that WaterData is able to load files
import JLD, Dierckx, VoronoiDelaunay, GeometricalPredicates

# add an override to call heat capacities with ValueSets
for HC in (WaterData.PTFuncHeatCapacity, WaterData.GridHeatCapacity,
           WaterData.TFuncHeatCapacity, WaterData.ConstantHeatCapacity)
    (cp::HC)(pv::PhysicalValues) = cp(pressure(pv), temperature(pv))
end


text:
@doc raw"""
    p = characteristic_polynomial(g_ol)

Determine the characteristic polynomial associated with open loop transfer function `g_ol`.

The characteristic polynomial is $1+g_{ol}(s)$. The roots of the characteristic polynomial determine the character of the response of the closed loop system to bounded inputs.

# Arguments
* `g_ol::TransferFunction`: open loop transfer function

# Returns
a polynomial of type `Polynomial`

# Example
```jldoctest
g_ol = 4 / (s + 3) / (s + 2) / (s + 1)
characteristic_polynomial(g_ol)
# output
Polynomials.Polynomial(10.0 + 11.0*s + 6.0*s^2 + 1.0*s^3)
```
"""
function characteristic_polynomial(g_ol::TransferFunction)
    if g_ol.time_delay != 0.0
        error("not implemented for time delay")
    end
    # TODO wut if time delay
    return g_ol.denominator + g_ol.numerator
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
#
# Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
# Copyright 2015 Oliver Heimlich (oheim@posteo.de)
#
# Original author: Marco Nehmeier (unit tests in libieeep1788)
# Converted into portable ITL format by Oliver Heimlich with minor corrections.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#Language imports

#Test library imports
using FactCheck

#Arithmetic library imports
using ValidatedNumerics

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_sqrRev_test") do

end

facts("minimal_sqrRevBin_test") do

end

facts("minimal_sqrRev_dec_test") do

end

facts("minimal_sqrRev_dec_bin_test") do

end

facts("minimal_absRev_test") do

end

facts("minimal_absRevBin_test") do

end

facts("minimal_absRev_dec_test") do

end

facts("minimal_absRev_dec_bin_test") do

end

facts("minimal_pownRev_test") do

end

facts("minimal_pownRevBin_test") do

end

facts("minimal_pownRev_dec_test") do

end

facts("minimal_pownRev_dec_bin_test") do

end

facts("minimal_sinRev_test") do

end

facts("minimal_sinRevBin_test") do

end

facts("minimal_sinRev_dec_test") do

end

facts("minimal_sinRev_dec_bin_test") do

end

facts("minimal_cosRev_test") do

end

facts("minimal_cosRevBin_test") do

end

facts("minimal_cosRev_dec_test") do

end

facts("minimal_cosRev_dec_bin_test") do

end

facts("minimal_tanRev_test") do

end

facts("minimal_tanRevBin_test") do

end

facts("minimal_tanRev_dec_test") do

end

facts("minimal_tanRev_dec_bin_test") do

end

facts("minimal_coshRev_test") do

end

facts("minimal_coshRevBin_test") do

end

facts("minimal_coshRev_dec_test") do

end

facts("minimal_coshRev_dec_bin_test") do

end

facts("minimal_mulRev_test") do

end

facts("minimal_mulRevTen_test") do

end

facts("minimal_mulRev_dec_test") do

end

facts("minimal_mulRev_dec_ten_test") do

end
# FactCheck.exitstatus()


text:
# Dantzig-Wolfe Decomposition
# Edward J. Xu
# 2019.5.28
# include("$(homedir())/Documents/Github/DantzigWolfeDecomposition/Case/TDSP/main.jl")
########################################################################################################################
push!(LOAD_PATH, "$(homedir())//Documents/Github/DantzigWolfeDecomposition/Case/SSCFLP")
cd("$(homedir())//Documents/Github/DantzigWolfeDecomposition/Case/SSCFLP")
using JuMP
# using CPLEX
using Gurobi
using LinearAlgebra
using Random
# using DantzigWolfeOptim
include("FuncSub.jl")
include("FuncMas.jl")
include("FuncOptim.jl")
include("FuncData.jl")
include("DantzigWolfeOptim.jl")
########################################################################################################################

function main()
    ## 1,  Data
    ## num, widthMax, heightMax, wCap, seed
    ## Some instances may be difficult, both for CPLEX/Gurobi and for the Dantzig-Wolfe approach.
    # (wCap, vec_w, vec_h) = getStripPack(3, 10, 10, 11, 2)
    (wCap, vec_w, vec_h) = getStripPack(4, 10, 10, 11, 2)
    # (wCap, vec_w, vec_h) = getStripPack(5, 10, 10, 12, 3)
    #(wCap, vec_w, vec_h) = getStripPack(6, 10, 10, 13, 4)
    #(wCap, vec_w, vec_h) = getStripPack(7, 10, 10, 20, 5)
    # (wCap, vec_w, vec_h) = getStripPack(8, 10, 10, 20, 6)
    #(wCap, vec_w, vec_h) = getStripPack(9, 10, 10, 20, 7)
    #(wCap, vec_w, vec_h) = getStripPack(10, 10, 10, 20, 8)
    #(wCap, vec_w, vec_h) = getStripPack(11, 10, 10, 20, 9)
    #(wCap, vec_w, vec_h) = getStripPack(12, 10, 10, 20, 1)
    #(wCap, vec_w, vec_h) = getStripPack(13, 10, 10, 20, 10)
    #(wCap, vec_w, vec_h) = getStripPack(14, 10, 10, 20, 11)
    # (wCap, vec_w, vec_h) = getStripPack(15, 10, 10, 20, 12)
    # println("wCap = $(wCap)")
    # println("vec_w = $(vec_w)")
    # println("vec_h = $(vec_h)")
    # solveStripPack(wCap, vec_w, vec_h)
    (vec_c, mat_a, vec_b, vecVecIndexSub, numSub, m_mat_h, vecRowMatD, vecVecIndexBinInSub) = getCoefMatrix(
        wCap, vec_w, vec_h
        )
    if numSub == 1
        solveLP(vec_c, mat_a, vec_b, vecVecIndexBinInSub)
    end
    ## Test the algorihtm using SSCFLP
    # (vec_c, mat_a, vec_b, vecVecIndexSub, numSub, m_mat_h, vecRowMatD, vecVecIndexBinInSub) = getData()
    ## 2,  Optim
    wheBranch = false
    vecWhiSub = [1]
    vecWhiVar = [1]
    vecWhiBranch = [1]
    doDantzigWolfeOptim(
        vec_c, mat_a, vec_b, vecVecIndexSub, numSub, m_mat_h, vecRowMatD, vecVecIndexBinInSub,
        wheBranch, vecWhiSub, vecWhiVar, vecWhiBranch
        )
end

main()


text:
import NMRSpectraSimulator

import Unicode
import JSON
import JSON3
import Destruct

output_folder = "/home/roy/Documents/repo/NMRData/input/fit_configs"
isdir(output_folder) || mkpath(output_folder)
save_path = joinpath(output_folder, "calibrate_700MHz_type1_select_compounds_default.json")

# get mapping from molecule names to their spin system info json files.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"
dict_compound_to_filename = JSON.parsefile("/home/roy/Documents/repo/NMRData/input/compound_mapping/select_compounds.json")

# get list of molecule names.
molecule_names = collect( key for (key, value) in dict_compound_to_filename )

# set up machine simulate and get
# machine values taken from the BMRB 700 MHz 20 mM glucose experiment.
fs = 14005.602240896402
SW = 20.0041938620844
ν_0ppm = 10656.011933076665

println("Timing: setupmixtureproxies()")
@time mixture_params = NMRSpectraSimulator.setupmixtureSH(molecule_names,
    H_params_path, dict_compound_to_filename, fs, SW,
    ν_0ppm)
As = mixture_params

println("Timing: mag equivalence")
@time MEs = NMRSpectraSimulator.getmageqinfomixture(molecule_names,
    H_params_path,
    dict_compound_to_filename;
    unique_cs_atol = 1e-6)

# corresponding to NMRSpectraSimulator.SpinSysFIDType1{T}
function populatedefaultvaluescompound1(N_spin_systems::Int, N_singlets::Int, fill_val::T) where T
    #
    systems = collect( fill_val for i = 1:N_spin_systems)
    singlets = collect( fill_val for i = 1:N_singlets)

    return [systems; singlets]
end

"""
Creates ordering of length N from equivalent indices.
`eq_inds` can only contain unique values from 1:N.
When used for a spin system, N is the number of nuclei in the system.

Returns the indices for the N nuclei, and the number of degrees of freedom.

Example:
import Destruct
tmp = collect( createorderingfromeqinds(ME[i], A.N_spins_sys[i]) for i = 1:length(A.N_spins_sys) )
κs_β_ordering, κs_β_DOF = Destruct.destruct(tmp)
"""
function createorderingfromeqinds(eq_inds::Vector{Vector{Int}}, N::Int)

    j = 0 # degrees of freedom counter.
    out = zeros(Int, N)

    if isempty(eq_inds)
        return collect(1:N), N
    end

    # check if `eq_inds` only contains unique values from 1:N.
    if all( all( (eq_inds[k][l] in 1:N) for l = 1:length(eq_inds[k])) for k = 1:length(eq_inds) )
        for k = 1:length(eq_inds)
            j += 1

            out[eq_inds[k]] .= j
        end

        # fill the rest.
        for n = 1:N
            if out[n] == 0
                # out[n] is unassigned.
                j += 1
                out[n] = j
            end
        end

        return out, j
    end

    println("Invalid ME, using default.")
    return collect(1:N), N
end

"""
based on the provided default values, and β ordering based on magnetic equivalence of nuclei.
"""
function createcalibrateconfigdict(As, MEs, molecule_names;
    max_cs_shift_default = 0.05,
    λ_lb_default = 0.7,
    λ_ub_default = 1.5)

    @assert length(As) == length(molecule_names) == length(MEs)

    dict_out = Dict()
    for n = 1:length(As)

        ## parse.
        A = As[n]
        name = molecule_names[n]
        ME = MEs[n]

        ## prepare.
        key_out = "$(name)"

        # spin systems and singlets.
        max_cs_shift = populatedefaultvaluescompound1(length(A.N_spins_sys), length(A.αs_singlets), max_cs_shift_default)
        λ_lb = populatedefaultvaluescompound1(length(A.N_spins_sys), length(A.αs_singlets), λ_lb_default)
        λ_ub = populatedefaultvaluescompound1(length(A.N_spins_sys), length(A.αs_singlets), λ_ub_default)

        # β for spin systems only. Local (to each spin system) indices.
        tmp = collect( createorderingfromeqinds(ME[i], A.N_spins_sys[i]) for i = 1:length(A.N_spins_sys) )
        κs_β_ordering, κs_β_DOF = Destruct.destruct(tmp)

        ## store.
        dict_out[key_out] = Dict( "λ_lb" => λ_lb,
        "λ_ub" => λ_ub,
        "maximum chemical shift" => max_cs_shift,
        "κs_β ordering" => κs_β_ordering,
        "κs_β degrees of freedom" => κs_β_DOF)
    end

    return dict_out
end

dict_out = createcalibrateconfigdict(As, MEs, molecule_names;
    max_cs_shift_default = 0.05,
    λ_lb_default = 0.7,
    λ_ub_default = 1.5)

stringdata = JSON.json(dict_out)

open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end


text:
module SyslogLoggingCompose

using SyslogLogging, LogCompose
import LogCompose: logcompose, log_min_level

const sysloglck = ReentrantLock()
function logcompose(::Type{SyslogLogging.SyslogLogger}, config::Dict{String,Any}, logger_config::Dict{String,Any})
    ident = logger_config["identity"]
    level = log_min_level(logger_config, "Info")

    kwargs = Dict{Symbol,Any}()

    kwargs[:facility] = Symbol(get(logger_config, "facility", "user"))

    if get(logger_config, "lock", false)
        kwargs[:lck] = sysloglck
    end

    server_type = get(logger_config, "server_type", "local")
    if (server_type == "tcp") || (server_type == "udp")
        kwargs[:tcp] = (server_type == "tcp")
        if haskey(logger_config, "server_host")
            kwargs[:host] = logger_config["server_host"]
        end
        if haskey(logger_config, "server_port")
            kwargs[:port] = logger_config["server_port"]
        end
    end
    SyslogLogging.SyslogLogger(ident, level; kwargs...)
end

end # module


text:
# # Executable Cells
#
# This file defines a custom CommonMark node type that provides executable code
# cells.

"""
A CommonMark rule used to define the "cell" parser. A `CellRule` holds a
`.cache` of the `Module`s that have been defined in a markdown document so that
cells can depend on definitions and values from previous cells.
"""
struct CellRule
    cache::Dict{String,Module}
    imports::Vector{Module}

    function CellRule(; cache = Dict{String,Module}(), imports = Module[])
        return new(cache, vcat(imports, Objects))
    end
end

struct Embedded <: CommonMark.AbstractBlock end

CommonMark.is_container(::Embedded) = true

CommonMark.write_html(::Embedded, w, n, ent) = nothing
CommonMark.write_latex(::Embedded, w, n, ent) = nothing
CommonMark.write_term(::Embedded, w, n, ent) = nothing
CommonMark.write_markdown(::Embedded, w, n, ent) = nothing

"""
    struct Cell

A custom node type for CommonMark.jl that holds an executable "cell" of code.
"""
struct Cell <: CommonMark.AbstractBlock
    node::CommonMark.Node
    value::Any
    output::String
end

CommonMark.block_modifier(c::CellRule) = CommonMark.Rule(100) do parser, node
    if isjuliacode(node) && iscell(node.meta)
        # Load the module for the current cell and evaluate the contents.
        sandbox = getmodule!(c, node)
        captured = IOCapture.capture(rethrow=InterruptException) do
            include_string(sandbox, node.literal)
        end
        # When the value is displayable as markdown then we reparse that
        # representation and include the resulting AST in it's place.
        # Otherwise we just capture it's value and output for display later as
        # a normal cell.
        if  get(node.meta, "markdown", "true") == "true" && showable(MIME("text/markdown"), captured.value)
            text = Base.invokelatest(() -> sprint(show, MIME("text/markdown"), captured.value))
            subparser = init_markdown_parser()
            ast = subparser(text)
            ast.t = Embedded()
            CommonMark.insert_after(node, ast)
            CommonMark.unlink(node)
        else
            cell = CommonMark.Node(Cell(node, captured.value, captured.output))
            CommonMark.insert_after(node, cell)
            if get(node.meta, "display", "true") == "false"
                cell.meta = node.meta
                CommonMark.unlink(node)
            end
        end
    end
    return nothing
end

struct EmbeddedInline <: CommonMark.AbstractInline end

CommonMark.is_container(::EmbeddedInline) = true

CommonMark.write_html(::EmbeddedInline, w, n, ent) = nothing
CommonMark.write_latex(::EmbeddedInline, w, n, ent) = nothing
CommonMark.write_term(::EmbeddedInline, w, n, ent) = nothing
CommonMark.write_markdown(::EmbeddedInline, w, n, ent) = nothing

CommonMark.inline_modifier(c::CellRule) = CommonMark.Rule(100) do parser, block
    for (node, ent) in block
        if ent && is_inline_code(node) && iscell(node.meta)
            sandbox = getmodule!(c, node)
            captured = IOCapture.capture(rethrow=InterruptException) do
                include_string(sandbox, node.literal)
            end
            if showable(MIME("text/markdown"), captured.value)
                text = Base.invokelatest(() -> sprint(show, MIME("text/markdown"), captured.value))
                subparser = init_markdown_parser()
                ast = subparser(text).first_child
                ast.t = EmbeddedInline()
                CommonMark.insert_after(node, ast)
                CommonMark.unlink(node)
            else
                node.literal = Base.invokelatest(() -> sprint(show, MIME("text/plain"), captured.value))
            end
        end
    end
    return nothing
end

function getmodule!(rule::CellRule, node::CommonMark.Node)
    id = get!(string ∘ gensym, node.meta, "cell")
    return get!(rule.cache, id) do
        sandbox = Module() # TODO: named.
        for each in rule.imports
            name = gensym()
            Core.eval(sandbox, :($name=$each; using .$name))
        end
        return sandbox
    end
end

isjuliacode(n::CommonMark.Node) = n.t isa CommonMark.CodeBlock && n.t.info == "julia"
is_inline_code(n::CommonMark.Node) = n.t isa CommonMark.Code
iscell(d::AbstractDict) = haskey(d, "cell") || get(d, "element", "") == "cell"

# ## Cell Evaluator

"""
    display_as(default, cell, writer, [mimes...])

Given a `cell` this function evaluates it and prints the output to `writer`
using the first available `MIME` from `mimes`. Uses the `default` printer
function to print any code blocks that are required in the output.
"""
function display_as(default, cell, w, mimes)
    ## Display options for cell:
    show_output = get(cell.node.meta, "output", "true")
    show_result = get(cell.node.meta, "result", "true")
    ## Evaluate the cell contents in a sandboxed module, possibly reusing one
    ## from an earlier cell if the names match.
    if !isempty(cell.output) && show_output == "true"
        ## There's been some output to the stream, put that in
        ## a verbatim block before the real output so long as
        ## `output=false` was not set for the cell.
        out = CommonMark.Node(CommonMark.CodeBlock())
        out.meta["class"] = ["plaintext", "cell-output", "cell-stream"]
        out.literal = cell.output
        default(out.t, w, out, true)
    end
    show_result == "true" || return nothing # Display result unless `result=false` was set.
    cell.value === nothing && return nothing # Skip `nothing` results.
    for mime in mimes
        if showable(mime, cell.value)
            ## We've found a suitable mimetype, display as that.
            limitedshow(w.buffer, default, mime, cell.value)
            return nothing
        end
    end
    ## Default output displays the result as in the REPL.
    code = CommonMark.Node(CommonMark.CodeBlock())
    code.t.info = "plaintext"
    code.meta["class"] = ["plaintext", "cell-output", "cell-result"]
    code.literal = limitedshow(default, MIME("text/plain"), cell.value)
    default(code.t, w, code, true)
    return nothing
end

"""
    limitedshow([io], mime, result)

Prints out a "limited" representation of `result` in the given `mime` to the
provided `io` stream, or returns a `String` of the output when no `io` is
given.
"""
function limitedshow end

limitedshow(io::IO, default, m, r) = Base.invokelatest(show, IOContext(io, :limit=>true), m, r)
limitedshow(default, m, r) = sprint(limitedshow, default, m, r)

# ## Supported image MIMES.

const SUPPORTED_MIMES = Dict{Symbol,Vector{MIME}}(
    :html  => map(MIME, [
        "image/svg+xml", # TODO: optimal ordering.
        "image/png",
        "image/jpeg",
        "image/gif",
        "text/html",
        "text/latex",
    ]),
    :latex => map(MIME, [
        "text/tikz", # TODO: optimal ordering.
        "image/png",
        "application/pdf",
        "text/latex",
    ]),
    :term  => MIME[],
)

const IMAGE_MIMES = Union{
    MIME"application/pdf",
    MIME"image/gif",
    MIME"image/jpeg",
    MIME"image/png",
    MIME"image/svg+xml",
    MIME"text/tikz",
}

function limitedshow(io::IO, fn, mime::IMAGE_MIMES, result)
    name = string(hash(result), _ext(mime))
    open(name, "w") do handle
        Base.invokelatest(show, handle, mime, result)
    end
    node = CommonMark.Node(CommonMark.Image())
    node.t.destination = _inline_image(fn, name)
    return cm_wrapper(fn)(io, node)
end

_inline_image(::typeof(CommonMark.write_html), name::AbstractString) = _base64resource(name)
_inline_image(::Any, name) = name

_ext(::MIME"application/pdf") = ".pdf"
_ext(::MIME"image/gif") = ".gif"
_ext(::MIME"image/jpeg") = ".jpeg"
_ext(::MIME"image/png") = ".png"
_ext(::MIME"image/svg+xml") = ".svg"
_ext(::MIME"text/tikz") = ".tikz"

# ## CommonMark Writers
#
# These definitions are needed by CommonMark to hook into it's display system.

function CommonMark.write_html(cell::Cell, w, n, ent)
    ent && display_as(CommonMark.write_html, cell, w, SUPPORTED_MIMES[:html])
    return nothing
end
cm_wrapper(::typeof(CommonMark.write_html)) = CommonMark.html # The wrapper function for write_html

function CommonMark.write_latex(cell::Cell, w, n, ent)
    ent && display_as(CommonMark.write_latex, cell, w, SUPPORTED_MIMES[:latex])
    return nothing
end
cm_wrapper(::typeof(CommonMark.write_latex)) = CommonMark.latex # The wrapper function for write_latex

# The following two definitions aren't really needed since Publish doesn't support
# output to terminal or markdown, but are defined to ensure the display system is
# complete for the [`Cell`](#) node type.

function CommonMark.write_term(cell::Cell, w, n, ent)
    if ent
        display_as(CommonMark.write_term, cell, w, SUPPORTED_MIMES[:term])
        ## Make sure to add a linebreak afterwards if needed.
        if !CommonMark.isnull(n.nxt)
            CommonMark.print_margin(w)
            CommonMark.print_literal(w, "\n")
        end
    end
    return nothing
end

## Markdown roundtrips, so shouldn't display cells.
CommonMark.write_markdown(cell::Cell, w, n, ent) = nothing

#
# Custom Tables and Figures
#

module Objects

export Figure, Table

# TODO: implement a table version of this as well.
Base.@kwdef struct Figure{T}
    object::T
    placement::Symbol = :h
    alignment::Symbol = :center
    maxwidth::String = "\\linewidth"
    landscape::Bool = false
    caption::String = ""
    desc::String = ""
end
Figure(object; options...) = Figure{typeof(object)}(; object=object, options...)

Base.@kwdef struct Table{T}
    object::T
    placement::Symbol = :h
    alignment::Symbol = :center
    landscape::Bool = false
    caption::String = ""
    desc::String = ""
    type::Symbol = :tabular
    pretty_table::NamedTuple = NamedTuple()
end
Table(object; options...) = Table{typeof(object)}(; object=object, options...)

end

function _format_caption(m::MIME, content::AbstractString)
    replacer(::MIME"text/latex", s) = replace(rstrip(s), r"\\par$" => "")
    replacer(::MIME, s) = rstrip(s)
    p = init_markdown_parser()
    return replacer(m, sprint(show, m, p(content)))
end

function _readable_unique_filename(desc::String, object)
    hash62 = string(hash(object); base = 62)
    return isempty(desc) ? hash62 : "$(CommonMark.slugify(desc))-$hash62"
end
_readable_unique_filename(f::Objects.Figure) = _readable_unique_filename(f.desc, f.object)

function Base.show(io::IO, m::MIME"text/latex", f::Objects.Figure)
    for mime in SUPPORTED_MIMES[:latex]
        if showable(mime, f.object)
            filename = string(_readable_unique_filename(f), _ext(mime))
            open(filename, "w") do handle
                Base.invokelatest(show, handle, mime, f.object)
            end
            f.landscape && println(io, "\\begin{landscape}")
            println(io, "\\begin{figure}[$(f.placement)]")
            println(io, "\\adjustimage{max height=\\textheight,max width=$(f.maxwidth),$(f.alignment)}{./$filename}")
            if !isempty(f.caption)
                desc = isempty(f.desc) ? "" : "[$(_format_caption(m, f.desc))]"
                println(io, "\\caption$(desc){$(_format_caption(m, f.caption))}")
            end
            println(io, "\\end{figure}")
            f.landscape && println(io, "\\end{landscape}")
            return nothing
        end
    end
    throw(ErrorException("cannot display type $(typeof(f.object)) as a figure."))
end

function Base.show(io::IO, m::MIME"text/latex", f::Objects.Figure{Matrix{T}}) where T
    objects = f.object
    for mime in SUPPORTED_MIMES[:latex]
        if all(object -> showable(mime, object), objects)
            f.landscape && println(io, "\\begin{landscape}")
            println(io, "\\begin{figure}[$(f.placement)]")
            # Write out the objects left to right.
            println(io, "\\begin{tabular}{$(repeat('c', size(objects, 2)))}")
            M, N = size(objects)
            max_height = round(1 / M; digits = 2)
            max_width = round(1 / N; digits = 2)
            for (row_num, row) in enumerate(eachrow(objects))
                for (ith, object) in enumerate(row)
                    filename = string(_readable_unique_filename(f.desc * "[$row_num,$ith]", object), _ext(mime))
                    open(filename, "w") do handle
                        Base.invokelatest(show, handle, mime, object)
                    end
                    print(io, "\\includegraphics[max width=$(max_width)\\textwidth,max height=$(max_height)\\textheight,valign=m]{./$filename}")
                    ith < N ? print(io, " & ") : println(io, "\\\\")
                end
            end
            println(io, "\\end{tabular}")
            if !isempty(f.caption)
                desc = isempty(f.desc) ? "" : "[$(_format_caption(m, f.desc))]"
                println(io, "\\caption$(desc){$(_format_caption(m, f.caption))}")
            end
            println(io, "\\end{figure}")
            f.landscape && println(io, "\\end{landscape}")
            return nothing
        end
    end
    throw(ErrorException("cannot display type $(typeof(f.object)) as a figure."))
end

function Base.show(io::IO, m::MIME"text/html", f::Objects.Figure)
    for mime in SUPPORTED_MIMES[:html]
        if showable(mime, f.object)
            println(io, "
")
            if isa(mime, MIME"text/html")
                # HTML mimes must be embedded directly into output.
                Base.invokelatest(show, io, mime, f.object)
            else
                # Other types get written to file then read back in.
                filename = string(_readable_unique_filename(f), _ext(mime))
                open(filename, "w") do handle
                    Base.invokelatest(show, handle, mime, f.object)
                end
                img = CommonMark.Node(CommonMark.Image())
                img.t.destination = _base64resource(filename)
                img.t.title = f.caption
                CommonMark.html(io, img)
                if !isempty(f.caption)
                    cap = _format_caption(m, "Figure: $(f.caption)")
                    println(io, "
$(cap)
")
                end
            end
            println(io, "
")
            return nothing
        end
    end
    throw(ErrorException("cannot display type $(typeof(f.object)) as a figure."))
end

const _LATEX_HORIZONTAL_ALIGNMENT_MAPPING = Dict(
    :center => "\\centering",
    :left => "\\raggedleft",
    :right => "\\raggedright",
)

function Base.show(io::IO, m::MIME"text/latex", t::Objects.Table)
    # We only wrap tabular environments, since longtable does all this for us.
    # It is nessecary to pass along the options to pretty_table for longtables
    # since it handles captions and the like internally.
    t.landscape && println(io, "\\begin{landscape}")
    if t.type == :tabular
        println(io, "\\begin{table}[$(t.placement)]")
        println(io, get(_LATEX_HORIZONTAL_ALIGNMENT_MAPPING, t.alignment, "\\centering"))
        if !isempty(t.caption)
            desc = isempty(t.desc) ? "" : "[$(_format_caption(m, t.desc))]"
            println(io, "\\caption$(desc){$(_format_caption(m, t.caption))}")
        end
        # Manually unwrap table until upstream deps are sorted.
        temp_io = IOBuffer()
        PrettyTables.pretty_table(
            temp_io, t.object;
            backend=Val(:latex),
            # wrap_table=false,
            table_type=:tabular,
            t.pretty_table... # pass through any extra options by user.
        )
        str = String(take!(temp_io))
        join(io, split(str, "\n")[2:end-2], "\n")
        println(io, "\\end{table}")
    else
        PrettyTables.pretty_table(
            io, t.object;
            backend=Val(:latex),
            # wrap_table=false,
            table_type=:longtable,
            title=_format_caption(m, t.caption),
            t.pretty_table..., # pass through any extra options by user.
        )
    end
    t.landscape && println(io, "\\end{landscape}")
    return nothing
end

function Base.show(io::IO, m::MIME"text/html", t::Objects.Table)
    println(io, "
")
    cap = _format_caption(m, "Table: $(t.caption)")
    println(io, "
$(cap)
")
    PrettyTables.pretty_table(
        io, t.object;
        backend=Val(:html),
    )
    println(io, "
")
end


text:
using FractionalDiffEq, Plots

function sys!(du, u, p, t)
    du[1] = -0.05*u[2] - 0.05*u[3] + 0.01*tanh(u[2])
    du[2] = 0.05*u[1] + 0.02*u[2] + 0.01*tanh(u[1])
    du[3] = 0.1 - 0.2*u[3] + 0.05*u[1]*u[3] + 0.01*tanh(u[3])
end
prob = FractionalDifferenceSystem(sys!, 0.98, [1, -1, 0])
result = solve(prob, 7, GL())

plot(result[1, :], result[2, :], result[3, :], seriestype=:scatter)


text:
module SearchSortedNearest

export searchsortednearest

# from: https://discourse.julialang.org/t/findnearest-function/4143/5
"""
    searchsortednearest(a, x; by=<transform>, lt=<comparison>, rev=false)

Find the index of (sorted) collection `a` that is nearest (smallest absolute difference) to `x`.
Ties go to the smallest index.

# Examples

    using SearchSortedNearest

    searchsortednearest(1:10, 1.1) == 1
    searchsortednearest(1:10, 1.9) == 2
"""
function searchsortednearest(a, x; by=identity, lt=isless, rev=false)
    idx = searchsortedfirst(a, x; by, lt, rev)
    idx == 1 && return idx
    idx > length(a) && return length(a)
    a[idx] == x && return idx
    idx = lt(abs(by(a[idx]) - by(x)), abs(by(a[idx - 1]) - by(x))) ? idx : idx - 1
    return rev ? length(a) - idx + 1 : idx
end

end


text: abstract type DDMArray{T,N} <: AbstractArray{T,N} end


text:
#=

# Error handling

Error handling is achieved via [ResultTypes.jl](https://github.com/iamed2/ResultTypes.jl) to avoid the large overhead introduced with `try`/`catch` blocks. All functions that need to perform an operation that returns a `VkResult` are wrapped into a `ResultTypes.Result` (distinct from `Vulkan.Result` which is the wrapped version of `VkResult`), which contains a [`VulkanError`](@ref) if a non-success code is encountered. Custom error handling can be performed using the following pattern

=#

using Vulkan

res = create_instance([], [])
if iserror(res) # handle the error
    err = unwrap_error(res)
    if err.code == ERROR_INCOMPATIBLE_DRIVER
        error(
            """
            No driver compatible with the requested API version could be found.
            Please make sure that a driver supporting Vulkan is installed, and
            that it is up to date with the requested version.
            """
        )
    elseif err.code == ERROR_LAYER_NOT_PRESENT
        error("One or more layers could not be found")
    else
        throw(err)
    end
else # get the instance
    unwrap(res)
end

# Note that calling `unwrap` directly on the result will throw any contained `VulkanError` if there is one. So, if you just want to throw an exception when encountering an error, you can do

unwrap(create_instance([], []))

# Because it may be tedious to unwrap everything by hand and explicitly set the create info structures, [convenience constructors](@ref expose-create-info-args) are defined for handle types so that you can just do

Instance([], [])

#=

However, note that exceptions are thrown whenever the result is an error with this shorter approach.

Furthermore, all functions that may return non-success (but non-error) codes return a `Vulkan.Result` type along with any other returned value, since the return code may still be of interest.

For more details on the `ResultTypes.Result` type and how to handle it, please consult the [ResultTypes.jl documentation](https://iamed2.github.io/ResultTypes.jl/stable/).

=#


text:
# VARIANT GETTERS_SETTERS

function com_variantset(var::Ref{Variant}, data::T) where {T}
    unsafe_store!(Ptr{T}(pointer_from_objref(var) + 8), data)
    return var
end

function com_variantset(pvar::Ptr{Variant}, data::T) where {T}
    unsafe_store!(Ptr{T}(pvar + 8), data)
end

function com_variantget(var::Ref{Variant}, ::Type{T}) where {T}
    return unsafe_load(Ptr{T}(pointer_from_objref(var) + 8))
end

function com_variantget(pvar::Ptr{Variant}, ::Type{T}) where {T}
    return unsafe_load(Ptr{T}(pvar + 8))
end

# VARIANT CONSTRUCTORS

function variant_free(var::Variant)
    if !(gettype(var.vt) <: Ptr)
        ccall((:VariantClear, "oleaut32.dll"), Cvoid, (Ref{Variant},), var)
    end
end

# This is not a copy constructor, just plumbing for for some algorithms.
# Variant structure is immutable therefore it can be copied just as v1 = v2.
#
Variant(pass::Variant) = pass

# Convert julia string to Variant with newly allocated BSTR.
#
function Variant(value::String)
    str = push!(transcode(UInt16, value), '\0')
    var = Variant(ccall((:SysAllocString, "oleaut32.dll"), Ptr{UInt16}, (Ptr{UInt16},), str))
    return var
end

# Allocate SafeArray of variants to be used wtih COM methods.
#
# vars : Array of variants.
# return: Variant with allocated SafeArray.
#
function com_createarray(vars::Array{Variant}) ::Variant

    abnds = SafeArrayBound.(collect(size(vars)), 0)
    array = ccall((:SafeArrayCreate, "oleaut32.dll"), Ptr{SafeArray}, (Cushort, Cuint, Ptr{SafeArrayBound}),
                VT_VARIANT, ndims(vars), pointer(abnds))
    if array != C_NULL
        pdata = Ref(Ptr{Variant}(0))
        ccall((:SafeArrayAccessData, "oleaut32.dll"), Cint, (Ptr{SafeArray}, Ref{Ptr{Variant}}), array, pdata)
        for n in 0:length(vars) - 1
            unsafe_store!(pdata[] + sizeof(Variant)*n, vars[n+1])
        end
        ccall((:SafeArrayUnaccessData, "oleaut32.dll"), Cint, (Ptr{SafeArray},), array)
    end
    return com_variantset(Ref(Variant(VT_VARIANT | VT_ARRAY, 0, 0, 0, 0)), array)[]
end

function Variant(value::Array{T}) where {T}
    return com_createarray(Variant.(value))
end

function Variant(value::T) where {T <: Union{Ptr, Number}}
    var = Ref(Variant(getindex(T), 0, 0, 0, 0))
    com_variantset(var, value)
    return var[]
end

# CONVERTERS

# Return the typed value of Variant, with special handling of arrays.
#
function value(var::Variant)
    if (var.vt & VT_ARRAY) > 0
        ret = convert(Array, var)
        if ret[1] === nothing
            popfirst!(ret)
            if length(ret) == 1
                return ret[1]
            end
        end
        return ret
    else
        return convert(get(indextype, var.vt & VT_TYPEMASK, Nothing), var)
    end
end

function Base.convert(::Type{Nothing}, var::Variant)
    return nothing
end

function Base.convert(::Type{T}, var::Variant) where {T <: Union{Ptr, Number}}
    value = com_variantget(Ref(var), T)
    if var.vt == 8
        return convert(String, value)
    else
        return value
    end
end

Base.String(var::Variant) = convert(String, var)
#
function Base.convert(::Type{String}, var::Variant)
    return String(convert(Ptr{UInt16}, var))
end

Base.String(str::Ptr{UInt16}) = convert(String, str)
#
function Base.convert(::Type{String}, str::Ptr{UInt16})
    buf = Vector{UInt16}()
    for n in 1:typemax(Int16)
        chr = unsafe_load(str, n)
        if chr == zero(UInt16)
            break
        end
        push!(buf, chr)
    end
    return transcode(String, buf)
end

function Base.convert(::Type{Array}, var::Variant)
    if (var.vt & VT_ARRAY) > 0
        if (var.vt & VT_BYREF) > 0
            apntr = unsafe_load(com_variantget(Ref(var), Ptr{Ptr{SafeArray}}))
        else
            apntr = com_variantget(Ref(var), Ptr{SafeArray})
        end
        array = unsafe_load(apntr)
        abnds = getfield.(unsafe_load.(Ptr{SafeArrayBound}(apntr + 24), 1:array.cDims), :cElements)
        if length(abnds) > 1
            abnds[1], abnds[2] = abnds[2], abnds[1]
        end
        vars = unsafe_load.(Ptr{Variant}(array.pvData), 1:prod(abnds))
        return reshape(convert.(gettype.(getfield.(vars, :vt)), vars), Tuple(abnds))
    end
    return nothing
end

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
using Documenter

makedocs(
    sitename="Peacock.jl Documentation",
    pages = [
        "index.md",
        "Tutorials" => [
            "tutorials/getting_started.md",
        ],
        "How-to guides" => [
            "how-tos/zoo.md",
            "how-tos/wilson_loops.md",
            "how-tos/gpu.md",
        ],
        "reference/index.md",
        "contributing.md",
    ]
)

deploydocs(
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
    repo="github.com/sp94/Peacock.jl.git",
)


text:
"""
    read_particles_in_geometry( filename::String, blocks::Vector{String},
                                geometry::AbstractGadgetGeometry;
                                parttype::Integer=0, verbose::Bool=true,
                                use_keys::Bool=true, do_shift_across_box_border::Bool=true)

Reads all particles within a space defined by an `AbstractGeometry` struct for a given particle type.
Returns a dictionary with all requested blocks.
If `shift_across_box_border` is `true`, the particles are moved beyond the borders of a periodic box,
if `false`, the periodicity of the box is still accounted for, but the particles' positions are not
shifted.
"""
function read_particles_in_geometry(filename::String, blocks::Vector{String},
                                    geometry::AbstractGadgetGeometry;
                                    parttype::Integer=0, verbose::Bool=true,
                                    use_keys::Bool=true, do_shift_across_box_border::Bool=true)

    h = read_header(filename * ".0")

    # make sure position is read - needed for determining if in geometry or not
    if "POS" ∉ blocks
        @info "POS not found in blocks, adding POS..."
        append!(blocks, "POS")
    end

    # get the corners of a box which contains the geometry
    corner_lowerleft, corner_upperright = get_geometry_box_corners(geometry)

    # read the box
    d = read_particles_in_box(filename, blocks, corner_lowerleft, corner_upperright; parttype, verbose, use_keys)

    # determine particles within geometry, shifting the particles across periodic boundaries
    r₀ = get_geometry_center(geometry)
    if do_shift_across_box_border
        d["POS"] .= shift_across_box_border.(d["POS"], r₀, h.boxsize, 1 // 2 * h.boxsize)
        mask = get_geometry_mask(geometry, d["POS"])
    else
        pos = shift_across_box_border.(d["POS"], r₀, h.boxsize, 1 // 2 * h.boxsize)
        mask = get_geometry_mask(geometry, pos)
    end

    # iterate blocks of dict to remove filtered out particles
    for block ∈ keys(d)
        # get all dimensions but the last (the last is masked)
        colons = repeat([:], ndims(d[block]) - 1)

        d[block] = d[block][colons..., mask]
    end

    return d
end

"""
    read_particles_in_geometry( filename::String, block::String,
                                geometry::AbstractGadgetGeometry;
                                parttype::Integer=0, verbose::Bool=true,
                                use_keys::Bool=true, do_shift_across_box_border::Bool=true)

Reads all particles within a space defined by an `AbstractGeometry` struct for a given particle type.
Returns a dictionary with the requested block.
"""
function read_particles_in_geometry(filename::String, block::String,
                                    geometry::AbstractGadgetGeometry;
                                    parttype::Integer=0, verbose::Bool=true,
                                    use_keys::Bool=true, do_shift_across_box_border::Bool=true)

    return read_particles_in_geometry(filename, [block], geometry; parttype, verbose, use_keys)
end


text:
export DataGenerator, UniformUnion, Sawtooth, BayesianConvNP, Mixture

"""
    DataGenerator

# Fields
- `process`: Something that can be called at inputs `x` and a noise level `noise` and gives
    back a distribution that can be fed to `randn` to sample values corresponding to those
    inputs `x` at observation noise `noise`.
- `batch_size::Integer=16`: Number of tasks in a batch.
- `x_context::Distribution=Uniform(-2, 2)`: Distribution to sample context inputs from.
- `x_target::Distribution=Uniform(-2, 2)`: Distribution to sample target inputs from.
- `num_context::Distribution=DiscreteUniform(10, 10)`: Distribution of number of context
    points in a task.
- `num_target::Distribution=DiscreteUniform(100, 100)`: Distribution of number of target
    points in a task.
- `σ²::Float64`: Noise variance to add to the data.
"""
struct DataGenerator
    process
    batch_size::Integer
    x_context::Distribution
    x_target::Distribution
    num_context::Distribution
    num_target::Distribution
    σ²::Float64

    function DataGenerator(
        process;
        batch_size::Integer=16,
        x_context::Distribution=Uniform(-2, 2),
        x_target::Distribution=Uniform(-2, 2),
        num_context::Distribution=DiscreteUniform(10, 10),
        num_target::Distribution=DiscreteUniform(100, 100),
        σ²::Float64=1e-8
    )
        return new(
            process,
            batch_size,
            x_context,
            x_target,
            num_context,
            num_target,
            σ²
        )
    end
end

"""
    struct UniformUnion{T<:Real} <: ContinuousUnivariateDistribution

A union of various `Uniform`s.

# Fields
- `uniforms::Vector{Uniform{T}}`: Underlying `Uniforms`s.
- `probs::ProbabilityWeights`: Probabilities associated to every `Uniform`.
"""
struct UniformUnion{T<:Real} <: ContinuousUnivariateDistribution
    uniforms::Vector{Uniform{T}}
    probs::ProbabilityWeights
end

"""
    UniformUnion(uniforms::Uniform...)

Construct a `UnionUniform` and weight the uniforms according to the size of their domain.

# Arguments
- `uniforms::Uniform...`: Underlying `Uniform`s.

# Returns
- `UniformUnion`: Associated `UniformUnion`.
"""
function UniformUnion(uniforms::Uniform...)
    lengths = [maximum(d) - minimum(d) for d in uniforms]
    return UniformUnion(collect(uniforms), pweights([x / sum(lengths) for x in lengths]))
end

Base.rand(u::UniformUnion) = rand(StatsBase.sample(u.uniforms, u.probs))
Base.rand(u::UniformUnion, n::Int64) = [rand(u) for _ = 1:n]

"""
    (generator::DataGenerator)(num_batches::Integer)

# Arguments
- `num_batches::Integer`: Number of batches to sample.

# Returns
- `Vector`: Vector of tasks.
"""
function (generator::DataGenerator)(num_batches::Integer)
    return [_make_batch(
        generator,
        rand(generator.num_context),
        rand(generator.num_target)
    ) for i in 1:num_batches]
end

_float32(x) = Float32.(x)

function _make_batch(generator::DataGenerator, num_context::Integer, num_target::Integer)
    # Sample tasks.
    tasks = []
    for i in 1:generator.batch_size
        # Determine context and target set.
        xc = rand(generator.x_context, num_context)
        xt = rand(generator.x_target, num_target)

        # Concatenate inputs and sample.
        x = vcat(xc, xt)
        y = rand(generator.process(x, generator.σ²))

        push!(tasks, _float32.((
            xc,
            y[1:num_context],
            xt,
            y[num_context + 1:end]
        )))
    end

    # Collect as a batch and return.
    return map(x -> cat(x...; dims=3), zip(tasks...))
end

"""
    FDD{T}

Finite-dimensional distribution of a process at particular inputs.

# Fields
- `x`: Inputs.
- `σ²`: Noise variance.
- `process::T`: Underlying process.
"""
struct FDD{T}
    x
    σ²
    process::T
end

"""
    Sawtooth

Random truncated Fourier expansion of a sawtooth wave.

# Fields
- `freq_dist::Distribution=Uniform(3, 5)`: Distribution of the frequency.
- `shift_dist::Distribution=Uniform(-5, 5)`: Distribution of the shift.
- `trunc_dist::Distribution=DiscreteUniform(10, 20)`: Distribution of the truncation of
    the Fourier expansion.
"""
struct Sawtooth
    freq_dist::Distribution
    shift_dist::Distribution
    trunc_dist::Distribution

    function Sawtooth(
        freq_dist::Distribution=Uniform(3, 5),
        shift_dist::Distribution=Uniform(-5, 5),
        trunc_dist::Distribution=DiscreteUniform(10, 20)
    )
        return new(freq_dist, shift_dist, trunc_dist)
    end
end

(s::Sawtooth)(x, noise) = FDD(x, noise, s)

function Base.rand(s::FDD{Sawtooth})
    # Sample parameters for particular sawtooth wave.
    amp = 1
    freq = rand(s.process.freq_dist)
    shift = rand(s.process.shift_dist)
    trunc = rand(s.process.trunc_dist)

    # Apply shift.
    x = s.x .- shift

    # Construct expansion.
    k = collect(range(1, trunc + 1, step=1))'
    f = 0.5amp .- (amp / pi) .* sum((-1).^k .* sin.(2pi .* k .* freq .* x) ./ k, dims=2)

    # Add noise and return.
    return f .+ sqrt(s.σ²) .* randn(eltype(x), size(x)...)
end

"""
    Bayesian ConvNP

Bayesian ConvNP.

# Fields
- `receptive_field::Float32`: Width of the receptive field.
- `num_layers::Integer`: Number of layers of the CNN, excluding an initial
    and final pointwise convolutional layer to change the number of channels
    appropriately.
- `num_channels::Integer`: Number of channels of the CNN.
- `points_per_unit::Float32=30f0`: Points per unit for the discretisation. See
     `UniformDiscretisation1D`.
"""
struct BayesianConvNP
    receptive_field::Float32
    num_layers::Integer
    num_channels::Integer
    points_per_unit::Float32

    function BayesianConvNP(;
        receptive_field::Float32=1f0,
        num_layers::Integer=4,
        num_channels::Integer=6,
        points_per_unit::Float32=30f0
    )
        new(receptive_field, num_layers, num_channels, points_per_unit)
    end
end

(convnp::BayesianConvNP)(x, noise) = FDD(x, noise, convnp)

function Base.rand(convnp::FDD{BayesianConvNP})
    # Contruct discretisation.
    disc = UniformDiscretisation1D(
        convnp.process.points_per_unit,
        convnp.process.receptive_field / 2,
        1
    )
    xz = reshape(disc(convnp.x), :, 1, 1)

    # Construct CNN with random initialisation.
    conv = build_conv(
        convnp.process.receptive_field,
        convnp.process.num_layers,
        convnp.process.num_channels;
        points_per_unit=convnp.process.points_per_unit,
        num_in_channels=1,
        num_out_channels=1,
        dimensionality=1,
        init_conv=_init_conv_random_bias,
        init_depthwiseconv=_init_depthwiseconv_random_bias
    )

    # Construct decoder.
    scale = 2 / convnp.process.points_per_unit
    decoder = set_conv(1, scale)

    # Draw random encoding.
    encoding = randn(Float32, length(xz))

    # Pass through CNN, which takes in images of height one.
    encoding = reshape(encoding, length(encoding), 1, 1, 1)
    latent = conv(encoding)

    # Perform decoding.
    xt = reshape(convnp.x, length(convnp.x), 1, 1)
    sample = code(decoder, xz, latent, xt)[2][:, 1, 1]  # Also returns inputs.

    # Normalise sample.
    sample = (sample .- mean(sample)) ./ std(sample)

    # Return with noise.
    return sample .+ sqrt(convnp.σ²) .* randn(Float32, size(sample)...)
end

"""
    struct Mixture

Mixture of processes.

# Fields
- `processes`: Processes in the mixture.
"""
struct Mixture
    processes
end

Mixture(processes...) = Mixture(processes)

(mixture::Mixture)(x, noise) = FDD(x, noise, mixture)

Base.rand(fmixture::FDD{Mixture}) =
    rand(rand(fmixture.process.processes)(fmixture.x, fmixture.σ²))


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Pattern <: TechnicalParams
        name::String
        multipliers::Union{Nothing, Vector{Float64}}
        internal::InfrastructureSystemsInternal
    end

Forecast pattern.

# Arguments
- `name::String`
- `multipliers::Union{Nothing, Vector{Float64}}`: Multiplier values of the pattern.
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct Pattern <: TechnicalParams
    name::String
    "Multiplier values of the pattern."
    multipliers::Union{Nothing, Vector{Float64}}
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Pattern(name, multipliers, )
    Pattern(name, multipliers, InfrastructureSystemsInternal(), )
end

function Pattern(; name, multipliers, )
    Pattern(name, multipliers, )
end

# Constructor for demo purposes; non-functional.
function Pattern(::Nothing)
    Pattern(;
        name="init",
        multipliers=nothing,
    )
end

"""Get Pattern name."""
get_name(value::Pattern) = value.name
"""Get Pattern multipliers."""
get_multipliers(value::Pattern) = value.multipliers
"""Get Pattern internal."""
get_internal(value::Pattern) = value.internal


text:
import Base.∈

export constrained_dimensions,
       tosimplehrep,
       remove_redundant_constraints,
       remove_redundant_constraints!,
       linear_map,
       chebyshev_center,
       an_element

# default LP solver for floating-point numbers
function default_lp_solver(N::Type{<:AbstractFloat})
    GLPKSolverLP(method=:Simplex)
end

# default LP solver for rational numbers
function default_lp_solver(N::Type{<:Rational})
    GLPKSolverLP(method=:Exact)
end

# fallback method
function default_polyhedra_backend(P, N)
    require(:Polyhedra; fun_name="default_polyhedra_backend")
    error("no default backend for numeric type $N")
end

"""
    ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N}) where {N<:Real}

Check whether a given point is contained in a polyhedron.

### Input

- `x` -- point/vector
- `P` -- polyhedron

### Output

`true` iff ``x ∈ P``.

### Algorithm

This implementation checks if the point lies inside each defining half-space.
"""
function ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    @assert length(x) == dim(P) "a $(length(x))-dimensional point cannot be " *
        "an element of a $(dim(P))-dimensional set"

    for c in constraints_list(P)
        if dot(c.a, x) > c.b
            return false
        end
    end
    return true
end

"""
    isuniversal(P::AbstractPolyhedron{N}, [witness]::Bool=false
               ) where {N<:Real}

Check whether a polyhedron is universal.

### Input

- `P`       -- polyhedron
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``P`` is universal
* If `witness` option is activated:
  * `(true, [])` iff ``P`` is universal
  * `(false, v)` iff ``P`` is not universal and ``v ∉ P``

### Algorithm

`P` is universal iff it has no constraints.

A witness is produced using `isuniversal(H)` where `H` is the first linear
constraint of `P`.
"""
function isuniversal(P::AbstractPolyhedron{N}, witness::Bool=false
                    ) where {N<:Real}
    constraints = constraints_list(P)
    if isempty(constraints)
        return witness ? (true, N[]) : true
    else
        return witness ? isuniversal(constraints[1], true) : false
    end
end

"""
    constrained_dimensions(P::AbstractPolyhedron) where {N<:Real}

Return the indices in which a polyhedron is constrained.

### Input

- `P` -- polyhedron

### Output

A vector of ascending indices `i` such that the polyhedron is constrained in
dimension `i`.

### Examples

A 2D polyhedron with constraint ``x1 ≥ 0`` is constrained in dimension 1 only.
"""
function constrained_dimensions(P::AbstractPolyhedron)
    constraints = constraints_list(P)
    if isempty(constraints)
        return Int[]
    end
    zero_indices = zeros(Int, dim(P))
    for constraint in constraints
        for i in constrained_dimensions(constraint)
            zero_indices[i] = i
        end
    end
    return filter(x -> x != 0, zero_indices)
end

"""
    tosimplehrep(constraints::AbstractVector{LC})
        where {N<:Real, LC<:LinearConstraint{N}}

Return the simple H-representation ``Ax ≤ b`` from a list of linear constraints.

### Input

- `constraints` -- a list of linear constraints

### Output

The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
vector of offsets.
"""
function tosimplehrep(constraints::AbstractVector{LC}
                     ) where {N<:Real, LC<:LinearConstraint{N}}
    n = length(constraints)
    if n == 0
        A = Matrix{N}(undef, 0, 0)
        b = Vector{N}(undef, 0)
        return (A, b)
    end
    A = zeros(N, n, dim(first(constraints)))
    b = zeros(N, n)
    @inbounds begin
        for (i, Pi) in enumerate(constraints)
            A[i, :] = Pi.a
            b[i] = Pi.b
        end
    end
    return (A, b)
end

"""
     remove_redundant_constraints!(
         constraints::AbstractVector{<:LinearConstraint{N}};
         [backend]=default_lp_solver(N)) where {N<:Real}

Remove the redundant constraints of a given list of linear constraints; the list
is updated in-place.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `default_lp_solver(N)`) the backend used
                   to solve the linear program
### Output

`true` if the function was successful and the list of constraints `constraints`
is modified by removing the redundant constraints, and `false` only if the
constraints are infeasible.

### Notes

Note that the result may be `true` even if the constraints are infeasible.
For example, ``x ≤ 0 && x ≥ 1`` will return `true` without removing any
constraint.
To check if the constraints are infeasible, use
`isempty(HPolyhedron(constraints)`.

### Algorithm

If there are `m` constraints in `n` dimensions, this function checks one by one
if each of the `m` constraints is implied by the remaining ones.

To check if the `k`-th constraint is redundant, an LP is formulated using the
constraints that have not yet being removed.
If, at an intermediate step, it is detected that a subgroup of the constraints
is infeasible, this function returns `false`.
If the calculation finished successfully, this function returns `true`.

For details, see [Fukuda's Polyhedra
FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html).
"""
function remove_redundant_constraints!(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N)) where {N<:Real}

    A, b = tosimplehrep(constraints)
    m, n = size(A)
    non_redundant_indices = 1:m

    i = 1 # counter over reduced constraints

    for j in 1:m    # loop over original constraints
        α = A[j, :]
        Ar = A[non_redundant_indices, :]
        br = b[non_redundant_indices]
        br[i] = b[j] + one(N)
        lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend)
        if lp.status == :Infeasible
            # the polyhedron is empty
            return false
        elseif lp.status == :Optimal
            objval = -lp.objval
            if _leq(objval, b[j])
                # the constraint is redundant
                non_redundant_indices = setdiff(non_redundant_indices, j)
            else
                # the constraint is not redundant
                i = i+1
            end
        else
            error("LP is not optimal; the status of the LP is $(lp.status)")
        end
    end

    deleteat!(constraints, setdiff(1:m, non_redundant_indices))
    return true
end

"""
    remove_redundant_constraints(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N)) where {N<:Real}

Remove the redundant constraints of a given list of linear constraints.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `default_lp_solver(N)`) the backend used
                   to solve the linear program
### Output

The list of constraints with the redundant ones removed, or an empty set if the
constraints are infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N)) where {N<:Real}
    constraints_copy = copy(constraints)
    if remove_redundant_constraints!(constraints_copy, backend=backend)
        return constraints_copy
    else  # the constraints are infeasible
        return EmptySet{N}(dim(constraints[1]))
    end
end

# =======================
# Concrete linear map
# =======================

struct LinearMapInverse{T, MT<:AbstractMatrix{T}} <: AbstractLinearMapAlgorithm
    inverse::MT
end

struct LinearMapInverseRight <: AbstractLinearMapAlgorithm
    #
end

struct LinearMapLift <: AbstractLinearMapAlgorithm
    #
end

struct LinearMapElimination{T, S} <: AbstractLinearMapAlgorithm
    backend::T
    method::S
end

struct LinearMapVRep <: AbstractLinearMapAlgorithm
    #
end

function _check_algorithm_applies(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                ::Type{LinearMapInverse}; cond_tol=DEFAULT_COND_TOL, throw_error=false) where {N}

    inv_condition = issquare(M) && isinvertible(M; cond_tol=cond_tol)
    if !inv_condition
        throw_error && throw(ArgumentError("algorithm \"inverse\" requires an " *
                            "invertible matrix"))
        return false
    end

    dense_condition = !issparse(M)
    if !dense_condition
        throw_error && throw(ArgumentError("the inverse of a sparse matrix is not " *
            "available; either convert your matrix to a dense matrix with `Matrix(M)`, " *
            "or try the \"inverse_right\" algorithm"))
        return false
    end
    return true
end

function _check_algorithm_applies(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
            ::Type{LinearMapInverseRight}; cond_tol=DEFAULT_COND_TOL, throw_error=false) where {N}
    inv_condition = issquare(M) && isinvertible(M; cond_tol=cond_tol)
    if !inv_condition
        throw_error && throw(ArgumentError("algorithm \"inverse_right\" requires an " *
                            "invertible matrix"))
        return false
    end
    return true
end

function _check_algorithm_applies(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                                   ::Type{LinearMapLift}; throw_error=false) where {N}

    m, n = size(M)
    size_condition = m > n
    if !size_condition
        throw_error && throw(ArgumentError("algorithm \"lift\" requires that the number " *
            "of rows of the linear map is greater than the number of columns, but " *
            "they are of size $m and $n respectively"))
        return false
    end

    # rank condition
    r = rank(M)
    if r != n
        throw_error && throw(ArgumentError("the rank of the given matrix is " *
            "$r, but the algorithm \"lift\" requires it to be $n"))
        return false
    end

    return true
end

function _check_algorithm_applies(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                                   ::Type{LinearMapVRep}; throw_error=false) where {N}

    # TODO: the second check should be !isbounded(P) but this may be expensive;
    # see also #998 and #1926
    is_polytopic = P isa AbstractPolytope
    if !is_polytopic
        throw_error && throw(ArgumentError("algorithm \"vrep\" requires that the " *
            "list of vertices of the polyhedron is available, but it is not for a polyhedron " *
            "of type $(typeof(P))"))
        return false
    end
    return true
end

function _get_elimination_instance(N, backend, elimination_method)
    require(:Polyhedra; fun_name="linear_map with elimination")
    if backend == nothing
        require(:CDDLib; fun_name="linear_map with elimination")
        backend = default_cddlib_backend(N)
    end
    if elimination_method == nothing
        elimination_method = Polyhedra.BlockElimination()
    end
    return LinearMapElimination(backend, elimination_method)
end

function _default_linear_map_algorithm(M::AbstractMatrix{N}, P::AbstractPolyhedron{N};
                cond_tol=DEFAULT_COND_TOL, backend=nothing, elimination_method=nothing) where {N}

    if _check_algorithm_applies(M, P, LinearMapInverse, cond_tol=cond_tol)
        algo = LinearMapInverse(inv(M))
    elseif _check_algorithm_applies(M, P, LinearMapLift)
        algo = LinearMapLift()
    else
        algo = _get_elimination_instance(N, backend, elimination_method)
    end
    return algo
end

"""
    linear_map(M::AbstractMatrix{N},
               P::AbstractPolyhedron{N};
               [algorithm]::Union{String, Nothing}=nothing,
               [check_invertibility]::Bool=true,
               [cond_tol]::Number=DEFAULT_COND_TOL,
               [inverse]::Union{AbstractMatrix{N}, Nothing}=nothing,
               [backend]=nothing,
               [elimination_method]=nothing) where {N<:Real}

Concrete linear map of a polyhedral set.

### Input

- `M`         -- matrix
- `P`         -- polyhedral set
- `algorithm` -- (optional; default: `nothing`) algorithm to be used; for the
                 description see the Algorithm section below; possible choices are:

    - `"inverse"`, alias: `"inv"`
    - `"inverse_right"`, alias: `"inv_right"`
    - `"elimination"`, alias: `"elim"`
    - `"lift"`
    - `"vrep"`

- `check_invertibility` -- (optional, default: `true`) if `true` check whether
                           given matrix `M` is invertible; set to `false` only
                           if you know that `M` is invertible
- `cond_tol`  -- (optional; default: `DEFAULT_COND_TOL`) tolerance of matrix
                 condition (used to check whether the matrix is invertible)
- `inverse`   -- (optional; default: `nothing`) matrix inverse `M⁻¹`; use this
                 option if you have already computed the inverse matrix of `M`
- `backend`   -- (optional: default: `nothing`) polyhedra backend
- `elimination_method`  -- (optional: default: `nothing`) elimination method for
                           the `"elimination"` algorithm

### Output

The type of the result is "as close as possible" to the the type of `P`.
Let `(m, n)` be the size of `M`, where `m ≠ n` is allowed for rectangular maps.

To fix the type of the output to something different than the default value,
consider post-processing the result of this function with a call to a suitable
`convert` method.

In particular, the output depends on the type of `P`, on `m`, and the algorithm
that was used:

- If the vertex-based approach was used:

    - If `P` is a `VPolygon` and `m = 2` then the output is a `VPolygon`.
    - If `P` is a `VPolytope` then the output is a `VPolytope`.
    - Otherwise, the output is an `Interval` if `m = 1`, a `VPolygon` if `m = 2`
      and a `VPolytope` in other cases.

- If the invertibility criterion was used:

    - The types of `HalfSpace`, `Hyperplane`, `Line` and `AbstractHPolygon` are
      preserved.
    - If `P` is an `AbstractPolytope`, then the output is an `Interval` if `m = 1`,
      an `HPolygon` if `m = 2` and an `HPolytope` in other cases.
    - Otherwise, the output is an `HPolyhedron`.

### Notes

Since the different linear map algorithms work at the level of constraints (not sets
representations), this function uses dispatch on two stages: once the algorithm
has been defined, first the helper functions `_linear_map_hrep_helper` (resp.
`_linear_map_vrep`) are invoked, which dispatch on the set type. Then, each helper
function calls the concrete implementation of `_linear_map_hrep`, which dispatches
on the algorithm, and returns a list of constraints.

To simplify working with different algorithms and options, the
types `<: AbstractLinearMapAlgorithm` are used. These types are singleton type
or types that carry only the key data for the given algorithm, such as the matrix
inverse or the polyhedra backend.

New subtypes of the `AbstractPolyhedron` interface may define their own helper functions
`_linear_map_vrep`, respectively `_linear_map_hrep_helper` for special handling
of the constraints returned by the implementations of `_linear_map_hrep`;
otherwise the fallback implementation for `AbstractPolyhedron` is used, which
instantiates an `HPolyhedron`.

### Algorithm

This function mainly implements several approaches for the linear map: inverse,
right inverse, transformation to the vertex representation, variable elimination,
and variable lifting. Depending on the properties of `M` and `P`, one algorithm
may be preferable over the other. Details on the algorithms are given in the
following subsections.

Otherwise, if the algorithm argument is not specified, a default option is chosen based
on heuristics on the types and values of `M` and `P`:

- If the `"inverse"` algorithm applies, it is used.
- If the `"inverse_right"` algorithm applies, it is used.
- Otherwise, if the `"lift"` algorithm applies, it is used.
- Otherwise, the `"elimination"` algorithm is used.

Note that `"inverse"` does not require the external library `Polyhedra`, and neither
does `"inverse_right"`. However, the fallback method `"elimination"` requires
`Polyhedra` as well as the library `CDDLib`.

The optional keyword arguments `inverse` and `check_invertibility`
modify the default behavior:

- If an inverse matrix is passed in `inverse`, the given algorithm is applied,
  and if none is given, either `"inverse"` or `"inverse_right"` is applied
  (in that order of preference).
- If `check_invertibility` is set to `false`, the given algorithm is applied,
  and if none is given, either `"inverse"` or `"inverse_right"` is applied
  (in that order of preference).

#### Inverse

This algorithm is invoked with the keyword argument `algorithm="inverse"`
(or `algorithm="inv"`). The algorithm requires that `M` is invertible, square,
and dense. If you know a priori that `M` is invertible, set the flag
`check_invertibility=false`, such that no extra checks are done within `linear_map`.
Otherwise, we check the sufficient condition that the condition number
of `M` is not too high. The threshold for the condition number can be modified
from its default value, `DEFAULT_COND_TOL`, by passing a custom `cond_tol`.

The algorithm is described next. Assuming that the matrix ``M`` is invertible
(which we check via a sufficient condition,), ``y = M x`` implies
``x = \\text{inv}(M) y`` and we can transform the polyhedron
``A x ≤ b`` to the polyhedron ``A \\text{inv}(M) y ≤ b``.

If the dense condition on `M` is not fullfilled, there are two suggested
workarounds: either transform to dense matrix, i.e. calling `linear_map` with
`Matrix(M)`, or use the `"inverse_right"` algorithm, which does not compute the
inverse matrix explicitly, but uses a polyalgorithm; see the documentation
of `?\` for details.

#### Inverse-right

This algorithm is invoked with the keyword argument `algorithm="inverse_right"`
(or `algorithm="inv_right"`). This algorithm applies to square and invertible matrices
`M`. The idea is essentially the same as for the `"inverse"` algorithm; the difference
is that in `"inverse"` the full matrix inverse is computed, and in `"inverse_right"`
only the left division on the normal vectors is used. In particular, `"inverse_right"`
is good as a workaround when `M` is sparse (since the `inv` function is not available
for sparse matrices).

### Elimination

This algorithm is invoked with the keyword argument `algorithm = "elimination"` or
`algorithm = "elim"`. The algorithm applies to any matrix `M` (invertible or not),
and any polyhedron `P` (bounded or not).

The idea is described next. If `P : Ax <= b` and `y = Mx` denote the polyhedron
and the linear map respectively, we consider the vector `z = [y, x]`, write the
given equalities and the inequalities, and then eliminate the last x variables
(there are `length(x)` in total) using a call to `Polyhedra.eliminate` to a backend
library that can do variable elimination, typically `CDDLib` with the
`BlockElimination()` algorithm. In this way we have eliminated the "old" variables
`x` and kept the "new" or transformed variables "y".

The default elimination method is block elimination. For possible options we refer
to the documentation of Polyhedra,
[projection/elimination](https://juliapolyhedra.github.io/Polyhedra.jl/latest/projection/).

### Lift

This algorithm is invoked with the keyword argument `algorithm="lift"`.
The algorithm applies if `M` is rectangular of size `m × n` with `m > n` and
full rank (i.e. of rank `n`).

The idea is to embed the polyhedron into the `m`-dimensional space by appending zeros,
i.e. extending all constraints of `P` to `m` dimensions, and constraining the last
`m - n` dimensions to `0`. The matrix resulting matrix is extended to an invertible
`m × m` matrix and the algorithm using the inverse of the linear map is applied.
For the technical details of the extension of `M` to a higher-dimensional
invertible matrix, see `LazySets.Arrays.extend`.

### Vertex representation

This algorithm is invoked with the keyword argument `algorithm="vrep"`.
The idea is to convert the polyhedron to its vertex representation and apply the
linear map to each vertex of `P`.

The returned set is a polytope in vertex representation. Note that conversion of
the result back to half-space representation is not computed by default, since this
may be costly. If you used this algorithm and still want to convert back to
half-space representation, apply `tohrep` to the result of this function.
Note that this method only works for bounded polyhedra.
"""
function linear_map(M::AbstractMatrix{N},
                    P::AbstractPolyhedron{N};
                    algorithm::Union{String, Nothing}=nothing,
                    check_invertibility::Bool=true,
                    cond_tol::Number=DEFAULT_COND_TOL,
                    inverse::Union{AbstractMatrix{N}, Nothing}=nothing,
                    backend=nothing,
                    elimination_method=nothing) where {N<:Real}

   size(M, 2) != dim(P) && throw(ArgumentError("a linear map of size $(size(M)) " *
                            "cannot be applied to a set of dimension $(dim(P))"))

    got_algorithm = algorithm != nothing
    got_inv = got_algorithm && (algorithm == "inv" || algorithm == "inverse")
    got_inv_right = got_algorithm && (algorithm == "inv_right" || algorithm == "inverse_right")

    if inverse != nothing
        if !got_algorithm || got_inv
            algo = LinearMapInverse(inverse)
        elseif got_inv_right
            algo = LinearMapInverseRight()
        else
            throw(ArgumentError("received an inverse matrix but only the algorithms " *
                                "\"inverse\" and \"inverse_right\" apply, got $algorithm"))
        end
        return _linear_map_hrep_helper(M, P, algo)
    elseif !check_invertibility
        if got_inv || !issparse(M)
            inverse = inv(M)
            algo = LinearMapInverse(inverse)
        else
            algo = LinearMapInverseRight()
        end
        return _linear_map_hrep_helper(M, P, algo)
    end

    if !got_algorithm
        algo = _default_linear_map_algorithm(M, P; cond_tol=cond_tol)
        return _linear_map_hrep_helper(M, P, algo)

    elseif algorithm == "vrep"
        return _linear_map_vrep(M, P)

    elseif got_inv
        check_invertibility && _check_algorithm_applies(M, P, LinearMapInverse;
                                cond_tol=cond_tol, throw_error=true)
        return _linear_map_hrep_helper(M, P, LinearMapInverse(inv(M)))

    elseif got_inv_right
        check_invertibility && _check_algorithm_applies(M, P, LinearMapInverseRight;
                                cond_tol=cond_tol, throw_error=true)
        return _linear_map_hrep_helper(M, P, LinearMapInverseRight())

    elseif algorithm == "elimination" || algorithm == "elim"
        algo = _get_elimination_instance(N, backend, elimination_method)
        return _linear_map_hrep_helper(M, P, algo)

    elseif algorithm == "lift"
        _check_algorithm_applies(M, P, LinearMapLift, throw_error=true)
        return _linear_map_hrep_helper(M, P, LinearMapLift())

    else
        throw(ArgumentError("got unknown algorithm \"$algorithm\"; " *
            "available choices: \"inverse\", \"inverse_right\", \"lift\", " *
            "\"elimination\", \"vrep\""))
    end
end

# handle different numeric types
function linear_map(M::AbstractMatrix{NM},
                    P::AbstractPolyhedron{NP};
                    kwargs...) where {NM<:Real, NP<:Real}
    N = promote_type(NM, NP)
    if N != NP
        error("conversion between numeric types of polyhedra not implemented " *
            "yet (see #1181)")
    else
        return linear_map(N.(M), P; kwargs...)
    end
end

# TODO: merge the preconditions into _check_algorithm_applies ?
# review this method after #998
function _linear_map_vrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    if !isbounded(P)
        throw(ArgumentError("the linear map in vertex representation for an " *
            "unbounded set is not defined"))
    end
    require(:Polyhedra; fun_name="linear_map",
            explanation="of a $(typeof(P)) by a non-invertible matrix")
    # since P is bounded, we pass an HPolytope and then convert it to vertex representation
    P = tovrep(HPolytope(constraints_list(P), check_boundedness=false))
    return _linear_map_vrep(M, P)
end

# generic function for the AbstractPolyhedron interface => returns an HPolyhedron
function _linear_map_hrep_helper(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                                 algo::AbstractLinearMapAlgorithm) where {N<:Real}
    constraints = _linear_map_hrep(M, P, algo)
    return HPolyhedron(constraints)
end

# preconditions should have been checked in the caller function
function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          algo::LinearMapInverse) where {N}
    inverse = algo.inverse
    constraints_P = constraints_list(P)
    constraints_MP = _preallocate_constraints(constraints_P)
    @inbounds for (i, c) in enumerate(constraints_P)
        cinv = vec(_At_mul_B(c.a, inverse))
        constraints_MP[i] = LinearConstraint(cinv, c.b)
    end
    return constraints_MP
end

# preconditions should have been checked in the caller function
function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          algo::LinearMapInverseRight) where {N}
    constraints_P = constraints_list(P)
    constraints_MP = _preallocate_constraints(constraints_P)
    @inbounds for (i, c) in enumerate(constraints_P)
        # take left division for each constraint c, transpose(M) \ c.a
        cinv = _At_ldiv_B(M, c.a)
        constraints_MP[i] = LinearConstraint(cinv, c.b)
    end
    return constraints_MP
end

# preconditions should have been checked in the caller function
function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          algo::LinearMapLift) where {N<:Real}
    m, n = size(M)

    # we extend M to an invertible m x m matrix by appending m-n columns
    # orthogonal to the column space of M
    Mext, inv_Mext = extend(M, check_rank=false)

    # append zeros to the existing constraints, in the last m-n coordinates
    # TODO: cast to common vector type instead of Vector(c.a), see #1942, #1952
    cext = [HalfSpace(vcat(Vector(c.a), zeros(N, m-n)), c.b) for c in constraints_list(P)]

    # now fix the last m-n coordinates to zero
    id_out = Matrix(one(N)*I, m-n, m-n)
    cext = vcat(cext, [HalfSpace(vcat(zeros(N, n), id_out[i, :]), zero(N)) for i in 1:(m-n)],
                      [HalfSpace(vcat(zeros(N, n), -id_out[i, :]), zero(N)) for i in 1:(m-n)])

    Pext = HPolyhedron(cext)

    # now Mext is invertible and we can apply the inverse algorithm
    return _linear_map_hrep(Mext, Pext, LinearMapInverse(inv_Mext))
end

# If P : Ax <= b and y = Mx, we consider the vector z = [y, x], write the
# equalities and the inequalities, and then eliminate the last x variables
# (there are length(x) in total) using Polyhedra.eliminate calls
# to a backend library that can do variable elimination, typically CDDLib,
# with the BlockElimination() algorithm.
function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          algo::LinearMapElimination) where {N<:Real}
    m, n = size(M)
    ₋Id_m = Matrix(-one(N)*I, m, m)
    backend = algo.backend
    method = algo.method

    # extend the polytope storing the y variables first
    # append zeros to the existing constraints, in the last m-n coordinates
    # TODO: cast to common vector type instead of hard-coding Vector(c.a), see #1942 and #1952
    Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(N, m), Vector(c.a)), c.b) for c in constraints_list(P)]
    y_eq_Mx = [Polyhedra.HyperPlane(vcat(₋Id_m[i, :], Vector(M[i, :])), zero(N)) for i in 1:m]

    Phrep = Polyhedra.hrep(y_eq_Mx, Ax_leq_b)
    Phrep = polyhedron(Phrep, backend) # define concrete subtype
    Peli_block = Polyhedra.eliminate(Phrep, (m+1):(m+n), method)
    Peli_block = Polyhedra.removeduplicates(Polyhedra.hrep(Peli_block))

    # TODO: take constraints directly -- see #1988
    return constraints_list(HPolyhedron(Peli_block))
end

@inline function _preallocate_constraints(constraints::Vector{LinearConstraint{N, VN}}) where {N, VN<:AbstractVector{N}}
    return Vector{LinearConstraint{N, Vector{N}}}(undef, length(constraints))
end

"""
    plot_recipe(P::AbstractPolyhedron{N}, [ε]::N=zero(N)) where {N<:Real}

Convert a (bounded) polyhedron to a pair `(x, y)` of points for plotting.

### Input

- `P` -- bounded polyhedron
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

A pair `(x, y)` of points that can be plotted.

### Algorithm

We first assert that `P` is bounded (i.e., that `P` is a polytope).

One-dimensional polytopes are converted to an `Interval`.
Three-dimensional or higher-dimensional polytopes are not supported.

For two-dimensional polytopes (i.e., polygons) we compute their set of vertices
using `vertices_list` and then plot the convex hull of these vertices.
"""
function plot_recipe(P::AbstractPolyhedron{N}, ε::N=zero(N)) where {N<:Real}
    @assert dim(P) <= 2 "cannot plot a $(dim(P))-dimensional $(typeof(P))"
    @assert isbounded(P) "cannot plot an unbounded $(typeof(P))"

    if dim(P) == 1
        return plot_recipe(convert(Interval, P), ε)
    else
        vlist = transpose(hcat(convex_hull(vertices_list(P))...))
        if isempty(vlist)
            @warn "received a polyhedron with no vertices during plotting"
            return plot_recipe(EmptySet{N}(2), ε)
        end
        x, y = vlist[:, 1], vlist[:, 2]

        if length(x) > 1
            # add first vertex to "close" the polygon
            push!(x, x[1])
            push!(y, y[1])
        end

        return x, y
    end
end

"""
    chebyshev_center(P::AbstractPolyhedron{N};
                     [get_radius]::Bool=false,
                     [backend]=default_polyhedra_backend(P, N),
                     [solver]=JuMP.with_optimizer(GLPK.Optimizer)
                     ) where {N<:AbstractFloat}

Compute the [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
of a polytope.

### Input

- `P`          -- polytope
- `get_radius` -- (optional; default: `false`) option to additionally return the
                  radius of the largest ball enclosed by `P` around the
                  Chebyshev center
- `backend`    -- (optional; default: `default_polyhedra_backend(P, N)`) the
                  backend for polyhedral computations
- `solver`     -- (optional; default: `JuMP.with_optimizer(GLPK.Optimizer)`) the
                  LP solver passed to `Polyhedra`

### Output

If `get_radius` is `false`, the result is the Chebyshev center of `P`.
If `get_radius` is `true`, the result is the pair `(c, r)` where `c` is the
Chebyshev center of `P` and `r` is the radius of the largest ball with center
`c` enclosed by `P`.

### Notes

The Chebyshev center is the center of a largest Euclidean ball enclosed by `P`.
In general, the center of such a ball is not unique (but the radius is).
"""
function chebyshev_center(P::AbstractPolyhedron{N};
                          get_radius::Bool=false,
                          backend=default_polyhedra_backend(P, N),
                          solver=JuMP.with_optimizer(GLPK.Optimizer)
                         ) where {N<:AbstractFloat}
    require(:Polyhedra; fun_name="chebyshev_center")
    # convert to HPolyhedron to ensure `polyhedron` is applicable (see #1505)
    Q = polyhedron(convert(HPolyhedron, P); backend=backend)
    c, r = Polyhedra.chebyshevcenter(Q, solver)

    if get_radius
        return c, r
    end
    return c
end

"""
    an_element(P::AbstractPolyhedron{N}) where {N<:Real}

Return some element of a convex set.

### Input

- `P` -- polyhedron

### Output

An element of a polyhedron.

### Algorithm

An element of the polyhedron is obtained by evaluating its support vector along
direction ``[1, 0, …, 0]``.
"""
function an_element(P::AbstractPolyhedron{N}) where {N<:Real}
    n = dim(P)
    if n == -1
        throw(ArgumentError("the dimension of this polyhedron is not defined, " *
                            "hence `an_element` is not available"))
    end
    e₁ = SingleEntryVector(1, n, one(N))
    return σ(e₁, P)
end


text:
module BenchDictGetExisting

using BenchmarkTools
using ConcurrentCollections

function generate(; datasize = 2^10)
    vs = UInt64.(1:datasize)
    ks = string.(vs; pad = 16)
    # ks = vs
    # ks = UInt32.(vs)
    dict = ConcurrentDict{eltype(ks),eltype(vs)}(zip(ks, vs))
    return dict, ks
end

function sumall(dict, ks)
    s = zero(valtype(dict))
    for k in ks
        s += dict[k]
    end
    return s
end

const CACHE = Ref{Any}()

function setup()
    cdict, ks = generate()
    dict = Dict(cdict)
    CACHE[] = (cdict = cdict, dict = dict, ks = ks)

    suite = BenchmarkGroup()
    suite["base-seq"] = @benchmarkable(
        sumall(dict, ks),
        setup = begin
            dict = CACHE[].dict::$(typeof(dict))
            ks = CACHE[].ks::$(typeof(ks))
        end,
        evals = 1,
    )
    suite["cdict-seq"] = @benchmarkable(
        sumall(dict, ks),
        setup = begin
            dict = CACHE[].cdict::$(typeof(cdict))
            ks = CACHE[].ks::$(typeof(ks))
        end,
        evals = 1,
    )
    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
"""
    local_clustering_coefficient(g, v)
    local_clustering_coefficient(g, vs)

Return the [local clustering coefficient](https://en.wikipedia.org/wiki/Clustering_coefficient)
for node `v` in graph `g`. If a list of vertices `vs` is specified, return a vector
of coefficients for each node in the list.
"""
function local_clustering_coefficient(g::AbstractGraph, v::Integer)
    ntriang, nalltriang = local_clustering(g, v)
    return nalltriang == 0 ? 0. : ntriang * 1.0 / nalltriang
end

function local_clustering_coefficient(g::AbstractGraph, vs = vertices(g))
    ntriang, nalltriang = local_clustering(g, vs)
    return map(p -> p[2] == 0 ? 0. : p[1] * 1.0 / p[2], zip(ntriang, nalltriang))
end

function local_clustering!(storage::AbstractVector{Bool}, g::AbstractGraph, v::Integer)
    k = degree(g, v)
    k <= 1 && return (0, 0)
    neighs = neighbors(g, v)
    tcount = 0
    storage[neighs] .= true

    @inbounds for i in neighs
        @inbounds for j in neighbors(g, i)
            if (i != j) && storage[j]
                tcount += 1
            end
        end
    end
    return is_directed(g) ? (tcount, k * (k - 1)) : (div(tcount, 2), div(k * (k - 1), 2))
end

function local_clustering!(storage::AbstractVector{Bool},
                           ntriang::AbstractVector{Int},
                           nalltriang::AbstractVector{Int},
                           g::AbstractGraph,
                           vs)
    i = 0
    for (i, v) in enumerate(vs)
        ntriang[i], nalltriang[i] = local_clustering!(storage, g, v)
        storage[neighbors(g, v)] .= false
    end
    return ntriang, nalltriang
end

"""
    local_clustering(g, v)
    local_clustering(g, vs)

Return a tuple `(a, b)`, where `a` is the number of triangles in the neighborhood
of `v` and `b` is the maximum number of possible triangles. If a list of vertices
`vs` is specified, return two vectors representing the number of triangles and
the maximum number of possible triangles, respectively, for each node in the list.

This function is related to the local clustering coefficient `r` by ``r=\\frac{a}{b}``.
"""
function local_clustering(g::AbstractGraph, v::Integer)
    storage = zeros(Bool, nv(g))
    return local_clustering!(storage, g, v)
end

function local_clustering(g::AbstractGraph, vs = vertices(g))
    storage = zeros(Bool, nv(g))
    ntriang = zeros(Int, length(vs))
    nalltriang = zeros(Int, length(vs))
    return local_clustering!(storage, ntriang, nalltriang, g, vs)
end

"""
    triangles(g[, v])
    triangles(g, vs)

Return the number of triangles in the neighborhood of node `v` in graph `g`.
If a list of vertices `vs` is specified, return a vector of number of triangles
for each node in the list. If no vertices are specified, return the number
of triangles for each node in the graph.
"""
triangles(g::AbstractGraph, v::Integer) = local_clustering(g, v)[1]
triangles(g::AbstractGraph, vs = vertices(g)) = local_clustering(g, vs)[1]

"""
    global_clustering_coefficient(g)

Return the [global clustering coefficient](https://en.wikipedia.org/wiki/Clustering_coefficient)
of graph `g`.
"""
function global_clustering_coefficient(g::AbstractGraph)
    c = 0
    ntriangles = 0
    for v in vertices(g)
        neighs = neighbors(g, v)
        for i in neighs, j in neighs
            i == j && continue
            if has_edge(g, i, j)
                c += 1
            end
        end
        k = degree(g, v)
        ntriangles += k * (k - 1)
    end
    ntriangles == 0 && return 1.
    return c / ntriangles
end


text:
# Types for problem classes and solvers

export PipeLayoutSolver, GroundStructureSolver, JunctionLocationSolver
export PipeDimensioningSolver, optimize

"Base type for all solvers"
abstract type PipeLayoutSolver end

"Base type for all solvers depending on ground structure"
abstract type GroundStructureSolver <: PipeLayoutSolver end

"Base type for all solvers of junction location with fixed topology"
abstract type JunctionLocationSolver <: PipeLayoutSolver end

"Base type for all solvers of pipe dimensioning with fixed geometry"
abstract type PipeDimensioningSolver <: PipeLayoutSolver end

function optimize(inst::Instance, topo::Topology, solver)
    throw(ArgumentError("No method defined for solver type $(typeof(solver))"))
end