text:
# blocking analysis of single time series
# `blocking_analysis()`
# `mean_and_se()`

"""
    BlockingResult(mean, err, err_err, p_cov, k, blocks)
Result of [`blocking_analysis()`](@ref).

### Fields:
- `mean`: sample mean
- `err`: standard error (estimated standard deviation of the mean)
- `err_err`: estimated uncertainty of `err`
- `p_cov`: estimated pseudo covariance of `mean`, relevant for complex time series
- `k::Int`: k-1 blocking steps were used to uncorrelate time series
- `blocks::Int`: number of uncorrelated values after blocking

Has methods for [`NamedTuple`](@ref), [`val_and_errs`](@ref), [`val`](@ref), [`errs`](@ref),
[`mean_and_se`](@ref), `Measurements.:±`,
`MonteCarloMeasurements.Particles`, and
`Statistics.cov` for `Complex` data.

**Example:**
```jldoctest; setup = :(Random.seed!(1234))
julia> blocking_analysis(smoothen(randn(2^10), 2^5))
BlockingResult{Float64}
  mean = -0.025 ± 0.025
  with uncertainty of ± 0.00311966837382259
  from 32 blocks after 5 transformations (k = 6).
```
"""
struct BlockingResult{T}
    mean::T
    err::Float64
    err_err::Float64
    p_cov::T # pseudo covariance for complex normal distribution
    k::Int
    blocks::Int
end
# constructor from NamedTuple of vectors
function BlockingResult(nt, k)
    T = eltype(nt.mean)
    if k < 0 # blocking failed
        return BlockingResult{T}(nt.mean[1], NaN, NaN, T(NaN), k, 0)
    end
    return BlockingResult{T}(
        nt.mean[1], nt.std_err[k], nt.std_err_err[k], nt.p_cov[k], k, nt.blocks[k]
    )
end

function Base.show(io::IO, r::BlockingResult{T}) where T
    println(io, "BlockingResult{$T}")
    println(io, "  mean = $(Measurements.measurement(r))")
    println(io, "  with uncertainty of ± $(r.err_err)")
    T<:Complex && println(io, "  cov = $(cov(r))")
    if r.k > 0
        println(io,"  from $(r.blocks) blocks after $(r.k-1) transformations (k = $(r.k)).")
    else
        println(io, "  Blocking unsuccessful. k = $(r.k). Try using more time steps!")
    end
end

"""
    cov(r::BlockingResult{<:Complex})
Return the covariance matrix of the multivariate normal distribution approximating
the uncertainty of the blocking result `r` of a complex time series.
See (https://en.wikipedia.org/wiki/Complex_normal_distribution).
"""
function Statistics.cov(r::BlockingResult{<:Complex})
    v_xx = real(r.err^2 + r.p_cov)/2
    v_xy = imag(r.p_cov)/2
    v_yy = real(r.err^2 - r.p_cov)/2
    return [v_xx v_xy; v_xy v_yy]
end

"""
    measurement(r::BlockingResult)
    Measurements.±(r::BlockingResult)
Convert a `BlockingResult` into a `Measurement` for linear error propagation with
[`Measurements`](@ref).

Limitation: Does not account for covariance in complex `BlockingResult`.
Consider using [`MonteCarloMeasurements.Particles(r)`](@ref)!
"""
function Measurements.measurement(r::BlockingResult{<:Real})
    return measurement(r.mean, r.err)
end
function Measurements.measurement(r::BlockingResult{<:Complex})
    Σ = cov(r) # real valued covariance matrix
    cm = complex(measurement(real(r.mean), √Σ[1,1]), measurement(imag(r.mean), √Σ[2,2]))
    return cm
end

"""
    MonteCarloMeasurements.Particles(r::BlockingResult; mc_samples = 2000)
    MonteCarloMeasurements.±(r::BlockingResult)
Convert a `BlockingResult` into a `Particles` object for nonlinear error propagation with
[`MonteCarloMeasurements`](@ref).
"""
function MonteCarloMeasurements.Particles(r::BlockingResult{<:Real}; mc_samples = 2000)
    return Particles(mc_samples, Normal(r.mean, r.err))
end
function MonteCarloMeasurements.Particles(r::BlockingResult{<:Complex}; mc_samples = 2000)
    Σ = cov(r) # real valued covariance matrix
    ps = Particles(mc_samples, MvNormal([real(r.mean), imag(r.mean)], Σ))
    return complex(ps[1], ps[2])
end
MonteCarloMeasurements.:±(r::BlockingResult) = Particles(r)

"""
    blocking_analysis(v::AbstractVector; α = 0.01, corrected = true, skip=0)
    -> BlockingResult(mean, err, err_err, p_cov, k, blocks)
Compute the sample mean `mean` and estimate the standard deviation of the mean
(standard error) `err` of a correlated time series using the blocking algorithm from
Flyvberg and Peterson [JCP (1989)](http://aip.scitation.org/doi/10.1063/1.457480)
and the M test of Jonsson
[PRE (2018)](https://link.aps.org/doi/10.1103/PhysRevE.98.043304) at significance level
``1-α``. `k` is the number of
blocking transformations required to pass the hypothesis test for an uncorrelated time
series and `err_err` the estimated standard error or `err`. Use `skip` to skip the first
`skip` elements in `v`.

If decorrelating the
time series fails according to the M test, `NaN` is returned as the standard error and `-1`
for `k`.
`corrected` controls whether
bias correction for variances is used.

See [`BlockingResult`](@ref), [`shift_estimator`](@ref), [`ratio_of_means`](ref).
"""
function blocking_analysis(v::AbstractVector; α = 0.01, corrected::Bool=true, skip=0)
    T = float(eltype(v))
    v = @view v[skip+1:end]
    if length(v) == 0
        @error "Attempted blocking on an empty vector"
        return BlockingResult(zero(T), NaN, NaN, T(NaN), -1, 0)
    elseif length(v) == 1 # treat like failed M test
        return BlockingResult(T(v[1]), NaN, NaN, T(NaN), -1, 1)
    end
    nt = blocks_with_m(v; corrected)
    k = mtest(nt.mj; α, warn=false)
    return BlockingResult(nt, k)
end

"""
    mean_and_se(v::AbstractVector; α = 0.01, corrected::Bool=true, skip=0) -> mean, err
    mean_and_se(r::BlockingResult) -> mean, err
Return the mean and standard error (as a tuple) of a time series obtained from
[`blocking_analysis`](@ref). See also [`BlockingResult`](@ref).
"""
mean_and_se(r::BlockingResult) = r.mean, r.err
function mean_and_se(v::AbstractVector; kwargs...)
    return mean_and_se(blocking_analysis(v; kwargs...))
end

"""
    blocker(v::Vector) -> new_v::Vector
Reblock the data by successively taking the mean of two adjacent data points to
form a new vector with a half of the `length(v)`. The last data point will be
discarded if `length(v)` is odd.
"""
function blocker(v::AbstractVector{T}) where T
    P = typeof(zero(T)/1)
    new_v = Array{P}(undef,(length(v)÷2))
    for i  in 1:length(v)÷2
        @inbounds new_v[i] = (v[2i-1]+v[2i])/2
    end
    return new_v
end

"""
    mtest(mj::AbstractVector; α = 0.01, warn = true) -> k
    mtest(table; α = 0.01, warn = true) -> k
Hypothesis test for decorrelation of a time series after blocking transformations
with significance level ``1-α`` after Jonson
[PRE (2018)](https://doi.org/10.1103/PhysRevE.98.043304).
`mj` or `table.mj` is expected to be a vector with relevant ``M_j`` values from a blocking analysis as
obtained from [`blocks_with_m()`](@ref).
Returns the row number `k` where the M-test is passed.
If the M-test has failed `mtest()` returns the value `-1` and optionally prints
a warning message.
"""
function mtest(mj::AbstractVector; α = 0.01, warn = true)
    m = reverse(cumsum(reverse(mj)))
    k = 1
    while k <= length(m)-1
       if m[k] < cquantile(Chisq(k), α) # compare to χ^2 distribution at quantile 1-α
           return k
       else
           k += 1
       end
    end
    if warn
        @warn "M test failed, more data needed"
    end
    return -1 # indicating the the M-test has failed
end
mtest(table) = mtest(table.mj)

"""
    blocks_with_m(v; corrected = true) -> (;blocks, mean, std_err, std_err_err, p_cov, mj)
Perform the blocking algorithm from Flyvberg and Peterson
[JCP (1989)](http://aip.scitation.org/doi/10.1063/1.457480).
Returns named tuple with the results from all blocking steps. See [`mtest()`](@ref).
"""
@inline function blocks_with_m(v; corrected::Bool=true)
    T = float(eltype(v))
    R = real(T)
    n_steps = floor(Int,log2(length(v)))

    # initialise arrays to be returned
    blocks = Vector{Int}(undef,n_steps)
    mean_arr = Vector{T}(undef,n_steps)
    std_err = Vector{R}(undef,n_steps)
    std_err_err = Vector{R}(undef,n_steps)
    p_cov = Vector{T}(undef,n_steps)
    mj = Vector{R}(undef,n_steps)

    for i in 1:n_steps
        n = length(v)
        blocks[i] = n
        mean_v = mean(v)
        mean_arr[i] = mean_v
        variance = var(v; corrected, mean=mean_v) # variance
        # sample covariance ŷ(1) Eq. (6) [Jonsson]
        gamma = real(autocovariance(v,1; corrected, mean=mean_v))
        # the M value Eq. (12) [Jonsson]
        mj[i] = n*((n-1)*variance/(n^2)+gamma)^2/(variance^2)
        stderr_v = sqrt(variance/n) # standard error
        std_err[i] = stderr_v
        std_err_err[i] = stderr_v/sqrt(2*(n-1)) # error on standard error Eq. (28) [F&P]
        p_cov[i] = pseudo_cov(v,v; xmean=mean_v, ymean=mean_v, corrected)/n
        v = blocker(v) # re-blocking the dataset
    end
    (length(v)≤ 0 || length(v)>2) && @error "Something went wrong in `blocks_with_m`."
    return (;blocks, mean=mean_arr, std_err, std_err_err, p_cov, mj)
