text:
using GraphPlot
using MetaGraphs
using Plots
using Compose
using Colors
using LightGraphs
using Cairo
using Fontconfig

"""
Splits the graph given symbol.
Return a list with the group each vertex is in. Example: you split
the graph by the :type property, so if it has 3 types ["a", "b", "c"]
it will return a list of being inside them like: [1, 1, 2, 3, 1, 2...]
"""
function divide(g::MetaDiGraph, s::Symbol, filterDivide::Bool=false)
    class = []
    for i ∈ 1:nv(g)
        p = props(g, i)
        if s in collect(keys(p))
            append!(class, [string(p[s])])
        else
            append!(class, ["other"])
        end
    end
    if filterDivide
        uniques = unique(class)
        for (origen, desti) in zip(uniques, 1:length(uniques))
            class = replace(class, origen => desti)
        end
    end
    return class
end

function divide_edges(g::MetaDiGraph, s::Symbol, filterDivide::Bool=false)
    class = []
    for ed in edges(g)
        p = props(g, Edge(ed.src, ed.dst))
        if s in collect(keys(p))
            append!(class, [string(p[s])])
        else
            append!(class, ["other"])
        end
    end
    if filterDivide
        uniques = unique(class)
        for (origen, desti) in zip(uniques, 1:length(uniques))
            class = replace(class, origen => desti)
        end
    end
    return class
end

""""""
function create_colors(g::MetaDiGraph, s::Symbol, pallete::AbstractVector)
    auxClass = divide(g, s, true)
    colors = distinguishable_colors(length(unique(auxClass)), pallete)
    return colors[auxClass]
end

function create_colors_edges(g::MetaDiGraph, s::Symbol, pallete::AbstractVector)
    auxClass = divide_edges(g, s, true)
    colors = distinguishable_colors(length(unique(auxClass)), pallete)
    return colors[auxClass]
end

"""
    plot(...)
A plot method for plotting Neo4j graphs

# Arguments
- `g::MetaDiGraph` :
- `layout` : random_layout, circular_layout, spectral_layout, shell_layout
"""
function plot(
    g::MetaDiGraph;
    # Node configuration:
    nodeLabel::Union{Nothing, Symbol} = nothing,
        nodeLabelColor::Union{RGB, Symbol} = colorant"black",
        nodeLabelOffsetDist::Real = 0,
        nodeLabelOffsetAngle::Real = π / 4.0,
    nodeSize::Union{Real, Symbol} = 1,
    nodeColor::Union{RGB, Symbol} = colorant"turquoise",

    # Edge configuration:
    edgeLabel::Union{Vector, Symbol} = [],
        edgeLabelColor::Union{RGB, Symbol} = colorant"black",
    edgeColor::Union{RGB, Symbol} = colorant"lightgray",
    edgeWidth::Real = 1,

    # Color pallete
    pallete::AbstractVector=distinguishable_colors(1),

    # Layout configuration
    layout=random_layout,

    #save config
    filename = nothing,
    )

    # Apply the configs
    if nodeLabel isa Symbol
        nodeLabel = divide(g, nodeLabel)

        # only color if label selected (for eficency)
        if nodeLabelColor isa Symbol
            nodeLabelColor = create_colors(g, nodeLabelColor, pallete)
        end
    end

    if nodeColor isa Symbol
        nodeColor = create_colors(g, nodeColor, pallete)
    end

    if nodeSize isa Symbol
        nodeSize = divide(g, nodeSize)
    end

    if edgeLabel isa Symbol
        edgeLabel = divide_edges(g, edgeLabel)

        if edgeLabelColor isa Symbol
            edgeLabelColor = create_colors_edges(g, edgeLabelColor, pallete)
        end
    end

    if edgeColor isa Symbol
        edgeColor = create_colors_edges(g, edgeColor, pallete)
    end

    p = gplot(
        g;
        nodelabel = nodeLabel,
        nodelabelc = nodeLabelColor,
        nodefillc = nodeColor,
        nodesize = nodeSize,
        nodelabeldist = nodeLabelOffsetDist,
        nodelabelangleoffset = nodeLabelOffsetAngle,
        nodelabelsize = nodeSize * 4,

        edgelabel = [],
        edgelabelc = edgeLabelColor,
        edgelabelsize = nodeSize * 4,
        edgelinewidth = edgeWidth,
        edgestrokec = edgeColor,

        layout = layout,
    )

    display(p)

    if filename !== nothing
        extension = filename[end-2:end]
        draw(eval(Meta.parse("$(uppercase(extension))(\"$(filename)\", 16cm, 16cm)")), p)
    end
end


text:
module Halide

# Write your package code here.
using Artifacts
halide = artifact"Halide"
root = joinpath(halide, readdir(halide)[1])

# 4 binaries
get_host_target() = run(` $(joinpath(root, "bin", "get_host_target"))`)
featurization_to_sample() = run(` $(joinpath(root, "bin", "featurization_to_sample"))`)
retrain_cost_model() = run(` $(joinpath(root, "bin", "retrain_cost_model"))`)
weightsdir_to_weightsfile() = run(` $(joinpath(root, "bin", "weightsdir_to_weightsfile"))`)

# 5 shared libs
libroot = joinpath(root, "lib")
libautoschedule_adams2019 = nothing ### ?
libautoschedule_mullapudi2016 = nothing ### ?
libautoschedule_li2018 = nothing ### ?
libHalide =  nothing ### libHalid.so

end

ext: jl


text:
function euler_step(q,p,dt,HS)
    Q = q+dt*HS.Minv*p
    P = p+dt*HS.F(q)
    return (dt,Q,P)
end

function i_euler_step(q,p,dt,HS;tol = 1e-6,max_iter=100,no_warning=false)
    Q=q;P=p
    for i in 1:max_iter
        Qn = q+dt*HS.Minv*P
        Pn = p+dt*HS.F(Q)
        if norm(Qn-Q)<=tol && norm(Pn-P)<=tol
            return (dt,Qn,Pn)
        end
        Q=Qn
        P=Pn
    end
    if !no_warning
        println("WARNING: Implicit euler step did not reach tol.")
    end
    return (dt,Q,P)
end

function midpoint_step(q,p,dt,HS)
    Q = q+dt*HS.Minv*p
    P = p+dt*HS.F(q)
    Q = q+dt*HS.Minv*((p+P)/2)
    P = p+dt*HS.F((q+Q)/2)
    return (dt,Q,P)
end

function i_midpoint_step(q,p,dt,HS;tol = 1e-6,max_iter=100,no_warning=false)
    Q=q;P=p
    for i in 1:max_iter
        Qn = q+dt*HS.Minv*((p+P)/2)
        Pn = p+dt*HS.F((q+Q)/2)
        if norm(Qn-Q)<=tol && norm(Pn-P)<=tol
            return (dt,Qn,Pn)
        end
        Q=Qn
        P=Pn
    end
    if !no_warning
        println("WARNING: Implicit midpoint step did not reach tol.")
    end
    return (dt,Q,P)
end

function theta_step(q,p,dt,HS;theta=0.5, kwargs...)
    _,Qe,Pe = euler_step(q,p,dt,HS)
    _,Qi,Pi = i_euler_step(q,p,dt,HS;kwargs...)
    Q,P = theta.*(Qe,Pe).+(1-theta).*(Qi,Pi)
    return (dt,Q,P)
end

function s_euler_step(q,p,dt,HS)
    P = p+dt*HS.F(q)
    Q = q+dt*HS.Minv*P
    return (dt,Q,P)
end

function sa_euler_step(q,p,dt,HS)
    Q = q+dt*HS.Minv*p
    P = p+dt*HS.F(Q)
    return (dt,Q,P)
end

function verlet_step(q,p,dt,HS)
    Q = q+dt*HS.Minv*p+dt^2/2*HS.Minv*HS.F(q)
    P = p+dt/2*(HS.F(q)+HS.F(Q))
    return (dt,Q,P)
end

function q_split_verlet_step(q,p,dt,HS)
    Q = q + dt/2*HS.Minv*p
    P = p + dt*(HS.F(Q))
    Q = Q + dt/2*HS.Minv*P
    return (dt,Q,P)
end

function p_split_verlet_step(q,p,dt,HS)
    P = p + dt/2*(HS.F(q))
    Q = q + dt*HS.Minv*P
    P = P + dt/2*(HS.F(Q))
    return (dt,Q,P)
end

