text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)

struct DataLoader
    data
    batchsize::Int
    nobs::Int
    partial::Bool
    imax::Int
    indices::Vector{Int}
    shuffle::Bool
end

"""
     DataLoader(data...; batchsize=1, shuffle=false, partial=true)

An object that iterates over mini-batches of `data`, each mini-batch containing `batchsize` observations
(except possibly the last one).

Takes as input one or more data tensors, e.g. X in unsupervised learning, X and Y in
supervised learning. The last dimension in each tensor is considered to be the observation
dimension.

If `shuffle=true`, shuffles the observations each time iterations are re-started.
If `partial=false`, drops the last mini-batch if it is smaller than the batchsize.

Example usage:

    Xtrain = rand(10, 100)
    dtrain = DataLoader(Xtrain, batchsize=2)
    # iterate over 50 mini-batches
    for x in dtrain:
        @assert size(x) == (10, 2)
        ...
    end

    Xtrain = rand(10, 100)
    Ytrain = rand(100)
    dtrain = DataLoader(Xtrain, Ytrain, batchsize=2, shuffle=true)
    for epoch in 1:100
        for (x, y) in dtrain:
            @assert size(x) == (10, 2)
            @assert size(y) == (2,)
            ...
        end
    end

    # train for 10 epochs
    using IterTools: ncycle
    Flux.train!(loss, ps, ncycle(dtrain, 10), opt)
"""
function DataLoader(data...; batchsize=1, shuffle=false, partial=true)
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))

    nx = size(data[1])[end]
    for i=2:length(data)
        nx != size(data[i])[end] && throw(DimensionMismatch("All data should contain same number of observations"))
    end
    if nx < batchsize
        @warn "Number of data points less than batchsize, decreasing the batchsize to $nx"
        batchsize = nx
    end
    imax = partial ? nx : nx - batchsize + 1
    ids = 1:min(nx, batchsize)
    DataLoader(data, batchsize, nx, partial, imax, [1:nx;], shuffle)
end

getdata(x::AbstractArray, ids) = x[(Base.Colon() for _=1:ndims(x)-1)..., ids]

@propagate_inbounds function Base.iterate(d::DataLoader, i=0)     # returns data in d.indices[i+1:i+batchsize]
    i >= d.imax && return nothing
    if d.shuffle && i == 0
        shuffle!(d.indices)
    end
    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i+1:nexti]
    if length(d.data) == 1
        batch = getdata(d.data[1], ids)
    else
        batch = ((getdata(x, ids) for x in d.data)...,)
    end
    return (batch, nexti)
end