end


text:
using LinearAlgebra: dot
using Plots
function init_position(c::Potentials.Configuration, dim)
    #=
    Initialise the positions for a series
    of particles.

    Parameters
    ----------
    N: int
        Number of particles in the
        system
    dim: int
        Number of dimensions

    Returns
    -------
    ndarray of floats
        Initial velocities for a series of
        particles (eVs/Åamu)
    =#
    N = c.num_atoms
    if dim == 2
        N = N-1
        r0 = 2 .*[cos.(2*pi/N .* (1:N))+1e-2*randn(N) sin.(2*pi/N.*(1:N))+1e-2*randn(N) 0.0 .*(1:N)]
        r0 = [r0; [0 0 0]]
    elseif dim == 3
        if N <= 7
            r0 = [-0.9523699364        0.0159052548       -0.0840802250;
            -0.2528949024        0.8598484003       -0.3332183372;
            -0.3357005052       -0.8500181306        0.2812584200;
            0.7960710440        0.5155096551       -0.1218613860;
            0.7448954577       -0.5412454414        0.2579017563;
            -0.0442130155        0.1953980103        0.5377653956;
            0.0442118576       -0.1953977484       -0.5377656238]
        elseif N == 38
            r0 = [0.163946667699999993 0.319538517000000022 1.71226674359999986;
            1.18570812759999988 -1.12890812809999996 -0.616849616600000039;
            1.47183871269999988 -0.0530362826999999995 0.944308599800000037;
            -0.156601128000000006 1.44806844500000009 -0.969234048100000023;
            1.43644861639999988 0.381943172100000006 -0.92280851779999995;
            -0.659303829199999969 -1.58099572070000005 -0.355825472500000017;
            -1.43644861639999988 -0.381943172100000006 0.92280851779999995;
            0.659303829199999969 1.58099572070000005 0.355825472600000026;
            -1.18570812749999988 1.12890812809999996 0.616849616600000039;
            0.156601128099999987 -1.44806844500000009 0.969234048100000023;
            -1.47183871259999988 0.0530362826999999995 -0.944308599700000029;
            -0.692522445699999989 0.622174532900000021 -1.48126712230000002;
            -0.903125893300000038 0.439946115799999982 1.43235882759999988;
            -0.943262934499999983 -0.888676767300000003 -1.17530822099999988;
            -0.163946667699999993 -0.319538517000000022 -1.71226674359999986;
            0.402594671899999978 -1.70081948270000005 -0.0772747840999999935;
            -1.71903085060000005 0.307018840199999976 0.107299306900000002;
            0.905297373100000025 1.32824468299999987 -0.690683359700000055;
            -0.402594671899999978 1.70081948270000005 0.0772747840999999935;
            1.71903085060000005 -0.307018840199999976 -0.107299306799999994;
            -0.905297373000000016 -1.32824468299999987 0.690683359700000055;
            0.943262934599999991 0.888676767300000003 1.17530822109999988;
            0.692522445799999997 -0.622174532900000021 1.48126712230000002;
            0.903125893400000046 -0.439946115799999982 -1.43235882749999988;
            -1.18439743100000006 -0.634064735699999993 -0.125267553799999987;
            0.120331581600000001 -1.00573849309999996 -0.891368452299999947;
            0.934263338099999952 -0.873132691499999947 0.430486517999999985;
            1.18439743110000006 0.634064735699999993 0.125267553799999987;
            -0.934263338000000054 0.873132691499999947 -0.430486517899999976;
            -0.370465674600000017 -0.501458934099999976 1.19658741639999988;
            0.370465674700000025 0.501458934099999976 -1.19658741639999988;
            -0.120331581600000001 1.00573849309999996 0.891368452299999947;
            0.649793600800000037 -0.185104513500000012 -0.381540884400000002;
            -0.405362064899999985 -0.0660416075000000019 -0.658322833699999999;
            0.124574169400000004 0.750628853799999995 -0.15200806280000001;
            -0.649793600800000037 0.185104513500000012 0.381540884400000002;
            0.405362064999999994 0.0660416075000000019 0.658322833699999999;
            -0.124574169299999996 -0.750628853799999995 0.15200806280000001]
        else
            r0 = rand(N, d)
        end

        r0 += 1e-2*randn(size(r0))
        com = mean(r0, dims = 1)
        r0 = r0 .- com
    else
        return false
    end

    for j = 1:c.num_atoms
        c.Positions[j] = Potentials.Position( r0[j, 1], r0[j, 2], r0[j, 3], c.Positions[j].type)
    end

    return c
end

function init_velocity(c::Potentials.Configuration, dim::Int, T::Float64)
    #=
    Initialise the velocities for a series
    of particles.

    Parameters
    ----------
    T: float
        Temperature of the system at
        initialisation (K)
    number_of_particles: int
        Number of particles in the
        system

    Returns
    -------
    ndarray of floats
        Initial velocities for a series of
        particles (eVs/Åamu)
    =#
    N = c.num_atoms
    r0 = randn(N, 3) .- 0.5
    if dim == 2
        r0[:, 3] .= 0.0
    end
    r0 = r0 * T
    for j = 1:c.num_atoms
        c.Velocities[j] = Potentials.Position( r0[j, 1], r0[j, 2], r0[j, 3], c.Positions[j].type)
    end
    return c
end

function update_velo(v::Vector, a::Vector, a1::Vector, dt::Float64)
    #=
    Update the particle velocities.

    Parameters
    ----------
    v: ndarray of floats
        The velocities of the particles in a
        single dimension (eVs/Åamu)
    a: ndarray of floats
        The accelerations of the particles in a
        single dimension at the previous
        timestep (eV/Åamu)
    a1: ndarray of floats
        The accelerations of the particles in a
        single dimension at the current
        timestep (eV/Åamu)
    dt: float
        The timestep length (s)

    Returns
    -------
    ndarray of floats:
        New velocities of the particles in a
        single dimension (eVs/Åamu)
    =#
    return v + 0.5 * (a + a1) * dt
end

function verlet(p::Potentials.ArbitraryPotential, c::Potentials.Configuration, dt::Float64, Temp:: Float64; γ = 1.0, k_b = 1.0/119.8)
   #=
    Update the particle positions and velocities by solving.
    dx/dt = dv/dt
    dv/dt = F(r)/m - γ * v + R/m = a

    Parameters
    ----------
    x : ndarray of floats
        The positions of the particles (Å)
    v: ndarray of floats
        The velocities of the particles
        single dimension (eVs/Åamu)
    dt: float
        The timestep length (s)

    Returns
    -------
    ndarray of floats:
        New velocities and positions of the particles in a
        single dimension (Å, eVs/Åamu)
    =#
    pos = vec(c.Positions)
    vels = vec(c.Velocities)
    R = sqrt(2*Temp*k_b*dt*γ)
    a = Potentials.force(c, p)
    a = a - γ .* vels
    noise = 0.5*R*randn(c.num_atoms, 3)
    if pos[1][3] == 0.0
        for j = 1:length(a)
            a[j][3] = 0.0
            noise[j, 3] = 0.0
        end
    end
    v_star = update_velo(vels, a, 0.0*a, dt)
    noise = noise .- mean( noise, dims = 1 )
    noise = [[noise[i, j] for j = 1:3] for i = 1:c.num_atoms]
    v_star += 0.5*noise
    # Second half-time step
    x_new = c.Positions + v_star*dt
    # println("Old Positions ")
    # show(stdout, "text/plain", c.Positions)
    # println(" ")
    # println("New positions ")
    # show(stdout, "text/plain", x_new)
    # println(" ")
    for (i, xi) in enumerate(x_new)
        c.Positions[i] = xi
    end
    a_new = Potentials.force(c, p)
    a_new = a_new - γ .* vels
    noise = 0.5*R*randn(c.num_atoms, 3)
    if pos[1][3] == 0.0
        for j = 1:length(a)
            a_new[j][3] = 0.0
            noise[j, 3] = 0.0
        end
    end
    noise = noise .- mean( noise, dims = 1 )
    noise = [[noise[i, j] for j = 1:3] for i = 1:c.num_atoms]
    v_new = update_velo(v_star, a, a_new, dt)
    v_new += 0.5*noise

    # Save
    for (i, (xi, vi)) in enumerate(zip(x_new, v_new))
        c.Positions[i] = xi
        c.Velocities[i] = Potentials.Position(vi[1], vi[2], vi[3], c.Positions[i].type)
    end
    return c