function puls_projected_step(q,p,dt,HS;step=verlet_step,kwargs...)
    (dt,Q,P) = verlet_step(q,p,dt,HS;kwargs...);
    E = getEnergy(HS,HS.q0,HS.p0)
    if P'*HS.Minv*P<0.001 || U(Q)>E
        gamma = 1
    else
        gamma = sqrt(2*(E-HS.U(Q))/(P'*HS.Minv*P))
    end
    return (dt,Q,gamma*P)
end

function rk4_step(q,p,dt,HS)
    k1q = dt*HS.Minv*p
    k1p = dt*HS.F(q)

    k2q = dt*HS.Minv*(p+k1p/2)
    k2p = dt*HS.F(q+k1q/2)

    k3q = dt*HS.Minv*(p+k2p/2)
    k3p = dt*HS.F(q+k2q/2)

    k4q = dt*HS.Minv*(p+k3p)
    k4p = dt*HS.F(q+k3q)

    Q = q + 1/6*(k1q + 2*k2q + 2*k3q + k4q)
    P = p + 1/6*(k1p + 2*k2p + 2*k3p + k4p)
    return (dt,Q,P)
end

function projected_step(q,p,dt,HS;tol=1e-5,sigma=0.5,max_iter=100,step=rk4_step,no_warning=false,kwargs...)
    dt,Q,P = step(q,p,dt,HS,kwargs...)
    E0 = getEnergy(HS,HS.q0,HS.p0)
    E=getEnergy(HS,Q,P)-E0
    i=0;
    while abs(E)>tol
        if i > max_iter
            if !no_warning
                println("WARNING: projected step did not reach tol")
            end
            return (dt,Q,P)
        end
        dq = HS.F(Q)*E
        dp = -HS.Minv*P*E
        s=1/sigma
        diff = 0
        E_new = E
        j = 0
        while diff>=0
            if j>=max_iter
                if !no_warning
                    println("WARNING: projected step could not find sigma small enough!")
                end
                return (dt,Q,P)
            end
            s*=sigma
            E_new = getEnergy(HS,(Q+s*dq),(P+s*dp))-E0
            diff = E_new^2-E^2
            j+=1
        end
        E=E_new
        Q=Q+s*dq
        P=P+s*dp
        i+=1
    end
    return (dt,Q,P)
end

println("Finished loading HamiltonSteps")

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
function vtk_xml_write_header(vtk::DatasetFile)
    xroot = create_root(vtk.xdoc, "VTKFile")
    set_attribute(xroot, "type", vtk.grid_type)
    set_attribute(xroot, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    set_attribute(xroot, "header_type", string(HeaderType))
    if vtk.compression_level > 0
        set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    end
    xroot::XMLElement
end

# Specifies the "extent" of a structured grid, e.g. (1:4, 1:3, 1:42)
const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}

# Switch to zero-based extent used by VTK.
to_vtk_extent(r::AbstractUnitRange) = r .- 1

function to_vtk_extent(extent::Tuple)
    ext = to_extent3(extent)
    map(to_vtk_extent, ext)
end

# Convert different extent specifications to Extent{3}.
to_extent3(extent::Extent{3}) = extent

function to_extent3(extent::Extent{N}) where {N}
    if N > 3
        throw(ArgumentError("dimensionalities N > 3 not supported (got N = $N)"))
    elseif N < 3
        M = 3 - N
        (extent..., ntuple(d -> Base.OneTo(1), Val(M))...) :: Extent{3}
    end
end

to_extent3(Ns::Dims) = to_extent3(map(Base.OneTo, Ns))

# Returns the "extent" attribute required for structured (including rectilinear)
# grids.
function extent_attribute(extent_in)
    ext = to_vtk_extent(extent_in)
    iter = Iterators.map(e -> string(first(e), " ", last(e)), ext)
    join(iter, " ")
end

# Number of cells in structured grids.
# In 3D, all cells are hexahedrons (i.e. VTK_HEXAHEDRON), and the number of
# cells is (Ni-1)*(Nj-1)*(Nk-1). In 2D, they are quadrilaterals (VTK_QUAD), and in
# 1D they are line segments (VTK_LINE).
function num_cells_structured(Ns::Dims)
    prod(N -> max(1, N - 1), Ns)
end

# Accept passing coordinates as a tuple.
vtk_grid(filename, coords::NTuple; kwargs...) =
    vtk_grid(filename, coords...; kwargs...)


text:
module CGS
    SpeedOfLight = 2.99792458e10;                           # cm / s
    GravitationalConstant = 6.673e-8;                       # cm^3 / g s^2
    PlancksConstantH = 6.62606896e-27;                      # g cm^2 / s
    PlancksConstantHbar = 1.05457162825e-27;                # g cm^2 / s
    AstronomicalUnit = 1.49597870691e13;                    # cm
    LightYear = 9.46053620707e17;                           # cm
    Parsec = 3.08567758135e18;                              # cm
    GravAccel = 9.80665e2;                                  # cm / s^2
    ElectronVolt = 1.602176487e-12;                         # g cm^2 / s^2
    MassElectron = 9.10938188e-28;                          # g
    MassMuon = 1.88353109e-25;                              # g
    MassProton = 1.67262158e-24;                            # g
    MassNeutron = 1.67492716e-24;                           # g
    Rydberg = 2.17987196968e-11;                            # g cm^2 / s^2
    Boltzmann = 1.3806504e-16;                              # g cm^2 / K s^2
    MolarGas = 8.314472e7;                                  # g cm^2 / K mol s^2
    StandardGasVolume = 2.2710981e4;                        # cm^3 / mol
    Minute = 6e1;                                           # s
    Hour = 3.6e3;                                           # s
    Day = 8.64e4;                                           # s
    Week = 6.048e5;                                         # s
    Inch = 2.54e0;                                          # cm
    Foot = 3.048e1;                                         # cm
    Yard = 9.144e1;                                         # cm
    Mile = 1.609344e5;                                      # cm
    NauticalMile = 1.852e5;                                 # cm
    Fathom = 1.8288e2;                                      # cm
    Mil = 2.54e-3;                                          # cm
    Point = 3.52777777778e-2;                               # cm
    Texpoint = 3.51459803515e-2;                            # cm
    Micron = 1e-4;                                          # cm
    Angstrom = 1e-8;                                        # cm
    Hectare = 1e8;                                          # cm^2
    Acre = 4.04685642241e7;                                 # cm^2
    Barn = 1e-24;                                           # cm^2
    Liter = 1e3;                                            # cm^3
    UsGallon = 3.78541178402e3;                             # cm^3
    Quart = 9.46352946004e2;                                # cm^3
    Pint = 4.73176473002e2;                                 # cm^3
    Cup = 2.36588236501e2;                                  # cm^3
    FluidOunce = 2.95735295626e1;                           # cm^3
    Tablespoon = 1.47867647813e1;                           # cm^3
    Teaspoon = 4.92892159375e0;                             # cm^3
    CanadianGallon = 4.54609e3;                             # cm^3
    UkGallon = 4.546092e3;                                  # cm^3
    MilesPerHour = 4.4704e1;                                # cm / s
    KilometersPerHour = 2.77777777778e1;                    # cm / s
    Knot = 5.14444444444e1;                                 # cm / s
    PoundMass = 4.5359237e2;                                # g
    OunceMass = 2.8349523125e1;                             # g
    Ton = 9.0718474e5;                                      # g
    MetricTon = 1e6;                                        # g
    UkTon = 1.0160469088e6;                                 # g
    TroyOunce = 3.1103475e1;                                # g
    Carat = 2e-1;                                           # g
    UnifiedAtomicMass = 1.660538782e-24;                    # g
    GramForce = 9.80665e2;                                  # cm g / s^2
    PoundForce = 4.44822161526e5;                           # cm g / s^2
    KilopoundForce = 4.44822161526e8;                       # cm g / s^2
    Poundal = 1.38255e4;                                    # cm g / s^2
    Calorie = 4.1868e7;                                     # g cm^2 / s^2
    Btu = 1.05505585262e10;                                 # g cm^2 / s^2
    Therm = 1.05506e15;                                     # g cm^2 / s^2
    Horsepower = 7.457e9;                                   # g cm^2 / s^3
    Bar = 1e6;                                              # g / cm s^2
    StdAtmosphere = 1.01325e6;                              # g / cm s^2
    Torr = 1.33322368421e3;                                 # g / cm s^2
    MeterOfMercury = 1.33322368421e6;                       # g / cm s^2
    InchOfMercury = 3.38638815789e4;                        # g / cm s^2
    InchOfWater = 2.490889e3;                               # g / cm s^2
    Psi = 6.89475729317e4;                                  # g / cm s^2
    Poise = 1e0;                                            # g / cm s
    Stokes = 1e0;                                           # cm^2 / s
    Stilb = 1e0;                                            # cd / cm^2
    Lumen = 1e0;                                            # cd sr
    Lux = 1e-4;                                             # cd sr / cm^2
    Phot = 1e0;                                             # cd sr / cm^2
    Footcandle = 1.076e-3;                                  # cd sr / cm^2
    Lambert = 1e0;                                          # cd sr / cm^2
    Footlambert = 1.07639104e-3;                            # cd sr / cm^2
    Curie = 3.7e10;                                         # 1 / s
    Roentgen = 2.58e-7;                                     # A s / g
    Rad = 1e2;                                              # cm^2 / s^2
    SolarMass = 1.98892e33;                                 # g
    BohrRadius = 5.291772083e-9;                            # cm
    Newton = 1e5;                                           # cm g / s^2
    Dyne = 1e0;                                             # cm g / s^2
    Joule = 1e7;                                            # g cm^2 / s^2
    Erg = 1e0;                                              # g cm^2 / s^2
    StefanBoltzmannConstant = 5.67040047374e-5;             # g / K^4 s^3
    ThomsonCrossSection = 6.65245893699e-25;                # cm^2
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 65edde2e-d03e-11eb-300b-897bdc66d875
begin
import Pkg;
	Pkg.activate(".")
	Pkg.add(url="https://github.com/unfoldtoolbox/Unfold.jl")
	Pkg.add("PlutoUI")
	Pkg.add("StatsPlots")
	Pkg.add("DataFrames")
	Pkg.add("StatsModels")
end

# ╔═╡ 0a821c04-1aff-4fca-aab8-07a4d450838c
let
	using Revise
	using Unfold
	using StatsModels,PlutoUI,StatsPlots,DataFrames,Random
end

# ╔═╡ bec07615-8ca3-4251-9e44-bf5044828274
include("../test/test_utilities.jl");

# ╔═╡ 72bffaed-d3d9-44d8-9a74-0e0a0d3b0882
md"""# Simple Unfold.jl Example with overlap correction"""

# ╔═╡ ee305b1e-0c62-4a9f-ab80-d937b1401347
md""" ### Setup
Setup might take a while because the libraries are very large. It might be better to install some of the packages to your julia-environment prior to starting Pluto
"""

# ╔═╡ 2ec25bc8-1154-4850-ba02-62ed634cb558
md"""### Loading Data & Setting up Unfold"""

# ╔═╡ 76afd758-ac50-4f3c-a484-7f268ab059ef
data,evts = loadtestdata("test_case_3b");

# ╔═╡ c7fb216e-d177-4623-a135-db4cd856d4ad
f = @formula 0~1+conditionA+continuousA

# ╔═╡ a55058c9-3b56-447b-8f1a-f7319a29f478
md"""### Fitting the models"""

# ╔═╡ 8bfb1b5c-1c6a-472c-9087-dd03b0884922
md"""### Plotting the results"""

# ╔═╡ 9a4c715a-7a51-411d-a354-f79ec2369cb7
let
	md"""change window size τ = (-0.3, $(@bind τ2 Slider(0:0.1:5,default=2, show_value=true)))
	"""
end

# ╔═╡ ba2f0c48-f5f8-40eb-a45b-441b8134ffb1
τ = (-0.3,τ2)

# ╔═╡ f99efde3-bd34-46e2-ad1c-1a97ad866b79
fir = firbasis(τ=τ,sfreq=20,name="basisA");

# ╔═╡ dcbdef13-675e-47f7-b484-d340f0e55934
f_dict = Dict("eventA"=>(f,fir)) # combine eventname, formula & basis

# ╔═╡ ed61f099-343c-47d4-9af4-c4486c93fd4a
let
	md"""change noise: σ = $(@bind σ Slider(0:0.1:5,default=0, show_value=true))
	"""
end

# ╔═╡ 09705728-c6e5-41a5-9859-2c7e28828d5b
# add some noise
data_noise = data .+ σ .* randn(size(data));

# ╔═╡ b31a80a5-cca2-4c69-82f5-2b477f5b62be
begin
	plot(data_noise[1:600],title="Simulated Data with overlap")
	vline!(evts.latency[1:20],legend=false)
end

# ╔═╡ 3f4dd059-e275-483d-af42-a09a2795fca4
data_e,times = Unfold.epoch(data=data_noise,tbl=evts,τ=τ,sfreq=10);

# ╔═╡ 9aa42667-2629-48e1-b7fe-979dabc28f96
# non-overlapping (mass univariate)
m_e,res_e = fit(UnfoldLinearModel,f,evts,data_e,times);

# ╔═╡ bce42119-33d7-4e10-ac5e-b37028381b09
# overlap-corrected
m,res = fit(UnfoldLinearModel,f_dict,evts,data_noise,eventcolumn="type");

# ╔═╡ de267143-c029-4103-9fd7-ada153588ab6
begin

p1 = @df res_e plot(:colname_basis,:estimate,group=:coefname,title="Without ...",ylims=(-3,5),legend=false,xlims=(-0.3,3))
	p2 = @df res plot(:colname_basis,:estimate,group=:coefname,title="... with overlap correction",ylims=(-3,5),xlims=(-0.3,3),legend=:bottomright)
	plot(p1,p2,layout=(1,2))
end

# ╔═╡ Cell order:
# ╟─72bffaed-d3d9-44d8-9a74-0e0a0d3b0882
# ╟─ee305b1e-0c62-4a9f-ab80-d937b1401347
# ╠═65edde2e-d03e-11eb-300b-897bdc66d875
# ╠═0a821c04-1aff-4fca-aab8-07a4d450838c
# ╠═bec07615-8ca3-4251-9e44-bf5044828274
# ╟─2ec25bc8-1154-4850-ba02-62ed634cb558
# ╠═76afd758-ac50-4f3c-a484-7f268ab059ef
# ╠═09705728-c6e5-41a5-9859-2c7e28828d5b
# ╠═b31a80a5-cca2-4c69-82f5-2b477f5b62be
# ╠═3f4dd059-e275-483d-af42-a09a2795fca4
# ╠═c7fb216e-d177-4623-a135-db4cd856d4ad
# ╠═f99efde3-bd34-46e2-ad1c-1a97ad866b79
# ╠═dcbdef13-675e-47f7-b484-d340f0e55934
# ╟─a55058c9-3b56-447b-8f1a-f7319a29f478
# ╠═bce42119-33d7-4e10-ac5e-b37028381b09
# ╠═9aa42667-2629-48e1-b7fe-979dabc28f96
# ╟─8bfb1b5c-1c6a-472c-9087-dd03b0884922
# ╠═ba2f0c48-f5f8-40eb-a45b-441b8134ffb1
# ╟─9a4c715a-7a51-411d-a354-f79ec2369cb7
# ╟─ed61f099-343c-47d4-9af4-c4486c93fd4a
# ╠═de267143-c029-4103-9fd7-ada153588ab6


text:
using Documenter, NeXLUncertainties

function addNISTHeaders(htmlfile::String)
    # read HTML
    html = transcode(String,read(htmlfile))
    # Find </head>
    i = findfirst(r"</[Hh][Ee][Aa][Dd]>", html)
    # Already added???
    j = findfirst("nist-header-footer", html)
    if isnothing(j) && (!isnothing(i))
        # Insert nist-pages links right before </head>
        res = html[1:i.start-1]*
            "<link rel=\"stylesheet\" href=\"https://pages.nist.gov/nist-header-footer/css/nist-combined.css\">\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/jquery-1.9.0.min.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/nist-header-footer.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            html[i.start:end]
        write(htmlfile, res)
        println("Inserting NIST header/footer into $htmlfile")
    end
    return htmlfile
end

rm("build", force=true, recursive=true)

include("../weave/buildweave.jl")

makedocs(
    modules = [NeXLUncertainties],
    sitename = "NeXLUncertainties.jl",
    pages = [ "Home" => "index.md",
              "Getting Started" => "gettingstarted.md",
              "Composing Multi-Step Models" => "composing.md",
              "Resistor Network Example"=>"resistors.md",
              "Methods" => "methods.md" ]
)

names = ( "gettingstarted", "composing", "resistors" )
addNISTHeaders(joinpath(@__DIR__, "build","index.html"))
addNISTHeaders.(map(name->joinpath(@__DIR__, "build", splitext(name)[1], "index.html"), names))

rm("src/gettingstarted.md", force=true)
rm("src/resistors.md", force=true)
rm("src/resistor.png", force=true)

# deploydocs(repo = "github.com/NeXLUncertainties.jl.git")


text:
#! format: off

using SIIPExamples
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder
using DataFrames

using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

print_tree(PowerSimulations.PM.AbstractPowerModel)

ed_template = template_economic_dispatch(network = StandardPTDFModel)

set_device_model!(ed_template, HydroEnergyReservoir, HydroDispatchRunOfRiver)

PTDF_matrix = PTDF(sys)

problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 1,
    optimizer = solver,
    balance_slack_variables = true,
    constraint_duals = [
        :CopperPlateBalance,
        :network_flow__Line,
        :network_flow__TapTransformer,
    ],
    PTDF = PTDF_matrix,
)
build!(problem, output_dir = mktempdir())

solve!(problem)

res = ProblemResults(problem)
duals = get_duals(res)
λ = convert(Array, duals[:CopperPlateBalance][:, :CopperPlateBalance])
flow_duals = outerjoin(
    [duals[k] for k in [:network_flow__Line, :network_flow__TapTransformer]]...,
    on = :DateTime,
)
μ = Matrix(flow_duals[:, PTDF_matrix.axes[1]])

buses = get_components(Bus, sys)
congestion_lmp = Dict()
for bus in buses
    congestion_lmp[get_name(bus)] = μ * PTDF_matrix[:, get_number(bus)]
end
congestion_lmp = DataFrame(congestion_lmp)

LMP = λ .- congestion_lmp

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#==============================================================================
This example generates two noiseless measurement vectors corresponding to one
binary matrix and two different frequency vectors.

For each measurement vector, we compute MAP estimates, all global minimizers of
the negative log-likelihood, posterior means, and posterior standard deviations.

This is the first set of examples presented in Section 4.2 of the article:

L. Mustonen, X. Gao, A. Santana, R. Mitchell, Y. Vigfusson, and L. Ruthotto,
A Bayesian framework for molecular strain identification from mixed diagnostic samples
https://arxiv.org/abs/1803.02916
==============================================================================#

using StrainRecon
srand(0) # initialize the random number generator for reproducibility

# binary matrix and frequency vectors (ground truths)
M = [0 1; 1 0; 1 1]
w1 = [0.6; 0.4]
w2 = [0.5; 0.5]

# measurement data (without noise)
d1 = M*w1
d2 = M*w2

# dimensions of the problem (the number of classes is implicitly p=2)
m = size(M,1)
n = size(M,2)

R = 100 # inverse square root of the diagonal entry in the matrix Gamma
param1 = getStrainReconParam(d1,n,R=R)
param2 = getStrainReconParam(d2,n,R=R)

# verify that the block coordinate descent solver finds a global minimum
Mreco1,wreco1 = solveMAP(param1)
Mreco2,wreco2 = solveMAP(param2)
misfit1 = misfit(param1,Mreco1,wreco1)[1]
misfit2 = misfit(param2,Mreco2,wreco2)[1]
println("Misfit for the first MAP estimate is $misfit1")
println("Misfit for the second MAP estimate is $misfit2")
println()

# generate all possible binary matrices of size (m,n)
Mbits = [[bit=='1' for bit in bin(k,m*n)] for k in 0:2^(m*n)-1]
Mall = true*reshape(hcat(Mbits...),(m,n,2^(m*n)))
N = size(Mall,3)

# compute minimizing frequency vectors and the corresponding misfits for
# each binary matrix
wmins1 = zeros(n,N)
misfits1 = zeros(N)
wmins2 = zeros(n,N)
misfits2 = zeros(N)
for idx = 1:N
    wmins1[:,idx] = findw(param1,Mall[:,:,idx])
    misfits1[idx] = misfit(param1,Mall[:,:,idx],wmins1[:,idx])[1]

    wmins2[:,idx] = findw(param2,Mall[:,:,idx])
    misfits2[idx] = misfit(param2,Mall[:,:,idx],wmins2[:,idx])[1]
end

# sort misfits in ascending order
srt1 = sortperm(misfits1)
srt2 = sortperm(misfits2)

# print all pairs (M,w) that correspond to the minimum misfit
println("Case 1 minimizers:")
for idx = 1:N
    if abs(misfits1[srt1[idx]] - misfits1[srt1[1]]) > 1e-8
        break # not a global minimum anymore
    end
    Base.showarray(STDOUT,1*Mall[:,:,srt1[idx]],false)
    println()
    println(wmins1[:,srt1[idx]])
end
println()

println("Case 2 minimizers:")
for idx = 1:N
    if abs(misfits2[srt2[idx]] - misfits2[srt2[1]]) > 1e-8
        break # not a global minimum anymore
    end
    Base.showarray(STDOUT,1*Mall[:,:,srt2[idx]],false)
    println()
    println(wmins2[:,srt2[idx]])
end

# compute posterior means and variances
Mcm1,wcm1,dummy1,Mvar1,wvar1 = wrapperIntegratePosterior(param1)
Mcm2,wcm2,dummy2,Mvar2,wvar2 = wrapperIntegratePosterior(param2)

println("Case 1 posterior means:")
Base.showarray(STDOUT,Mcm1,false)
println()
println(wcm1)
println()

println("Case 1 posterior standard deviations:")
Base.showarray(STDOUT,sqrt.(Mvar1),false)
println()
println(sqrt.(wvar1))
println()

println("Case 2 posterior means:")
Base.showarray(STDOUT,Mcm2,false)
println()
println(wcm2)
println()

println("Case 2 posterior standard deviations:")
Base.showarray(STDOUT,sqrt.(Mvar2),false)
println()
println(sqrt.(wvar2))
println()


text:
module ExaTronAMDGPUKernels
    macro amdlocalmem(T, n)
        id = gensym("static_shmem")
        ex = quote
            ROCDeviceArray(($n,),alloc_local($(QuoteNode(id)), $T, $n))
        end
        esc(ex)
    end

    macro amdlocalmem(T, n1, n2)
        id = gensym("static_shmem")
        ex = quote
            ROCDeviceArray(($n1,$n2),alloc_local($(QuoteNode(id)), $T, $n1 * $n2))
        end
        esc(ex)
    end
    export @amdlocalmem
    using ..ExaTronKernels
    using ..AMDGPU
    include("daxpy.jl")
    include("dcopy.jl")
    include("ddot.jl")
    include("dmid.jl")
    include("dnrm2.jl")
    include("dgpnorm.jl")
    include("dscal.jl")
    include("dssyax.jl")
    include("dnsol.jl")
    include("dtsol.jl")
    include("dtrqsol.jl")
    include("dbreakpt.jl")
    include("dgpstep.jl")
    include("dicf.jl")
    include("dicfs.jl")
    include("dprsrch.jl")
    include("dcauchy.jl")
    include("dtrpcg.jl")
    include("dspcg.jl")
    include("dtron.jl")
    include("TronMatrix.jl")
    ExaTronKernels.supports(::ExaTronKernels.AMDGPUKernels) = true
end


text:
@inline ∂(a,I::CartesianIndex{d},f::AbstractArray{T,d}) where {T,d} = @inbounds f[I]-f[I-δ(a,I)]
@inline ∂(a,I::CartesianIndex{m},u::AbstractArray{T,n}) where {T,n,m} = @inbounds u[I+δ(a,I),a]-u[I,a]
@inline ϕ(a,I,f) = @inbounds (f[I]+f[I-δ(a,I)])*0.5
@fastmath quick(u,c,d) = median((5c+2d-u)/6,c,median(10c-9u,c,d))
@fastmath vanLeer(u,c,d) = (c≤min(u,d) || c≥max(u,d)) ? c : c+(d-c)*(c-u)/(d-u)
@inline ϕu(a,I,f,u,λ=quick) = @inbounds u>0 ? u*λ(f[I-2δ(a,I)],f[I-δ(a,I)],f[I]) : u*λ(f[I+δ(a,I)],f[I],f[I-δ(a,I)])
@fastmath @inline div(I::CartesianIndex{m},u) where {m} = sum(@inbounds ∂(i,I,u) for i ∈ 1:m)

@fastmath function tracer_transport!(r,f,u;Pe=0.1)
    N,n = size_u(u)
    for j ∈ 1:n
        @loop r[I] += ϕ(j,I,f)*u[I,j]-Pe*∂(j,I,f) over I ∈ slice(N,2,j,2)
        @loop (Φ = ϕu(j,I,f,u[I,j])-Pe*∂(j,I,f);
                r[I] += Φ; r[I-δ(j,I)] -= Φ) over I ∈ inside_u(N,j)
        @loop r[I-δ(j,I)] -= ϕ(j,I,f)*u[I,j]-Pe*∂(j,I,f) over I ∈ slice(N,N[j],j,2)
    end
end

@fastmath function conv_diff!(r,u;ν=0.1)
    r .= 0.
    N,n = size_u(u)
    for i ∈ 1:n, j ∈ 1:n
        @loop r[I,i] += ϕ(j,CI(I,i),u)*ϕ(i,CI(I,j),u)-ν*∂(j,CI(I,i),u) over I ∈ slice(N,2,j,2)
        @loop (Φ = ϕu(j,CI(I,i),u,ϕ(i,CI(I,j),u))-ν*∂(j,CI(I,i),u);
                r[I,i] += Φ; r[I-δ(j,I),i] -= Φ) over I ∈ inside_u(N,j)
        @loop r[I-δ(j,I),i] -= ϕ(j,CI(I,i),u)*ϕ(i,CI(I,j),u)-ν*∂(j,CI(I,i),u) over I ∈ slice(N,N[j],j,2)
    end
end

"""
    Flow{N,M,P}

Composite type for a multidimensional immersed boundary flow simulation.

Flow solves the unsteady incompressible [Navier-Stokes equations](https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations) on a Cartesian grid.
Solid boundaries are modelled using the [Boundary Data Immersion Method](https://eprints.soton.ac.uk/369635/).
The primary variables are the scalar pressure `p` (an array of dimension `N`)
and the velocity vector field `u` (an array of dimension `M=N+1`).
"""
struct Flow{N,M,P,T}
    # Fluid fields
    u :: Array{T,M} # velocity vector
    u⁰:: Array{T,M} # previous velocity
    f :: Array{T,M} # force vector
    p :: Array{T,N} # pressure scalar
    σ :: Array{T,N} # divergence scalar
    # BDIM fields
    V :: Array{T,M} # body velocity vector
    σᵥ:: Array{T,N} # body velocity divergence
    μ₀:: Array{T,M} # zeroth-moment on faces
    μ₁:: Array{T,P} # first-moment vector on faces
    # Non-fields
    U :: Vector{T}  # domain boundary values
    Δt:: Vector{T}  # time step
    ν :: T          # kinematic viscosity
    function Flow(N::Tuple,U::Vector;Δt=0.25,ν=0.,uλ::Function=(i,x)->0.,T=Float64)
        d = length(N); Nd = (N...,d)
        @assert length(U)==d
        u = Array{T}(undef,Nd...); apply!(uλ,u); BC!(u,U)
        u⁰ = copy(u)
        f,p,σ = zeros(T,Nd),zeros(T,N),zeros(T,N)
        V,σᵥ = zeros(T,Nd),zeros(T,N)
        μ₀ = ones(T,Nd); BC!(μ₀,zeros(T,d))
        μ₁ = zeros(T,N...,d,d)
        new{d,d+1,d+2,T}(u,u⁰,f,p,σ,V,σᵥ,μ₀,μ₁,U,[Δt],ν)
    end
end

@fastmath function BDIM!(a::Flow{n}) where n
    @. a.f = a.u⁰+a.Δt[end]*a.f-a.V
    for j ∈ 1:n, i ∈ 1:n
        @loop a.u[I,i] += μddn(j,CI(I,i),a.μ₁,a.f) over I ∈ inside(a.p)
    end
    @. a.u += a.V+a.μ₀*a.f
end
@inline μddn(j,I::CartesianIndex,μ,f) = @inbounds 0.5μ[I,j]*(f[I+δ(j,I)]-f[I-δ(j,I)])

@fastmath function project!(a::Flow{n},b::AbstractPoisson{n},w=1) where n
    @inside a.σ[I] = (div(I,a.u)+w*a.σᵥ[I])/a.Δt[end]
    solver!(a.p,b,a.σ)
    for i ∈ 1:n
        @loop a.u[I,i] -= a.Δt[end]*a.μ₀[I,i]*∂(i,I,a.p) over I ∈ inside(a.σ)
    end
end

"""
    mom_step!(a::Flow,b::AbstractPoisson)

Integrate the `Flow` one time step using the [Boundary Data Immersion Method](https://eprints.soton.ac.uk/369635/)
and the `AbstractPoisson` pressure solver to project the velocity onto an incompressible flow.
"""
@fastmath function mom_step!(a::Flow,b::AbstractPoisson)
    a.u⁰ .= a.u; a.u .= 0
    # predictor u → u'
    conv_diff!(a.f,a.u⁰,ν=a.ν)
    BDIM!(a); BC!(a.u,a.U)
    project!(a,b); BC!(a.u,a.U)
    # corrector u → u¹
    conv_diff!(a.f,a.u,ν=a.ν)
    BDIM!(a); BC!(a.u,a.U,2)
    project!(a,b,2); a.u ./= 2; BC!(a.u,a.U)
    push!(a.Δt,CFL(a))
end

function CFL(a::Flow{n}) where n
    mx = maximum(fout(I,a.u) for I ∈ inside(a.p))
    min(10.,inv(mx+5a.ν))
end
@fastmath @inline fout(I::CartesianIndex{d},u) where {d} =
    sum(@inbounds(max(0.,u[I+δ(a,I),a])+max(0.,-u[I,a])) for a ∈ 1:d)


text:
using StanBase
using Test

@testset "Basic HelpModel" begin

  @test CMDSTAN_HOME == get_cmdstan_home()

end


text:
using FastTransforms, Compat
using Compat.Test

import FastTransforms: normalizecolumns!, maxcolnorm

if VERSION ≥ v"0.7-"
    vecnorm(A, p...) = norm(A, p...)
end

@testset "Synthesis and analysis" begin

    # Starting with normalized spherical harmonic coefficients,

    n = 50
    F = sphrandn(Float64, n, n);
    normalizecolumns!(F);

    # we convert to bivariate Fourier series.

    G = sph2fourier(F);

    # At equispaced points in angle,

    θ = (0.5:n-0.5)*π/n
    φ = (0:2n-2)*2π/(2n-1)

    # the Fourier series evaluates to:

    SF = [sum(G[ℓ,1]*cos((ℓ-1)*θ)/sqrt(2π) for ℓ in 1:n) + sum(G[ℓ,2m]*sin(ℓ*θ)*sin(m*φ)/sqrt(π) for ℓ in 1:n, m in 1:2:n-1) + sum(G[ℓ,2m+1]*sin(ℓ*θ)*cos(m*φ)/sqrt(π) for ℓ in 1:n, m in 1:2:n-1) + sum(G[ℓ,2m]*cos((ℓ-1)*θ)*sin(m*φ)/sqrt(π) for ℓ in 1:n, m in 2:2:n-1) + sum(G[ℓ,2m+1]*cos((ℓ-1)*θ)*cos(m*φ)/sqrt(π) for ℓ in 1:n, m in 2:2:n-1) for θ in θ, φ in φ]

    # but that was slow, so we accelerate it via in-place FFTW technology:

    Ps = FastTransforms.plan_synthesis(G);

    Y = zero(G);

    mul!(Y, Ps, G)

    @test maxcolnorm(SF - Y) < 10000eps()

    # Retracing our steps, function values on the sphere are converted to Fourier coefficients:

    Pa = FastTransforms.plan_analysis(Y);

    Z = zero(Y);

    mul!(Z, Pa, Y)

    @test maxcolnorm(Z - G) < 10eps()

    # And Fourier coefficients are converted back to spherical harmonic coefficients:

    H = fourier2sph(Z)

    @test maxcolnorm(F - H) < 100eps()
end

@testset "Test for dropping last column" begin
    for f in ((θ,φ)->1/(3 + cospi(φ) + sinpi(θ)),
              (θ,φ)->cos(50*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)),
              (θ,φ)->cos(50*cospi(φ)*sinpi(θ)+80*sinpi(φ)*sinpi(θ)),
              (θ,φ)->sqrt(5+cospi(φ)*sinpi(θ)+exp(sinpi(φ)*sinpi(θ))+sin(cospi(θ))))
        n = 20
        UO = sphones(Float64, n, n)
        UE = sphones(Float64, n, n)
        while norm(UO[:,end]) > vecnorm(UO)*eps()
            θ = (0.5:n-0.5)/n
            φ = (0:2n-2)*2/(2n-1)
            F = [f(θ,φ) for θ in θ, φ in φ]
            V = zero(F)
            mul!(V, FastTransforms.plan_analysis(F), F)
            UO = fourier2sph(V)

            φ = (0:2n-3)*2/(2n-2)
            F = [f(θ,φ) for θ in θ, φ in φ]
            V = zero(F)
            mul!(V, FastTransforms.plan_analysis(F), F)
            UE = fourier2sph(V)

            n *= 2
        end
        @test vecnorm(UO[:,1:end-1] - UE) < n*vecnorm(UO)*eps()
    end