function Base.length(d::DataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end


text:
using Pkg
using RegistryTools.Compress: compress_versions

@testset "compress_versions()" begin
    # Test exact version matching
    vs = [v"1.1.0", v"1.1.1", v"1.1.2"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test holes
    vs = [v"1.1.0", v"1.1.1", v"1.1.4"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test patch variation with length(subset) > 1
    vs = [v"1.1.0", v"1.1.1", v"1.1.2", v"1.1.3", v"1.2.0"]
    @test compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1.1.2")

    # Test minor variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1-1.1")

    # Test major variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0", v"2.0.0"]
    @test compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1")

    # Test build numbers and prerelease values are ignored
    vs = [v"1.1.0-alpha", v"1.1.0+0", v"1.1.0+1"]
    @test compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1")
end

ext: jl


text:
module Gomah

using PyCall
using PyCall: @pywith
export @pywith

export chainer
# defined in mnist.jl
export train, predict
export debug, get_model
export reversedims
export ch2dense, ch2conv, ch2dwconv, ch2bn
export Conv2DBNActiv, BottleNeckA, BottleNeckB, ResBlock, ResNet

# chainer object
const chainer = PyNULL()
const chainercv = PyNULL()
const math = PyNULL()
const np = PyNULL()
const F = PyNULL()
const L = PyNULL()
const SerialIterator = PyNULL()
const training = PyNULL()
const extensions = PyNULL()
const optimizers = PyNULL()

# module initializer
function __init__()
    copy!(chainer, pyimport("chainer"))
    copy!(chainercv, pyimport("chainercv"))
    copy!(math, pyimport("math"))
    copy!(np, pyimport("numpy"))
    copy!(F, chainer.functions)
    copy!(L, chainer.links)
    copy!(SerialIterator, chainer.iterators.SerialIterator)
    copy!(training, chainer.training)
    copy!(extensions, training.extensions)
    copy!(optimizers, chainer.optimizers)
end

include("mnist.jl")
include("chainer2flux.jl")
include("resnet.jl")

end # module


text:
nk = (2,2,2)
R  = DFW.Vec3(2, 0, 0)
fermi = 13.129
n_ωh = 300
n_ωv = 50
ωh = -30.0
ωv = 0.5

i, structure = DFW.DFControl.wan_read_input(assetfile("wanup.win"))
hami = readhami(DFW.read_chk(assetfile("wanup.chk")), assetfile("wanup.eig"))

wbands = wannierbands(hami, uniform_kgrid(nk...))
@test isapprox(sum(wbands[3].eigvals), 26.24675680998617)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_edim_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_edim")
JLLWrappers.@generate_main_file("GAP_pkg_edim", UUID("aa0ab9cb-918b-5ca1-812e-3cecb1c413b6"))
end  # module GAP_pkg_edim_jll


text:
function save(f::Stream{format"PLY_BINARY"}, msh::AbstractMesh)
    io = stream(f)
    vts = decompose(Point{3, Float32}, msh)
    fcs = decompose(Face{3, Int32, -1}, msh)

    nV = length(vts)
    nF = length(fcs)

    # write the header
    write(io, "ply\n")
    write(io, "format binary_little_endian 1.0\n")
    write(io, "element vertex $nV\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    write(io, "element face $nF\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    write(io, vts)

    for f in fcs
        write(io, convert(UInt8, 3))
        write(io, f...)
    end
    close(io)
end

function save(f::Stream{format"PLY_ASCII"}, msh::AbstractMesh)
    io = stream(f)
    vts = vertices(msh)
    fcs = faces(msh)

    nV = length(vts)
    nF = length(fcs)

    # write the header
    write(io, "ply\n")
    write(io, "format ascii 1.0\n")
    write(io, "element vertex $nV\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    write(io, "element face $nF\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    for v in vts
        println(io, join(Point{3, Float32}(v), " "))
    end
    for f in fcs
        println(io, length(f), " ", join(Face{3, Cuint, -1}(f), " "))
    end
    close(io)
end

function load(fs::Stream{format"PLY_ASCII"}, MeshType=GLNormalMesh)
    io = stream(fs)
    nV = 0
    nF = 0

    properties = Compat.UTF8String[]

    # read the header
    line = readline(io)

    while !startswith(line, "end_header")
        if startswith(line, "element vertex")
            nV = parse(Int, split(line)[3])
        elseif startswith(line, "element face")
            nF = parse(Int, split(line)[3])
        elseif startswith(line, "property")
            push!(properties, line)
        end
        line = readline(io)
    end
    VertexType  = vertextype(MeshType)
    FaceType    = facetype(MeshType)
    FaceEltype  = eltype(FaceType)

    vts         = Array(VertexType, nV)
    #fcs         = Array(FaceType, nF)
    fcs         = FaceType[]

    # read the data
    for i = 1:nV
        vts[i] = VertexType(split(readline(io))) # line looks like: "-0.018 0.038 0.086"
    end

    for i = 1:nF
        line    = split(readline(io))
        len     = parse(Int, shift!(line))
        if len == 3
            push!(fcs, Face{len, FaceEltype, -1}(line)) # line looks like: "3 0 1 3"
        elseif len == 4
            push!(fcs, decompose(FaceType, Face{len, FaceEltype, -1}(line))...) # line looks like: "4 0 1 2 3"
        end
    end

    return MeshType(vts, fcs)
end


text:
# patches to metrics.jl of Distances.jl

const metrics = (SqEuclidean, Euclidean, Cityblock, Minkowski, Hamming, TotalVariation)
const UnionMetrics = Distances.UnionMetrics

# Before evaluation, unwrap the AbstractGray colorant and promote storage type
#
# (Gray{N0f8}, Gray{N0f8}) -> (Float32, Float32)
# (Gray{N0f8}, N0f8      ) -> (Float32, Float32)
# (N0f8      , Gray{N0f8}) -> (Float32, Float32)
# (N0f8      , N0f8      ) -> (Float32, Float32)
#
# For Color3, only promote storage type; basic operation such as `abs2`
# returns a `Number` for `Color3`(e.g., `RGB`)
#
# (RGB{N0f8}, RGB{N0f8}) -> (RGB{Float32}, RGB{Float32} )
#
# we don't extend `evaluate` here because it makes the dispatching rules too complicated
# we only need to extend `eval_op`, since other operations
# (e.g., `eval_start`, `eval_reduce`, and `eval_end`) will automatically get promoted
for M in metrics
    for (Ta, Tb) in ((AbstractGray, AbstractGray),
                     (AbstractGray, Number),
                     (Number, AbstractGray),
                     (PromoteType, PromoteType))
        @eval function Distances.eval_op(d::$M, a::$Ta, b::$Tb)
            T1 = eltype(floattype(typeof(a)))
            T2 = eltype(floattype(typeof(a)))
            Distances.eval_op(d, T1(a), T2(b))
        end
    end
    @eval function Distances.eval_op(d::$M,
                                     a::Color3{<:PromoteType},
                                     b::Color3{<:PromoteType})
        CT1 = floattype(typeof(a))
        CT2 = floattype(typeof(b))
        Distances.eval_op(d, CT1(a), CT2(b))
    end
end

# ambiguities
for (ATa, ATb) in ((AbstractGray, AbstractGray),
                   (AbstractGray, Number),
                   (Number, AbstractGray),
                   (PromoteType, PromoteType),
                   (Color3, Color3))
    @eval function result_type(dist::UnionMetrics, ::Type{Ta}, ::Type{Tb}) where {Ta <: $ATa,Tb <: $ATb}
        T1 = eltype(floattype(Ta))
        T2 = eltype(floattype(Tb))
        result_type(dist, T1, T2)
    end
end


text:
function add_born_sources!(issp::Int64, pac, pap)

	δx24I=pac.fc[:δx24I]; δz24I=pac.fc[:δz24I];
	δxI=pac.fc[:δxI]; δzI=pac.fc[:δzI];
	δt=pac.fc[:δt]
	δtI=pac.fc[:δtI]
	δmodKI=pac.δmod[:KI]; modK=pac.mod[:K];
	δmodrhovxI=pac.δmod[:rhovxI]; δmodrhovzI=pac.δmod[:rhovzI]
	nx=pac.ic[:nx]; nz=pac.ic[:nz]

	born_svalue_stack=pap[1].born_svalue_stack

	p1=pap[1].w1[:t][:p]
	p2=pap[2].w1[:t][:p]
	p1p=pap[1].w1[:tp][:p]
	p1pp=pap[1].w1[:tpp][:p]
	dpdx1=pap[1].w1[:dx][:p]
	dpdx2=pap[2].w1[:dx][:p]
	dpdz1=pap[1].w1[:dz][:p]
	dpdz2=pap[2].w1[:dz][:p]

	# secondary sources for Born modeling
	# adding born sources from pressure(:,:,1) to pressure(:,:,2)
	# upto until [it-2]
	# lambdaI scatterrer source term at [it-1]
	# p is at [it], pp is at [it-1], ppp is at [it-2]
	# dpdx is at [it-1] and dpdz is at [it-1]
	# modrhovxI scatterrer source term at [it-1]
	# modrhovzI scatterrer source term at [it-1]

	# compute strength of secondary sources to δmod
	born_stackKI!(born_svalue_stack,p1,p1p,p1pp,δmodKI,nx,nz,δtI,δt)
	born_stackrhovxI!(born_svalue_stack,dpdx1,δmodrhovxI,nx,nz,δx24I,δt)
	born_stackrhovzI!(born_svalue_stack,dpdz1,δmodrhovzI,nx,nz,δz24I,δt)
	born_stack!(p2,born_svalue_stack,modK,nx,nz,δt)
end
@inbounds @fastmath function born_stackKI!(born_svalue_stack,p1,p1p,p1pp,δmodKI,nx,nz,δtI,δt)
	for ix=1:nx
		@simd for iz=1:nz
			born_svalue_stack[iz,ix] += δt * ((-1.0 * (p1pp[iz,ix] + p1[iz,ix] - 2.0 * p1p[iz,ix]) * δmodKI[iz,ix] * δtI * δtI))
		end
	end
end
@inbounds @fastmath function born_stackrhovxI!(born_svalue_stack,dpdx1,δmodrhovxI,nx,nz,δx24I,δt)
	for ix=3:nx-1
		@simd for iz=1:nz
			born_svalue_stack[iz,ix] +=
				δt * ((27.e0*dpdx1[iz,ix] * δmodrhovxI[iz,ix] -27.e0*dpdx1[iz,ix-1] * δmodrhovxI[iz,ix-1] -dpdx1[iz,ix+1] * δmodrhovxI[iz,ix+1] +dpdx1[iz,ix-2] * δmodrhovxI[iz,ix-2] ) * (δx24I))
		end
	end

end
@inbounds @fastmath function born_stackrhovzI!(born_svalue_stack,dpdz1,δmodrhovzI,nx,nz,δz24I,δt)
	for ix=1:nx
		@simd for iz=3:nz-1
			born_svalue_stack[iz,ix] +=
				δt * ((27.e0*dpdz1[iz,ix] * δmodrhovzI[iz,ix] -27.e0*dpdz1[iz-1,ix] * δmodrhovzI[iz-1,ix] -dpdz1[iz+1,ix] * δmodrhovzI[iz+1,ix] +dpdz1[iz-2,ix] * δmodrhovzI[iz-2,ix] ) * (δz24I))
		end
	end
end
# add secondary sources to second pw
@inbounds @fastmath function born_stack!(p2,born_svalue_stack,modK,nx,nz,δt)
	for ix=1:nx
		@simd for iz=1:nz
			p2[iz,ix] += born_svalue_stack[iz,ix] * modK[iz,ix] * δt  #* δxI * δzI
		end
	end
end


text:
"""
Simplex in 2D.
PyDEC: Software and Algorithms for Discretization of Exterior Calculus
by Bell and Hirani.
"""
struct Simplex2D
    # vertex locations
    v::Array{Float64,2}
    # vertices that form triangles
    s::Array{Int64,2}
end

function buildsimplex2d(v, s)
    # The simplex has to be normalized.
    # The hilbert order increases locality of the geometry.
    reorder = hilbert_order(v)
    v = v[:, reorder]
    remap = sortperm(reorder)
    for j in 1:size(s, 2)
        for i in 1:size(s, 1)
            s[i, j] = remap[s[i, j]]
        end
    end
    # We will rely on simplex vertices being sorted.
    # circular permuations perserve orientation.
    for shift_idx in 1:size(s, 2)
        verts = s[:, shiftidx]
        s[:, shiftidx] = circshift(verts, findmin(verts))
    end
    s = sort(s, dims = 1)
    Simplex2D(v, s)
end

"""
Assume s is a list of triangles.
Pydec paper page 13. But transposed b/c this is column first ordering.
"""
function boundary_faces(s)
    d, n = size(s)  # dimension (3) and number of simplices
    splus = zeros(Int64, d + 2, n)
    splus[1:d, :] = s
    # need orientation of simplex for edges.
    for orient_idx in 1:d
        orient = iseven(countedinsertionsort!(splus[1:d, orient_idx]))
        splus[d + 1, orient_idx] = orient ? 1 : -1
        splus[d + 2, orient_idx] = orient_idx
    end
    resultant = zeros(Int64, d + 1, d * n)
    for vert in 1:d
        for tri in 1:n
            for copy_idx in 1:vert
                cnt = 1
                if copy_idx != vert
                    resultant[cnt, tri] = splus[copy_idx, tri]
                    cnt += 1
                end
                resultant[d, tri] = splus[d + 1, tri] * (-1)^(vert - 1)
                resultant[d + 1, tri] = splus[d + 2, tri]
            end
        end
    end
    ordered_resultant = sort(resultant, dims = 1)

    faces_pad = zeros(Int64, d - 1, size(resultant, 2))
    face_cnt = 1
    faces_pad[:, 1] = resultant[1:(d - 1), 1]
    for fidx in 2:size(resultant, 2)
        if faces_pad[1:(d - 1), fidx] != faces_pad[:, face_cnt]
            face_cnt += 1
            faces_pad[:, face_cnt] = resultant[1:(d - 1), fidx]
        end
    end
    faces = faces_pad[:, 1:face_cnt]

    edges = SparseMatrixCSC{Int64,Int64}(
        n,  # row count is number of simplices
        size(resultant, 2),  # column count is number of edges
        1:size(resultant, 2),  # one column per edge
        resultant[d + 1, :],  # row indices are the simplex index
        resultant[d, :]  # values are the orientation of the edge
        )
    edges, faces
end


text:
"""
  extractReads(reader::XAM.BAM.Reader, chrom::String, chromrange::UnitRange)

Extracts reads from a region as a FASTQ file. Note: This will NOT keep
track of pair information!
"""
function extractReads(reader::XAM.BAM.Reader, chrom::String,
  chromrange::UnitRange, outfile::String)
  open(FASTX.FASTQ.Writer, outfile) do w
    for record in XAM.BAM.GenomicFeatures.eachoverlap(reader, chrom, chromrange)
    	rec = FASTX.FASTQ.Record(XAM.BAM.tempname(record), XAM.BAM.sequence(record),
        XAM.BAM.quality(record))
        write(w, rec)
    end
  end
end
#TODO: outfile as gzipped file

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
baremodule StaticStorages

import UUIDs

struct BucketKey
    value::UUIDs.UUID
end

struct StorageKey
    value::UUIDs.UUID
end

function put! end
function get end
function getbucket end

module Implementations

using UUIDs: UUID

if VERSION >= v"1.6"
    using UUIDs: uuid4
else
    using Random: RandomDevice
    using UUIDs: UUIDs
    uuid4() = UUIDs.uuid4(RandomDevice())
end

using ..StaticStorages: StaticStorages, BucketKey, StorageKey

include("utils.jl")
include("implementations.jl")

end  # module Implementations

Implementations.define_docstrings()

end  # baremodule StaticStorages


text:
using BenchmarkTools

# Get the matrix dimensions N from the command line argument.
N = parse(Int, ARGS[1])

A = ones(N, N)
for i = 1:N
    A[i, i] = 6
end

println("--------------------------")
println(@sprintf "Matric square root: %d" N)
println("--------------------------")

# Take the square root of matrix A.
@btime begin
    B = sqrtm(A)
end


text:
function gather_check_dims(X::AbstractArray{Tx,Nx},
                           Y::AbstractArray{Ty,Ny},
                           idx::AbstractArray{Tidx,Nidx}) where
                           {Tx,Ty,Tidx<:IntOrIntTuple,Nx,Ny,Nidx}
    M = NNlib.typelength(Tidx)
    dims = gather_check_dims(Nx, Ny, M, Nidx)
    size(X)[1:dims] == size(Y)[1:dims] || throw(ArgumentError("Incompatible input shapes."))
    size(Y)[dims+1:end] == size(idx) || throw(ArgumentError("Incompatible input shapes."))
    return dims
end

function gather_check_dims(X::AbstractArray{Tx,Nx},
                           Y::AbstractArray{Ty,Ny},
                           idx::AbstractArray{CartesianIndex{M},Nidx}) where
                           {Tx,Ty,Nx,Ny,M,Nidx}
    dims = gather_check_dims(Nx, Ny, M, Nidx)
    size(X)[1:dims] == size(Y)[1:dims] || throw(ArgumentError("Incompatible input shapes."))
    size(Y)[dims+1:end] == size(idx) || throw(ArgumentError("Incompatible input shapes."))
    return dims
end

function gather_check_dims(Nx, Ny, M, Nidx)
    @assert Nx - M == Ny - Nidx "Incompatible input shapes of (dst, src, idx) = ($Nx, $Ny, $Nidx)."
    dims = Nx - M
    dims < 0 && throw(ArgumentError("dims must be non-negative but got dims=$dims."))
    return dims
end

function gather_kernel!(dst, src, idx, max_idx, max_dims_idx, dims_size)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        j, k = divrem(index-1, max_dims_idx)
        dims_i = CartesianIndices(dims_size)[k+1]
        dst[index] = src[dims_i, idx[j+1]...]
    end
    return nothing
end

function NNlib.gather!(dst::AnyCuArray, src::AnyCuArray, idx::AnyCuArray)
    dims = gather_check_dims(src, dst, idx)
    dims_size = size(src)[1:dims]
    max_dims_idx = prod(dims_size)
    max_idx = max_dims_idx * length(idx)
    args = dst, src, idx, max_idx, max_dims_idx, dims_size

    kernel = @cuda launch=false gather_kernel!(args...)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = min(max_idx, config.threads)
    blocks = cld(max_idx, threads)
    kernel(args...; threads=threads, blocks=blocks)
    return dst
end


text:
module TestStack

using Julio
using Julio.Internal: maybepop_nowait!
using Test

function test_pushpop()
    s = Julio.Internal.TreiberStack{Int}()
    push!(s, 111)
    push!(s, 222)
    @test something(maybepop_nowait!(s)) == 222
    @test something(maybepop_nowait!(s)) == 111
    @test maybepop_nowait!(s) === nothing
end

function test_iter()
    s = Julio.Internal.TreiberStack{Int}()
    push!(s, 111)
    push!(s, 222)
    @test collect(s) == [222, 111]
end

end  # module


text:
# AbstractOperators binding
# special cases
import Base: *, reshape

"""
`reshape(x::AbstractExpression, dims...)`

Reshapes the codomain of the expression.

# Example

```julia
julia> A,b = randn(10,3), randn(10);

julia> reshape(A*x-b,2,5)

```

"""
function reshape(a::AbstractExpression, dims...)
	A = convert(Expression,a)
	op = Reshape(A.L, dims...)
	return Expression{length(A.x)}(A.x,op)
end
#Reshape

imported = [
            :getindex :GetIndex;
            :exp      :Exp;
            :cos      :Cos;
            :sin      :Sin;
            :atan     :Atan;
            :tanh     :Tanh;
           ]

importedFFTW = [
            :fft      :(AbstractOperators.DFT);
            :rfft     :RDFT;
            :irfft    :IRDFT;
            :ifft     :IDFT;
            :dct      :DCT;
            :idct     :IDCT;
           ]

importedDSP = [
            :conv     :Conv;
            :xcorr    :Xcorr;
            :filt     :Filt;
           ]

exported = [
            :finitediff :FiniteDiff;
            :variation  :Variation;
            :mimofilt   :MIMOFilt;
            :zeropad    :ZeroPad;
            :sigmoid    :Sigmoid;
            :σ          :Sigmoid; #alias
            :pow        :Pow; #alias
           ]

#importing functions from Base
for f in  imported[:,1]
	@eval begin
		import Base: $f
	end
end
#importing functions from FFTW
for f in  importedFFTW[:,1]
	@eval begin
		import FFTW: $f
        export $f
	end
end
#importing functions from DSP
for f in  importedDSP[:,1]
	@eval begin
		import DSP: $f
		export $f
	end
end
#exporting functions
for f in  exported[:,1]
	@eval begin
		export $f
	end
end

fun = [imported; importedFFTW; importedDSP; exported]
for i = 1:size(fun,1)
	f,fAbsOp = fun[i,1],fun[i,2]
	@eval begin
		function $f(a::AbstractExpression, args...)
			A = convert(Expression,a)
			op = $fAbsOp(codomainType(operator(A)),size(operator(A),1), args...)
			return op*A
		end
	end
end

## docs

"""
`getindex(x::AbstractExpression, dims...)`

Slices the codomain of `ex`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> getindex(x,1:4)

julia> X = Variable(10,4)
Variable(Float64, (10, 4))

julia> X[:,1:3]

```
"""
getindex

"""
`fft(x::AbstractExpression)`

Applies the Fourier transform.
See documentation for `Base.fft` and `AbstractOperator.DFT`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = fft(x)

julia> operator(ex)
ℱ  ℝ^10 -> ℂ^10

```
"""
fft

"""
`rfft(x::AbstractExpression, [, dims] )`

Applies the Fourier transform exploiting the conjugate symmetry.
See documentation for `Base.rfft` and `AbstractOperator.RDFT`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = rfft(x)

julia> operator(ex)
ℱ  ℝ^10 -> ℂ^6

```
"""
rfft

"""
`ifft(x::AbstractExpression)`

Applies the inverse Fourier transform.
See documentation for `Base.ifft` and `AbstractOperator.IDFT`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = ifft(x)

julia> operator(ex)
ℱ⁻¹  ℝ^10 -> ℂ^10

```
"""
ifft

"""
`irfft(x::AbstractExpression, d, [, dims] )`

Applies the inverse Fourier transform exploiting the conjugate symmetry.

`d` must indicate the length of the real codomain.

See documentation for `Base.irfft` and `AbstractOperator.IRDFT`.

# Example

```julia
julia> x = Variable(Complex{Float64}, 10)
Variable(Float64, (10,))

julia> ex = irfft(x,19)

julia> operator(ex)
ℱ⁻¹  ℂ^10 -> ℝ^19

```
"""
irfft

"""
`dct(x::AbstractExpression)`

Applies the discrete cosine transform.
See documentation for `Base.dct` and `AbstractOperator.DCT`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = dct(x)

julia> operator(ex)
ℱc  ℝ^10 -> ℝ^10

```
"""
dct

"""
`idct(x::AbstractExpression)`

Applies the inverse discrete cosine transform.
See documentation for `Base.idct` and `AbstractOperator.IDCT`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = idct(x)

julia> operator(ex)
ℱc⁻¹  ℝ^10 -> ℝ^10

```
"""
idct

"""
`conv(x::AbstractExpression, h::AbstractVector)`

Perform discrete convolution with `h`.
See documentation for `Base.conv` and `AbstractOperator.Conv`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = conv(x, randn(5))

julia> operator(ex)
★  ℝ^10 -> ℝ^14

```
"""
conv

"""
`xcorr(x::AbstractExpression, h::AbstractVector)`

Performs the cross correlation of the codomain of `ex` with `h`.
See documentation for `Base.xcorr` and `AbstractOperator.Xcorr`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = xcorr(x, randn(5))

julia> operator(ex)
◎  ℝ^10 -> ℝ^19

```
"""
xcorr

"""
`filt(x::AbstractExpression, b::AbstractVector, [a::AbstractVector])`

Filter with the finite impulse response (FIR) `b`.
Alternatively infinite impulse responses filters (IIR) can be used as well by specifying the coefficients `a`.

See documentation for `Base.filt` and `AbstractOperator.Filt`.

# Example

```julia
julia> x = Variable(10)
Variable(Float64, (10,))

julia> ex = filt(x, [1.;0.;0.],[1.;0.;1.])

julia> operator(ex)
IIR  ℝ^10 -> ℝ^10

julia> ex = filt(x, [1.;0.;0.])

julia> operator(ex)
FIR  ℝ^10 -> ℝ^10

```
"""
filt

"""
`finitediff(X::AbstractExpression, dir = 1)`

Performs the discretized gradient over the specified direction `dir` obtained using forward finite differences.

# Example

```julia
julia> X = Variable(10,10)
Variable(Float64, (10,10))

julia> ex = finitediff(X)

julia> operator(ex)
δx  ℝ^(10, 10) -> ℝ^(9, 10)

julia> ex = finitediff(X,2)

julia> operator(ex)
δy  ℝ^(10, 10) -> ℝ^(10, 9)

```
"""
finitediff

"""
`variation(X::AbstractExpression)`

Returns the variation of `X` using forward finite differences.
Specifically if `X` is in `ℝ^(n, m)` the codomain of `variation(X)` will
consist of `ℝ^(n*m,2)`, having in the ``i``-th column the vectorized gradient over the ``i``-th direction.

See documentation of `AbstractOperator.Variation`.

# Example

```julia
julia> X = Variable(4,5,6)
Variable(Float64, (4,5,6))

julia> ex = variation(X)

julia> operator(ex)
Ʋ  ℝ^(4, 5, 6) -> ℝ^(120, 3)

```
"""
variation

"""
`mimofilt(X::AbstractExpression, B::Vector{AbstractVector}, [A::Vector{AbstractVector}])`

Multiple-input multiple-output filter: like `filt` but allows for multipule inputs.

```math
\\mathbf{y}_i = \\sum_{j = 1}^{M} \\mathbf{h}_{i,j} * \\mathbf{x}_j
```
where ``\\mathbf{y}_i`` and ``\\mathbf{x}_j`` are the ``i``-th and ``j``-th columns of the output `Y` and input `X` matrices respectively.

The filters ``\\mathbf{h}_{i,j}`` can be represented either by providing coefficients `B` and `A` (IIR) or `B` alone (FIR). These coefficients must be given in a `Vector` of `Vector`s.

For example for a `3` by `2` MIMO system (i.e. `size(X,2) == 3` inputs and `size(Y,2) == 2` outputs) `B` must be:

`B = [b11, b12, b13, b21, b22, b23]`

where `bij` are vector containing the filter coeffients of `h_{i,j}`.

See documentation of `AbstractOperator.MIMOFilt`.
"""
mimofilt

"""
`zeropad(x::AbstractExpression, zp::Tuple)`

Performs zeropadding:
```math
zp(\\mathbf{x}) = [\\mathbf{x}, 0, \\dots, 0  ]^T
```
where `zp` is a tuple containing the number of nonzero elements that are added in each dimension.

See documentation of `AbstractOperator.ZeroPad`.
"""
zeropad

"""
`sigmoid(x::AbstractExpression, γ = 1.0)`
`σ(x::AbstractExpression, γ = 1.0)`

Sigmoid function:
```math
\\sigma(\\mathbf{x}) = \\frac{1}{1+e^{-\\gamma \\mathbf{x} } }
```

See documentation of `AbstractOperator.Sigmoid`.
"""
sigmoid
σ

"""
`exp(x::AbstractExpression)`

Exponential function:
```math
e^{ \\mathbf{x}  }
```

See documentation of `AbstractOperator.Exp`.
"""
exp

"""
`sin(x::AbstractExpression)`

Sine function:
```math
\\sin( \\mathbf{x}  )
```

See documentation of `AbstractOperator.Sin`.
"""
sin

"""
`cos(x::AbstractExpression)`

Cosine function:
```math
\\cos( \\mathbf{x}  )
```

See documentation of `AbstractOperator.Cos`.
"""
cos

"""
`atan(x::AbstractExpression)`

Inverse tangent function:
```math
\\tan^{-1}( \\mathbf{x}  )
```

See documentation of `AbstractOperator.Atan`.
"""
atan

"""
`tanh(x::AbstractExpression)`

Hyperbolic tangent function:
```math
\\tanh ( \\mathbf{x}  )
```

See documentation of `AbstractOperator.Tanh`.
"""
tanh

"""
`pow(x::AbstractExpression, n)`

Elementwise power 'n' of 'x':
```math
x_i^{n} \\ \\forall \\  i = 0,1, \\dots
```

See documentation of `AbstractOperator.Pow`.
"""
pow


text:
## definition of the different types and constant

DEG2RAD =  π / 180.
Grav    = 4.30091e-3          ## Gravitational constant in km/s, pc, Msun

## Gaia data struct.
mutable struct Df
    ndata::Int32
    data::Array{Float64}
    raw::Array{Float64}
    err::Array{Float64}
    sourceid::Array{Int64}
end

## dbcan parameters
mutable struct model
    eps::Float64
    min_nei::Int
    min_cl::Int
end

## dbcan and weighting parameters
mutable struct modelfull
    eps::Float64
    min_nei::Int
    min_cl::Int
    w3d::Float64
    wvel::Float64
    whrd::Float64
end

#### mutable Markov Chains to get the stationary df for the ABC-MCMC/dbscan opt,
mutable struct mc
    eps::Array{Float64}
    mne::Array{Int32}
    mcl::Array{Int32}
    qc::Array{Float64}
    qn::Array{Int32}
end

#### mutable Markov Chains to get the stationary df for the ABC-MCMC/dbscan opt,
mutable struct mcfull
    eps::Array{Float64}
    mne::Array{Int32}
    mcl::Array{Int32}
    w3d::Array{Float64}
    wvel::Array{Float64}
    whrd::Array{Float64}
    qc::Array{Float64}
    qn::Array{Int32}
end

## parameters for the ABC MCMC optimization of the dbscan
struct abc
    minQ::Float64
    minstars::Int
    epsmean::Float64
    epsdisp::Float64
    min_nei::Int
    min_cl::Int
    ncoredisp::Int
    nburnout::Int
    niter::Int
end

## parameters for the ABC MCMC optimization of the dbscan
struct abcfull
    minQ::Float64
    minstars::Int
    forcedminstars::Int
    epsmean::Float64
    epsdisp::Float64
    min_nei::Int
    min_cl::Int
    ncoredisp::Int
    w3dmean::Float64
    w3ddisp::Float64
    wvelmean::Float64
    wveldisp::Float64
    whrdmean::Float64
    whrddisp::Float64
    nburnout::Int
    niter::Int
end

## Basic properties of a SC
struct SCproperties
    nstars::Int
    distance::Float64
    ra::Float64
    dec::Float64
    l::Float64
    b::Float64
    parallax::Float64
    pmra::Float64
    pmdec::Float64
    pml::Float64
    pmb::Float64
    vl::Float64
    vb::Float64
    vrad::Float64
    xdisp::Float64
    ydisp::Float64
    zdisp::Float64
    vldisp::Float64
    vbdisp::Float64
    vraddisp::Float64
end

## 2MASS  type
mutable struct twomass
    J::Array{Float64}
    Jerr::Array{Float64}
    H::Array{Float64}
    Herr::Array{Float64}
    K::Array{Float64}
    Kerr::Array{Float64}
end

### Spatial Structure Types
##
## position from the gravity centered
struct sc2dcentered
    nxy::Int
    xx::Array{Float64}
    yy::Array{Float64}
    radius::Array{Float64}
    dens::Array{Float64}
    nbin::Int
    radbin::Array{Float64}
    densbin::Array{Float64}
    densbinerr::Array{Float64}
end

## prior model / mcmc
mutable struct mcmcCauchy
    Cmean::Float64
    Cdisp::Float64
    smean::Float64
    sdisp::Float64
    mmean::Float64
    mdisp::Float64
    nburnin::Int
    niter::Int
end

mutable struct modelCauchy
    C::Float64
    s::Float64
    m::Float64
end

mutable struct mcCauchy
    C::Array{Float64}
    s::Array{Float64}
    m::Array{Float64}
end

### METADATA for standalone, setings, optimization, etc
mutable struct meta
    rootdir::String
    wdir::String
    votdir::String
    plotdir::String
    ocdir::String
    votname::String
    prefile::String
    uuid::UUID

    optim::String
    w3d::Float64
    wvel::Float64
    whrd::Float64
    eps::Float64
    mcl::Int
    mnei::Int
    maxdist::Float64
    mindist::Float64

    minQc::Float64
    minQn::Int
    maxQn::Int
    forcedminstars::Int
    mingoodsolution::Int
    niterqminq::Int
    nburnout::Int
    nchain::Int
    maxiter::Int

    cyclemax::Int
    minstarselection::Int
    minstarstop::Int
    minchainreached::Int
    qcminstop::Float64
    wratiominstop::Float64

    epsmean::Float64
    epsdisp::Float64
    min_nei::Int
    min_cl::Int
    ncoredisp::Int
    w3dmean::Float64
    w3ddisp::Float64
    wvelmean::Float64
    wveldisp::Float64
    whrdmean::Float64
    whrddisp::Float64

    eps_rw::Float64
    minnei_rw::Int
    mincl_rw::Int
    w3d_rw::Float64
    wvel_rw::Float64
    whrd_rw::Float64

    clustermax::Int
    aperture2d::Float64
    maxaperture2d::Float64
    aperture3d::Float64
    maxaperture3d::Float64
    aperturev::Float64
    maxaperturev::Float64
    nboot::Int
    labels::String

    meta()= new()
end

## Basic properties of a SC v2
mutable struct SCproperties2
    nstars::Int
    distance::Float64
    ra::Float64
    dec::Float64
    l::Float64
    b::Float64
    parallax::Float64
    pmra::Float64
    pmdec::Float64
    pml::Float64
    pmb::Float64
    vl::Float64
    vb::Float64
    vrad::Float64
    xdisp::Float64
    ydisp::Float64
    zdisp::Float64
    vldisp::Float64
    vbdisp::Float64
    vraddisp::Float64
    offdeg::Float64          ## offset in degree wrt. field center
    edgratg::Float64         ## edge ratio of the center
    edgratm::Float64         ## maximum edge ratio

    SCproperties2()= new()
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
using Behavior
using CoffeeMachine

@given("a machine without coffee") do context
    context[:machine] = Machine(coffee=0.0, milk=5.0)
end

@given("a machine without milk") do context
    context[:machine] = Machine(coffee=5.0, milk=0.0)
end

@then("the machine displays {String}") do context, message
    m = context[:machine]
    @expect readdisplay(m) == message
end


text:
module IntrusiveCollections

export TaggedIntrusiveList, IntrusiveList, IntrusiveListNode
export TaggedIntrusiveSList, IntrusiveSList, IntrusiveSListNode
export deleteafter!, insertafter!

include("api.jl")
include("utility.jl")
include("list_node.jl")
include("list.jl")
include("slist_node.jl")
include("slist.jl")

end


text:
using PlanarMaps: PlanarMap, EmbeddedMap
using NetworkLayout: Spring, Stress, SFDP

function generate_locs(g::PlanarMultigraph; he_start = nothing)
    g_simp, hes_splitted = simple_connected_planar_graph(g)
    g_norm, v_map, he_map, f_map = normalize(g_simp)
    he0 = (he_start === nothing ? surrounding_half_edge(g_norm, f_map[0]) : he_map[he_start])

    adjlist = [[dst(g_norm, he) for he in trace_vertex(g_norm, v)] for v in 1:nv(g_norm)]
    p = PlanarMap(adjlist)
    e = EmbeddedMap(p; outeredge = (src(g_norm, he0), dst(g_norm, he0)))
    locs_norm = e.locs
    v_locs = Dict(v => locs_norm[v_norm] for (v, v_norm) in v_map)
    he_locs = Dict(he => locs_norm[v_map[v_he]] for (he, v_he) in hes_splitted)

    N = nv(g_simp)
    hes_map = Dict(he => v_map[v] for (he, v) in hes_splitted)
    vs_map = Dict(v => v_map[v] for v in vertices(g))
    new_vs = collect(1:N)
    new_es = []
    for he in half_edges(g)
        if haskey(hes_map, he)
            push!(new_es, (vs_map[src(g, he)], hes_map[he]))
            if src(g, he) == dst(g, he) # self-loop
                (v1, v2) = (hes_map[he], hes_map[twin(g, he)])
                v1 < v2 && push!(new_es, (v1, v2))
            end
        else
            s = src(g, he)
            d = dst(g, he)
            s <= d && (push!(new_es, (vs_map[s], vs_map[d])))
        end
    end

    return v_locs, he_locs, vs_map, hes_map, new_vs, new_es
end

function plot_graph(vs, es, locs;
        vlabel = nothing, vsize = nothing, vfill = nothing,
        vstroke = nothing, vrots = nothing,
        tfill = nothing, tsize = nothing, kwargs...)
    N = length(vs)
    v_map = Dict(zip(vs, 1:N))
    v_map_inv = Dict(zip(1:N, vs))
    adjmat = zeros(Bool, N, N)
    for e in es
        adjmat[v_map[e[1]], v_map[e[2]]] = true
        adjmat[v_map[e[2]], v_map[e[1]]] = true
    end
    vlabel !== nothing && (vlabel = [vlabel[v_map_inv[v]] for v in 1:N])
    vlabel === nothing && (vlabel = ["$(v_map_inv[v])" for v in 1:N])
    vsize !== nothing && (vsize = [vsize[v_map_inv[v]] for v in 1:N])
    vfill !== nothing && (vfill = [vfill[v_map_inv[v]] for v in 1:N])
    vstroke !== nothing && (vstroke = [vstroke[v_map_inv[v]] for v in 1:N])
    tfill !== nothing && (tfill = [tfill[v_map_inv[v]] for v in 1:N])
    tsize !== nothing && (tsize = [tsize[v_map_inv[v]] for v in 1:N])
    vrots !== nothing && (vrots = [vrots[v_map_inv[v]] for v in 1:N])
    return plot_graph(adjmat, [locs[v_map[v]] for v in vs];
        vlabel = vlabel, vsize = vsize, vfill = vfill, vstroke = vstroke, vrots = vrots,
        tfill = tfill, tsize = tsize, kwargs...)
end

function plot_graph(adjmat, locs; padding = 0.2,
        vlabel = nothing, vsize = nothing, vfill = nothing,
        vstroke = nothing, vrots = nothing,
        tfill = nothing, tsize = nothing,
        layout = nothing)
    layout === :spring && (locs = Spring(; initialpos = locs, C = 0.2)(adjmat))
    layout === :stress && (locs = Stress(; initialpos = locs)(adjmat))
    layout === :sfdp && (locs = SFDP(; initialpos = locs, K = 0.01)(adjmat))

    N = size(adjmat, 1)
    x_locs = [loc[1] for loc in locs]
    y_locs = -[loc[2] for loc in locs]
    x_min = minimum(x_locs)
    x_max = maximum(x_locs)
    y_min = minimum(y_locs)
    y_max = maximum(y_locs)
    vlabel === nothing && (vlabel = ["$v" for v in 1:N])
    vsize === nothing && (vsize = [0.04 for _ in 1:N])
    vfill === nothing && (vfill = ["white" for _ in 1:N])
    vstroke === nothing && (vstroke = ["black" for _ in 1:N])
    tfill === nothing && (tfill = ["red" for _ = 1:N])
    tsize === nothing && (tsize = [8pt for _ = 1:N])
    vrots === nothing && (vrots = [0.0 for v in 1:N])
    line_locs = []
    for i = 1:N
        for j = i:N
            adjmat[i, j] && push!(line_locs, [(x_locs[i], y_locs[i]), (x_locs[j], y_locs[j])])
        end
    end

    return compose(context(units=UnitBox(x_min-padding, y_min-padding, x_max-x_min+2padding, y_max-y_min+2padding)),
        ((context(), text(x_locs[v], y_locs[v], vlabel[v], hcenter, vcenter, Rotation(vrots[v], x_locs[v], y_locs[v])),
            fill(tfill[v]), fontsize(tsize[v])) for v in 1:N)...,
        (context(), circle(x_locs, y_locs, vsize),
            fill(vfill), stroke(vstroke)),
        (context(), line(line_locs), stroke("gray")),
        (context(), rectangle(), fill("white"))
    )
end

function plot_planar(g::PlanarMultigraph; kwargs...)
    v_locs, he_locs, v_map, hes_map, new_vs, new_es = generate_locs(g)
    locs = Dict()
    for v in vertices(g)
        locs[v_map[v]] = v_locs[v]
    end
    for (he, v) in hes_map
        locs[v] = he_locs[he]
    end
    vlabel = Dict(v1 => "$v0" for (v0, v1) in v_map)
    return plot_graph(new_vs, new_es, locs; vlabel = vlabel, kwargs...)
end

function plot_planar(q::Tait; scale = 1, kwargs...)
    if length(q.inputs) > 0
        hes0 = trace_vertex(q, q.inputs[1])
        he_start = hes0[findfirst(he -> (face(q, he) == 0), hes0)]
        he_start = prev(q, he_start)
    else
        he_start = nothing
    end
    v_locs, he_locs, v_map, hes_map, new_vs, new_es = generate_locs(q.g; he_start = he_start)
    locs = Dict()
    vsize = Dict()
    vfill = Dict()
    vstroke = Dict()
    vrots = Dict()
    tfill = Dict()
    tsize = Dict()
    for v in vertices(q)
        locs[v_map[v]] = (v_locs[v][1] * scale, v_locs[v][2] * scale)
        vsize[v_map[v]] = is_open_vertex(q, v) ? 0.06 : 0.04
        vfill[v_map[v]] = is_open_vertex(q, v) ? "white" : (is_genus(q, v) ? "white" : "royalblue")
        vstroke[v_map[v]] = is_open_vertex(q, v) ? "white" : "royalblue"
        vrots[v_map[v]] = 0.0
        tfill[v_map[v]] = "red"
        tsize[v_map[v]] = 8/scale*pt
    end
    vlabel = Dict(v1 => "$v0" for (v0, v1) in v_map)
    for v in q.inputs
        vlabel[v_map[v]] = "$(v)ᵢ"
    end
    for v in q.outputs
        vlabel[v_map[v]] = "$(v)ₒ"
    end
    self_loops = Set{Int}()
    for (he, v) in hes_map
        if src(q, he) != dst(q, he) # not a self-loop
            locs[v] = (he_locs[he][1]*scale, he_locs[he][2]*scale)
            v_src = v_map[src(q, he)]
            v_dst = v_map[dst(q, he)]
            vrots[v] = angle(-(locs[v_src][1] - locs[v_dst][1]) + im * (locs[v_src][2] - locs[v_dst][2]))
            if is_open_half_edge(q, he)
                vfill[v] = "transparent"
                vstroke[v] = "transparent"
                vsize[v] = 0.04/scale
                vlabel[v] = "-$he→\n←$(twin(q, he))-"
            else
                vfill[v] = "gray"
                vstroke[v] = "gray"
                vsize[v] = 0.04/scale
                vlabel[v] = "-$he→\n$(quon_param(q, he))\n←$(twin(q, he))-"
            end
            tfill[v] = "black"
            tsize[v] = 5/scale*pt
        else
            he1 = he
            he2 = twin(q, he)
            if !(he1 in self_loops)
                push!(self_loops, he1, he2)
                continue
            end
            v1 = v
            v2 = hes_map[he2]
            locs[v1] = (he_locs[he1][1]*scale, he_locs[he1][2]*scale)
            locs[v2] = (he_locs[he2][1]*scale, he_locs[he2][2]*scale)
            vfill[v1] = "gray"
            vfill[v2] = "gray"
            vstroke[v1] = "gray"
            vstroke[v2] = "gray"
            vrots[v1] = angle(-(locs[v1][1] - locs[v2][1]) + im * (locs[v1][2] - locs[v2][2]))
            vrots[v2] = 0.0
            if is_open_half_edge(q, he1)
                vfill[v1] = "transparent"
                vfill[v2] = "transparent"
                vstroke[v1] = "transparent"
                vstroke[v2] = "transparent"
                vsize[v1] = 0.04/scale
                vsize[v2] = 0.04/scale
                vlabel[v1] = "-$he→\n←$(twin(q, he))-"
                vlabel[v2] = ""
            else
                vfill[v1] = "gray"
                vfill[v2] = "gray"
                vstroke[v1] = "gray"
                vstroke[v2] = "gray"
                vsize[v1] = 0.04/scale
                vlabel[v1] = "-$he→\n$(quon_param(q, he))\n←$(twin(q, he))-"
                vsize[v2] = 0
                vlabel[v2] = ""
            end
            tfill[v1] = "black"
            tfill[v2] = "black"
            tsize[v1] = 5/scale*pt
            tsize[v2] = 5/scale*pt
        end
    end

    hes_recorded = Set(keys(hes_map))
    v_max = length(new_vs)
    for he in half_edges(q)
        (he in hes_recorded || twin(q, he) in hes_recorded) && continue
        push!(hes_recorded, he)
        v_max += 1
        push!(new_vs, v_max)
        v1 = v_map[src(q, he)]
        v2 = v_map[dst(q, he)]
        push!(new_es, (v_max, v1))
        push!(new_es, (v_max, v2))
        locs[v_max] = ((locs[v1][1] + locs[v2][1])/2, (locs[v1][2] + locs[v2][2])/2)
        tfill[v_max] = "black"
        tsize[v_max] = 5/scale*pt
        vrots[v_max] = angle(-(locs[v1][1] - locs[v2][1]) + im * (locs[v1][2] - locs[v2][2]))
        if is_open_half_edge(q, he)
            vfill[v_max] = "transparent"
            vstroke[v_max] = "transparent"
            vsize[v_max] = 0.04/scale
            vlabel[v_max] = "-$(he)→\n←$(twin(q, he))-"
        else
            vfill[v_max] = "gray"
            vstroke[v_max] = "gray"
            vsize[v_max] = 0.04/scale
            vlabel[v_max] = "-$(he)→\n$(quon_param(q, he))\n←$(twin(q, he))-"
        end
    end

    return plot_graph(new_vs, new_es, locs;
        vlabel = vlabel, vsize = vsize, vfill = vfill,
        vstroke = vstroke, vrots = vrots,
        tfill = tfill, tsize = tsize, kwargs...)
end


text:
"""
    Equation(R, K)

Construct an Equation that represents a PDE of the form

    ∂ϕ/∂t = ∂/∂z (K ∂ϕ/∂z) + R ,

so that K is a diffusivity and R captures 'remaining' terms.

In OceanTurb, K and R are tuples of functions.
The functions K[j](m, i) and R[j](m, i) calculate the
diffusivity K and remaining terms R for the model `m` at
grid point `i`.
"""
struct Equation{TA, TL, TK, TR, U}
    A       :: TA
    L       :: TL
    K       :: TK
    R       :: TR
    update! :: U
end

funkynothing(args...) = nothing

function Equation(;
         K,
         A = Tuple(minuszero for i=1:length(K)),
         L = Tuple(minuszero for i=1:length(K)),
         R = Tuple(minuszero for i=1:length(K)),
    update = funkynothing
    )
    return Equation(A, L, K, R, update)
end


text:
"""
   Abstract type representing F-Localizations.
"""
abstract type FLocalization end

"""
   Abstract type representing a fitted F-Localization
   (i.e., wherein the F-localization has already been determined by data).
"""
abstract type FittedFLocalization end

# Holy trait
abstract type FLocalizationVexity end
struct LinearVexity <: FLocalizationVexity end
struct ConvexVexity <: FLocalizationVexity end

StatsBase.fit(floc::FittedFLocalization, args...; kwargs...) = floc

function nominal_alpha(floc::FLocalization)
    floc.α
end

function flocalization_constraint!(model, floc, prior::PriorVariable)
    model
end

"""
    DvoretzkyKieferWolfowitz(;α = 0.05, max_constraints = 1000) <: FLocalization

The Dvoretzky-Kiefer-Wolfowitz band (based on the Kolmogorov-Smirnov distance)
at confidence level `1-α` that bounds the distance of the true distribution function
to the ECDF ``\\widehat{F}_n`` based on ``n`` samples. The constant of the band is the sharp
constant derived by Massart:

```math
F \\text{ distribution}:  \\sup_{t \\in \\mathbb R}\\lvert F(t) - \\widehat{F}_n(t) \\rvert  \\leq  \\sqrt{\\log(2/\\alpha)/(2n)}
```
The supremum above is enforced discretely on at most `max_constraints` number of points.
"""
Base.@kwdef struct DvoretzkyKieferWolfowitz{T,N} <: FLocalization
    α::T = 0.05
    max_constraints::N = 1000
end

# for backwards compatibility
DvoretzkyKieferWolfowitz(α) = DvoretzkyKieferWolfowitz(;α=α)

vexity(::DvoretzkyKieferWolfowitz) = LinearVexity()

struct FittedDvoretzkyKieferWolfowitz{T,S,D<:AbstractDict{T,S},DKW} <:
        FittedFLocalization
    summary::D
    band::S
    dkw::DKW
    homoskedastic::Bool
end

vexity(dkw::FittedDvoretzkyKieferWolfowitz) = vexity(dkw.dkw)

function nominal_alpha(dkw::FittedDvoretzkyKieferWolfowitz)
    nominal_alpha(dkw.dkw)
end

# TODO: Allow this to work more broadly.
function StatsBase.fit(dkw::DvoretzkyKieferWolfowitz, Zs::AbstractVector{<:EBayesSample})
    StatsBase.fit(dkw, summarize(Zs))
end

function StatsBase.fit(dkw::DvoretzkyKieferWolfowitz, Zs_summary::MultinomialSummary{T}) where T
    α = nominal_alpha(dkw)
    n = nobs(Zs_summary)

    if skedasticity(Zs_summary) === Heteroskedastic()
        homoskedastic = false
        multiplier = exp(1)
        Zs_summary = compound(Zs_summary)
    elseif skedasticity(Zs_summary) === Homoskedastic()
        if T<:CompoundSample
            homoskedastic = false
            multiplier = exp(1)
        else
            homoskedastic = true
            multiplier = 1
        end
    end
    n_constraints = length(Zs_summary)

    band = sqrt(log(2 *multiplier / α) / (2n))

    cdf_probs = cumsum([v for (k, v) in Zs_summary.store]) #SORTED important here
    cdf_probs /= cdf_probs[end]
    _Zs = collect(keys(Zs_summary.store))

    issorted(_Zs) || error("MultinomialSummary not sorted.")

    max_constraints = dkw.max_constraints

    if max_constraints < n_constraints - 10
        _step = div(n_constraints-2, max_constraints)
        idxs = [1; 2:_step:(n_constraints-1); n_constraints]

        _Zs = _Zs[idxs]
        cdf_probs = cdf_probs[idxs]
    end
    # TODO: report issue with SortedDict upstream
    _dict = SortedDict(Dict(_Zs .=> cdf_probs))

    FittedDvoretzkyKieferWolfowitz(_dict, band, dkw, homoskedastic)
end

function flocalization_constraint!(
    model,
    dkw::FittedDvoretzkyKieferWolfowitz,
    prior::PriorVariable,
)
    band = dkw.band
    for (Z, cdf_value) in dkw.summary
        marginal_cdf = cdf(prior, Z::EBayesSample)
        if cdf_value + band < 1
            @constraint(model, marginal_cdf <= cdf_value + band)
        end
        if cdf_value - band > 0
            @constraint(model, marginal_cdf >= cdf_value - band)
        end
    end
    model
end

# add discrete version of this?
@recipe function f(fitted_dkw::FittedDvoretzkyKieferWolfowitz)
    x_dkw = response.(collect(keys(fitted_dkw.summary)))
    F_dkw = collect(values(fitted_dkw.summary))

    band = fitted_dkw.band
    lower = max.(F_dkw .- band, 0.0)
    upper = min.(F_dkw .+ band, 1.0)
    cis_ribbon  = F_dkw .- lower, upper .- F_dkw
    fillalpha --> 0.36
    legend --> :topleft
    seriescolor --> "#018AC4"
    ribbon --> cis_ribbon

    x_dkw, F_dkw
end


text:
using Printf, FileIO, JLD2

export ParticleGroup

abstract type AbstractParticleGroup end

"""
    ParticleGroup{D,V}( n_particles, charge, mass, q, weights)

- `n_particles` : number of particles
- `charge`      : charge of the particle species
- `mass`        : mass of the particle species
- `n_weights`   : number of differents weights
"""
struct ParticleGroup{D,V} <: AbstractParticleGroup
    dims::Tuple{Int,Int}
    n_particles::Int
    array::Array{Float64,2}
    common_weight::Float64
    charge::Float64
    mass::Float64
    n_weights::Int
    q_over_m::Float64
    n_spin::Int

    function ParticleGroup{D,V}(
        n_particles; charge=1.0, mass=1.0, n_weights=1, common_weight=0.0, n_spin=0
    ) where {D,V}
        dims = (D, V)
        array = zeros(Float64, (sum(dims) + n_weights + n_spin, n_particles))
        if common_weight == 0.0
            common_weight = 1.0 / n_particles
        end
        q_over_m = charge / mass

        return new(
            dims,
            n_particles,
            array,
            common_weight,
            charge,
            mass,
            n_weights,
            q_over_m,
            n_spin,
        )
    end
end

"""
    get_x( p, i )

Get position of ith particle of p
"""
@inline get_x(p::ParticleGroup{D,V}, i::Int) where {D,V} = p.array[1:D, i]

"""
    get_v( p, i )

Get velocity of ith particle of p
"""
@inline get_v(p::ParticleGroup{D,V}, i::Int) where {D,V} = p.array[(D + 1):(D + V), i]

"""
    get_charge( p, i; i_wi=1)

Get charge of ith particle of p (q * particle_weight)
"""
@inline function get_charge(p::ParticleGroup{D,V}, i::Int; i_wi=1) where {D,V}
    return p.charge * p.array[D + V + i_wi, i] * p.common_weight
end

"""
    get_mass( p, i; i_wi=1)

Get mass of ith particle of p (m * particle_weight)
"""
@inline function get_mass(p::ParticleGroup{D,V}, i::Int; i_wi=1) where {D,V}
    return p.mass * p.array[D + V + i_wi, i] * p.common_weight
end

"""
    get_weights( p, i)

Get ith particle weights of group p
"""
@inline function get_weights(p::ParticleGroup{D,V}, i::Int) where {D,V}
    return p.array[(D + V + 1):(D + V + p.n_weights), i]
end

"""
    get_spin( p, i, j)

Get the jth weight of the ith particle weights of group p
"""
@inline function get_spin(p::ParticleGroup{D,V}, i::Int, j::Int) where {D,V}
    return p.array[D + V + p.n_weights + j, i]
end

"""
    set_spin!( p, i, j)

Set the jth weight of the ith particle weights of group p
"""
@inline function set_spin!(p::ParticleGroup{D,V}, i::Int, j::Int, s) where {D,V}
    return p.array[D + V + p.n_weights + j, i] = s
end

"""
    set_x!( p, i, x )

Set position of ith particle of p to x
"""
@inline function set_x!(p::ParticleGroup{D,V}, i::Int, x::Vector{Float64}) where {D,V}
    for j in 1:D
        p.array[j, i] = x[j]
    end
end

"""
    set_x!( p, i, x)

Set position of ith particle of p to x

!!! note
    if `x` is a scalar value, only the first x dimension will be set.
"""
@inline function set_x!(p::ParticleGroup{D,V}, i::Int, x::Float64) where {D,V}
    return p.array[1, i] = x
end

"""
    set_v!( p, i, v)

Set velocity of ith particle of p to v
"""
@inline function set_v!(p::ParticleGroup{D,V}, i::Int, v::Vector{Float64}) where {D,V}
    for j in 1:V
        p.array[D + j, i] = v[j]
    end
end

"""
    set_v!( p, i, v)

Set velocity of ith particle of p to v
"""
@inline function set_v!(p::ParticleGroup{D,V}, i::Int, v::Float64) where {D,V}
    return p.array[D + 1, i] = v
end

"""
    set_weights!( p, i, w)

Set weights of ith particle of p to w
"""
function set_weights!(p::ParticleGroup{D,V}, i::Int, w::Vector{Float64}) where {D,V}
    for j in 1:(p.n_weights)
        p.array[D + V + j, i] = w[j]
    end
end

"""
    set_weights!( p, i, w)

Set weights of particle @ i
"""
function set_weights!(p::ParticleGroup{D,V}, i::Int, w::Float64) where {D,V}
    return p.array[D + V + 1, i] = w
end

function save(file, step, p::ParticleGroup{D,V}) where {D,V}
    datafile = @sprintf("%s-%06d.jld2", file, step)

    if p.n_spin > 0
        FileIO.save(
            datafile,
            Dict(
                "x" => p.array[1:D, :],
                "v" => p.array[(D + 1):(D + V), :],
                "w" => p.array[(D + V + 1):(D + V + p.n_weights), :],
                "s" => p.array[(D + V + p.n_weights + 1):end, :],
            ),
        )

    else
        FileIO.save(
            datafile,
            Dict(
                "x" => p.array[1:D, :],
                "v" => p.array[(D + 1):(D + V), :],
                "w" => p.array[(D + V + 1):end, :],
            ),
        )
    end
end


text:
# © Przemyslaw Szufel, 2018, run under Julia 1.0.0

struct ModelConfig
    # computational grid size
    GRID_SIZE :: Tuple{Int,Int}
    #data matrix size on a single process (all matrices are squares)
    dim_i::Int
    dim_j::Int
    #probabilities for empty, red and blue cells
    probs_w::Vector{Float64}
    #minimal required number of neighbours
    min_neighbours::Int
    slaves::Vector{Int64}
    master_id::Int64
end

mutable struct NodeData
    grid_i :: Int
    grid_j :: Int
    field :: Matrix{Int8}
    free_cells::Vector{Tuple{Int,Int}}
    rng::MersenneTwister
end

struct BorderData
    left :: Vector{Int8}
    right :: Vector{Int8}
    top :: Vector{Int8}
    bottom :: Vector{Int8}
end

struct CornerData
    topleft :: Int8
    topright :: Int8
    bottomleft :: Int8
    bottomright :: Int8
end

#The below const statements could be an alternative to barrier functions used
#const modelC  = ModelConfig((0,0),0,0,Vector{Float64}(),0)
#const node = NodeData(0,0,Matrix{Int8}(0,0),Vector{Tuple{Int,Int}}(),MersenneTwister(0),(0,0))
#const worker_grid = Matrix{Int}(0,0)
function get_field_interior():: Matrix{Int8}
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    field = zeros(Int8,(dim_i,dim_j))
    for i = (1+1):(dim_i+1), j = (1+1):(dim_j+1)
        field[i-1,j-1]=node.field[i,j]
    end
    return field
end

function sample(rng::AbstractRNG, weights::Vector{Float64})::Int
    t = rand(rng) * sum(weights)
    n = length(weights)
    i = 1
    cw = weights[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += weights[i]
    end
    return i
end

function init_node()::Bool
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    values = Int8(0):Int8(length(modelC.probs_w)-1)
    node.field = zeros(Int8,(1+dim_i+1,1+dim_j+1))
    node.free_cells=Vector{Tuple{Int,Int}}()
    for i = (1+1):(dim_i+1), j = (1+1):(dim_j+1)
       val = values[sample(node.rng,modelC.probs_w)]
       if val == 0
           push!(node.free_cells,(i,j))
       end
       node.field[i,j] = val
    end
    return true
end

function get_border()::BorderData
    return get_border(modelC, node)
end
function get_border(modelC :: ModelConfig, node ::NodeData)::BorderData
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    #Please note that 4 bytes are redundant, however this increases code readibility
    return BorderData(
        node.field[(1+1):(1+dim_i),1+1],
        node.field[(1+1):(1+dim_i),1+dim_j],
        node.field[1+1,(1+1):(1+dim_j)],
        node.field[1+dim_i,(1+1):(1+dim_j)] )
end

function get_border_top()::Vector{Int8}
    return get_border_top(modelC, node)
end
function get_border_bottom()::Vector{Int8}
    return get_border_bottom(modelC, node)
end
function get_border_left()::Vector{Int8}
    return get_border_left(modelC, node)
end
function get_border_right()::Vector{Int8}
    return get_border_right(modelC, node)
end

function get_border_top(modelC :: ModelConfig, node ::NodeData)::Vector{Int8}
    return node.field[(1+1),(1+1):(1+modelC.dim_j)]
end
function get_border_bottom(modelC :: ModelConfig, node ::NodeData)::Vector{Int8}
    return node.field[(1+modelC.dim_i),(1+1):(1+modelC.dim_j)]
end
function get_border_left(modelC :: ModelConfig, node ::NodeData)::Vector{Int8}
    return node.field[(1+1):(1+modelC.dim_i),(1+1)]
end
function get_border_right(modelC :: ModelConfig, node ::NodeData)::Vector{Int8}
    return node.field[(1+1):(1+modelC.dim_i),(1+modelC.dim_j)]
end

function get_corner_top_left()::Int8
    return get_corner_top_left(modelC, node)
end
function get_corner_top_right()::Int8
    return get_corner_top_right(modelC, node)
end
function get_corner_bottom_left()::Int8
    return get_corner_bottom_left(modelC, node)
end
function get_corner_bottom_right()::Int8
    return get_corner_bottom_right(modelC, node)
end

function get_corner_top_left(modelC :: ModelConfig, node ::NodeData)::Int8
    return node.field[1+1,1+1]
end
function get_corner_top_right(modelC :: ModelConfig, node ::NodeData)::Int8
    return node.field[1+1,1+modelC.dim_j]
end
function get_corner_bottom_left(modelC :: ModelConfig, node ::NodeData)::Int8
    return node.field[1+modelC.dim_i,1+1]
end
function get_corner_bottom_right(modelC :: ModelConfig, node ::NodeData)::Int8
    return node.field[1+modelC.dim_i,1+modelC.dim_j]
end

function step_setborder_nodes!()::Float64
    return step_setborder_nodes!(modelC, node)
end
function step_setborder_nodes!(modelC :: ModelConfig, node ::NodeData)::Bool
    start_t = time()
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    #borders
    for j = 1:dim_j
        node.field[1,j+1]=border_data.top[j]
        node.field[1+dim_i+1,j+1]=border_data.bottom[j]
    end
    for i = 1:dim_i
        node.field[i+1,1]=border_data.left[i]
        node.field[i+1,1+dim_j+1]=border_data.right[i]
    end
    node.field[1,1] = corner_data.topleft
    node.field[1,1+dim_j+1] = corner_data.topright
    node.field[1+dim_i+1,1] = corner_data.bottomleft
    node.field[1+dim_i+1,1+dim_j+1] = corner_data.bottomright
    return time()-start_t
end

function step_setborder_nodes_p2p!()::Float64
    return step_setborder_nodes_p2p!(modelC, node, worker_grid)
end
function step_setborder_nodes_p2p!(modelC :: ModelConfig, node ::NodeData, worker_grid::Matrix{Int})::Float64
    start_t = time()
    j_left = node.grid_j-1
    if j_left < 1; j_left = modelC.GRID_SIZE[2]; end
    j_right = node.grid_j+1
    if j_right > modelC.GRID_SIZE[2]; j_right = 1; end
    i_top = node.grid_i-1
    if i_top < 1; i_top = modelC.GRID_SIZE[1]; end
    i_bottom = node.grid_i+1
    if i_bottom > modelC.GRID_SIZE[1]; i_bottom = 1; end
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    #borders
    @sync begin
        #top
        @async node.field[1,(1+1):(dim_j+1)] = remotecall_fetch(get_border_bottom,worker_grid[i_top,node.grid_j])
        #bottom
        @async node.field[(1+dim_i+1),(1+1):(dim_j+1)] = remotecall_fetch(get_border_top,worker_grid[i_bottom,node.grid_j])
        #left
        @async node.field[(1+1):(dim_i+1),1] = remotecall_fetch(get_border_right,worker_grid[node.grid_i,j_left])
        #right
        @async node.field[(1+1):(dim_i+1),(1+dim_j+1)] = remotecall_fetch(get_border_left,worker_grid[node.grid_i,j_right])
        #top-left
        @async node.field[1,1] = remotecall_fetch(get_corner_bottom_right,worker_grid[i_top,j_left])
        #top-right
        @async node.field[1,(1+dim_j+1)] = remotecall_fetch(get_corner_bottom_left,worker_grid[i_top,j_right])
        #bottom-left
        @async node.field[(1+dim_i+1),1] = remotecall_fetch(get_corner_top_right,worker_grid[i_bottom,j_left])
        #bottom-right
        @async node.field[(1+dim_i+1),(1+dim_j+1)] = remotecall_fetch(get_corner_top_left,worker_grid[i_bottom,j_right])
    end
    return time()-start_t
end

struct LeaveLocData
    relocate_counts :: Vector{Int}
    free_cells_count :: Int
    elapsed_time :: Float64
end

function step_leave_location!()::LeaveLocData
    return step_leave_location!(modelC, node)::LeaveLocData
end
function step_leave_location!(modelC :: ModelConfig, node ::NodeData)::LeaveLocData
    start_t = time()
    field = node.field
    free_count = length(node.free_cells)
    dim_i = modelC.dim_i
    dim_j = modelC.dim_j
    ds = ((1,0),  (0,-1), (0,1), (-1,0), (1,1),  (1,-1), (-1,1), (-1,-1))
    for i = (1+1):(dim_i+1), j = (1+1):(dim_j+1)
        if field[i,j] > 0
            c = 0
            for (d1, d2) ∈ ds
                c += field[i,j] == field[i+d1,j+d2]
            end
            #print("c=",c," i=",i," j=",j," f=",field[i-1:i+1,j-1:j+1],"\r\n")
            if c < modelC.min_neighbours
                push!(node.free_cells,(i,j))
            end
        end
    end
    relocate_counts::Vector{Int} = zeros(Int,length(modelC.probs_w)-1)
    for (i, j) in node.free_cells[free_count+1:length(node.free_cells)]
        relocate_counts[field[i,j]] +=1
        field[i,j]=0
    end
    return LeaveLocData(relocate_counts,length(node.free_cells),time()-start_t)
end

function step_set_new_location!()::Float64
    return step_set_new_location!(modelC, node)
end
function step_set_new_location!(modelC :: ModelConfig, node ::NodeData)::Float64
    start_t = time()
    shuffle!(node.rng,node.free_cells)
    ix = length(node.free_cells)
    for z = 1:length(new_citizens)
        for a = 1:new_citizens[z]
            node.field[node.free_cells[ix]...]=z
            ix -=1
            pop!(node.free_cells)
        end
    end
    return time()-start_t
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "4 Trig functions" begin
    include("4.1 Sine/.jl")
    include("4.2 Cosine/.jl")
    include("4.3 Tangent/.jl")
    include("4.4 Cotangent/.jl")
    include("4.5 Secant/.jl")
    include("4.6 Cosecant/.jl")
    include("4.7 Miscellaneous/.jl")
end

max_stars_repo_path: test/4 Trig functions/.jl

max_issues_repo_path: test/4 Trig functions/.jl

max_forks_repo_path: test/4 Trig functions/.jl


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/print.jl
# Testing $fa pretty-printing-related functionality
#############################################################################
using JuMP, FactCheck
import JuMP.REPLMode, JuMP.IJuliaMode
import JuMP.repl, JuMP.ijulia

# Helper function to test IO methods work correctly
fname = VERSION < v"0.5.0-dev+4340" ? :writemime : :show
@eval function io_test(mode, obj, exp_str; repl=:both)
    if mode == REPLMode
        repl != :show  && @fact sprint(print, obj) --> exp_str
        repl != :print && @fact sprint(show,  obj) --> exp_str
    else
        @fact sprint($fname, "text/latex", obj) --> string("\$\$ ",exp_str," \$\$")
    end
end

facts("[print] JuMPContainer{Variable}") do
    le, ge, fa = repl[:leq], repl[:geq], repl[:for_all]
    inset, dots = repl[:in], repl[:dots]
    infty, union = repl[:infty], repl[:union]

    m = Model()

    #------------------------------------------------------------------
    # Test bound printing
    context("bound printing") do
    @variable(m,      bnd_free[2:5])
    @variable(m,      bnd_lowb[2:5] >= 2)
    @variable(m,      bnd_high[2:5] <= 5)
    @variable(m, 2 <= bnd_both[2:5] <= 5)
    @variable(m,      bnd_difflo[i=2:5] >= i)
    @variable(m,      bnd_diffup[i=2:5] <= i)
    @variable(m, i <= bnd_diffbo[i=2:5] <= 2i)
    @variable(m, i <= bnd_difflo_with_up[i=2:5] <= 5)
    @variable(m, 2 <= bnd_diffup_with_lo[i=2:5] <= i)

    io_test(REPLMode, bnd_free, "bnd_free[i] free $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_lowb, "bnd_lowb[i] $ge 2 $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_high, "bnd_high[i] $le 5 $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_both, "2 $le bnd_both[i] $le 5 $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_difflo, "bnd_difflo[i] $ge $dots $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_diffup, "bnd_diffup[i] $le $dots $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_diffbo, "$dots $le bnd_diffbo[i] $le $dots $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_difflo_with_up, "$dots $le bnd_difflo_with_up[i] $le 5 $fa i $inset {2,3,4,5}")
    io_test(REPLMode, bnd_diffup_with_lo, "2 $le bnd_diffup_with_lo[i] $le $dots $fa i $inset {2,3,4,5}")

    io_test(IJuliaMode, bnd_free, "bnd_free_{i} free \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_lowb, "bnd_lowb_{i} \\geq 2 \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_high, "bnd_high_{i} \\leq 5 \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_both, "2 \\leq bnd_both_{i} \\leq 5 \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_difflo, "bnd_difflo_{i} \\geq \\dots \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_diffup, "bnd_diffup_{i} \\leq \\dots \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_diffbo, "\\dots \\leq bnd_diffbo_{i} \\leq \\dots \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_difflo_with_up, "\\dots \\leq bnd_difflo_with_up_{i} \\leq 5 \\quad\\forall i \\in \\{2,3,4,5\\}")
    io_test(IJuliaMode, bnd_diffup_with_lo, "2 \\leq bnd_diffup_with_lo_{i} \\leq \\dots \\quad\\forall i \\in \\{2,3,4,5\\}")
    end

    #------------------------------------------------------------------
    # Test index set printing
    context("index set printing") do
    @variable(m, rng_unit1[1:10])  # Array{Variable}
    @variable(m, rng_unit2[-2:3])  # JuMPArray
    @variable(m, rng_unit3[[1:10;]])  # JuMPDict
    @variable(m, rng_step1[1:2:10])
    @variable(m, rng_step2[-2:5:10])
    @variable(m, rng_step3[1:5:3])
    @variable(m, rng_step4[0:2:2])
    @variable(m, arr_1[[:a,:b,:c]])
    @variable(m, arr_2[[:a,1,"test"]])
    @variable(m, arr_3[[:apple,:banana,:carrot,:diamonds]])
    @variable(m, rng2_1[1:10,[:a,:b,:c]])
    @variable(m, tri_1[i=1:3,j=i:3])
    @variable(m, tri_2[i=1:3,j=-i])
    @variable(m, tri_3[(i,j)=[(i,i+2) for i in 1:5],k=i:j])
    @variable(m, iter_1[keys(Dict(:a => 1))])

    io_test(REPLMode, rng_unit1, "rng_unit1[i] free $fa i $inset {1,2,$dots,9,10}")
    io_test(REPLMode, rng_unit2, "rng_unit2[i] free $fa i $inset {-2,-1,$dots,2,3}")
    io_test(REPLMode, rng_unit3, "rng_unit3[i] free $fa i $inset {1,2,$dots,9,10}")
    io_test(REPLMode, rng_step1, "rng_step1[i] free $fa i $inset {1,3,$dots,7,9}")
    io_test(REPLMode, rng_step2, "rng_step2[i] free $fa i $inset {-2,3,8}")
    io_test(REPLMode, rng_step3, "rng_step3[i] free $fa i $inset {1}")
    io_test(REPLMode, rng_step4, "rng_step4[i] free $fa i $inset {0,2}")
    io_test(REPLMode, arr_1, "arr_1[i] free $fa i $inset {a,b,c}")
    io_test(REPLMode, arr_2, "arr_2[i] free $fa i $inset {a,1,test}")
    io_test(REPLMode, arr_3, "arr_3[i] free $fa i $inset {apple,banana,carrot,diamonds}")
    io_test(REPLMode, rng2_1, "rng2_1[i,j] free $fa i $inset {1,2,$dots,9,10}, j $inset {a,b,c}")
    io_test(REPLMode, tri_1, "tri_1[i,j] free $fa i $inset {1,2,3}, j $inset {$dots}")
    io_test(REPLMode, tri_2, "tri_2[i,j] free $fa i $inset {1,2,3}, j $inset {$dots}")
    io_test(REPLMode, tri_3, "tri_3[(i,j),k] free $fa (i,j) $inset {(1,3),(2,4),$dots,(4,6),(5,7)}, k $inset {$dots}")
    io_test(REPLMode, iter_1, "iter_1[i] free $fa i $inset {a}")

    io_test(IJuliaMode, rng_unit1, "rng_unit1_{i} free \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}")
    io_test(IJuliaMode, rng_unit2, "rng_unit2_{i} free \\quad\\forall i \\in \\{-2,-1,\\dots,2,3\\}")
    io_test(IJuliaMode, rng_unit3, "rng_unit3_{i} free \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}")
    io_test(IJuliaMode, rng_step1, "rng_step1_{i} free \\quad\\forall i \\in \\{1,3,\\dots,7,9\\}")
    io_test(IJuliaMode, rng_step2, "rng_step2_{i} free \\quad\\forall i \\in \\{-2,3,8\\}")
    io_test(IJuliaMode, rng_step3, "rng_step3_{i} free \\quad\\forall i \\in \\{1\\}")
    io_test(IJuliaMode, rng_step4, "rng_step4_{i} free \\quad\\forall i \\in \\{0,2\\}")
    io_test(IJuliaMode, arr_1, "arr_1_{i} free \\quad\\forall i \\in \\{a,b,c\\}")
    io_test(IJuliaMode, arr_2, "arr_2_{i} free \\quad\\forall i \\in \\{a,1,test\\}")
    io_test(IJuliaMode, arr_3, "arr_3_{i} free \\quad\\forall i \\in \\{apple,banana,carrot,diamonds\\}")
    io_test(IJuliaMode, rng2_1, "rng2_1_{i,j} free \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}, j \\in \\{a,b,c\\}")
    io_test(IJuliaMode, tri_1, "tri_1_{i,j} free \\quad\\forall i \\in \\{1,2,3\\}, j \\in \\{\\dots\\}")
    io_test(IJuliaMode, tri_2, "tri_2_{i,j} free \\quad\\forall i \\in \\{1,2,3\\}, j \\in \\{\\dots\\}")
    io_test(IJuliaMode, tri_3, "tri_3_{(i,j),k} free \\quad\\forall (i,j) \\in \\{(1,3),(2,4),\\dots,(4,6),(5,7)\\}, k \\in \\{\\dots\\}")
    io_test(IJuliaMode, iter_1, "iter_1_{i} free \\quad\\forall i \\in \\{a\\}")
    end

    #------------------------------------------------------------------
    # Test category printing
    context("category printing") do
    @variable(m, cat_bin[1:3], Bin)
    @variable(m, 2 <= cat_int[1:3] <= 5, Int)
    @variable(m, cat_semiint_both[2:3] >= 2, SemiInt)
    @variable(m, i <= cat_semiint_difflow[i=2:3] <= 4, SemiInt)
    @variable(m, 2 <= cat_semiint_diffup[i=2:3] <= i, SemiInt)
    @variable(m, i <= cat_semiint_none[i=2:3] <= 2i, SemiInt)
    @variable(m, cat_semicont_both[2:3] >= 2, SemiCont)
    @variable(m, i <= cat_semicont_difflow[i=2:3] <= 4, SemiCont)
    @variable(m, 2 <= cat_semicont_diffup[i=2:3] <= i, SemiCont)
    @variable(m, i <= cat_semicont_none[i=2:3] <= 2i, SemiCont)
    @variable(m, fixed_var[i=2:3] == i)

    io_test(REPLMode, cat_bin, "cat_bin[i] $inset {0,1} $fa i $inset {1,2,3}")
    io_test(REPLMode, cat_int, "2 $le cat_int[i] $le 5, integer, $fa i $inset {1,2,3}")
    io_test(REPLMode, cat_semiint_both, "cat_semiint_both[i] $inset {2,$dots,$infty} $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semiint_difflow, "cat_semiint_difflow[i] $inset {$dots,$dots,4} $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semiint_diffup, "cat_semiint_diffup[i] $inset {2,$dots,$dots} $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semiint_none, "cat_semiint_none[i] $inset {$dots,$dots,$dots} $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semicont_both, "cat_semicont_both[i] $inset [2,$infty] $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semicont_difflow, "cat_semicont_difflow[i] $inset [$dots,4] $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semicont_diffup, "cat_semicont_diffup[i] $inset [2,$dots] $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, cat_semicont_none, "cat_semicont_none[i] $inset [$dots,$dots] $union {0} $fa i $inset {2,3}")
    io_test(REPLMode, fixed_var, "fixed_var[i] = $dots $fa i $inset {2,3}")

    io_test(IJuliaMode, cat_bin, "cat_bin_{i} \\in \\{0,1\\} \\quad\\forall i \\in \\{1,2,3\\}")
    io_test(IJuliaMode, cat_int, "2 \\leq cat_int_{i} \\leq 5, \\in \\mathbb{Z}, \\quad\\forall i \\in \\{1,2,3\\}")
    io_test(IJuliaMode, cat_semiint_both, "cat_semiint_both_{i} \\in \\{2,\\dots,\\infty\\} \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semiint_difflow, "cat_semiint_difflow_{i} \\in \\{\\dots,\\dots,4\\} \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semiint_diffup, "cat_semiint_diffup_{i} \\in \\{2,\\dots,\\dots\\} \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semiint_none, "cat_semiint_none_{i} \\in \\{\\dots,\\dots,\\dots\\} \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semicont_both, "cat_semicont_both_{i} \\in \\[2,\\infty\\] \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semicont_difflow, "cat_semicont_difflow_{i} \\in \\[\\dots,4\\] \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semicont_diffup, "cat_semicont_diffup_{i} \\in \\[2,\\dots\\] \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, cat_semicont_none, "cat_semicont_none_{i} \\in \\[\\dots,\\dots\\] \\cup \\{0\\} \\quad\\forall i \\in \\{2,3\\}")
    io_test(IJuliaMode, fixed_var, "fixed_var_{i} = \\dots \\quad\\forall i \\in \\{2,3\\}")
    end

    #------------------------------------------------------------------
    # Tests for particular issues
    context("Empty JuMPContainer printing (#124)") do
    @variable(m, empty_free[1:0])
    io_test(REPLMode, empty_free, "Empty Array{Variable} (no indices)")
    io_test(IJuliaMode, empty_free, "Empty Array{Variable} (no indices)")
    @variable(m, empty_set[[]])
    io_test(REPLMode, empty_set, "empty_set (no indices)")
    io_test(IJuliaMode, empty_set, "empty_set (no indices)")
    end
end

facts("[print] JuMPContainer{Number}") do
    # The same output for REPL and IJulia, so only testing one
    mod = Model()
    @variable(mod, i*j <= w[i=9:10, [:Apple,5,:Banana], j=-1:+1] <= i*j)
    @variable(mod, i*j*k <= x[i=9:11,j=99:101,k=3:4] <= i*j*k)
    @variable(mod, i*j <= y[i=9:11,j=i:11] <= i*j)
    @variable(mod, j <= z[i=[:a,'b'],j=1:3] <= j)
    solve(mod)

    # Deal with hashing variations
    if hash(5) < hash(:Apple)
        io_test(REPLMode, getvalue(w), """
w: 3 dimensions:
[ 9,:,:]
  [ 9, Apple,:]
    [ 9, Apple,-1] = -9.0
    [ 9, Apple, 0] = 0.0
    [ 9, Apple, 1] = 9.0
  [ 9,     5,:]
    [ 9,     5,-1] = -9.0
    [ 9,     5, 0] = 0.0
    [ 9,     5, 1] = 9.0
  [ 9,Banana,:]
    [ 9,Banana,-1] = -9.0
    [ 9,Banana, 0] = 0.0
    [ 9,Banana, 1] = 9.0
[10,:,:]
  [10, Apple,:]
    [10, Apple,-1] = -10.0
    [10, Apple, 0] = 0.0
    [10, Apple, 1] = 10.0
  [10,     5,:]
    [10,     5,-1] = -10.0
    [10,     5, 0] = 0.0
    [10,     5, 1] = 10.0
  [10,Banana,:]
    [10,Banana,-1] = -10.0
    [10,Banana, 0] = 0.0
    [10,Banana, 1] = 10.0
""", repl=:print)
    else
        io_test(REPLMode, getvalue(w), """
w: 3 dimensions:
[ 9,:,:]
  [ 9, Apple,:]
    [ 9, Apple,-1] = -9.0
    [ 9, Apple, 0] = 0.0
    [ 9, Apple, 1] = 9.0
  [ 9,     5,:]
    [ 9,     5,-1] = -9.0
    [ 9,     5, 0] = 0.0
    [ 9,     5, 1] = 9.0
  [ 9,Banana,:]
    [ 9,Banana,-1] = -9.0
    [ 9,Banana, 0] = 0.0
    [ 9,Banana, 1] = 9.0
[10,:,:]
  [10, Apple,:]
    [10, Apple,-1] = -10.0
    [10, Apple, 0] = 0.0
    [10, Apple, 1] = 10.0
  [10,     5,:]
    [10,     5,-1] = -10.0
    [10,     5, 0] = 0.0
    [10,     5, 1] = 10.0
  [10,Banana,:]
    [10,Banana,-1] = -10.0
    [10,Banana, 0] = 0.0
    [10,Banana, 1] = 10.0
""", repl=:print)
    end

    io_test(REPLMode, getvalue(x), """
x: 3 dimensions:
[ 9,:,:]
  [ 9, 99,:]
    [ 9, 99,3] = 2673.0
    [ 9, 99,4] = 3564.0
  [ 9,100,:]
    [ 9,100,3] = 2700.0
    [ 9,100,4] = 3600.0
  [ 9,101,:]
    [ 9,101,3] = 2727.0
    [ 9,101,4] = 3636.0
[10,:,:]
  [10, 99,:]
    [10, 99,3] = 2970.0
    [10, 99,4] = 3960.0
  [10,100,:]
    [10,100,3] = 3000.0
    [10,100,4] = 4000.0
  [10,101,:]
    [10,101,3] = 3030.0
    [10,101,4] = 4040.0
[11,:,:]
  [11, 99,:]
    [11, 99,3] = 3267.0
    [11, 99,4] = 4356.0
  [11,100,:]
    [11,100,3] = 3300.0
    [11,100,4] = 4400.0
  [11,101,:]
    [11,101,3] = 3333.0
    [11,101,4] = 4444.0
""", repl=:print)

    io_test(REPLMode, getvalue(y), """
y: 2 dimensions, 6 entries:
 [ 9, 9] = 81.0
 [ 9,10] = 90.0
 [ 9,11] = 99.0
 [10,10] = 100.0
 [10,11] = 110.0
 [11,11] = 121.0""")

    io_test(REPLMode, getvalue(z), """
z: 2 dimensions:
[a,:]
  [a,1] = 1.0
  [a,2] = 2.0
  [a,3] = 3.0
[b,:]
  [b,1] = 1.0
  [b,2] = 2.0
  [b,3] = 3.0
""")

end

facts("[print] SOS constraints") do
    modS = Model()
    a = [1,2,3]
    @variable(modS, x[1:3], Bin)
    addSOS1(modS, [a[i]x[i] for i in 1:3])
    s1 = JuMP.SOSConstraint([x[i] for i in 1:3],
                            [a[i] for i in 1:3], :SOS1)
    io_test(REPLMode, s1, "SOS1: {1 x[1], 2 x[2], 3 x[3]}")
    io_test(IJuliaMode, s1, "SOS1: \\{1 x[1], 2 x[2], 3 x[3]\\}")

    b = [5,4,7,2,1]
    @variable(modS, y[1:5], Bin)
    s2 = JuMP.SOSConstraint([y[i] for i in 1:5],
                            [b[i] for i in 1:5], :SOS2)
    io_test(REPLMode, s2, "SOS2: {5 y[1], 4 y[2], 7 y[3], 2 y[4], 1 y[5]}")
    io_test(IJuliaMode, s2, "SOS2: \\{5 y[1], 4 y[2], 7 y[3], 2 y[4], 1 y[5]\\}")
end

facts("[print] Model") do
    le, ge, fa = repl[:leq], repl[:geq], repl[:for_all]
    inset, dots = repl[:in], repl[:dots]
    infty, union = repl[:infty], repl[:union]
    Vert, sub2 = repl[:Vert], repl[:sub2]

    #------------------------------------------------------------------

    mod_1 = Model()
    @variable(mod_1, a>=1)
    @variable(mod_1, b<=1)
    @variable(mod_1, -1<=c<=1)
    @variable(mod_1, a1>=1,Int)
    @variable(mod_1, b1<=1,Int)
    @variable(mod_1, -1<=c1<=1,Int)
    @variable(mod_1, x, Bin)
    @variable(mod_1, y)
    @variable(mod_1, z, Int)
    @variable(mod_1, sos[1:3], Bin)
    @variable(mod_1, 2 <= si <= 3, SemiInt)
    @variable(mod_1, 2 <= sc <= 3, SemiCont)
    @variable(mod_1, fi == 9)
    @objective(mod_1, Max, a - b + 2a1 - 10x)
    @constraint(mod_1, a + b - 10c - 2x + c1 <= 1)
    @constraint(mod_1, a*b <= 2)
    addSOS1(mod_1, [i*sos[i] for i in 1:3])
    @constraint(mod_1, norm(sos) + a <= 1)

    io_test(REPLMode, mod_1, """
Max a - b + 2 a1 - 10 x
Subject to
 a + b - 10 c - 2 x + c1 $le 1
 a*b - 2 $le 0
 SOS1: {1 sos[1], 2 sos[2], 3 sos[3]}
 $(Vert)[sos[1],sos[2],sos[3]]$(Vert)$(sub2) $le -a + 1
 sos[i] $inset {0,1} $fa i $inset {1,2,3}
 a $ge 1
 b $le 1
 -1 $le c $le 1
 a1 $ge 1, integer
 b1 $le 1, integer
 -1 $le c1 $le 1, integer
 x $inset {0,1}
 y free
 z free, integer
 si $inset {2,$dots,3} $union {0}
 sc $inset [2,3] $union {0}
 fi = 9
""", repl=:print)

    io_test(REPLMode, mod_1, """
Maximization problem with:
 * 1 linear constraint
 * 1 quadratic constraint
 * 1 SOS constraint
 * 1 SOC constraint
 * 15 variables: 4 binary, 4 integer, 1 semicontinuous, 1 semi-integer
Solver is default solver""", repl=:show)

    io_test(IJuliaMode, mod_1, """
\\begin{alignat*}{1}\\max\\quad & a - b + 2 a1 - 10 x\\\\
\\text{Subject to} \\quad & a + b - 10 c - 2 x + c1 \\leq 1\\\\
 & a\\times b - 2 \\leq 0\\\\
 & SOS1: \\{1 sos[1], 2 sos[2], 3 sos[3]\\}\\\\
 & \\Vert[sos_{1},sos_{2},sos_{3}]\\Vert_2 $le -a + 1\\\\
 & sos_{i} \\in \\{0,1\\} \\quad\\forall i \\in \\{1,2,3\\}\\\\
 & a \\geq 1\\\\
 & b \\leq 1\\\\
 & -1 \\leq c \\leq 1\\\\
 & a1 \\geq 1, \\in \\mathbb{Z}\\\\
 & b1 \\leq 1, \\in \\mathbb{Z}\\\\
 & -1 \\leq c1 \\leq 1, \\in \\mathbb{Z}\\\\
 & x \\in \\{0,1\\}\\\\
 & y free\\\\
 & z free, \\in \\mathbb{Z}\\\\
 & si \\in \\{2,\\dots,3\\} \\cup \\{0\\}\\\\
 & sc \\in \\[2,3\\] \\cup \\{0\\}\\\\
 & fi = 9\\\\
\\end{alignat*}
""")

    #------------------------------------------------------------------

    mod_2 = Model()
    @variable(mod_2, x, Bin)
    @variable(mod_2, y, Int)
    @constraint(mod_2, x*y <= 1)

    io_test(REPLMode, mod_2, """
Feasibility problem with:
 * 0 linear constraints
 * 1 quadratic constraint
 * 2 variables: 1 binary, 1 integer
Solver is default solver""", repl=:show)

    mod_2 = Model()
    @variable(mod_2, x)
    @constraint(mod_2, x <= 3)

    io_test(REPLMode, mod_2, """
Feasibility problem with:
 * 1 linear constraint
 * 1 variable
Solver is default solver""", repl=:show)

    #------------------------------------------------------------------

    mod_3 = Model()

    @variable(mod_3, x[1:5])
    @NLconstraint(mod_3, x[1]*x[2] == 1)
    @NLconstraint(mod_3, x[3]*x[4] == 1)
    @NLconstraint(mod_3, x[5]*x[1] == 1)
    @NLobjective(mod_3, Min, x[1]*x[3])

    io_test(REPLMode, mod_3, """
Min (nonlinear expression)
Subject to
 3 nonlinear constraints
 x[i] free $fa i $inset {1,2,$dots,4,5}
""", repl=:print)
    io_test(REPLMode, mod_3, """
Minimization problem with:
 * 0 linear constraints
 * 3 nonlinear constraints
 * 5 variables
Solver is default solver""", repl=:show)
    io_test(IJuliaMode, mod_3, """
\\begin{alignat*}{1}\\min\\quad & (nonlinear expression)\\\\
\\text{Subject to} \\quad & 3 nonlinear constraints\\\\
 & x_{i} free \\quad\\forall i \\in \\{1,2,\\dots,4,5\\}\\\\
\\end{alignat*}
""", repl=:print)
end

facts("[print] changing variable categories") do
    le, ge, fa = repl[:leq], repl[:geq], repl[:for_all]
    inset, dots = repl[:in], repl[:dots]
    infty, union = repl[:infty], repl[:union]

    mod = Model()
    @variable(mod, x[1:3])
    @variable(mod, y[i=1:3,i:3])
    setcategory(x[3], :SemiCont)
    setcategory(y[1,3], :Int)

    io_test(REPLMode, mod, """
Min 0
Subject to
 x[1] free
 x[2] free
 x[3] $inset [-$infty,$infty] $union {0}
 y[1,1] free
 y[1,2] free
 y[1,3] free, integer
 y[2,2] free
 y[2,3] free
 y[3,3] free
""", repl=:print)

    io_test(IJuliaMode, mod, """
\\begin{alignat*}{1}\\min\\quad & 0\\\\
\\text{Subject to} \\quad & x_{1} free\\\\
 & x_{2} free\\\\
 & x_{3} \\in \\[-\\infty,\\infty\\] \\cup \\{0\\}\\\\
 & y_{1,1} free\\\\
 & y_{1,2} free\\\\
 & y_{1,3} free, \\in \\mathbb{Z}\\\\
 & y_{2,2} free\\\\
 & y_{2,3} free\\\\
 & y_{3,3} free\\\\
\\end{alignat*}
""")
end

facts("[print] expressions") do
    # Most of the expression logic is well covered by test/operator.jl
    # This is really just to check IJulia printing for expressions
    le, ge = repl[:leq], repl[:geq]

    #------------------------------------------------------------------
    mod = Model()
    @variable(mod, x[1:5])
    @variable(mod, y[i=2:4,j=i:5])
    @variable(mod, z)

    @constraint(mod, x[1] + 2*y[2,3] <= 3)
    io_test(REPLMode, mod.linconstr[end], "x[1] + 2 y[2,3] $le 3")
    io_test(IJuliaMode, mod.linconstr[end], "x_{1} + 2 y_{2,3} \\leq 3")

    @constraint(mod, (x[1]+x[2])*(y[2,2]+3.0) <= 1)
    io_test(REPLMode, mod.quadconstr[end], "x[1]*y[2,2] + x[2]*y[2,2] + 3 x[1] + 3 x[2] - 1 $le 0")
    io_test(IJuliaMode, mod.quadconstr[end], "x_{1}\\times y_{2,2} + x_{2}\\times y_{2,2} + 3 x_{1} + 3 x_{2} - 1 \\leq 0")

    @constraint(mod, (y[2,2]+3.0)*(x[1]+x[2]) <= 1)
    io_test(REPLMode, mod.quadconstr[end], "x[1]*y[2,2] + x[2]*y[2,2] + 3 x[1] + 3 x[2] - 1 $le 0")
    io_test(IJuliaMode, mod.quadconstr[end], "x_{1}\\times y_{2,2} + x_{2}\\times y_{2,2} + 3 x_{1} + 3 x_{2} - 1 \\leq 0")
end

facts("[print] Variable") do
    m = Model()
    @variable(m, 0 <= x <= 2, inconstraints=ConstraintRef{Model,LinearConstraint}[], objective=0.0, coefficients=Float64[] )

    @fact    getname(x) --> "x"
    io_test(REPLMode,   x, "x")
    io_test(IJuliaMode, x, "x")

    setname(x, "x2")
    @fact    getname(x) --> "x2"
    io_test(REPLMode,   x, "x2")
    io_test(IJuliaMode, x, "x2")

    setname(x, "")
    @fact    getname(x) --> "col_1"
    io_test(REPLMode,   x, "col_1")
    io_test(IJuliaMode, x, "col_1")

    @variable(m, z[1:2,3:5])
    @fact       getname(z[1,3]) --> "z[1,3]"
    io_test(REPLMode,   z[1,3],    "z[1,3]")
    io_test(IJuliaMode, z[1,3],    "z_{1,3}")
    @fact       getname(z[2,4]) --> "z[2,4]"
    io_test(REPLMode,   z[2,4],    "z[2,4]")
    io_test(IJuliaMode, z[2,4],    "z_{2,4}")
    @fact       getname(z[2,5]) --> "z[2,5]"
    io_test(REPLMode,   z[2,5],    "z[2,5]")
    io_test(IJuliaMode, z[2,5],    "z_{2,5}")

    @variable(m, w[3:9,["red","blue","green"]])
    @fact    getname(w[7,"green"]) --> "w[7,green]"
    io_test(REPLMode,   w[7,"green"], "w[7,green]")
    io_test(IJuliaMode, w[7,"green"], "w_{7,green}")

    rng = 2:5
    @variable(m, v[rng,rng,rng,rng,rng,rng,rng])
    a_v = v[4,5,2,3,2,2,4]
    @fact    getname(a_v) --> "v[4,5,2,3,2,2,4]"
    io_test(REPLMode,   a_v, "v[4,5,2,3,2,2,4]")
    io_test(IJuliaMode, a_v, "v_{4,5,2,3,2,2,4}")
end

facts("[print] User-created Array{Variable}") do
    m = Model()
    @variable(m, x)
    @variable(m, y)

    v = [x,y,x]
    A = [x y; y x]
    if VERSION >= v"0.5.0-"
        io_test(REPLMode,   v, "JuMP.Variable[x,y,x]")
        io_test(IJuliaMode, v, "JuMP.Variable[x,y,x]")
    else
        io_test(REPLMode,   v, "[x,y,x]")
        io_test(IJuliaMode, v, "[x,y,x]")
    end

    if VERSION >= v"0.5.0-dev+3642"
        io_test(REPLMode,   A, "JuMP.Variable[x y; y x]")
        io_test(IJuliaMode, A, "JuMP.Variable[x y; y x]")
    else
        io_test(REPLMode,   A, """
2x2 Array{JuMP.Variable,2}:
 x  y
 y  x""")
        io_test(IJuliaMode, A, """
2x2 Array{JuMP.Variable,2}:
 x  y
 y  x""")
    end
end

facts("[print] basename keyword argument") do
    m = Model()
    @variable(m, x, basename="foo")
    @variable(m, y[1:3], basename=:bar)
    num = 123
    @variable(m, z[[:red,:blue]], basename="color_$num")
    @variable(m, v[1:2,1:2], SDP, basename=string("i","$num",num))
    @variable(m, w[1:3,1:3], Symmetric, basename="symm")

    io_test(REPLMode,   x, "foo")
    io_test(IJuliaMode, x, "foo")
    io_test(REPLMode,   y[2], "bar[2]")
    io_test(IJuliaMode, y[2], "bar_{2}")
    io_test(REPLMode,   z[:red], "color_123[red]")
    io_test(IJuliaMode, z[:red], "color_123_{red}")
    io_test(REPLMode,   v[2,1], "i123123[1,2]")
    io_test(IJuliaMode, v[2,1], "i123123_{1,2}")
    io_test(REPLMode,   w[1,3], "symm[1,3]")
    io_test(IJuliaMode, w[1,3], "symm_{1,3}")
end


text:
# fallbacks for inbands_getindex and inbands_setindex!
@inline function inbands_getindex(x::AbstractArray, i::Integer, j::Integer)
    @inbounds r = getindex(x, i, j)
    r
end
@inline function inbands_setindex!(x::AbstractArray, v, i::Integer, j::Integer)
    @inbounds r = setindex!(x, v, i, j)
    r
end

inbands_getindex(x::Adjoint, i::Integer, j::Integer) =
    inbands_getindex(parent(x), j, i)'
inbands_getindex(x::Transpose, i::Integer, j::Integer) =
    transpose(inbands_getindex(parent(x), j, i))
inbands_setindex!(x::Adjoint, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), v', j, i)
inbands_setindex!(x::Transpose, v, i::Integer, j::Integer) =
    inbands_setindex!(parent(x), transpose(v), j, i)

###
# Lazy getindex
# this uses a lazy-materialize idiom to construct a matrix based
# on the memory layout
###

sub_materialize(::AbstractBandedLayout, V) = BandedMatrix(V)

@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::Colon) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::Colon) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractVector, jr::AbstractVector) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::AbstractVector) = lazy_getindex(A, kr, jr)
@inline getindex(A::AbstractBandedMatrix, kr::AbstractVector, jr::Colon) = lazy_getindex(A, kr, jr)

@inline getindex(A::AbstractMatrix, b::Band) = lazy_getindex(A, b)
@inline getindex(A::AbstractMatrix, kr::BandRangeType, j::Integer) = lazy_getindex(A, kr, j)
@inline getindex(A::AbstractMatrix, k::Integer, jr::BandRangeType) = lazy_getindex(A, k, jr)


text:
const _gr_mimeformats = Dict(
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
)

const _gr_wstype_default = @static if is_linux()
    "x11"
    # "cairox11"
elseif is_apple()
    "quartz"
else
    "use_default"
end

const _gr_wstype = Ref(get(ENV, "GKS_WSTYPE", _gr_wstype_default))
gr_set_output(wstype::String) = (_gr_wstype[] = wstype)

for (mime, fmt) in _gr_mimeformats
    @eval function Plots._show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{GRBackend})
        GR.emergencyclosegks()
        filepath = tempname() * "." * $fmt
        ENV["GKS_WSTYPE"] = $fmt
        ENV["GKS_FILEPATH"] = filepath
        gr_display(plt)
        GR.emergencyclosegks()
        write(io, readstring(filepath))
        rm(filepath)
    end
end

function Plots._display(plt::Plot{GRBackend})
    if plt[:display_type] == :inline
        GR.emergencyclosegks()
        filepath = tempname() * ".pdf"
        ENV["GKS_WSTYPE"] = "pdf"
        ENV["GKS_FILEPATH"] = filepath
        gr_display(plt)
        GR.emergencyclosegks()
        content = string("\033]1337;File=inline=1;preserveAspectRatio=0:", base64encode(open(readbytes, filepath)), "\a")
        println(content)
        rm(filepath)
    else
        ENV["GKS_DOUBLE_BUF"] = true
        if _gr_wstype[] != "use_default"
            ENV["GKS_WSTYPE"] = _gr_wstype[]
        end
        gr_display(plt)
    end
end

Plots.closeall(::GRBackend) = GR.emergencyclosegks()


text:
#=
Each complex should have collection of cells per dimension:
- cells::Dict{Int,Vector{C}} or Vector{Vector{C}}
=#
abstract type AbstractComplex end

function boundary(cplx::AbstractComplex, ch::Chain{PID, IX}) where {PID, IX<:Integer}
    d = dim(ch)
    cc = Chain(d-1, PID, IX)
    for (coef,elem) in ch
        append!(cc, coef * boundary(cplx, elem, d, PID))
    end
    return simplify(cc)
end

function coboundary(cplx::AbstractComplex, ch::Chain{PID, IX}) where {PID, IX<:Integer}
    d = dim(ch)
    cc = Chain(d+1, PID, IX)
    for (coef,elem) in ch
        append!(cc, coef * coboundary(cplx, elem, d, PID))
    end
    return simplify(cc)
end

#
# AbstractComplex Public Interface
#
"""Return a complex boundary given element and dimension"""
boundary(cplx::AbstractComplex, i::Integer, d::Int, ::Type{PID}) where {PID} = throw(MethodError(boundary, (typeof(cplx),Int,Int,PID)))
boundary(cplx::AbstractComplex, i::Integer, d::Int) = boundary(cplx, i, d, Int)

"""Return a complex coboundary given element and dimension"""
coboundary(cplx::AbstractComplex, i::Integer, d::Int, ::Type{PID}) where {PID} = throw(MethodError(coboundary, (typeof(cplx),Int,Int,PID)))
coboundary(cplx::AbstractComplex, i::Integer, d::Int) = coboundary(cplx, i, d, Int)

"""Return a complex cell type"""
celltype(cplx::AbstractComplex) = throw(MethodError(celltype, (typeof(cplx),)))

"""Return a cell collection per dimension (increasing)"""
cells(cplx::AbstractComplex) = throw(MethodError(cells, (typeof(cplx),)))

"""Return a nullable cell collection per dimension (increasing)"""
cells(cplx::AbstractComplex, d::Int) = throw(MethodError(cells, (typeof(cplx),Int)))

"""
    push!(cplx::AbstractComplex, cell::AbstractCell; recursive=false) -> Vector{AbstractCell}

Insert a `cell` to a complex `cplx`, and returns an array of inserted cell(s). If `recursive=true` is passed then all faces of the `cell` are also added to the complex `cplx`.
"""
Base.push!(cplx::AbstractComplex, c::AbstractCell; recursive=false) = throw(MethodError(push!, (typeof(cplx),typeof(c))))

#
# Public Methods
#
"""Return a number of cells per dimension"""
Base.size(cplx::AbstractComplex) = (map(length, cells(cplx))...,)

"""Return a dimension of the complex"""
dim(cplx::AbstractComplex) = length(size(cplx))-1

"""Return a total number of cells in the complex"""
Base.length(cplx::AbstractComplex) = sum(size(cplx))

"""Return a number of the cell in the complex of a dimension `d` (0-based)"""
function Base.size(cplx::AbstractComplex, d::Int)
    sz = size(cplx)
    szlen = length(sz)
    (d < 0 || d >= szlen) && return 0
    return sz[d+1]
end

"""
    position(complex, index, dimenion)

Return a position of the cell in an order of cells of the same dimenion of the `complex` given its `index` and `dimenion`.
"""
function Base.position(cplx::AbstractComplex, idx::Integer, d::Int)
    dcells = cells(cplx, d)
    dcells === nothing && return nothing
    cidx = findfirst(c->hash(c) == idx, dcells)
    return cidx
end

"""
    position(complex, cell)

Return a position of the `cell` in an order of cells of the same dimenion of the `complex`.
"""
function Base.position(cplx::AbstractComplex, c::C) where {C<:AbstractCell}
    @assert celltype(cplx) == C "Incorrect cell type"
    return position(cplx, hash(c), dim(c))
end

"""
    cplx[cell]

Return an identifier of the k-dimensional `cell` in the complex `cplx`. If the cell is not in the complex, `nothing` is returned.
"""
function Base.getindex(cplx::AbstractComplex, c::C) where {C <: AbstractCell}
    cidx = position(cplx, c)
    cidx === nothing && return nothing #size(cplx, d)+1
    return hash(cells(cplx, dim(c))[cidx])
end

"""
    cplx[idx, d]

Return a `d`-dimensional cell given its index `idx` and dimenion `d`.
"""
function Base.getindex(cplx::AbstractComplex, idx::Integer, d::Int)
    cidx = position(cplx, idx, d)
    cidx === nothing && return nothing
    return cells(cplx, d)[cidx]
end

"""
    boundary(complex, d, PID)

Generate a boundary matrix from the cell `complex` of the dimension `d` in using coefficients of `PID`.
"""
function boundary(cplx::AbstractComplex, d::Int, ::Type{PID}) where {PID}
    csize = size(cplx)
    rows = d > 0 ? csize[d] : 0
    cols = d <= dim(cplx) ? csize[d+1] : 0
    bm = spzeros(PID, rows, cols)
    if d>=0 && d <= dim(cplx)
        for c in cells(cplx, d)
            idx = hash(c)
            i = position(cplx, idx, d)
            for (coef,elem) in boundary(cplx, idx, d, PID)
                j = position(cplx, elem, d-1)
                bm[j, i] = coef
            end
        end
    end
    return bm
end

"""
    cell in cplx

Checks if the `cell` is in the complex `cplx`
"""
function Base.in(c::C, cplx::AbstractComplex) where {C <: AbstractCell}
    return position(cplx, c) !== nothing
end

"""
    cochain(complex, d, coefficients)

Return a cochain of the dimension `d` for the `complex` with PID `coefficients`.
"""
function cochain(cplx::AbstractComplex, d::Int, coefs::Vector{PID}) where {PID}
    cs = cells(cplx, d)
    @assert length(cs) == length(coefs) "Number of coefficients must match number of cells of dimension $d"
    return Chain(d, coefs, map(c->hash(c), cs))
end


text:
@testset "Narp.jl" begin
    filename = "data\\test.narp"
    df = DDR2import.Narp.read(filename)
    @test df[1, :AIRPORT_ID] == "AGGA"
    @test df[2, :AIRPORT_NAME] == "BALALAE"
    @test df[3, :LAT_DEG] == -9.4316666667
    @test df[3, :LON_DEG] == 160.0533333333
    @test df[9, :TIS] == 0
    @test df[9, :TRS] == 0
    @test df[9, :TAXITIME] == 0
    @test df[9, :ALTITUDE_FL] == 82
    #LAST
    @test df[13, :AIRPORT_ID] == "AGGV"
end


text:
@testset "Katz" begin
    g5 = DiGraph(4)
    add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
    for g in testdigraphs(g5)
        z = @inferred(katz_centrality(g, 0.4))
        @test round.(z, 2) == [0.32, 0.44, 0.62, 0.56]
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
module TestSQLite
using Test
using SQLite
using StructDatabaseMapping
using Dates

include("../includes/basic_test.jl")

DB_FILE = "test_db"

function test()
    test_create_tables()
    test_sqlite()
end
function cleanup()
    try
        rm(DB_FILE)
    catch
    end
end
function test_create_tables()
    mapper = DBMapper(()->SQLite.DB(DB_FILE))

    register!(mapper, Author)
    register!(mapper, Book)
    configure_relation(mapper, Book, :author, on_delete=Cascade())
    @test (StructDatabaseMapping.create_table_query(mapper, Author)
    == "CREATE TABLE IF NOT EXISTS author (" *
       "age INTEGER NOT NULL, " *
       "country VARCHAR NOT NULL, " *
       "date DATETIME NOT NULL, " *
       "id INTEGER PRIMARY KEY, " *
       "name VARCHAR NOT NULL)")
    @test (StructDatabaseMapping.create_table_query(mapper, Book)
     == "CREATE TABLE IF NOT EXISTS book (" *
         "author_id INTEGER NOT NULL, " *
         "data JSON NOT NULL, " *
         "id VARCHAR PRIMARY KEY, " *
         "title VARCHAR NOT NULL, "*
         "FOREIGN KEY(author_id) REFERENCES author(id) ON DELETE CASCADE ON UPDATE NO ACTION)")
    configure_relation(mapper, Book, :author, on_delete=Restrict(), on_update=Cascade())
    @test (StructDatabaseMapping.create_table_query(mapper, Book)
     == "CREATE TABLE IF NOT EXISTS book (" *
         "author_id INTEGER NOT NULL, " *
         "data JSON NOT NULL, " *
         "id VARCHAR PRIMARY KEY, " *
         "title VARCHAR NOT NULL, " *
         "FOREIGN KEY(author_id) REFERENCES author(id) ON DELETE RESTRICT ON UPDATE CASCADE)")
end
function test_sqlite()

    _test_basic_functionalities(()->SQLite.DB(DB_FILE))

end
function test_cleanup()
    mapper = DBMapper(()->SQLite.DB(DB_FILE))

    register!(mapper, Author)
    register!(mapper, Book)

    clean_table!(mapper, Author)
    clean_table!(mapper, Book)

    drop_table!(mapper, Author)
    drop_table!(mapper, Book)
end
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
import Base.Test # Even if using was used here full qualification of the functions in Base.Test is needed, or if any other module was used.
import JSON
function main()
    s = "{\"number\" : 5, \"arr\" : [3, 2, 5, 4]}"
    j = JSON.parse(s)
    Base.Test.@test  maximum(j["arr"]) == j["number"]
    println(join(ARGS,","))
    println("It works")
end