end

function solve(p::Potentials.ArbitraryPotential, c0::Potentials.Configuration, dt::Float64, Nt::Int, Temp ::Array{Float64}; save_dt = 1 :: Int)
    N = Int(Nt / save_dt)
    r = Vector{Potentials.Configuration}(undef, N)
    r = fill!(r, c0)

    t = zeros( (N, 1));

    τ = 0.0
    count = 1
    c_temp = deepcopy(c0)
    for j = 2:(Nt)
        τ +=  dt
        # println("Time ", τ)
        c_temp = verlet(p, c_temp, dt, Temp[j-1])
        if (j-1) % save_dt == 0
            count += 1
            r[count] = deepcopy(c_temp)
            t[count] = τ
        end
    end

    return r, t
end

function get_positions(r::Vector{PotentialUQ.Potentials.Configuration})
    N = length(r)
    num_atoms = r[1].num_atoms

    pos = Vector{Vector{Float64}}(undef, N)
    for (i, ri) in enumerate(r)
        pos[i] = PotentialUQ.Potentials.norm.(ri.Positions)
    end
    return pos
end

function get_velocities(r::Vector{Potentials.Configuration})
    N = length(r)
    num_atoms = r[1].num_atoms

    pos = Vector{Vector{Float64}}(undef, N)
    for (i, ri) in enumerate(r)
        pos[i] = Potentials.norm.(ri.Velocities)
    end
    return pos
end

function calculate_rdf(r::Vector{Potentials.Configuration}; L = 2.5, maxbin = 1000)

    bins = vcat(range(0, L; length = maxbin+1)...)
    rdf = zeros(maxbin)
    bin_centers = 0.5 .* (bins[1:end-1] + bins[2:end])
    n = length(r)
    num_atoms = r[1].num_atoms
    dr = bin_centers[2] - bin_centers[1]
    hist = zeros(maxbin)
    for c ∈ r
        distances = Potentials.get_interparticle_distance(c)
        for d in distances
            bin = ceil(Int, d / dr)
            if (bin > 1) && (bin < maxbin)
                hist[bin] += 1
            end
        end
    end

    volume = 4 / 3 * π * L^3
    density = num_atoms / volume
    rdf = zeros(maxbin)
    for bin = 1:maxbin
        rlower = bins[bin]
        rupper = bins[bin+1]
        volume_of_shell = 4/3 * π * (rupper^3 - rlower^3)
        rdf[bin] = (hist[bin] / n) / volume_of_shell / density
    end

    return (bin_centers, rdf)
end

function calculate_mean_squared_displacement(r::Vector{Potentials.Configuration})
    n = length(r)
    num_atoms = r[1].num_atoms
    p0 = vec(r[1].Positions)
    disp = zeros(n)
    for i = 1:n
        for j = 1:num_atoms
            p = vec(r[i].Positions)
            dr = p - p0
            disp[i] += dot(dr, dr) / num_atoms
        end
    end
    return disp
end

function calculate_mass_diffusivity(r::Vector{Potentials.Configuration}, t)
    disp = calculate_mean_squared_displacement(r)
    D = sum( disp[end-100:end] ./ (6 .* t[end-100:end] ) ) / 100
    return D
end

function animate_atoms(r::Vector{Potentials.Configuration}, t; dim = 2, fps = 60, dT = 20, title = "lennard_jones_cluster.gif")
    n = length(r)
    num_atoms = r[1].num_atoms

    if dim == 3

        anim = @animate for i = 1:dT:n
            pos = hcat(vec(r[i].Positions)...)
            scatter(pos[1, :], pos[2, :], pos[3, :], xaxis = ("x", (-2, 2), -2:0.5:2),
                            yaxis = ("y", (-2, 2), -2:0.5:2),
                            zaxis = ("z", (-2, 2), -2:0.5:2),
                            markersize = 20,
                            c=colormap("Blues",num_atoms),
                            legend = false,
                            title = string("τ = ", floor(t[i]) ))
        end

    elseif dim==2
        anim = @animate for i = 1:dT:n
            pos = hcat(vec(r[i].Positions)...)
            scatter(pos[1, :], pos[2, :], xaxis = ("x", (-4, 4), -4:0.5:4),
                            yaxis = ("y", (-2, 2), -4:0.5:4),
                            c=colormap("Blues",num_atoms),
                            legend = false,
                            title = string("τ = ", floor(t[i]) ) )

        end
    end

    gif(anim, title, fps = fps)
end

ext: jl


text:
using GLFW, ModernGL

# whenever the window size changed (by OS or user resize) this callback function executes
function framebuffer_size_callback(window::GLFW.Window, width::Cint, height::Cint)
    # make sure the viewport matches the new window dimensions; note that width and
    # height will be significantly larger than specified on retina displays.
	glViewport(0, 0, width, height)
end

# process all input: query GLFW whether relevant keys are pressed/released this frame and react accordingly
function processInput(window::GLFW.Window)
    GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# misc. config
const SCR_WIDTH = 800
const SCR_HEIGHT = 600

GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

# create window
window = GLFW.CreateWindow(SCR_WIDTH, SCR_HEIGHT, "LearnOpenGL")
window == C_NULL && error("Failed to create GLFW window.")
GLFW.MakeContextCurrent(window)
GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)

# render loop
while !GLFW.WindowShouldClose(window)
    processInput(window)
    # swap buffers and poll IO events
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using PyPlot
PyPlot.svg(true)

x = [1:100]
y = [i^2 for i in x]
p = plot(x,y)
xlabel("X")
ylabel("Y")
title("Basic plot")
grid("on")

x = linspace(0, 3pi, 1000)
y = cos(2*x + 3*sin(3*x));
plot(x, y, color="orange", linewidth=2.0, linestyle="--");
title("Another plot using sine and cosine");

X = linspace(0, 2 * pi, 100)
Ya = sin(X)
Yb = cos(X)

plot(X, Ya)
plot(X, Yb)
show()

Pkg.add("UnicodePlots")

using UnicodePlots

myPlot = lineplot([1, 2, 3, 7], [1, 2, -5, 7], title="My Plot", border=:dotted)

myPlot = densityplot([1:100], randn(100), border=:dotted)

using ASCIIPlots

x = collect(-pi:0.2:pi)
lineplot(x, sin(x))

using Vega

x = [0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9]
y = [28, 43, 81, 19, 52, 24, 87, 17, 68, 49, 55, 91, 53, 87, 48, 49, 66, 27, 16, 15]
g = [0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]

a = areaplot(x = x, y = y, group = g, stacked = true)

fruit = ["peaches", "plums", "blueberries", "strawberries", "bananas"];
bushels = [100, 32, 180, 46, 21];
piechart(x = fruit, y = bushels, holesize = 125)

using Bokeh
using Dates # for version 0.3

start = Date(2015, 6, 21)
days = 120
x = map(d -> start + Dates.Day(d), 1:days)
y = 15 + randn(days) * 4
plot(x, y, title="A typical British Summer", legends=["Temperature"])

max_stars_repo_path: Chapter 5/ch05-Pyplot.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 5/ch05-Pyplot.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 5/ch05-Pyplot.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
export Simulation

struct Simulation{A, C, D}
    x::Array{Float64, 2}
    v::Array{Float64, 2}
    Δ_no::MVector{3, Float64}
    nn_arr::Array{Int64, 2}
    dhp_neutral::Array{Float64, 2}
    dhp_ion::Array{Float64, 2}
    dhp_coup::Array{Float64, 2}
    hp::MVector{3, Float64}
    F::Array{Float64, 2}
    surfp::MVector{A, Int16}
    occnum::MVector{C, Int16}
    surfh::MVector{D, Int16}
    surfpinit::MVector{A, Int16}
    surfpnew::MVector{A, Int16}
    trajzmin::MVector{1, Float64}
    trajtheta::MVector{1, Float64}
    H::Array{Float64, 2}
    ψ::Array{ComplexF64, 2}
    ϕ::Array{ComplexF64, 2}
    Γ::Array{Float64, 2}
    λ::MVector{C, Float64}
    dhdea::Array{Float64, 2}
    dhdv::Array{Float64, 2}
    akl::MVector{1, ComplexF64}
    akk::MVector{1, Float64}
    dm::Array{Float64, 2}
    blk::Array{Float64, 2}
    blk2::Array{Float64, 2}
    Pb::Array{Float64, 1}
    storage_aop::Array{Float64, 1}
    storage_op::Array{Float64, 2}
    storage_e::Array{Float64, 2}
    storage_xno::Array{Float64, 2}
    storage_vno::Array{Float64, 2}
    storage_temp::Array{Float64, 1}
    storage_phop::Array{Float64, 1}
    phipsi::MVector{1, ComplexF64}
    Pbmaxest::MVector{1, Float64}
    vdot::Array{Float64, 2}
    vtemp::Array{Float64, 2}
    vscale::MVector{1, Float64}
    storage_K::Array{Float64, 1}
    storage_P::Array{Float64, 1}
    storage_state::Array{Float64, 1}
    storage_deltaKNO::Array{Float64, 1}
    storage_deltaKAu::Array{Float64, 1}
    storage_hoptimes::Array{Int64, 1}
    exnum::MVector{1, Int64}
    attnum::MVector{1, Int64}
    nf::MVector{1, Int64}
end

#standard type for arrays:

const float_array = Array{Float64, 2}


text:
input = joinpath(@__DIR__, "input")
hexes = readlines(input)

function chart(hexes)
    hs = Dict{CartesianIndex,Bool}()
    for hex in hexes
        i = 1
        dest = CartesianIndex(0, 0)
        while i <= length(hex)
            if hex[i:i] == "e"
                dest += CartesianIndex(1, 0)
                i += 1
            elseif hex[i:i] == "w"
                dest += CartesianIndex(-1, 0)
                i += 1
            elseif hex[i:i + 1] == "nw"
                dest += CartesianIndex(0, -1)
                i += 2
            elseif hex[i:i + 1] == "sw"
                dest += CartesianIndex(-1, 1)
                i += 2
            elseif hex[i:i + 1] == "ne"
                dest += CartesianIndex(1, -1)
                i += 2
            elseif hex[i:i + 1] == "se"
                dest += CartesianIndex(0, 1)
                i += 2
            end
        end
        flip!(hs, dest)
    end
    return hs
end

function flip!(hexes, dest)
    if haskey(hexes, dest)
        hexes[dest] = !hexes[dest]
    else
        hexes[dest] = true
    end
end

function adjacent(hex)
     adj = [
         CartesianIndex(1, 0),
         CartesianIndex(-1, 0),
         CartesianIndex(0, -1),
         CartesianIndex(-1, 1),
         CartesianIndex(1, -1),
         CartesianIndex(0, 1)
     ]
     return map(a -> hex + a, adj)
end

function game(hexes, rounds)
    next = copy(hexes)
    for _ in 1:rounds
        floor = next
        neighbors = Dict{CartesianIndex, Int}()
        for hex in keys(floor)
            if floor[hex]
                adj = adjacent(hex)
                for a in adj
                    haskey(neighbors, a) ? neighbors[a] = neighbors[a] + 1 : neighbors[a] = 1
                end
            end
        end
        next = Dict{CartesianIndex, Bool}()
        for (coord, count) in neighbors
            if haskey(floor, coord) && haskey(neighbors, coord) && floor[coord] && neighbors[coord] in [1, 2]
                next[coord] = true
            elseif neighbors[coord] == 2
                next[coord] = true
            end
        end
    end
    return next
end

initial = chart(hexes)

p1 = length(filter(v -> v, collect(values(initial))))
p2 = length(filter(v -> v, collect(values(game(initial, 100)))))

println("-----------------------------------------------------------------------")
println("lobby layout -- part one :: $p1")
println("lobby layout -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 420)
@assert(p2 == 4206)


text:
export PeriodicSegment

"""
	PeriodicSegment(a,b)

represents a periodic interval from `a` to `b`, that is, the point
`b` is identified with `a`.
"""
struct PeriodicSegment{T} <: PeriodicDomain{T}
    a::T
    b::T
    PeriodicSegment{T}() where {T} = new{T}(0,2convert(T,π))
    PeriodicSegment{T}(a,b) where {T} = new{T}(a,b)
end

PeriodicSegment()=PeriodicSegment{Float64}()
PeriodicSegment(a::Int,b::Int) = PeriodicSegment(Float64(a),Float64(b)) #convenience method
PeriodicSegment(a,b) = PeriodicSegment{promote_type(typeof(a),typeof(b))}(a,b)
PeriodicSegment(a::Tuple,b::Tuple) = Interval(Vec(a...),Vec(b...))

function convert(::Type{PeriodicSegment}, d::ClosedInterval)
	a,b = d.left,d.right
    @assert isfinite(a) && isfinite(b)
    PeriodicSegment(a,b)
end

Segment(d::PeriodicSegment) = Segment(leftendpoint(d),rightendpoint(d))
Interval(d::PeriodicSegment) = Interval(leftendpoint(d),rightendpoint(d))
PeriodicSegment(d::Segment) = PeriodicSegment(leftendpoint(d),rightendpoint(d))

convert(::Type{PeriodicSegment{T}}, d::PeriodicSegment) where {T<:Number} = PeriodicSegment{T}(leftendpoint(d),rightendpoint(d))

isambiguous(d::PeriodicSegment) = all(isnan(leftendpoint(d))) && all(isnan(rightendpoint(d)))
convert(::Type{PeriodicSegment{T}},::AnyDomain) where {T<:Number} = PeriodicSegment{T}(NaN,NaN)
convert(::Type{PeriodicSegment{Vec{d,T}}},::AnyDomain) where {d,T} = PeriodicSegment(Vec(fill(NaN,d)...),Vec(fill(NaN,d)...))
convert(::Type{PeriodicSegment{T}},::AnyDomain) where {T} = PeriodicSegment(nan(T),nan(T))
PeriodicSegment{T}(d) where T = convert(PeriodicSegment{T}, d)
PeriodicSegment(d) = convert(PeriodicSegment, d)

## Information
leftendpoint(d::PeriodicSegment) = d.a
rightendpoint(d::PeriodicSegment) = d.b

first(d::PeriodicSegment) = leftendpoint(d)

issubset(a::PeriodicSegment,b::PeriodicSegment) = first(a)∈b && (a.b∈b || a.b == b.b)

# we disable last since the domain is "periodic"

## Map periodic interval

tocanonical(d::PeriodicSegment{T},x) where {T} = π*(tocanonical(Segment(d),x)+1)
tocanonicalD(d::PeriodicSegment{T},x) where {T} = π*tocanonicalD(Segment(d),x)
fromcanonical(d::PeriodicSegment,θ) = fromcanonical(Segment(d),θ/π-1)
fromcanonicalD(d::PeriodicSegment,θ) = fromcanonicalD(Segment(d),θ/π-1)/π

arclength(d::PeriodicSegment) = norm(complexlength(d))
angle(d::PeriodicSegment) = angle(complexlength(d))
complexlength(d::PeriodicSegment) = rightendpoint(d) - leftendpoint(d)
reverseorientation(d::PeriodicSegment) = PeriodicSegment(rightendpoint(d), leftendpoint(d))

==(d::PeriodicSegment,m::PeriodicSegment) = leftendpoint(d) == m.a && rightendpoint(d) == m.b

## algebra

for op in (:*,:+,:-)
    @eval begin
        $op(c::Number,d::PeriodicSegment) = PeriodicSegment($op(c,leftendpoint(d)),$op(c,rightendpoint(d)))
        $op(d::PeriodicSegment,c::Number) = PeriodicSegment($op(leftendpoint(d),c),$op(rightendpoint(d),c))
    end
end

@eval /(d::PeriodicSegment,c::Number) = PeriodicSegment(/(leftendpoint(d),c),/(rightendpoint(d),c))

+(d1::PeriodicSegment,d2::PeriodicSegment) = PeriodicSegment(d1.a+d2.a,d1.b+d2.b)


text:
module Inequality

    using Statistics
    using StatsBase: mean, AbstractWeights, weights

    export
        atkinson,
        watkinson,
        gini,
        wgini,
        lorenz_curve,
        wlorenz_curve,
        mld,
        wmld,
        watts,
        wwatts,
        theil,
        wtheil,
        gen_entropy,
        wgen_entropy,
        fgt,
        wfgt,
        headcount,
        wheadcount,
        poverty_gap,
        wpoverty_gap

    include("Atkinson.jl")
    include("Gini.jl")
    include("Lorenz.jl")
    include("Mld.jl")
    include("Watts.jl")
    include("Theil.jl")
    include("Entropy.jl")
    include("Fgt.jl")
    include("Headcount.jl")
    include("Poverty_gap.jl")
    include("Misc.jl")

end # module


text:
@testset "fsm_ack_rst_relaxed.jl" begin
    function ack_listen_test(seqno::WrappingInt32, ackno::WrappingInt32;
        cap::Int=64000,
        retx_timeout::UInt16=UInt16(1000),
        fixed_isn=nothing)
        conn = fsm_in_listen(;cap, retx_timeout, fixed_isn)
        send_ack!(conn, seqno, ackno)
        expect_state(conn, JLSponge.LISTEN)
        expect_no_seg(conn)
    end

    function ack_rst_syn_sent_test(base_seq::WrappingInt32,
        seqno::WrappingInt32,
        ackno::WrappingInt32;
        cap::Int=64000,
        retx_timeout::UInt16=UInt16(1000),
        fixed_isn=nothing)
        conn = fsm_in_syn_sent(base_seq)
        send_ack!(conn, seqno, ackno)
        expect_state(conn, JLSponge.SYN_SENT)
        expect_no_seg(conn)
    end

    base_seq = WrappingInt32(1 << 31)
    DEFAULT_CAPACITY = 64000
    @testset "in ESTABLISHED, send unacceptable segments and ACKs" begin
        conn = fsm_in_established(base_seq - 1, base_seq - 1)
        send_ack!(conn, base_seq, base_seq)
        expect_no_seg(conn)
        send_ack!(conn, base_seq, base_seq - 1)
        expect_no_seg(conn)

        send_byte!(conn, base_seq - 1, base_seq, Char(1))

        @test conn |> unassembled_bytes == 0
        expect_one_seg(conn; ack=true, ackno=base_seq)
        send_byte!(conn, base_seq + DEFAULT_CAPACITY, base_seq, Char(1))
        @test conn |> unassembled_bytes == 0
        expect_one_seg(conn; ack=true, ackno=base_seq)
        send_byte!(conn, base_seq, base_seq, Char(1))
        expect_data!(conn)
        expect_one_seg(conn; ack=true, ackno=base_seq+1)
        send_rst!(conn, base_seq+1)
        expect_state(conn, JLSponge.RESET)
    end
    @testset "in LISTEN, send RSTs" begin
        conn = fsm_in_listen()
        send_rst!(conn, base_seq)
        send_rst!(conn, base_seq - 1)
        send_rst!(conn, base_seq - DEFAULT_CAPACITY)
        expect_no_seg(conn)
    end
    @testset "ACKs in LISTEN" begin
       ack_listen_test(base_seq, base_seq)
       ack_listen_test(base_seq - 1, base_seq)
       ack_listen_test(base_seq, base_seq - 1)
       ack_listen_test(base_seq - 1, base_seq)
       ack_listen_test(base_seq - 1, base_seq - 1)
       ack_listen_test(base_seq + DEFAULT_CAPACITY, base_seq)
       ack_listen_test(base_seq, base_seq + DEFAULT_CAPACITY)
       ack_listen_test(base_seq + DEFAULT_CAPACITY, base_seq)
       ack_listen_test(base_seq + DEFAULT_CAPACITY, base_seq + DEFAULT_CAPACITY)
    end

    @testset "ACK and RST in SYN_SENT" begin
        conn = fsm_in_syn_sent(base_seq)
        send_rst!(conn, base_seq + 1)
        expect_state(conn, JLSponge.RESET)
        expect_no_seg(conn)
    end

    @testset "ack/rst in SYN_SENT" begin
        ack_rst_syn_sent_test(base_seq, base_seq, base_seq)
        ack_rst_syn_sent_test(base_seq, base_seq, base_seq + 2)
    end
end

"""

"""

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
# -----------------------------------------------------------------------------
# - MILP and direct generalisations (like strict inequalities or indicators)
# -----------------------------------------------------------------------------

MOIU.@model(
    MILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

const StrictlyGreaterThan{T} = CP.Strictly{MOI.GreaterThan{T}, T}
const StrictlyLessThan{T} = CP.Strictly{MOI.LessThan{T}, T}

MOIU.@model(
    PseudoMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval, StrictlyGreaterThan, StrictlyLessThan),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    StrictlyLessThanMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval, StrictlyLessThan),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    StrictlyGreaterThanMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval, StrictlyGreaterThan),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

const EqualToIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}
const EqualToIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.EqualTo{T}}
const GreaterThanIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{T}}
const GreaterThanIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.GreaterThan{T}}
const LessThanIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}
const LessThanIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{T}}
const StrictlyGreaterThanIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE, CP.Strictly{MOI.GreaterThan{T}, T}}
const StrictlyGreaterThanIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.Strictly{MOI.GreaterThan{T}, T}}
const StrictlyLessThanIndicatorOne{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ONE, CP.Strictly{MOI.LessThan{T}, T}}
const StrictlyLessThanIndicatorZero{T} =
    MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.Strictly{MOI.LessThan{T}, T}}

MOIU.@model(
    IndicatorMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        EqualToIndicatorOne,
        EqualToIndicatorZero,
        GreaterThanIndicatorOne,
        GreaterThanIndicatorZero,
        LessThanIndicatorOne,
        LessThanIndicatorZero
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    AbsoluteValueIndicatorPseudoMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval, StrictlyGreaterThan, StrictlyLessThan),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.AbsoluteValue,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        EqualToIndicatorOne,
        EqualToIndicatorZero,
        GreaterThanIndicatorOne,
        GreaterThanIndicatorZero,
        LessThanIndicatorOne,
        LessThanIndicatorZero,
        StrictlyGreaterThanIndicatorOne,
        StrictlyGreaterThanIndicatorZero,
        StrictlyLessThanIndicatorOne,
        StrictlyLessThanIndicatorZero,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    IndicatorPseudoMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval, StrictlyGreaterThan, StrictlyLessThan),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        EqualToIndicatorOne,
        EqualToIndicatorZero,
        GreaterThanIndicatorOne,
        GreaterThanIndicatorZero,
        LessThanIndicatorOne,
        LessThanIndicatorZero,
        StrictlyGreaterThanIndicatorOne,
        StrictlyGreaterThanIndicatorZero,
        StrictlyLessThanIndicatorOne,
        StrictlyLessThanIndicatorZero,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# -----------------------------------------------------------------------------
# - Generic CP constraints.
# -----------------------------------------------------------------------------

# All different.