end

@testset "Test for sampling through the poles" begin
    for f in ((θ,φ)->cos(50*cospi(φ)*sinpi(θ)*sinpi(φ)*sinpi(θ)),
              (θ,φ)->cos(50*cospi(φ)*sinpi(θ)+80*sinpi(φ)*sinpi(θ)),
              (θ,φ)->sqrt(5+cospi(φ)*sinpi(θ)+exp(sinpi(φ)*sinpi(θ))+sin(cospi(θ))))
        n = 200

        θ = (0.5:n-0.5)/n
        φ = (0:2n-2)*2/(2n-1)
        F = [f(θ,φ) for θ in θ, φ in φ]
        V = zero(F)
        mul!(V, FastTransforms.plan_analysis(F), F)
        G = zero(V)
        mul!(G, FastTransforms.plan_synthesis(V), V)

        θ2 = (0.0:n-1)/(n-1)
        F2 = [f(θ,φ) for θ in θ2, φ in φ]
        V2 = zero(F2)
        mul!(V2, FastTransforms.plan_analysis2(F2), F2)
        G2 = zero(V2)
        mul!(G2, FastTransforms.plan_synthesis2(V2), V2)

        @test vecnorm(V-V2) < n*vecnorm(V)*eps()
        @test vecnorm(F-G) < n*vecnorm(F)*eps()
        @test vecnorm(F2-G2) < n*vecnorm(F)*eps()
    end
end

@testset "Degree-3 polynomial" begin
    # This test confirms numerically that [P_4(z⋅y) - P_4(x⋅y)]/(z⋅y - x⋅y) is actually a degree-3 polynomial on 𝕊²
    x = [0,0,1]
    y = normalize!([.123,.456,.789])

    z = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]

    P4 = x -> (35*x^4-30*x^2+3)/8

    n = 5
    θ = (0.5:n-0.5)/n
    φ = (0:2n-2)*2/(2n-1)
    F = [(P4(z(θ,φ)⋅y) - P4(x⋅y))/(z(θ,φ)⋅y - x⋅y) for θ in θ, φ in φ]
    V = zero(F)
    mul!(V, FastTransforms.plan_analysis(F), F)
    U3 = fourier2sph(V)

    # U3 is degree-3

    F = [P4(z(θ,φ)⋅y) for θ in θ, φ in φ]
    V = zero(F)
    mul!(V, FastTransforms.plan_analysis(F), F)
    U4 = fourier2sph(V)

    # U4 is degree-4
end


text:
# state_dict = Dict{String,String}(
#     "India"=>"UN",
#     "Andhra Pradesh"=>"AP",
#     "Bihar"=>"BR",
#     "Karnataka"=>"KA",
#     "Assam"=>"AS",
#     "Chandigarh"=>"CH",
#     "West Bengal"=>"WB",
#     "Jharkhand"=>"JH",
#     "Madhya Pradesh"=>"MP",
#     "Delhi"=>"DL",
#     "Dadra and Nagar Haveli"=>"DN",
#     "Tripura"=>"TR",
#     "Uttarakhand"=>"UT",
#     "Haryana"=>"HR",
#     "Andaman and Nicobar Islands"=>"AN",
#     "Gujarat" => "GJ",
#     "Manipur" => "MN",
#     "Telangana" => "TG",
#     "Jammu and Kashmir" => "JK",
#     "Tamil Nadu" => "TN",
#     "Goa" => "GA",
#     "Nagaland" => "NL",
#     "TT" => "TT",
#     "Uttar Pradesh" => "UP",
#     "Puducherry" => "PY",
#     "Chhattisgarh" => "CT",
#     "Mizoram" => "MZ",
#     "Punjab" => "PB",
#     "Himachal Pradesh" => "HP",
#     "Arunachal Pradesh" => "AR",
#     "Ladakh"=>"LA",
#     "Meghalaya"=>"ML",
#     "Kerala"=>"KL",
#     "Maharashtra"=>"MH",
#     "Lakshadweep"=>"LD",
#     "Odisha"=>"OD",
#     "Rajasthan"=>"RJ",
#     "Sikkim"=>"SK"

# )


text:
module ADjson

using Statistics
using GZip
using JSON
using TensorCast
using Dates
using TimeZones
using ADerrors