MOIU.@model(
    AllDifferentIndexingModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.AllDifferent,
        CP.ElementVariableArray,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    AllDifferentInverseModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.AllDifferent,
        CP.Inverse,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Count.

MOIU.@model(
    CountModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Count{MOI.EqualTo{Int}},
        CP.Count{MOI.EqualTo{Float64}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Different from.

MOIU.@model(
    DifferentFromModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        CP.DifferentFrom,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Domain.

MOIU.@model(
    DomainModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        CP.Domain,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Lexicographically greater than.

MOIU.@model(
    LexicographicallyGreaterThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.LexicographicallyGreaterThan,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    IntStrictlyLexicographicallyGreaterThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Strictly{CP.LexicographicallyGreaterThan, Int},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    FloatStrictlyLexicographicallyGreaterThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Strictly{CP.LexicographicallyGreaterThan, Float64},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Lexicographically less than.

MOIU.@model(
    LexicographicallyLessThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.LexicographicallyLessThan,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    IntStrictlyLexicographicallyLessThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Strictly{CP.LexicographicallyLessThan, Int},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    FloatStrictlyLexicographicallyLessThanModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Strictly{CP.LexicographicallyLessThan, Float64},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Global cardinality.

MOIU.@model(
    GlobalCardinalityModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        CP.Domain,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        CP.GlobalCardinalityFixedOpen,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    GlobalCardinalityVariableModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Membership,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        CP.GlobalCardinalityVariableOpen,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Implication.

MOIU.@model(
    ImplicationModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Implication{MOI.LessThan{Int}, MOI.LessThan{Int}},
        CP.Implication{MOI.LessThan{Float64}, MOI.LessThan{Float64}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Indicators about different from.

MOIU.@model(
    FloatDifferentFromIndicatorMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, CP.DifferentFrom{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.EqualTo{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.GreaterThan{Float64}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{Float64}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    IntDifferentFromIndicatorMILPModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.LessThan{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.GreaterThan{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ONE, CP.DifferentFrom{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.EqualTo{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, MOI.GreaterThan{Int}},
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{Int}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Reification of equal-to.

MOIU.@model(
    FloatReificationEqualToModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        CP.DifferentFrom,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Reification{MOI.EqualTo{Float64}},
        CP.Reification{CP.DifferentFrom{Float64}},
        CP.Reification{MOI.LessThan{Float64}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    IntReificationEqualToModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        CP.DifferentFrom,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Reification{MOI.EqualTo{Int}},
        CP.Reification{CP.DifferentFrom{Int}},
        CP.Reification{MOI.LessThan{Int}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    DisjunctionModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.Disjunction{NTuple{2, MOI.LessThan{Int}}},
        CP.Disjunction{NTuple{2, MOI.LessThan{Float64}}},
        CP.Disjunction{NTuple{4, MOI.LessThan{Int}}},
        CP.Disjunction{NTuple{4, MOI.LessThan{Float64}}},
        CP.Disjunction{NTuple{6, MOI.LessThan{Int}}},
        CP.Disjunction{NTuple{6, MOI.LessThan{Float64}}},
        CP.Disjunction{NTuple{2, MOI.EqualTo{Int}}},
        CP.Disjunction{NTuple{2, MOI.EqualTo{Float64}}},
        CP.Disjunction{NTuple{4, MOI.EqualTo{Int}}},
        CP.Disjunction{NTuple{4, MOI.EqualTo{Float64}}},
        CP.Disjunction{NTuple{6, MOI.EqualTo{Int}}},
        CP.Disjunction{NTuple{6, MOI.EqualTo{Float64}}},
        CP.Disjunction{Tuple{CP.Domain{Int}, CP.Domain{Int}, CP.DifferentFrom{Int}}},
        CP.Disjunction{Tuple{CP.Domain{Float64}, CP.Domain{Float64}, CP.DifferentFrom{Float64}}},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# -----------------------------------------------------------------------------
# - Higher-level CP functions.
# -----------------------------------------------------------------------------

# Absolute value.

MOIU.@model(
    AbsoluteValueModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.AbsoluteValue
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# -----------------------------------------------------------------------------
# - Higher-level CP constraints.
# -----------------------------------------------------------------------------

# Bin packing.

MOIU.@model(
    BinPackingModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, CP.BinPackingNoCapacity),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    VariableCapacityBinPackingModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, CP.BinPackingVariableCapacity),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Knapsack.

MOIU.@model(
    KnapsackModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, CP.KnapsackFixedCapacityUnvalued),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

MOIU.@model(
    VariableCapacityKnapsackModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
    ),
    (MOI.PowerCone, MOI.DualPowerCone, CP.KnapsackVariableCapacityUnvalued),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Non-overlapping orthotopes.

MOIU.@model(
    ConditionallyNonOverlappingOrthotopesModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.NonOverlappingOrthotopes{CP.CONDITIONAL_NONVERLAPPING_ORTHOTOPES},
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

# Sorting.

MOIU.@model(
    SortPermutationModel,
    (),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
    ),
    (
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.ExponentialCone,
        CP.SortPermutation,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)


text:
function eigenvector(index::Integer,n::Integer)
    res = 1
    for i in bin(index,n)
        if i=='0'
            res = kron(res,1/sqrt(2)*[1,1])
        elseif i=='1'
            res = kron(res,1/sqrt(2)*[1,-1])
        end
    end
    return res
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
# toplevel
# --------

abstract type ToplevelErrorReport end

# `ToplevelErrorReport` interface
function Base.getproperty(er::ToplevelErrorReport, sym::Symbol)
    return if sym === :file
        getfield(er, sym)::String
    elseif sym === :line
        getfield(er, sym)::Int
    else
        getfield(er, sym) # fallback
    end
end

struct SyntaxErrorReport <: ToplevelErrorReport
    err::ErrorException
    file::String
    line::Int
    SyntaxErrorReport(msg::AbstractString, file, line) = new(ErrorException(msg), file, line)
end

struct RecursiveIncludeErrorReport <: ToplevelErrorReport
    duplicated_file::String
    files::Set{String}
    file::String
    line::Int
end

# wraps general errors from actual execution
struct ActualErrorWrapped <: ToplevelErrorReport
    err
    st::Vector{Base.StackTraces.StackFrame}
    file::String
    line::Int

    # default constructor
    ActualErrorWrapped(err, st, file, line) = new(err, st, file, line)

    # bypass syntax error
    function ActualErrorWrapped(err::ErrorException, st, file, line)
        return if startswith(err.msg, "syntax: ")
            SyntaxErrorReport(err.msg, file, line)
        else
            new(err, st, file, line)
        end
    end
end

# wraps an error that might happen because of inappropriate top-level code abstraction
struct MissingConcretization <: ToplevelErrorReport
    err
    st::Vector{Base.StackTraces.StackFrame}
    file::String
    line::Int
end

# inference
# ---------

# TODO: maybe we want to use https://github.com/aviatesk/Mixin.jl
abstract type InferenceErrorReport end

# to help inference
function Base.getproperty(er::InferenceErrorReport, sym::Symbol)
    return if sym === :st
        getfield(er, sym)::VirtualStackTrace
    elseif sym === :msg
        getfield(er, sym)::String
    elseif sym === :sig
        getfield(er, sym)::Vector{Any}
    elseif sym === :lin # only needed for ExceptionReport
        getfield(er, sym)::LineInfoNode
    else
        getfield(er, sym) # fallback
    end
end

function Base.show(io::IO, report::T) where {T<:InferenceErrorReport}
    print(io, T.name.name, '(')
    for a in report.sig
        _print_signature(io, a, (; annotate_types = true); bold = true)
    end
    print(io, ')')
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", report::T) where {T<:InferenceErrorReport} =
    return report

@withmixedhash struct VirtualFrame
    file::Symbol
    line::Int
    sig::Vector{Any}
    linfo::MethodInstance
end

# `InferenceErrorReport` is supposed to keep its virtual stack trace in order of
# "from entry call site to error point"
const VirtualStackTrace = Vector{VirtualFrame}

struct InferenceErrorReportCache
    T::Type{<:InferenceErrorReport}
    st::VirtualStackTrace
    msg::String
    sig::Vector{Any}
    spec_args::NTuple{N,Any} where N
end

function cache_report!(cache, report::T) where {T<:InferenceErrorReport}
    st = copy(report.st)
    new = InferenceErrorReportCache(T, st, report.msg, report.sig, spec_args(report))
    return push!(cache, new)
end

function restore_cached_report!(cache::InferenceErrorReportCache,
                                interp#=::JETInterpreter=#,
                                )
    report = restore_cached_report(cache)
    if isa(report, UncaughtExceptionReport)
        stash_uncaught_exception!(interp, report)
    else
        report!(interp, report)
    end
    return report
end

function restore_cached_report(cache::InferenceErrorReportCache)
    T = cache.T
    st = copy(cache.st)
    return T(st, cache.msg, cache.sig, cache.spec_args)::InferenceErrorReport
end

macro reportdef(ex, kwargs...)
    T = esc(first(ex.args))
    args = map(ex.args) do x
        # unwrap @nospecialize
        if @isexpr(x, :macrocall) && first(x.args) === Symbol("@nospecialize")
            x = esc(last(x.args)) # `esc` is needed because `@nospecialize` escapes its argument anyway
        end
        return x
    end
    spec_args = args[4:end] # keep those additional, specific fields

    local track_from_frame = false
    local supertype = InferenceErrorReport
    for ex in kwargs
        @assert @isexpr(ex, :(=))
        kw, val = ex.args
        if kw === :track_from_frame
            track_from_frame = val
        elseif kw === :supertype
            supertype = val
        end
    end

    args′ = strip_type_decls.(args)
    spec_args′ = strip_type_decls.(spec_args)
    interp_constructor = Expr(:function, ex, Expr(:block, __source__, quote
        interp = $(args′[2])
        sv = $(args′[3])

        msg = get_msg(#= T, interp, sv, ... =# $(args′...))
        sig = get_sig(#= T, interp, sv, ... =# $(args′...))

        $(if track_from_frame quote
            # when report is constructed _after_ the inference on `sv` has been done,
            # collect location information from `sv.linfo`
            st = VirtualFrame[get_virtual_frame(interp, sv.linfo)]
        end else quote
            st = VirtualFrame[get_virtual_frame(interp, sv)]
        end end)

        return new(st, msg, sig, $(spec_args′...))
    end))

    spec_types = extract_type_decls.(spec_args)

    cache_constructor_sig = :($(T)(st::VirtualStackTrace,
                                   msg::AbstractString,
                                   sig::AbstractVector,
                                   @nospecialize(spec_args),
                                   ))
    cache_constructor_call = :(new(st, msg, sig))
    for (i, spec_type) in enumerate(spec_types)
        push!(cache_constructor_call.args,
              :($(esc(:spec_args))[$(i)]::$(spec_type)), # `esc` is needed because `@nospecialize` escapes its argument anyway
              )
    end
    cache_constructor = Expr(:function, cache_constructor_sig, Expr(:block, __source__,
        :(return @inbounds $(cache_constructor_call))
    ))

    spec_getter_sig = :($(esc(:spec_args))(report::$(T)))
    spec_getter_tuple = Expr(:tuple)
    for spec_arg in QuoteNode.(strip_escape.(spec_args′))
        push!(spec_getter_tuple.args,
              Expr(:call, GlobalRef(Base, :getproperty), :report, spec_arg),
              )
    end
    spec_args_getter = Expr(:function, spec_getter_sig, Expr(:block, __source__,
        :(return $(spec_getter_tuple)::Tuple{$(spec_types...)})
    ))

    return Expr(:block, __source__, quote
        struct $(T) <: $(supertype)
            st::VirtualStackTrace
            msg::String
            sig::Vector{Any}
            $(spec_args...)

            # constructor from abstract interpretation process by `JETInterpreter`
            $(interp_constructor)

            # constructor from cache
            $(cache_constructor)
        end

        $(spec_args_getter)
    end)
end

function strip_type_decls(x)
    @isexpr(x, :escape) && return Expr(:escape, strip_type_decls(first(x.args))) # keep escape
    return @isexpr(x, :(::)) ? first(x.args) : x
end

strip_escape(x) = @isexpr(x, :escape) ? first(x.args) : x
extract_type_decls(x) = @isexpr(x, :(::)) ? last(x.args) : Any

@reportdef NoMethodErrorReport(interp, sv, unionsplit::Bool, @nospecialize(atype::Type))

@reportdef InvalidBuiltinCallErrorReport(interp, sv, argtypes::Vector{Any})

@reportdef NoFieldErrorReport(interp, sv, @nospecialize(typ::Type), name::Symbol)

@reportdef GlobalUndefVarErrorReport(interp, sv, mod::Module, name::Symbol)

@reportdef LocalUndefVarErrorReport(interp, sv, name::Symbol) track_from_frame = true

@reportdef NonBooleanCondErrorReport(interp, sv, @nospecialize(t::Union{Type,Vector{Type}}))

@reportdef DivideErrorReport(interp, sv)

# TODO we may want to hoist `InvalidConstXXX` errors into top-level errors

@reportdef InvalidConstantRedefinition(interp, sv, mod::Module, name::Symbol, @nospecialize(t′), @nospecialize(t))

@reportdef InvalidConstantDeclaration(interp, sv, mod::Module, name::Symbol)

"""
    ExceptionReport <: InferenceErrorReport

The abstract type for "serious" errors that are invoked by `throw` calls but should be
    reported even if they may be caught in actual execution.
In order to avoid duplicated reports for the `throw` call, any subtype of `ExceptionReport`
    should keep `lin::LineInfoNode` field, which represents where the report gets collected.
"""
abstract type ExceptionReport <: InferenceErrorReport end

# # NOTE: this mixin implementation is cleaner but doesn't help inference,
# # because this `getproperty` interface relies on constant prop' and it won't happen when
# # there're multiple applicable methods.
# function Base.getproperty(er::ExceptionReport, sym::Symbol)
#     sym === :lin && return getfield(er, :lin)::LineInfoNode
#     return @invoke getproperty(er::InferenceErrorReport, sym::Symbol)
# end

@reportdef UndefKeywordErrorReport(interp, sv, err::UndefKeywordError, lin::LineInfoNode) supertype = ExceptionReport

"""
    UncaughtExceptionReport <: InferenceErrorReport

Represents general `throw` calls traced during inference.
They are reported only when they're not caught by any control flow.
"""
:(UncaughtExceptionReport)
@reportdef UncaughtExceptionReport(interp, sv, throw_calls::Vector{Expr}) track_from_frame = true

"""
    NativeRemark <: InferenceErrorReport

This special `InferenceErrorReport` is just for wrapping remarks from `NativeInterpreter`.

!!! note
    Currently JET.jl doesn't make any use of `NativeRemark`.
"""
:(NativeRemark)
@reportdef NativeRemark(interp, sv, s::String)

function get_virtual_frame(interp#=::JETInterpreter=#, loc::Union{InferenceState,MethodInstance})
    sig = get_sig(interp, loc)
    file, line = get_file_line(loc)
    linfo = isa(loc, InferenceState) ? loc.linfo : loc
    return VirtualFrame(file, line, sig, linfo)
end

get_file_line(frame::InferenceState) = get_file_line(get_lin(frame))
get_file_line(linfo::LineInfoNode)   = linfo.file, linfo.line
# this location is not exact, but this is whay we know at best
function get_file_line(linfo::MethodInstance)
    def = linfo.def

    isa(def, Method) && return def.file, Int(def.line)

    # toplevel
    src = linfo.uninferred::CodeInfo
    return get_file_line(first(src.linetable::Vector)::LineInfoNode)
end

# adapted from https://github.com/JuliaLang/julia/blob/0f11a7bb07d2d0d8413da05dadd47441705bf0dd/base/show.jl#L989-L1011
function get_sig(interp#=::JETInterpreter=#, l::MethodInstance)
    def = l.def
    ret = if isa(def, Method)
        if isdefined(def, :generator) && l === def.generator
            # print(io, "MethodInstance generator for ")
            # show(io, def)
            sprint(show, def)
        else
            # print(io, "MethodInstance for ")
            @static if hasmethod(Base.show_tuple_as_call, (IO, Symbol, Type), (:demangle, :kwargs, :argnames, :qualified))
                # show_tuple_as_call(io, def.name, l.specTypes; qualified=true)
                sprint((args...)->Base.show_tuple_as_call(args...; qualified=true), def.name, l.specTypes)
            else
                # show_tuple_as_call(io, def.name, l.specTypes, false, nothing, nothing, true)
                sprint(Base.show_tuple_as_call, def.name, l.specTypes, false, nothing, nothing, true)
            end
        end
    else
        # print(io, "Toplevel MethodInstance thunk")
        # # `thunk` is not very much information to go on. If this
        # # MethodInstance is part of a stacktrace, it gets location info
        # # added by other means.  But if it isn't, then we should try
        # # to print a little more identifying information.
        # if !from_stackframe
        #     linetable = l.uninferred.linetable
        #     line = isempty(linetable) ? "unknown" : (lt = linetable[1]; string(lt.file) * ':' * string(lt.line))
        #     print(io, " from ", def, " starting at ", line)
        # end
        "toplevel"
    end
    return Any[ret]
end

# entry
get_sig(::Type{<:InferenceErrorReport}, interp#=::JETInterpreter=#, sv, @nospecialize(args...)) = get_sig(interp, sv)
get_sig(interp#=::JETInterpreter=#, sv::InferenceState) = _get_sig(interp, sv, get_stmt(sv))

# special cased entries
get_sig(::Type{LocalUndefVarErrorReport}, interp#=::JETInterpreter=#, sv, name) = Any[""] # TODO
function get_sig(::Type{UncaughtExceptionReport}, interp#=::JETInterpreter=#, sv, throw_calls)
    sig = Any[]
    ncalls = length(throw_calls)
    for (i, call) in enumerate(throw_calls)
        call_sig = _get_sig(interp, sv, call)
        append!(sig, call_sig)
        i ≠ ncalls && push!(sig, ", ")
    end
    return sig
end

_get_sig(args...) = first(_get_sig_type(args...))::Vector{Any}

function _get_sig_type(interp#=::JETInterpreter=#, sv::InferenceState, expr::Expr)
    head = expr.head
    return if head === :call
        f = first(expr.args)
        args = expr.args[2:end]

        # special case splat call signature
        if isa(f, GlobalRef) && f.name === :_apply_iterate && begin
                itf = first(args)
                isa(itf, GlobalRef) && itf.name === :iterate
            end
            f = args[2]
            args = args[3:end]

            sig = _get_sig(interp, sv, f)
            push!(sig, '(')

            nargs = length(args)
            for (i, arg) in enumerate(args)
                arg_sig = _get_sig(interp, sv, arg)
                append!(sig, arg_sig)
                i ≠ nargs ? push!(sig, ", ") : push!(sig, "...)")
            end
        else
            sig = _get_sig(interp, sv, f)
            push!(sig, '(')

            nargs = length(args)
            for (i, arg) in enumerate(args)
                arg_sig = _get_sig(interp, sv, arg)
                append!(sig, arg_sig)
                i ≠ nargs && push!(sig, ", ")
            end
            push!(sig, ')')
        end

        sig, nothing
    elseif head === :(=)
        _get_sig_type(interp, sv, last(expr.args))
    elseif head === :static_parameter
        typ = widenconst(sv.sptypes[first(expr.args)])
        Any['_', typ], typ
    else
        Any[string(expr)], nothing
    end
end
function _get_sig_type(interp#=::JETInterpreter=#, sv::InferenceState, ssa::SSAValue)
    sig, sig_typ = _get_sig_type(interp, sv, sv.src.code[ssa.id])
    typ = widenconst(ignorelimited(ignorenotfound(sv.src.ssavaluetypes[ssa.id])))
    sig_typ == typ || push!(sig, typ)
    return sig, typ
end
function _get_sig_type(interp#=::JETInterpreter=#, sv::InferenceState, slot::SlotNumber)
    name = get_slotname(sv, slot)
    sig = string(name)
    if isempty(sig)
        sig = string(slot) # fallback if no explicit slotname
    end
    if istoplevel(interp, sv)
        # this is a abstract global variable, form the global reference
        return _get_sig_type(interp, sv, GlobalRef(interp.toplevelmod, name))
    else
        typ = widenconst(ignorelimited(get_slottype(sv, slot)))
        return Any[sig, typ], typ
    end
end
_get_sig_type(interp#=::JETInterpreter=#, ::InferenceState, gr::GlobalRef) = Any[string(gr.mod, '.', gr.name)], nothing
function _get_sig_type(interp#=::JETInterpreter=#, sv::InferenceState, s::Symbol)
    if istoplevel(interp, sv)
        # this is concrete global variable, form the global reference
        return _get_sig_type(interp, sv, GlobalRef(interp.toplevelmod, s))
    else
        return Any[repr(s; context = :compact => true)]
    end
end
function _get_sig_type(interp, sv::InferenceState, gotoifnot::GotoIfNot)
    sig  = Any[string("goto %", gotoifnot.dest, " if not "), _get_sig(interp, sv, gotoifnot.cond)...]
    return sig, nothing
end
function _get_sig_type(interp#=::JETInterpreter=#, sv::InferenceState, rn::ReturnNode)
    sig = is_unreachable(rn) ? Any["unreachable"] : Any["return ", _get_sig(interp, sv, rn.val)...]
    return sig, nothing
end
function _get_sig_type(interp#=::JETInterpreter=#, ::InferenceState, qn::QuoteNode)
    typ = typeof(qn.value)
    return Any[string(qn), typ], typ
end
_get_sig_type(interp#=::JETInterpreter=#, ::InferenceState, @nospecialize(x)) = Any[repr(x; context = :compact => true)], nothing

# TODO count invalid unon split case
get_msg(::Type{NoMethodErrorReport}, interp, sv, unionsplit, @nospecialize(args...)) = unionsplit ?
    "for any of the union split cases, no matching method found for call signature" :
    "no matching method found for call signature"
get_msg(::Type{InvalidBuiltinCallErrorReport}, interp, sv, @nospecialize(args...)) =
    "invalid builtin function call"
get_msg(::Type{NoFieldErrorReport}, interp, sv, @nospecialize(typ), name) =
    "type $(typ) has no field $(name)"
get_msg(::Type{GlobalUndefVarErrorReport}, interp, sv, mod, name) =
    "variable $(mod).$(name) is not defined"
get_msg(::Type{LocalUndefVarErrorReport}, interp, sv, name) =
    "local variable $(name) is not defined"
get_msg(::Type{NonBooleanCondErrorReport}, interp, sv, @nospecialize(t::Type)) =
    "non-boolean ($t) used in boolean context"
get_msg(::Type{NonBooleanCondErrorReport}, interp, sv, ts::Vector{Type}) =
    "for $(length(ts)) of union split cases, non-boolean ($(join(ts, ','))) used in boolean context"
@eval get_msg(::Type{DivideErrorReport}, interp, sv) = $(let
    io = IOBuffer()
    showerror(io, DivideError())
    String(take!(io))
end)
get_msg(::Type{InvalidConstantRedefinition}, interp, sv, mod, name, @nospecialize(t′), @nospecialize(t)) =
    "invalid redefinition of constant $(mod).$(name) (from $(t′) to $(t))"
get_msg(::Type{InvalidConstantDeclaration}, interp, sv, mod, name) =
    "cannot declare a constant $(mod).$(name); it already has a value"
get_msg(::Type{UndefKeywordErrorReport}, interp, sv, err, lin) = sprint(showerror, err)
get_msg(::Type{UncaughtExceptionReport}, interp, sv, throw_blocks) = isone(length(throw_blocks)) ?
    "may throw" :
    "may throw either of"
get_msg(::Type{NativeRemark}, interp, sv, s) = s


text:
## Intro

using ApproxFun,Base.Test
x = Fun(identity,[0.,10.])
f = sin(x^2)
g = cos(x)

@test_approx_eq_eps f(.1) sin(.1^2) 1000eps()

h = f + g^2
r = roots(h)
rp = roots(differentiate(h))

@test norm(h(r))<1000eps()

@test norm(h'(rp))<100000eps()

## Differentiation and Integration
f = Fun(x->exp(x),[-1.,1.])
@test norm(f-f')<1000eps()

g = cumsum(f)
g = g + f(-1)
@test norm(f-g)<100eps()

x = Fun(identity)
f = exp(x)
g = f/sqrt(1-x^2)

space(f),domain(f)
space(g),domain(g)

f = Fun(x->cospi(5x))
g = abs(f)
space(f)
space(g)

x = Fun()
f = erf(x)
g = besselj(3,exp(f))
h = airyai(10asin(f)+2g)

## Solving ODEs

x = Fun(identity,[-1000.,200.])
d = domain(x)
D = Derivative(d)
B = dirichlet(d)
L = D^2 - x
u = [B;L] \ [airyai(d.a);airyai(d.b)]

@test_approx_eq_eps u(0.) airyai(0.) 10000eps()

## Nonlinear BVPs
x=Fun()
u0=0.x

N=u->[u(-1.)-1.,u(1.)+0.5,0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)

@test norm(N(u)[end]) ≤ 1000eps()

## Periodic Functions

f = Fun(cos,Fourier([-π,π]))
@test norm(differentiate(f) + Fun(sin,Fourier([-π,π])))<100eps()

s = Chebyshev([-π,π])
a = Fun(t-> 1+sin(cos(2t)),s)
L = Derivative() + a
f = Fun(t->exp(sin(10t)),s)
B = periodic(s,0)
uChebyshev = [B;L]\[0.,f]

s = Fourier([-π,π])
a = Fun(t-> 1+sin(cos(2t)),s)
L = Derivative() + a
f = Fun(t->exp(sin(10t)),s)
uFourier = L\f

@ test_approx_eq uChebyshev(0.) uFourier(0.)

## Sampling

f = abs(Fun(sin,[-5,5]))
x = ApproxFun.sample(f,10000)

## PDEs

d = Interval()^2                            # Defines a rectangle

u = [dirichlet(d);lap(d)+100I]\ones(4)      # First four entries of rhs are

d = Interval()^2
u0 = Fun((x,y)->exp(-40(x-.1)^2-40(y+.2)^2),d)
B = dirichlet(d)
D = Derivative(Interval())
L = (0.01D^2-4D)⊗I + I⊗(0.01D^2-3D)
h = 0.002

with_bigfloat_precision(1000) do
    d=Interval{BigFloat}(0,1)
    D=Derivative(d)
    u=[ldirichlet();D-I]\[1]
    @test_approx_eq u(1) exp(BigFloat(1))
end


text:
module NASTRAN_import_export_examples
using FinEtools
using FinEtools.MeshImportModule: import_NASTRAN
using FinEtools.MeshExportModule

function NASTRAN_import_export()
    output = import_NASTRAN("$(@__DIR__)" * "/Slot-coarser.nas")
    File = "Slot-coarser.vtk"
    MeshExportModule.vtkexportmesh(File, output["fens"], output["fesets"][1])
    @async run(`"paraview.exe" $File`)
end # NASTRAN_import_export

function allrun()
    println("#####################################################")
    println("# NASTRAN_import_export ")
    NASTRAN_import_export()
    return true
end # function allrun

end # module NASTRAN_import_export_examples


text:
using Test

@testset "Collision: String" begin

    # Parameters
    origin = Origin{Float64}()
    pbody = Sphere(0.5, 1.0)
    cbody = Sphere(0.5, 1.0)
    joint = JointConstraint(Fixed(origin, pbody))

    bodies = [pbody, cbody]
    joints = [joint]

    collision = StringCollision{Float64,2,3,6}(
            szeros(3),
            szeros(3),
            2.0)
    friction_parameterization = SA{Float64}[
        1.0  0.0
        0.0  1.0
    ]
    body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

    contacts = [ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)]

    mech = Mechanism(origin, bodies, joints, contacts,
                gravity=1.0 * -9.81,
                timestep=0.1)

    mech.bodies[1].state.x2 = [0.0, 0.0, 0.0]
    mech.bodies[2].state.x2 = [0.0, 0.0, 1.0]

    storage = simulate!(mech, 2.0,
        verbose=false,
        record=true)

    # vis = Visualizer()
    # render(vis)
    # visualize(mech, storage,
    #     vis=vis)

    @test norm(mech.bodies[2].state.x2 - mech.bodies[1].state.x2) < 2.0
end

# Parameters
origin = Origin{Float64}()
pbody = Sphere(0.5, 1.0)
cbody = Sphere(0.5, 1.0)
joint = JointConstraint(Fixed(origin, pbody))

bodies = [pbody, cbody]
joints = [joint]

# collision = SphereSphereCollision{Float64,2,3,6}(
#         szeros(3),
#         szeros(3),
#         pbody.shape.r,
#         cbody.shape.r)
# friction_parameterization = SA{Float64}[
#     1.0  0.0
#     0.0  1.0
# ]
# body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

# contacts = [ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)]

collision = StringCollision{Float64,0,3,0}(
        1.0 * [0.0; 0.0; -0.5],
        1.0 * [0.1; 0.0; 0.0],
        2.0)
friction_parameterization = szeros(Float64, 0, 2)

#SA{Float64}[
#     1.0  0.0
#     0.0  1.0
# ]
body_body_contact = ImpactContact{Float64,2}(friction_parameterization, collision)

contacts = [ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)]

mech = Mechanism(origin, bodies, joints, contacts,
            gravity=1.0 * -9.81,
            timestep=0.1)

mech.bodies[1].state.x2 = [0.0, 0.0, 0.0]
mech.bodies[2].state.x2 = [0.25, 0.5, 1.5]

storage = simulate!(mech, 5.0,
    verbose=true,
    record=true)

vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis)

@test norm(mech.bodies[2].state.x2 - mech.bodies[1].state.x2) < 2.0

function test_solmat(;
    ϵ=1.0e-6,
    tsim=0.1,
    ctrl=(m, k)->nothing,
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    verbose=false,
    T=Float64,
    kwargs...)

    origin = Origin{Float64}()
    pbody = Sphere(0.5, 1.0)
    cbody = Sphere(0.5, 1.0)
    joint = JointConstraint(Fixed(origin, pbody))

    bodies = [pbody, cbody]
    joints = [joint]

    collision = StringCollision{Float64,2,3,6}(
            szeros(3),
            szeros(3),
            2.0)
    friction_parameterization = SA{Float64}[
        1.0  0.0
        0.0  1.0
    ]
    body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

    contacts = [ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)]

    mechanism = Mechanism(origin, bodies, joints, contacts,
                gravity=1.0 * -9.81,
                timestep=0.1)

    mechanism.bodies[1].state.x2 = [0.0, 0.0, 0.0]
    mechanism.bodies[2].state.x2 = [0.0, 0.0, 1.0]

    # simulate
    storage = simulate!(mechanism, tsim, ctrl,
        record=true,
        verbose=verbose,
        opts=SolverOptions(rtol=ϵ, btol=ϵ))

    # Set data
    Nb = length(mechanism.bodies)
    data = Dojo.get_data(mechanism)
    Dojo.set_data!(mechanism, data)
    sol = Dojo.get_solution(mechanism)
    attjac = Dojo.attitude_jacobian(data, Nb)

    # IFT
    solmat = Dojo.full_matrix(mechanism.system)
    # finite diff
    fd_solmat = finite_difference_solution_matrix(mechanism, data, sol,
        δ=1.0e-5,
        verbose=verbose)
    # @test norm(fd_solmat + solmat, Inf) < ϵ

    return fd_solmat, solmat
end

get_data(model::NonlinearContact) = [model.friction_coefficient; model.collision.origin_parent; model.collision.origin_child;]
# Contact
function set_data!(model::NonlinearContact, data::AbstractVector)
	model.friction_coefficient = data[1]
    model.collision.origin_parent = data[2:4]
    model.collision.origin_child = data[5:7]
    return nothing
end
data_dim(model::NonlinearContact) = 7 # [friction_coefficient, contact_radius, p]

f1, f2 = test_solmat()
Gray.(f1)
Gray.(f2)
Gray.(f1 + f2)
norm(f1 + f2, Inf)
Gray.(1e6 * abs.(f1 + f2))

f1[end-1:end, :]
-f2[end-1:end, :]