"""
    load_json(fname::String)

Load data from a json.gz file

# Currently supported datatypes
- `Obs`
- `List`

# Arguments
- `fname::string`: File name. Suffix '.json.gz' is automatically added if not specified explicitly.

# Notes
- ADerrors always assumes that replica are labeled as r0, r1, etc. Custom replica names get lost in the conversion.
- Irregular Monte Carlo chains are currently not supported as not all cases pyerrors can output can be safely converted to ADerrors
"""
function load_json(fname::String)

    if !endswith(fname, ".json.gz")
        fname *= ".json.gz"
    end

    if !(isfile(fname) || islink(fname))
        error("File '" * fname * "' does not exist.")
    end

    df = GZip.open(fname, "r") do io
        JSON.parse(io)
    end

    println("Data has been written using ", df["program"])
    println("Format version ", df["version"])
    println("Written by ", df["who"], " on ", df["date"], " on host ", df["host"])
    if haskey(df, "description")
        println("\nDescription: ", df["description"])
    end

    res_list = Vector{Any}()

    for entry in df["obsdata"]
        if !(entry["type"] in ["Obs", "List", "Array", "Corr"])
            error("Type '" * entry["type"] * "' is not implemented." )
        end

        res = Vector{uwreal}()
        out_length = prod([parse(Int, o) for o in split(entry["layout"], ", ")])

        if haskey(entry, "data")
            for element in entry["data"]
                conc_deltas = [Float64[] for b in 1:out_length]
                int_cnfg_numbers = Int[]
                rep_lengths = Int[]

                for rep in element["replica"]
                    my_arr = rep["deltas"]
                    @cast deltas[i][j] := my_arr[j][i]
                    cnfg_numbers = convert(Array{Int,1}, deltas[1])

                    if length(element["replica"]) > 1
                        for ch in 1:length(cnfg_numbers)
                            if ch != cnfg_numbers[ch]
                                error("Irregular Monte Carlo chains for multiple replica cannot be safely read into ADerrors.")
                            end
                        end
                    end

                    append!(int_cnfg_numbers, cnfg_numbers)
                    append!(rep_lengths, length(cnfg_numbers))
                    for i in 1:out_length
                        append!(conc_deltas[i], Vector{Float64}(deltas[i + 1]) .- mean(deltas[i + 1]))
                    end
                end

                if length(res) == 0
                    for i in 1:out_length
                        if length(element["replica"]) > 1
                            push!(res, uwreal(Vector{Float64}(conc_deltas[i]), element["id"], rep_lengths))
                        else
                            push!(res, uwreal(Vector{Float64}(conc_deltas[i]), element["id"], int_cnfg_numbers, int_cnfg_numbers[end]))
                        end
                    end
                else
                    for i in 1:out_length
                        if length(element["replica"]) > 1
                            res[i] += uwreal(Vector{Float64}(conc_deltas[i]), element["id"], rep_lengths)
                        else
                            res[i] += uwreal(Vector{Float64}(conc_deltas[i]), element["id"], int_cnfg_numbers, sum(rep_lengths))
                        end
                    end
                end
            end
        end

        if haskey(entry, "cdata")
            for element in entry["cdata"]
                if element["layout"] != "1, 1"
                   error("ADerror does not support non-scalar constants at the moment.")
                end

                if length(res) == 0
                    for i in 1:out_length
                        push!(res, uwreal([0.0, sqrt(element["grad"][1][i] * element["cov"][1] * element["grad"][1][i])], element["id"]))
                    end
                else
                    for i in 1:out_length
                        res[i] += uwreal([0.0, sqrt(element["grad"][1][i] * element["cov"][1] * element["grad"][1][i])], element["id"])
                    end
                end
            end
        end

        for i in 1:out_length
            res[i] += entry["value"][i]
        end

        if length(res) == 1
            push!(res_list, res[1])
        else
            push!(res_list, reshape(res, Tuple(Int(x) for x in [parse(Int, o) for o in split(entry["layout"], ", ")])))
        end
    end

    if length(res_list) == 1
        return res_list[1]
    else
        return res_list
    end
end

"""
    dump_to_json(p, fname::String, description="", indent=1)

Dump data to a json.gz file

# Currently supported datatypes
- `Obs`

# Arguments
- `p::uwreal` or `p::Vector{uwreal}`: Data to be written to the file
- `fname::String`: File name. Suffix '.json.gz' is automatically added if not specified explicitly.
- `description::String`: Optional description to be added to the file
- `indent::Int`: Indent of the json structure. Use `ident=nothing` for smallest file size but less human readability.
"""
function dump_to_json(data, fname::String, description="", indent=1)
    jsonstring = Dict{String, Any}()
    jsonstring["program"] = "ADjson 1.1"
    jsonstring["who"] = ENV["USER"]
    jsonstring["host"] = gethostname() * ", " * Sys.MACHINE
    jsonstring["date"] = Dates.format(now(localzone()), "Y-mm-dd HH:MM:SS zzzz")
    jsonstring["version"] = "1.1"
    jsonstring["obsdata"] = []
    if !isempty(description)
        jsonstring["description"] = description
    end

    ws = ADerrors.wsg

    vec_data = Vector{uwreal}()
    if data isa uwreal
        push!(vec_data, data)
    elseif data isa Vector{uwreal}
        append!(vec_data, data)
    else
        error("Unkown data type.")
    end

    for (index, p) in enumerate(vec_data)
        push!(jsonstring["obsdata"], Dict{String, Any}())
        jsonstring["obsdata"][index]["layout"] = "1"
        jsonstring["obsdata"][index]["type"] = "Obs"
        jsonstring["obsdata"][index]["value"] = [p.mean]

        i_data = 0
        i_cdata = 0
        for i in 1:convert(Int32, ADerrors.unique_ids!(p, ws))
            # Collect and rename replica names
            my_rep_names = Vector{String}()
            for el in ADerrors.get_repnames_from_id(p.ids[i], ws)
                loc = findlast("_", el)
                push!(my_rep_names, el[1:Vector{Int}(loc)[1] - 1] * "|" * el[Vector{Int}(loc)[1] + 1:end])
            end

            # Get deltas
            dt = zeros(Float64, ws.fluc[ws.map_ids[p.ids[i]]].nd)
            for j in 1:length(p.prop)
                if (p.prop[j] && (ws.map_nob[j] == p.ids[i]))
                    dt .= dt .+ p.der[j] .* ws.fluc[j].delta
                end
            end

            idc = convert(Vector{Int32}, ADerrors.get_repidc_from_id(p.ids[i], ws))
            if length(idc) == 2
                # Write covobs data
                i_cdata += 1
                if !haskey(jsonstring["obsdata"][index], "cdata")
                    jsonstring["obsdata"][index]["cdata"] = []
                end
                push!(jsonstring["obsdata"][index]["cdata"], Dict{String, Any}())
                jsonstring["obsdata"][index]["cdata"][i_cdata]["id"] = ADerrors.get_name_from_id(p.ids[i], ws)
                jsonstring["obsdata"][index]["cdata"][i_cdata]["layout"] = "1, 1"
                jsonstring["obsdata"][index]["cdata"][i_cdata]["cov"] = [dt[1] ^ 2]
                jsonstring["obsdata"][index]["cdata"][i_cdata]["grad"] = [[1.0]]
            else
                # Write Monte Carlo data
                i_data += 1
                if !haskey(jsonstring["obsdata"][index], "data")
                    jsonstring["obsdata"][index]["data"] = []
                end
                push!(jsonstring["obsdata"][index]["data"], Dict{String, Any}())
                jsonstring["obsdata"][index]["data"][i_data]["id"] = ADerrors.get_name_from_id(p.ids[i], ws)
                jsonstring["obsdata"][index]["data"][i_data]["replica"] = []
                rep_indices = prepend!(cumsum(convert(Vector{Int32}, ws.fluc[ws.map_ids[p.ids[i]]].ivrep)), 0)
                for j in 1:length(my_rep_names)
                    for delta_entry in dt[rep_indices[j] + 1:rep_indices[j + 1]]
                        if isapprox(delta_entry, 0.0; atol=1e-14)
                            error("Irregular Monte Carlo chain cannot be safely written to json.gz file.")
                        end
                    end
                    tmp = []
                    push!(jsonstring["obsdata"][index]["data"][i_data]["replica"], Dict{String, Any}())
                    jsonstring["obsdata"][index]["data"][i_data]["replica"][j]["name"] = my_rep_names[j]
                    push!(tmp, idc[rep_indices[j] + 1:rep_indices[j + 1]])
                    push!(tmp, dt[rep_indices[j] + 1:rep_indices[j + 1]])
                    @cast deltas[i][j] := tmp[j][i]
                    jsonstring["obsdata"][index]["data"][i_data]["replica"][j]["deltas"] = deltas
                end
            end
        end
    end
    if !endswith(fname, ".json.gz")
        fname *= ".json.gz"
    end

    GZip.open(fname, "w") do f
        JSON.print(f, jsonstring, indent)
    end
end

export load_json, dump_to_json

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Result (ModiaMath/Result/_module.jl)
#

#---------------------------- Acessing the time series results (various options)
appendUnit(name, unit) = unit == "" ? string(name) : string(name, " [", unit, "]")

"""
    indexOftrailingDot(name::AbstractString)::Int

Return the index of the trailing "." of name. If there is no ".xxx", zero is returned
"""
function indexOftrailingDot(name::AbstractString)::Int
   i = first(something(findlast(".", name), 0:-1))
   return i > 0 && i < length(name) ? i : 0
end

"""
    getStringDictResult(model, res)

Return dictionary of result, from raw result data structure `raw`.
"""
function getStringDictResult(model, res::RawResult)
    nt     = res.nt
    data   = res.data
    result = Dict{AbstractString,Any}()
    for (i, name) in enumerate(res.names)
        result[name] = @view data[1:nt,i]
    end

    return result, nt
end

const StringDictAnyResult = Dict{AbstractString,Any}
resultHeading(result::StringDictAnyResult) = ""
variablesDependingOnStruct(state::Any) = 0

function getSignal(seriesDict::StringDictAnyResult, name::AbstractString)
    hasSignal = false
    if haskey(seriesDict, name)
        sig       = seriesDict[name]
        keyName   = name
        hasSignal = true
    else
        # Handle signal arrays, auch as  a.b.c[3]
        if name[end] == ']'
            indexRange = something(findlast("[", name), 0:-1)
            i = indexRange[1]
            @assert(i >= 2)
            keyName = name[1:i-1]

            if haskey(seriesDict, keyName)
                sig2 = seriesDict[keyName]

                # Build AST and evaluate: sig = sig2[:,...]
                ex1 = Meta.parse(name[i:end])
                ex2 = :( $sig2[:] )
                append!(ex2.args, ex1.args)
                sig = @eval $ex2
                if ndims(sig) != 1
                    @warn "ModiaMath.plot: argument name (= $name) does not characterize one array element\nIndex ranges are not yet supported."
                    return (nothing, keyName, name)
                end
                hasSignal = true
            end

        elseif (i = indexOftrailingDot(name)) > 0
            # Handle vector of structs, such as a.b.c and a.b is the key and c is the field in struct a.b[i]
            keyName   = name[1:i-1]
            fieldName = Symbol(name[i+1:end])
            if haskey(seriesDict, keyName)
                sig3 = seriesDict[keyName]
                if typeof(sig3) <: AbstractVector && length(sig3) > 0 && isstructtype(typeof(sig3[1]))
                    # Signal is a vector of structs and has at least one element
                    if isdefined(sig3[1], fieldName) && fieldtype(typeof(sig3[1]), fieldName) <: Number
                        # Struct has fieldName and fieldName is a number
                        sig = zeros(length(sig3))
                        for i in 1:length(sig3)
                           sig[i] = getfield(sig3[i], fieldName)
                        end
                        hasSignal = true
                    else
                        # Check whether dependent variables are defined
                        vdict = variablesDependingOnStruct(sig3[1])
                        key   = name[i+1:end]
                        if typeof(vdict) != Int
                            if haskey(vdict, key)
                                # Dependent variable fieldName is defined
                                sig = zeros(length(sig3))
                                fc  = vdict[key]
                                for i in 1:length(sig3)
                                    sig[i] = fc(sig3[i])
                                end
                                hasSignal = true
                            end
                        end
                    end
                end
            end
        end
    end

    if hasSignal
        return (sig, keyName, name)

    else
        @warn "ModiaMath.plot: argument name (= $name) is not correct or does not identify a signal in the result."
        return (nothing, name, name)
    end
end

"""
    (xsig, xsigLegend, ysig, ysigLegend) =
           ModiaMath.resultTimeSeries(result, name, xLabel::Bool, xAxis)

For a desired `result` data structure, this function has to be provided
to return the x-vector (`xsig`), the y-vector (`ysig`) and the legend
of the x-vector (`xsigLegend`), and of the y-vector (`ysigLegend`) as Strings, given
the key of the y-vector (`name`) and the key of the x-vector (`xAxis`).
If `xLabel=false` the legend of the x-vector should be an empty string (`""`).
"""
function resultTimeSeries(result::StringDictAnyResult, name, xLabel::Bool, xAxis)
    # Get x-axis signal
    (xsig, xkeyName, xName) = getSignal(result, string(xAxis))
    if xsig == nothing
        return (nothing, nothing, nothing, nothing)
    end
    if ndims(xsig) != 1
        @warn "ModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector."
        return (nothing, nothing, nothing, nothing)
    end
    xsigLegend = xLabel ? appendUnit(string(xAxis), string(unit(xsig[1]))) : ""
    xsig       = convert.(Float64, ustrip.(xsig))

    # Get y-axis signals
    (ysig, ykeyName, yName) = getSignal(result, string(name))

    if ysig == nothing
        return (nothing, nothing, nothing, nothing)
    end

    if ndims(ysig) == 0
        # ysig is a scalar
        # Construct a constant time series with two points at the first and the last value of the time vector
        ysigLegend = [appendUnit(yName, string(unit(ysig)))]
        xsig = [xsig[1], xsig[end]]
        ysig = [ysig   , ysig     ]
    elseif ndims(ysig) == 1
        # ysig is a vector
        ysigLegend = [appendUnit(yName, string(unit(ysig[1])))]
    elseif ndims(ysig) == 2
        # ysig is a matrix
        nLegend = size(ysig,2)
        ysigLegend = Array{String}(undef, nLegend)

        for i = 1:nLegend
            ysigLegend[i] = appendUnit(ykeyName*"["*string(i)*"]", string(unit(ysig[1,i])))
        end
    else
        @warn "ModiaMath.plot: Variable $ykeyName is not plotted because variables with more as one dimension not yet supported."
        return (nothing, nothing, nothing, nothing)
    end

    ysig = convert.(Float64, ustrip.(ysig))

    return (xsig, xsigLegend, ysig, ysigLegend)
end

sizeToString(value::Any)    = string( size(value) )
sizeToString(value::Number) = ""

"""
    table = resultTable(result)

Return the variables stored in `result` in form of a DataFrames table
(which can then be printed/showed in various forms).

`Base.show(io, result)` is defined to print `resultTable(result)`,
in case result is of type ModiaMath.ResultWithVariables.

# Examples
```julia
import ModiaMath
using  Unitful

t = range(0.0, stop=10.0, length=100)
result = Dict{AbstractString,Any}()
result["time"] = t * u"s";
result["phi"]  = sin.(t)u"rad";

# Print table of the variables that are stored in result
println("result variables = ", ModiaMath.resultTable(result))

# Results in
result variables =
│ Row │ name   │ elType  │ sizeOrValue   │ unit   │
│     │ String │ String  │ String        │ String │
├─────┼────────┼─────────┼───────────────┼────────┤
│ 1   │ phi    │ Float64 │ (100,)        │ rad    │
│ 2   │ time   │ Float64 │ (100,)        │ s      │

```

"""
function resultTable(result::StringDictAnyResult)
    resultTable = DataFrames.DataFrame(name=String[], elType=String[], sizeOrValue=String[], unit=String[])

    for key in sort( collect( keys(result) ) )
        value = result[key]

        # Determine unit as string (if columns have different units, provide unit per column)
        strippedValue = ustrip.(value) # Strip units from value
        tvalue = typeof( strippedValue )
        tsize  = ndims(value) > 0 ? sizeToString( strippedValue ) : string( strippedValue )
        if tvalue <: Number
            tunit = string( unit(value) )
        elseif tvalue <: AbstractVector || (tvalue <: AbstractMatrix && size(value,2) == 1)
            tunit = string( unit( value[1] ) )
        elseif tvalue <: AbstractMatrix
            tunit = string( unit( value[1] ) )
            columnsHaveSameUnit = true
            for i in 2:size(value,2)
                if string( unit( value[1,i] ) ) != tunit
                    columnsHaveSameUnit = false
                    break
                end
            end
            if !columnsHaveSameUnit
                tunit = "[" * tunit
                for i in 2:size(value,2)
                     tunit = tunit * ", " * string( unit( value[1,i] ) )
                end
                tunit = tunit * "]"
            end
        else
            tunit = "???"
        end

        push!(resultTable, [key, string( typeof(strippedValue[1]) ), tsize, tunit] )
    end
    return resultTable
end

"""
    mutable struct ResultWithVariables

Struct that is generated as a result of ModiaMath.simulate!(...),
if the model is a Modia3D.AbstractComponentWithVariables struct.
"""
mutable struct ResultWithVariables
    "Dictionary of time series (elements can be Float64/Int64/Bool Vectors or Matrices)"
    series::StringDictAnyResult

    "Dictionary of variables (elements are <: ModiaMath.AbstractVariable with required fields value and unit)"
    var::Dict{Symbol,Any}

    "String used as optional heading of a plot window"
    resultHeading::String

    ResultWithVariables(series, var, resultHeading="") = new(series, var, resultHeading)
end

resultHeading(result::ResultWithVariables) = result.resultHeading

function resultTimeSeries(result::ResultWithVariables, name, xLabel::Bool, xAxis)
    seriesDict = result.series

    # Get x-axis signal
    (xsig, xkeyName, xName) = getSignal(seriesDict, string(xAxis))
    if xsig == nothing
        return (nothing, nothing, nothing, nothing)
    end
    if ndims(xsig) != 1
        @warn "ModiaMath.plot: argument xAxis (= $xAxis) does not characterize a signal vector."
        return (nothing, nothing, nothing, nothing)
    end
    xsigLegend = xLabel ? appendUnit(xName, result.var[Symbol(xkeyName)].unit) : ""

    # Get y-axis signal(s)
    (ysig, ykeyName, yNameAsString) = getSignal(seriesDict, string(name))
    symbol_ykeyName = Symbol(ykeyName)

    if ysig == nothing
        return (nothing, nothing, nothing, nothing)
    end

    yvar = result.var[symbol_ykeyName]
    if ndims(ysig) == 0
        # ysig is a scalar
        # Construct a constant time series with two points at the first and the last value of the time vector
        ysigLegend = [appendUnit(yNameAsString, yvar.unit)]
        xsig = [xsig[1], xsig[end]]
        ysig = [ysig   , ysig     ]

    elseif ndims(ysig) == 1
        # ysig is a vector
        ysigLegend = [appendUnit(yNameAsString, yvar.unit)]

    else
        # sig has more as one dimension
        ysigLegend = Array{String}(undef, length(yvar.value))
        for i in eachindex(yvar.value)
            ysigLegend[i] = appendUnit(ModiaMath.indexToString(symbol_ykeyName, yvar.value, i), yvar.unit)
        end
    end

    return (xsig, xsigLegend, ysig, ysigLegend)
end

function resultTable(result::ResultWithVariables)
    resultTable = DataFrames.DataFrame(name=String[], elType=String[], sizeOrValue=String[], unit=String[], info=String[])
    series      = result.series
    vars        = result.var

    for key in sort( collect( keys(series) ) )
        value = series[key]
        var   = vars[Symbol(key)]
        tsize = ndims(value) > 0 ? sizeToString( value ) : string( value )
        push!(resultTable, [key, string( typeof(value[1]) ), tsize, var.unit, var.info] )
    end
    return resultTable
end

function Base.show(io::IO, result::ResultWithVariables)
    show(io, resultTable(result), summary=false, splitcols=true)
end

getHeading(result, heading) = heading != "" ? heading : resultHeading(result)


text:
using Test
using MPI

MPI.Init()
try
    comm = MPI.COMM_WORLD
    size = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)

    source = Cint[rank]
    dest = Cint.(rank:(size-1))
    degree = Cint[length(dest)]
    graph_comm = MPI.Dist_graph_create(comm, source, degree, dest)

    send = [rank*rank]
    recv = collect(1:(rank+1))
    fill!(recv,-1)
    MPI.Neighbor_allgather!(send, recv, graph_comm)

    @test sort(recv) == collect(0:rank).^2
catch e
    if isa(e, MPI.FeatureLevelError)
        @test_broken e.min_version <= MPI.Get_version()
    else
        rethrow(e)
    end
end

MPI.Finalize()
@test MPI.Finalized()


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

## types ##

typealias Dims{N} NTuple{N,Int}
typealias DimsInteger{N} NTuple{N,Integer}
typealias Indices{N} NTuple{N,AbstractUnitRange}

"""
    <:(T1, T2)

Subtype operator, equivalent to `issubtype(T1,T2)`.
"""
const (<:) = issubtype

"""
    supertype(T::DataType)

Return the supertype of DataType `T`.

```jldoctest
julia> supertype(Int32)
Signed
```
"""
supertype(T::DataType) = T.super

## generic comparison ##

==(x, y) = x === y

"""
    isequal(x, y)

Similar to `==`, except treats all floating-point `NaN` values as equal to each other, and
treats `-0.0` as unequal to `0.0`. The default implementation of `isequal` calls `==`, so if
you have a type that doesn't have these floating-point subtleties then you probably only
need to define `==`.

`isequal` is the comparison function used by hash tables (`Dict`). `isequal(x,y)` must imply
that `hash(x) == hash(y)`.

This typically means that if you define your own `==` function then you must define a
corresponding `hash` (and vice versa). Collections typically implement `isequal` by calling
`isequal` recursively on all contents.

Scalar types generally do not need to implement `isequal` separate from `==`, unless they
represent floating-point numbers amenable to a more efficient implementation than that
provided as a generic fallback (based on `isnan`, `signbit`, and `==`).
"""
isequal(x, y) = x == y

# TODO: these can be deleted once the deprecations of ==(x::Char, y::Integer) and
# ==(x::Integer, y::Char) are gone and the above returns false anyway
isequal(x::Char, y::Integer) = false
isequal(x::Integer, y::Char) = false

## minimally-invasive changes to test == causing NotComparableError
# export NotComparableError
# =={T}(x::T, y::T) = x === y
# immutable NotComparableError <: Exception end
# const NotComparable = NotComparableError()
# ==(x::ANY, y::ANY) = NotComparable
# !(e::NotComparableError) = throw(e)
# isequal(x, y) = (x == y) === true

## alternative NotComparableError which captures context
# immutable NotComparableError; a; b; end
# ==(x::ANY, y::ANY) = NotComparableError(x, y)

isequal(x::AbstractFloat, y::AbstractFloat) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::Real,          y::AbstractFloat) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::AbstractFloat, y::Real         ) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)

isless(x::AbstractFloat, y::AbstractFloat) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::Real,          y::AbstractFloat) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::AbstractFloat, y::Real         ) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)

function ==(T::Type, S::Type)
    @_pure_meta
    typeseq(T, S)
end
function !=(T::Type, S::Type)
    @_pure_meta
    !(T == S)
end
==(T::TypeVar, S::Type) = false
==(T::Type, S::TypeVar) = false

## comparison fallbacks ##

"""
    !=(x, y)
    ≠(x,y)

Not-equals comparison operator. Always gives the opposite answer as `==`. New types should
generally not implement this, and rely on the fallback definition `!=(x,y) = !(x==y)` instead.
"""
!=(x,y) = !(x==y)
const ≠ = !=

"""
    ===(x,y) -> Bool
    ≡(x,y) -> Bool

Determine whether `x` and `y` are identical, in the sense that no program could distinguish
them. Compares mutable objects by address in memory, and compares immutable objects (such as
numbers) by contents at the bit level. This function is sometimes called `egal`.
"""
===
const ≡ = ===

"""
    !==(x, y)
    ≢(x,y)

Equivalent to `!(x === y)`.
"""
!==(x,y) = !(x===y)
const ≢ = !==

"""
    <(x, y)

Less-than comparison operator. New numeric types should implement this function for two
arguments of the new type. Because of the behavior of floating-point NaN values, `<`
implements a partial order. Types with a canonical partial order should implement `<`, and
types with a canonical total order should implement `isless`.
"""
<(x,y) = isless(x,y)

"""
    >(x, y)

Greater-than comparison operator. Generally, new types should implement `<` instead of this
function, and rely on the fallback definition `>(x,y) = y<x`.
"""
>(x,y) = y < x

"""
    <=(x, y)
    ≤(x,y)

Less-than-or-equals comparison operator.
"""
<=(x,y) = !(y < x)
const ≤ = <=

"""
    >=(x, y)
    ≥(x,y)

Greater-than-or-equals comparison operator.
"""
>=(x,y) = (y <= x)
const ≥ = >=

"""
    .>(x, y)

Element-wise greater-than comparison operator.
"""
.>(x,y) = y .< x

"""
    .>=(x, y)
    .≥(x,y)

Element-wise greater-than-or-equals comparison operator.
"""
.>=(x,y) = y .<= x
const .≥ = .>=

# this definition allows Number types to implement < instead of isless,
# which is more idiomatic:
isless(x::Real, y::Real) = x<y
lexcmp(x::Real, y::Real) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)

"""
    ifelse(condition::Bool, x, y)

Return `x` if `condition` is `true`, otherwise return `y`. This differs from `?` or `if` in
that it is an ordinary function, so all the arguments are evaluated first. In some cases,
using `ifelse` instead of an `if` statement can eliminate the branch in generated code and
provide higher performance in tight loops.
"""
ifelse(c::Bool, x, y) = select_value(c, x, y)

"""
    cmp(x,y)

Return -1, 0, or 1 depending on whether `x` is less than, equal to, or greater than `y`,
respectively. Uses the total order implemented by `isless`. For floating-point numbers, uses `<`
but throws an error for unordered arguments.
"""
cmp(x,y) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)

"""
    lexcmp(x, y)

Compare `x` and `y` lexicographically and return -1, 0, or 1 depending on whether `x` is
less than, equal to, or greater than `y`, respectively. This function should be defined for
lexicographically comparable types, and `lexless` will call `lexcmp` by default.
"""
lexcmp(x,y) = cmp(x,y)

"""
    lexless(x, y)

Determine whether `x` is lexicographically less than `y`.
"""
lexless(x,y) = lexcmp(x,y)<0

# cmp returns -1, 0, +1 indicating ordering
cmp(x::Integer, y::Integer) = ifelse(isless(x,y), -1, ifelse(isless(y,x), 1, 0))

"""
    max(x, y, ...)

Return the maximum of the arguments. See also the [`maximum`](:func:`maximum`) function
to take the maximum element from a collection.
"""
max(x,y) = ifelse(y < x, x, y)

"""
    min(x, y, ...)

Return the minimum of the arguments. See also the [`minimum`](:func:`minimum`) function
to take the minimum element from a collection.
"""
min(x,y) = ifelse(y < x, y, x)

"""
    minmax(x, y)

Return `(min(x,y), max(x,y))`. See also: [`extrema`](:func:`extrema`) that returns `(minimum(x), maximum(x))`.

```jldoctest
julia> minmax('c','b')
('b','c')
```
"""
minmax(x,y) = y < x ? (y, x) : (x, y)

scalarmax(x,y) = max(x,y)
scalarmax(x::AbstractArray, y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmax(x               , y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmax(x::AbstractArray, y               ) = throw(ArgumentError("ordering is not well-defined for arrays"))

scalarmin(x,y) = min(x,y)
scalarmin(x::AbstractArray, y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmin(x               , y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmin(x::AbstractArray, y               ) = throw(ArgumentError("ordering is not well-defined for arrays"))

## definitions providing basic traits of arithmetic operators ##

"""
    identity(x)

The identity function. Returns its argument.
"""
identity(x) = x

+(x::Number) = x
*(x::Number) = x
(&)(x::Integer) = x
(|)(x::Integer) = x
xor(x::Integer) = x

const ⊻ = xor

# foldl for argument lists. expand recursively up to a point, then
# switch to a loop. this allows small cases like `a+b+c+d` to be inlined
# efficiently, without a major slowdown for `+(x...)` when `x` is big.
afoldl(op,a) = a
afoldl(op,a,b) = op(a,b)
afoldl(op,a,b,c...) = afoldl(op, op(a,b), c...)
function afoldl(op,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,qs...)
    y = op(op(op(op(op(op(op(op(op(op(op(op(op(op(op(a,b),c),d),e),f),g),h),i),j),k),l),m),n),o),p)
    for x in qs; y = op(y,x); end
    y
end

for op in (:+, :*, :&, :|, :xor, :min, :max, :kron)
    @eval begin
        # note: these definitions must not cause a dispatch loop when +(a,b) is
        # not defined, and must only try to call 2-argument definitions, so
        # that defining +(a,b) is sufficient for full functionality.
        ($op)(a, b, c, xs...) = afoldl($op, ($op)(($op)(a,b),c), xs...)
        # a further concern is that it's easy for a type like (Int,Int...)
        # to match many definitions, so we need to keep the number of
        # definitions down to avoid losing type information.
    end
end

"""
    \\(x, y)

Left division operator: multiplication of `y` by the inverse of `x` on the left. Gives
floating-point results for integer arguments.
"""
\(x,y) = (y'/x')'

# .<op> defaults to <op>
./(x::Number,y::Number) = x/y
.\(x::Number,y::Number) = y./x
.*(x::Number,y::Number) = x*y
.^(x::Number,y::Number) = x^y
.+(x::Number,y::Number) = x+y
.-(x::Number,y::Number) = x-y
.<<(x::Integer,y::Integer) = x<<y
.>>(x::Integer,y::Integer) = x>>y

.==(x::Number,y::Number) = x == y

"""
    .!=(x, y)
    .≠(x,y)

Element-wise not-equals comparison operator.
"""
.!=(x::Number,y::Number) = x != y
.<( x::Real,y::Real) = x < y

"""
    .<=(x, y)
    .≤(x,y)

Element-wise less-than-or-equals comparison operator.
"""
.<=(x::Real,y::Real) = x <= y

const .≤ = .<=
const .≠ = .!=

# Core <<, >>, and >>> take either Int or UInt as second arg. Signed shift
# counts can shift in either direction, and are translated here to unsigned
# counts. Integer datatypes only need to implement the unsigned version.

"""
    <<(x, n)

Left bit shift operator, `x << n`. For `n >= 0`, the result is `x` shifted left
by `n` bits, filling with `0`s. This is equivalent to `x * 2^n`. For `n < 0`,
this is equivalent to `x >> -n`.

```jldoctest
julia> Int8(3) << 2
12

julia> bits(Int8(3))
"00000011"

julia> bits(Int8(12))
"00001100"
```
See also [`>>`](:func:`>>`), [`>>>`](:func:`>>>`).
"""
function <<(x::Integer, c::Integer)
    typemin(Int) <= c <= typemax(Int) && return x << (c % Int)
    (x >= 0 || c >= 0) && return zero(x)
    oftype(x, -1)
end
<<(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x << (c % UInt) : zero(x)
<<(x::Integer, c::Int) = c >= 0 ? x << unsigned(c) : x >> unsigned(-c)

"""
    >>(x, n)

Right bit shift operator, `x >> n`. For `n >= 0`, the result is `x` shifted
right by `n` bits, where `n >= 0`, filling with `0`s if `x >= 0`, `1`s if `x <
0`, preserving the sign of `x`. This is equivalent to `fld(x, 2^n)`. For `n <
0`, this is equivalent to `x << -n`.

```jldoctest
julia> Int8(13) >> 2
3

julia> bits(Int8(13))
"00001101"

julia> bits(Int8(3))
"00000011"

julia> Int8(-14) >> 2
-4

julia> bits(Int8(-14))
"11110010"

julia> bits(Int8(-4))
"11111100"
```
See also [`>>>`](:func:`>>>`), [`<<`](:func:`<<`).
"""
function >>(x::Integer, c::Integer)
    typemin(Int) <= c <= typemax(Int) && return x >> (c % Int)
    (x >= 0 || c < 0) && return zero(x)
    oftype(x, -1)
end
>>(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x >> (c % UInt) : zero(x)
>>(x::Integer, c::Int) = c >= 0 ? x >> unsigned(c) : x << unsigned(-c)

"""
    >>>(x, n)

Unsigned right bit shift operator, `x >>> n`. For `n >= 0`, the result is `x`
shifted right by `n` bits, where `n >= 0`, filling with `0`s. For `n < 0`, this
is equivalent to `x << -n`.

For `Unsigned` integer types, this is equivalent to [`>>`](:func:`>>`). For
`Signed` integer types, this is equivalent to `signed(unsigned(x) >> n)`.

```jldoctest
julia> Int8(-14) >>> 2
60

julia> bits(Int8(-14))
"11110010"

julia> bits(Int8(60))
"00111100"
```
`BigInt`s are treated as if having infinite size, so no filling is required and this
is equivalent to [`>>`](:func:`>>`).

See also [`>>`](:func:`>>`), [`<<`](:func:`<<`).
"""
>>>(x::Integer, c::Integer) =
    typemin(Int) <= c <= typemax(Int) ? x >>> (c % Int) : zero(x)
>>>(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x >>> (c % UInt) : zero(x)
>>>(x::Integer, c::Int) = c >= 0 ? x >>> unsigned(c) : x << unsigned(-c)

# fallback div, fld, and cld implementations
# NOTE: C89 fmod() and x87 FPREM implicitly provide truncating float division,
# so it is used here as the basis of float div().
div{T<:Real}(x::T, y::T) = convert(T,round((x-rem(x,y))/y))

"""
    fld(x, y)

Largest integer less than or equal to `x/y`.

```jldoctest
julia> fld(7.3,5.5)
1.0
```
"""
fld{T<:Real}(x::T, y::T) = convert(T,round((x-mod(x,y))/y))

"""
    cld(x, y)

Smallest integer larger than or equal to `x/y`.
```jldoctest
julia> cld(5.5,2.2)
3.0
```
"""
cld{T<:Real}(x::T, y::T) = convert(T,round((x-modCeil(x,y))/y))
#rem{T<:Real}(x::T, y::T) = convert(T,x-y*trunc(x/y))
#mod{T<:Real}(x::T, y::T) = convert(T,x-y*floor(x/y))
modCeil{T<:Real}(x::T, y::T) = convert(T,x-y*ceil(x/y))

# operator alias

"""
    rem(x, y)
    %(x, y)

Remainder from Euclidean division, returning a value of the same sign as `x`, and smaller in
magnitude than `y`. This value is always exact.

```julia
x == div(x,y)*y + rem(x,y)
```
"""
rem
const % = rem
.%(x::Real, y::Real) = x%y

"""
    div(x, y)
    ÷(x, y)

The quotient from Euclidean division. Computes `x/y`, truncated to an integer.
"""
div
const ÷ = div
.÷(x::Real, y::Real) = x÷y

"""
    mod1(x, y)

Modulus after flooring division, returning a value `r` such that `mod(r, y) == mod(x, y)`
 in the range ``(0, y]`` for positive `y` and in the range ``[y,0)`` for negative `y`.
"""
mod1{T<:Real}(x::T, y::T) = (m=mod(x,y); ifelse(m==0, y, m))
# efficient version for integers
mod1{T<:Integer}(x::T, y::T) = mod(x+y-T(1),y)+T(1)

"""
    fld1(x, y)

Flooring division, returning a value consistent with `mod1(x,y)`

```julia
x == fld(x,y)*y + mod(x,y)
x == (fld1(x,y)-1)*y + mod1(x,y)
```
"""
fld1{T<:Real}(x::T, y::T) = (m=mod(x,y); fld(x-m,y))
# efficient version for integers
fld1{T<:Integer}(x::T, y::T) = fld(x+y-T(1),y)

"""
    fldmod1(x, y)

Return `(fld1(x,y), mod1(x,y))`.
"""
fldmod1{T<:Real}(x::T, y::T) = (fld1(x,y), mod1(x,y))
# efficient version for integers
fldmod1{T<:Integer}(x::T, y::T) = (fld1(x,y), mod1(x,y))

# transpose

"""
    ctranspose(A)

The conjugate transposition operator (`'`).
"""
ctranspose(x) = conj(transpose(x))
conj(x) = x

# transposed multiply

"""
    Ac_mul_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ⋅B``.
"""
Ac_mul_B(a,b)  = ctranspose(a)*b

"""
    A_mul_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A⋅Bᴴ``.
"""
A_mul_Bc(a,b)  = a*ctranspose(b)

"""
    Ac_mul_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ Bᴴ``.
"""
Ac_mul_Bc(a,b) = ctranspose(a)*ctranspose(b)

"""
    At_mul_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ⋅B``.
"""
At_mul_B(a,b)  = transpose(a)*b

"""
    A_mul_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A⋅Bᵀ``.
"""
A_mul_Bt(a,b)  = a*transpose(b)

"""
    At_mul_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ⋅Bᵀ``.
"""
At_mul_Bt(a,b) = transpose(a)*transpose(b)

# transposed divide

"""
    Ac_rdiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ / B``.
"""
Ac_rdiv_B(a,b)  = ctranspose(a)/b

"""
    A_rdiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A / Bᴴ``.
"""
A_rdiv_Bc(a,b)  = a/ctranspose(b)

"""
    Ac_rdiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ / Bᴴ``.
"""
Ac_rdiv_Bc(a,b) = ctranspose(a)/ctranspose(b)

"""
    At_rdiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ / B``.
"""
At_rdiv_B(a,b)  = transpose(a)/b

"""
    A_rdiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A / Bᵀ``.
"""
A_rdiv_Bt(a,b)  = a/transpose(b)

"""
    At_rdiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ / Bᵀ``.
"""
At_rdiv_Bt(a,b) = transpose(a)/transpose(b)

"""
    Ac_ldiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``B``.
"""
Ac_ldiv_B(a,b)  = ctranspose(a)\b

"""
    A_ldiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A`` \\ ``Bᴴ``.
"""
A_ldiv_Bc(a,b)  = a\ctranspose(b)

"""
    Ac_ldiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``Bᴴ``.
"""
Ac_ldiv_Bc(a,b) = ctranspose(a)\ctranspose(b)

"""
    At_ldiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ`` \\ ``B``.
"""
At_ldiv_B(a,b)  = transpose(a)\b

"""
    A_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A`` \\ ``Bᵀ``.
"""
A_ldiv_Bt(a,b)  = a\transpose(b)

"""
    At_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ`` \\ ``Bᵀ``.
"""
At_ldiv_Bt(a,b) = At_ldiv_B(a,transpose(b))

"""
    Ac_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``Bᵀ``.
"""
Ac_ldiv_Bt(a,b) = Ac_ldiv_B(a,transpose(b))

widen{T<:Number}(x::T) = convert(widen(T), x)

"""
    eltype(type)

Determine the type of the elements generated by iterating a collection of the given `type`.
For associative collection types, this will be a `Pair{KeyType,ValType}`. The definition
`eltype(x) = eltype(typeof(x))` is provided for convenience so that instances can be passed
instead of types. However the form that accepts a type argument should be defined for new
types.
"""
eltype(::Type) = Any
eltype(::Type{Any}) = Any
eltype(t::DataType) = eltype(supertype(t))
eltype(x) = eltype(typeof(x))

# function pipelining

"""
    |>(x, f)

Applies a function to the preceding argument. This allows for easy function chaining.

```jldoctest
julia> [1:5;] |> x->x.^2 |> sum |> inv
0.01818181818181818
```
"""
|>(x, f) = f(x)

# array shape rules

promote_shape(::Tuple{}, ::Tuple{}) = ()

function promote_shape(a::Tuple{Int,}, b::Tuple{Int,})
    if a[1] != b[1]
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

function promote_shape(a::Tuple{Int,Int}, b::Tuple{Int,})
    if a[1] != b[1] || a[2] != 1
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

promote_shape(a::Tuple{Int,}, b::Tuple{Int,Int}) = promote_shape(b, a)

function promote_shape(a::Tuple{Int, Int}, b::Tuple{Int, Int})
    if a[1] != b[1] || a[2] != b[2]
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

"""
    promote_shape(s1, s2)

Check two array shapes for compatibility, allowing trailing singleton dimensions, and return
whichever shape has more dimensions.

```jldoctest
julia> a = ones(3,4,1,1,1);

julia> b = ones(3,4);

julia> promote_shape(a,b)
(Base.OneTo(3),Base.OneTo(4),Base.OneTo(1),Base.OneTo(1),Base.OneTo(1))

julia> promote_shape((2,3,1,4), (2,3,1,4,1))
(2,3,1,4,1)
```
"""
function promote_shape(a::Dims, b::Dims)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end

function promote_shape(a::AbstractArray, b::AbstractArray)
    promote_shape(indices(a), indices(b))
end

function promote_shape(a::Indices, b::Indices)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1:1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end

function throw_setindex_mismatch(X, I)
    if length(I) == 1
        throw(DimensionMismatch("tried to assign $(length(X)) elements to $(I[1]) destinations"))
    else
        throw(DimensionMismatch("tried to assign $(dims2string(size(X))) array to $(dims2string(I)) destination"))
    end
end

# check for valid sizes in A[I...] = X where X <: AbstractArray
# we want to allow dimensions that are equal up to permutation, but only
# for permutations that leave array elements in the same linear order.
# those are the permutations that preserve the order of the non-singleton
# dimensions.
function setindex_shape_check(X::AbstractArray, I::Integer...)
    li = ndims(X)
    lj = length(I)
    i = j = 1
    while true
        ii = length(indices(X,i))
        jj = I[j]
        if i == li || j == lj
            while i < li
                i += 1
                ii *= length(indices(X,i))
            end
            while j < lj
                j += 1
                jj *= I[j]
            end
            if ii != jj
                throw_setindex_mismatch(X, I)
            end
            return
        end
        if ii == jj
            i += 1
            j += 1
        elseif ii == 1
            i += 1
        elseif jj == 1
            j += 1
        else
            throw_setindex_mismatch(X, I)
        end
    end
end

setindex_shape_check(X::AbstractArray) =
    (_length(X)==1 || throw_setindex_mismatch(X,()))

setindex_shape_check(X::AbstractArray, i::Integer) =
    (_length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer) =
    (_length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer, j::Integer) =
    (_length(X)==i*j || throw_setindex_mismatch(X, (i,j)))

function setindex_shape_check{T}(X::AbstractArray{T,2}, i::Integer, j::Integer)
    if length(X) != i*j
        throw_setindex_mismatch(X, (i,j))
    end
    sx1 = length(indices(X,1))
    if !(i == 1 || i == sx1 || sx1 == 1)
        throw_setindex_mismatch(X, (i,j))
    end
end
setindex_shape_check(X, I...) = nothing # Non-arrays broadcast to all idxs

# convert to a supported index type (Array, Colon, or Int)
to_index(i::Int) = i
to_index(i::Integer) = convert(Int,i)::Int
to_index(c::Colon) = c
to_index(I::AbstractArray{Bool}) = find(I)
to_index(A::AbstractArray) = A
to_index{T<:AbstractArray}(A::AbstractArray{T}) = throw(ArgumentError("invalid index: $A"))
to_index(A::AbstractArray{Colon}) = throw(ArgumentError("invalid index: $A"))
to_index(i) = throw(ArgumentError("invalid index: $i"))

to_indexes() = ()
to_indexes(i1) = (to_index(i1),)
to_indexes(i1, I...) = (@_inline_meta; (to_index(i1), to_indexes(I...)...))

# Addition/subtraction of ranges
for f in (:+, :-)
    @eval begin
        function $f(r1::OrdinalRange, r2::OrdinalRange)
            r1l = length(r1)
            (r1l == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            range($f(first(r1),first(r2)), $f(step(r1),step(r2)), r1l)
        end

        function $f{T<:AbstractFloat}(r1::FloatRange{T}, r2::FloatRange{T})
            len = r1.len
            (len == r2.len ||
             throw(DimensionMismatch("argument dimensions must match")))
            divisor1, divisor2 = r1.divisor, r2.divisor
            if divisor1 == divisor2
                FloatRange{T}($f(r1.start,r2.start), $f(r1.step,r2.step),
                              len, divisor1)
            else
                d1 = Int(divisor1)
                d2 = Int(divisor2)
                d = lcm(d1,d2)
                s1 = div(d,d1)
                s2 = div(d,d2)
                FloatRange{T}($f(r1.start*s1, r2.start*s2),
                              $f(r1.step*s1, r2.step*s2),  len, d)
            end
        end

        function $f{T<:AbstractFloat}(r1::LinSpace{T}, r2::LinSpace{T})
            len = r1.len
            (len == r2.len ||
             throw(DimensionMismatch("argument dimensions must match")))
            divisor1, divisor2 = r1.divisor, r2.divisor
            if divisor1 == divisor2
                LinSpace{T}($f(r1.start, r2.start), $f(r1.stop, r2.stop),
                            len, divisor1)
            else
                linspace(convert(T, $f(first(r1), first(r2))),
                         convert(T, $f(last(r1), last(r2))), len)
            end
        end

        $f(r1::Union{FloatRange, OrdinalRange, LinSpace},
           r2::Union{FloatRange, OrdinalRange, LinSpace}) =
               $f(promote(r1, r2)...)
    end
end

# vectorized ifelse

function ifelse(c::AbstractArray{Bool}, x, y)
    [ifelse(ci, x, y) for ci in c]
end

function ifelse(c::AbstractArray{Bool}, x::AbstractArray, y::AbstractArray)
    [ifelse(c_elem, x_elem, y_elem) for (c_elem, x_elem, y_elem) in zip(c, x, y)]
end

function ifelse(c::AbstractArray{Bool}, x::AbstractArray, y)
    [ifelse(c_elem, x_elem, y) for (c_elem, x_elem) in zip(c, x)]
end

function ifelse(c::AbstractArray{Bool}, x, y::AbstractArray)
    [ifelse(c_elem, x, y_elem) for (c_elem, y_elem) in zip(c, y)]
end

# Pair

immutable Pair{A,B}
    first::A
    second::B
end

const => = Pair

start(p::Pair) = 1
done(p::Pair, i) = i>2
next(p::Pair, i) = (getfield(p,i), i+1)

indexed_next(p::Pair, i::Int, state) = (getfield(p,i), i+1)

hash(p::Pair, h::UInt) = hash(p.second, hash(p.first, h))

==(p::Pair, q::Pair) = (p.first==q.first) & (p.second==q.second)
isequal(p::Pair, q::Pair) = isequal(p.first,q.first) & isequal(p.second,q.second)

isless(p::Pair, q::Pair) = ifelse(!isequal(p.first,q.first), isless(p.first,q.first),
                                                             isless(p.second,q.second))
getindex(p::Pair,i::Int) = getfield(p,i)
getindex(p::Pair,i::Real) = getfield(p, convert(Int, i))
reverse{A,B}(p::Pair{A,B}) = Pair{B,A}(p.second, p.first)

endof(p::Pair) = 2
length(p::Pair) = 2

convert{A,B}(::Type{Pair{A,B}}, x::Pair{A,B}) = x
function convert{A,B}(::Type{Pair{A,B}}, x::Pair)
    Pair{A, B}(convert(A, x[1]), convert(B, x[2]))
end

promote_rule{A1, B1, A2, B2}(::Type{Pair{A1, B1}}, ::Type{Pair{A2, B2}}) =
    Pair{promote_type(A1, A2), promote_type(B1, B2)}

# some operators not defined yet
global //, >:, <|, hcat, hvcat, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛

this_module = current_module()
baremodule Operators

export
    !,
    !=,
    !==,
    ===,
    xor,
    %,
    .%,
    ÷,
    .÷,
    &,
    *,
    +,
    -,
    .!=,
    .+,
    .-,
    .*,
    ./,
    .<,
    .<=,
    .==,
    .>,
    .>=,
    .\,
    .^,
    /,
    //,
    <,
    <:,
    >:,
    <<,
    <=,
    ==,
    >,
    >=,
    ≥,
    ≤,
    ≠,
    .≥,
    .≤,
    .≠,
    >>,
    .>>,
    .<<,
    >>>,
    \,
    ^,
    |,
    |>,
    <|,
    ~,
    ⋅,
    ×,
    ∈,
    ∉,
    ∋,
    ∌,
    ⊆,
    ⊈,
    ⊊,
    ∩,
    ∪,
    √,
    ∛,
    ⊻,
    colon,
    hcat,
    vcat,
    hvcat,
    getindex,
    setindex!,
    transpose,
    ctranspose

import ..this_module: !, !=, xor, %, .%, ÷, .÷, &, *, +, -, .!=, .+, .-, .*, ./, .<, .<=, .==, .>,
    .>=, .\, .^, /, //, <, <:, <<, <=, ==, >, >=, >>, .>>, .<<, >>>,
    <|, |>, \, ^, |, ~, !==, ===, >:, colon, hcat, vcat, hvcat, getindex, setindex!,
    transpose, ctranspose,
    ≥, ≤, ≠, .≥, .≤, .≠, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛, ⊻

end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia