text:
export TensorMesh3D, getTensorMesh3D
export getCellCenteredGrid, getNodalGrid, getFaceGrids, getEdgeGrids
export getCellCenteredAxes, getNodalAxes,getBoundaryNodes
export getVolume, getVolumeInv, getFaceArea, getFaceAreaInv, getLength, getLengthInv

"""
	mutable struct jInv.Mesh.TensorMesh3D <: AbstractTensorMesh

	Fields:

		h1::Vector{Float64}     - cell size in x1 direction
		h2::Vector{Float64}     - cell size in x2 direction
		h3::Vector{Float64}     - cell size in x3 direction
		x0::Vector{Float64}     - origin
		dim::Int                - dimension (dim=3)
		n::Vector{Int64}        - number of cells in each direction
		nc::Int                 - nc total number of cells (nc=prod(n))
		nf::Vector{Int64}       - number of faces
		ne::Vector{Int64}       - number of edges

	Persistent Operators:

	Operators should not be accessed directly. They will be built, if needed,
	when accessing them using specified method. clear!(M) will release all
	memory.

		Div::SparseMatrixCSC    - divergence (faces -> cell-centers)
								  Access via: getDivergenceMatrix(M)
		Grad::SparseMatrixCSC   - gradient (nodal -> edges)
								  Access via: getNodalGradientMatrix(M)
		Curl::SparseMatrixCSC   - curl (edges -> faces)
								  Access via: getCurlMatrix(M)
		Af::SparseMatrixCSC     - face average (faces -> cell-centers)
								  Access via: getFaceAverageMatrix(M)
		Ae::SparseMatrixCSC     - edge average (edges -> cell-centers)
								  Access via: getEdgeAverageMatrix(M)
		An::SparseMatrixCSC     - nodal average (nodes -> cell-centers)
								  Access via: getNodalAverageMatrix(M)
		V::SparseMatrixCSC      - cell volumes (diagonal matrix)
								  Access via: getVolume(M)
		F::SparseMatrixCSC      - face area (diagonal matrix)
								  Access via: getFaceArea(M)
		L::SparseMatrixCSC      - edge length (diagonal matrix)
								  Access via: getLength(M)
		Vi::SparseMatrixCSC     - inverse cell volumes (diagonal matrix)
								  Access via: getVolumeInv(M)
		Fi::SparseMatrixCSC     - inverse face area (diagonal matrix)
								  Access via: getFaceAreaInv(M)
		Li::SparseMatrixCSC     - inverse edge length (diagonal matrix)
								  Access via: getLengthAreaInv(M)
		nLap::SparseMatrixCSC   - nodal Laplacian
								  Access via: getNodalLaplacian(M)

	Example:

	h1 = rand(4); h2 = rand(6); h3 = rand(5);
	M  = getTensorMesh2D(h1,h2,h3)

"""
mutable struct TensorMesh3D <: AbstractTensorMesh
	h1::Vector{Float64}
	h2::Vector{Float64}
	h3::Vector{Float64}
	x0::Vector{Float64}
	dim::Int
	n::Vector{Int64}
	nc::Int
	nf::Vector{Int64}
	ne::Vector{Int64}
	Div::SparseMatrixCSC
	Grad::SparseMatrixCSC
	Curl::SparseMatrixCSC
	Af::SparseMatrixCSC
	Ae::SparseMatrixCSC
	An::SparseMatrixCSC
	V::SparseMatrixCSC
	F::SparseMatrixCSC
	L::SparseMatrixCSC
	Vi::SparseMatrixCSC
	Fi::SparseMatrixCSC
	Li::SparseMatrixCSC
	nLap::SparseMatrixCSC
end

"""
	function jInv.Mesh.getTensorMesh3D

	constructs TensorMesh3D

	Required Input:

	   h1::Array - cell-sizes in x1 direction
	   h2::Array - cell-sizes in x2 direction
	   h3::Array - cell-sizes in x3 direction

	Optional Input
	   x0::Array - origin (default = zeros(3))

"""
function getTensorMesh3D(h1::Array,h2::Array,h3::Array,x0=zeros(3))
	n = [length(h1); length(h2);  length(h3)]
	nc = prod(n)
	nf = [(n[1]+1)*n[2]*n[3]; n[1]*(n[2]+1)*n[3]; n[1]*n[2]*(n[3]+1) ]
	ne = [n[1]*(n[2]+1)*(n[3]+1); (n[1]+1)*n[2]*(n[3]+1); (n[1]+1)*(n[2]+1)*n[3] ]
	empt = spzeros(0,0);
	dim = 3
return  TensorMesh3D(h1,h2,h3,x0,dim,n,nc,nf,ne,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt)
end

import Base.==
function ==(M1::TensorMesh3D,M2::TensorMesh3D)
	isEqual = fill(true,20)

	# check mandatory fields
	isEqual[1] =  (M1.h1 == M2.h1) & (M1.h2 == M2.h2) & (M1.h3 == M2.h3)
	isEqual[3] =  (M1.x0    == M2.x0)
	isEqual[4] =  (M1.dim   == M2.dim)
	isEqual[5] =  (M1.n     == M2.n)
	isEqual[6] =  (M1.nc    == M2.nc)
	isEqual[7] =  (M1.nf    == M2.nf)
	isEqual[8] =  (M1.ne    == M2.ne)

	# check fields that might be empty
	if !(isempty(M1.Div)) && !(isempty(M2.Div))
		isEqual[9] = (M1.Div == M2.Div)
	end
	if !(isempty(M1.Grad)) && !(isempty(M2.Grad))
		isEqual[10] = (M1.Grad == M2.Grad)
	end
	if !(isempty(M1.Curl)) && !(isempty(M2.Curl))
		isEqual[11] = (M1.Curl == M2.Curl)
	end
	if !(isempty(M1.Af)) && !(isempty(M2.Af))
		isEqual[12] = (M1.Af == M2.Af)
	end
	if !(isempty(M1.Ae)) && !(isempty(M2.Ae))
		isEqual[13] = (M1.Ae == M2.Ae)
	end
	if !(isempty(M1.An)) && !(isempty(M2.An))
		isEqual[14] = (M1.An == M2.An)
	end
	if !(isempty(M1.V)) && !(isempty(M2.V))
		isEqual[15] = (M1.V == M2.V)
	end
	if !(isempty(M1.F)) && !(isempty(M2.F))
		isEqual[16] = (M1.F == M2.F)
	end
	if !(isempty(M1.L)) && !(isempty(M2.L))
		isEqual[17] = (M1.L == M2.L)
	end
	if !(isempty(M1.Vi)) && !(isempty(M2.Vi))
		isEqual[18] = (M1.Vi == M2.Vi)
	end
	if !(isempty(M1.Fi)) && !(isempty(M2.Fi))
		isEqual[19] = (M1.Fi == M2.Fi)
	end
	if !(isempty(M1.Li)) && !(isempty(M2.Li))
		isEqual[20] = (M1.Li == M2.Li)
	end
	return all(isEqual)
end

function getNodalAxes(Mesh::TensorMesh3D)
	xn,yn,zn = getNodalAxes(Mesh.h1,Mesh.h2,Mesh.h3)
	return xn.+Mesh.x0[1],yn.+Mesh.x0[2],zn.+Mesh.x0[3]
end

function getCellCenteredAxes(Mesh::TensorMesh3D)
	xc,yc,zc = getCellCenteredAxes(Mesh.h1,Mesh.h2,Mesh.h3)
	return xc.+Mesh.x0[1],yc.+Mesh.x0[2],zc.+Mesh.x0[3]
end

function getNodalGrid(Mesh::TensorMesh3D)
# X = getNodalGrid(Mesh::TensorMesh3D)
	xn,yn,zn = getNodalAxes(Mesh)
	Xn,Yn,Zn = ndgrid(xn,yn,zn)
	return [vec(Xn) vec(Yn) vec(Zn)]
end

function getCellCenteredGrid(Mesh::TensorMesh3D)
# X = getCellCenteredGrid(Mesh::TensorMesh3D)
	xc,yc,zc = getCellCenteredAxes(Mesh)
	Xc,Yc,Zc = ndgrid(xc,yc,zc);
	return [vec(Xc) vec(Yc) vec(Zc)]
end

function getEdgeGrids(Mesh::TensorMesh3D)
# [Xe1, Xe2, Xe3] = getEdgeGrids(Mesh::TensorMesh3D)

	xn,yn,zn    = getNodalAxes(Mesh)
	xc,yc,zc    = getCellCenteredAxes(Mesh)

	Xe1,Ye1,Ze1 = ndgrid(xc,yn,zn)
	Xe2,Ye2,Ze2 = ndgrid(xn,yc,zn)
	Xe3,Ye3,Ze3 = ndgrid(xn,yn,zc)

	return  [vec(Xe1) vec(Ye1) vec(Ze1)], [vec(Xe2) vec(Ye2) vec(Ze2)], [vec(Xe3) vec(Ye3) vec(Ze3)]

end

function getFaceGrids(Mesh::TensorMesh3D)
# [Xf1, Xf2, Xf3] = getFaceGrids(Mesh::TensorMesh3D)

	xn,yn,zn    = getNodalAxes(Mesh)
	xc,yc,zc    = getCellCenteredAxes(Mesh)

	Xf1,Yf1,Zf1 = ndgrid(xn,yc,zc)
	Xf2,Yf2,Zf2 = ndgrid(xc,yn,zc)
	Xf3,Yf3,Zf3 = ndgrid(xc,yc,zn)

	return [vec(Xf1) vec(Yf1) vec(Zf1)], [ vec(Xf2) vec(Yf2) vec(Zf2)], [vec(Xf3) vec(Yf3) vec(Zf3)]

end

function getNodalAxes(h1,h2,h3)
	nc = [length(h1); length(h2);  length(h3)]

	x = zeros(nc[1]+1); for i=1:nc[1], x[i+1] = x[i] + h1[i]; end
	y = zeros(nc[2]+1); for i=1:nc[2], y[i+1] = y[i] + h2[i]; end
	z = zeros(nc[3]+1); for i=1:nc[3], z[i+1] = z[i] + h3[i]; end
	return x,y,z
end

function getCellCenteredAxes(h1,h2,h3)
	x,y,z = getNodalAxes(h1,h2,h3)
	# cell centered grids
	xc = x[1:end-1] + h1/2;
	yc = y[1:end-1] + h2/2;
	zc = z[1:end-1] + h3/2;
	return xc,yc,zc
end

function getNodalGrid(h1,h2,h3)
# X = getNodalGrid(h1,h2,h3)
	xn,yn,zn = getNodalAxes(h1,h2,h3)
	Xn,Yn,Zn = ndgrid(xn,yn,zn)
	return [vec(Xn) vec(Yn) vec(Zn)]

end

function getCellCenteredGrid(h1,h2,h3)
# X = getCellCenteredGrid(h1,h2,h3)
	xc,yc,zc = getCellCenteredAxes(h1,h2,h3)
	Xc,Yc,Zc = ndgrid(xc,yc,zc);
	return [vec(Xc) vec(Yc) vec(Zc)]
end

function getEdgeGrids(h1,h2,h3)
# [Xe1,Xe2,Xe3] = getEdgeGrids(h1,h2,h3)
	xn,yn,zn = getNodalAxes(h1,h2,h3)
	xc,yc,zc = getCellCenteredAxes(h1,h2,h3)

	Xe1,Ye1,Ze1 = ndgrid(xc,yn,zn)
	Xe2,Ye2,Ze2 = ndgrid(xn,yc,zn)
	Xe3,Ye3,Ze3 = ndgrid(xn,yn,zc)
	return  [vec(Xe1) vec(Ye1) vec(Ze1)], [vec(Xe2) vec(Ye2) vec(Ze2)], [vec(Xe3) vec(Ye3) vec(Ze3)]

end

function getFaceGrids(h1,h2,h3)
# [Xf1, Xf2, Xf3] = getFaceGrids(h1,h2,h3)
	xn,yn,zn = getNodalAxes(h1,h2,h3)
	xc,yc,zc = getCellCenteredAxes(h1,h2,h3)

	Xf1,Yf1,Zf1 = ndgrid(xn,yc,zc)
	Xf2,Yf2,Zf2 = ndgrid(xc,yn,zc)
	Xf3,Yf3,Zf3 = ndgrid(xc,yc,zn)

	return [vec(Xf1) vec(Yf1) vec(Zf1)], [ vec(Xf2) vec(Yf2) vec(Zf2)], [vec(Xf3) vec(Yf3) vec(Zf3)]
end

# --- linear operators for tensor mesh
function getVolume(Mesh::TensorMesh3D;saveMat::Bool=true)
# Mesh.V = getVolume(Mesh::TensorMesh3D) computes volumes v, returns diag(v)
	if isempty(Mesh.Vi)
		V = kron(sdiag(Mesh.h3),kron(sdiag(Mesh.h2),sdiag(Mesh.h1)))
		if saveMat
			Mesh.V = V;
		end
		return V;
	else
		return Mesh.V
	end
end
function getVolumeInv(Mesh::TensorMesh3D;saveMat::Bool=true)
# Mesh.Vi = getVolumeInv(Mesh::TensorMesh3D) returns sdiag(1 ./v)
	if isempty(Mesh.Vi)
		Vi = kron(sdiag(1 ./Mesh.h3),kron(sdiag(1 ./Mesh.h2),sdiag(1 ./Mesh.h1)))
		if saveMat
			Mesh.Vi = Vi;
		end
		return Vi;
	else
		return Mesh.Vi
	end
end
function getFaceArea(Mesh::TensorMesh3D)
# Mesh.F = getFaceArea(Mesh::TensorMesh3D) computes face areas a, returns  sdiag(a)
	if isempty(Mesh.F)
		f1  = kron(sdiag(Mesh.h3)      ,kron(sdiag(Mesh.h2)      ,sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		f2  = kron(sdiag(Mesh.h3)      ,kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1) ,sdiag(Mesh.h1)))
		f3  = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1) ,kron(sdiag(Mesh.h2)      ,sdiag(Mesh.h1)))
		Mesh.F = blockdiag(blockdiag(f1,f2),f3)
	end
return Mesh.F
end
function getFaceAreaInv(Mesh::TensorMesh3D)
# Mesh.Fi = getFaceAreaInv(Mesh::TensorMesh3D) computes inverse of face areas, returns sdiag(1 ./a)
	if isempty(Mesh.Fi)
		f1i  = kron(sdiag(1 ./Mesh.h3)   ,kron(sdiag(1 ./Mesh.h2)   ,sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		f2i  = kron(sdiag(1 ./Mesh.h3)   ,kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1) ,sdiag(1 ./Mesh.h1)))
		f3i  = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1) ,kron(sdiag(1 ./Mesh.h2)   ,sdiag(1 ./Mesh.h1)))
		Mesh.Fi = blockdiag(blockdiag(f1i,f2i),f3i)
	end
return Mesh.Fi
end

function getLength(Mesh::TensorMesh3D)
# Mesh.L = getLength(Mesh::TensorMesh3D) computes edge lengths l, returns sdiag(l)
	if isempty(Mesh.L)
		l1  = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1),kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1),sdiag(Mesh.h1)))
		l2  = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1),kron(sdiag(Mesh.h2),sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		l3  = kron(sdiag(Mesh.h3),kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1),sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		Mesh.L   = blockdiag(blockdiag(l1,l2),l3);
	end
return Mesh.L
end

function getLengthInv(Mesh::TensorMesh3D)
# Mesh.L = getLength(Mesh::TensorMesh3D) computes inverse of edge lengths l, returns sdiag(1 ./l)
	if isempty(Mesh.Li)
		l1i = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1),kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1),sdiag(1 ./Mesh.h1)))
		l2i = kron(sparse(1.0I,Mesh.n[3]+1,Mesh.n[3]+1),kron(sdiag(1 ./Mesh.h2),sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		l3i = kron(sdiag(1 ./Mesh.h3),kron(sparse(1.0I,Mesh.n[2]+1,Mesh.n[2]+1),sparse(1.0I,Mesh.n[1]+1,Mesh.n[1]+1)))
		Mesh.Li  = blockdiag(blockdiag(l1i,l2i),l3i);
	end
return Mesh.Li
end

"""
function jInv.Mesh.getBoundaryNodes(Mesh)

	Returns an array tuple (boundary node indices, inner node indices)

	Input:
		Mesh::Abstract Mesh

	Output:
		Tuple{Array{Int64,1},Array{Int64,1}}
"""
function getBoundaryNodes(M::AbstractTensorMesh)

    # number of cells in each direction
    N = M.n

    # nodal matrix
    nodal_mat = 1:prod(N.+1)

    # check dimension
    if M.dim==2

        nodal_mat = reshape(nodal_mat, N[1]+1, N[2]+1)

        iin = nodal_mat[2:end-1,2:end-1]
        ib = setdiff(nodal_mat, iin)

    elseif M.dim==3

        nodal_mat = reshape(nodal_mat, N[1]+1, N[2]+1, N[3]+1)

        iin = nodal_mat[2:end-1,2:end-1,2:end-1]
        ib = setdiff(nodal_mat, iin)

    end

    return (vec(ib), vec(iin))
end

include("getEdgeIntegralOfPolygonalChain.jl")


text:
export FileSource, GitSource, DirectorySource

"""
An `AbstractSource` is something used as source to build the package.  Sources
are installed to `\${WORKSPACE}/srcdir` in the build environment.

Concrete subtypes of `AbstractSource` are:

* [`FileSource`](@ref): a remote file to download from the Internet;
* [`GitSource`](@ref): a remote Git repository to clone;
* [`DirectorySource`](@ref): a local directory to mount.
"""
abstract type AbstractSource end

"""
    FileSource(url::String, hash::String; unpack_target::String = "")

Specify a remote file to be downloaded from the Internet from `url`.  `hash` is
the 64-character SHA256 checksum of the file.

If the file is an archive in one of the supported archive formats (e.g., TAR or
ZIP balls), it will be automatically unpacked to `\${WORKSPACE}/srcdir`, or in
its subdirectory pointed to by the optional keyword `unpack_target`, if
provided.
"""
struct FileSource <: AbstractSource
    url::String
    hash::String
    unpack_target::String
end
FileSource(url::String, hash::String; unpack_target::String = "") =
    FileSource(url, hash, unpack_target)

"""
    GitSource(url::String, hash::String; unpack_target::String = "")

Specify a remote Git repository to clone form `url`.  `hash` is the 40-character
SHA1 revision to checkout after cloning.

The repository will be cloned in `\${WORKSPACE}/srcdir`, or in its subdirectory
pointed to by the optional keyword `unpack_target`, if provided.
"""
struct GitSource <: AbstractSource
    url::String
    hash::String
    unpack_target::String
end
GitSource(url::String, hash::String; unpack_target::String = "") =
    GitSource(url, hash, unpack_target)

"""
    DirectorySource(path::String; unpack_target::String = "")

Specify a local directory to mount from `path`.

The content of the directory will be mounted in `\${WORKSPACE}/srcdir`, or in its
subdirectory pointed to by the optional keyword `unpack_target`, if provided.
"""
struct DirectorySource <: AbstractSource
    path::String
    unpack_target::String
end
DirectorySource(path::String; unpack_target::String = "") =
    DirectorySource(path, unpack_target)

# This is not meant to be used as source in the `build_tarballs.jl` scripts but
# only to set up the source in the workspace.
struct SetupSource
    path::String
    hash::String
    unpack_target::String
end
SetupSource(path::String, hash::String; unpack_target::String = "") =
    SetupSource(path, hash, unpack_target)

function download_source(source::FileSource; verbose::Bool = false)
    if isfile(source.url)
        # Immediately abspath() a src_url so we don't lose track of
        # sources given to us with a relative path
        src_path = abspath(source.url)

        # And if this is a locally-sourced tarball, just verify
        verify(src_path, source.hash; verbose=verbose)
    else
        # Otherwise, download and verify
        src_path = storage_dir("downloads", string(source.hash, "-", basename(source.url)))
        download_verify(source.url, source.hash, src_path; verbose=verbose)
    end
    return SetupSource(src_path, source.hash, source.unpack_target)
end

function download_source(source::GitSource; verbose::Bool = false)
    src_path = storage_dir("downloads", basename(source.url))

    # If this git repository already exists, ensure that its origin remote actually matches
    if isdir(src_path)
        origin_url = LibGit2.with(LibGit2.GitRepo(src_path)) do repo
            LibGit2.url(LibGit2.get(LibGit2.GitRemote, repo, "origin"))
        end

        # If the origin url doesn't match, wipe out this git repo.  We'd rather have a
        # thrashed cache than an incorrect cache.
        if origin_url != source.url
            rm(src_path; recursive=true, force=true)
        end
    end

    if isdir(src_path)
        # If we didn't just mercilessly obliterate the cached git repo, use it!
        LibGit2.with(LibGit2.GitRepo(src_path)) do repo
            LibGit2.fetch(repo)
        end
    else
        # If there is no src_path yet, clone it down.
        repo = LibGit2.clone(source.url, src_path; isbare=true)
    end
    return SetupSource(src_path, source.hash, source.unpack_target)
end

function download_source(source::DirectorySource; verbose::Bool = false)
    if !isdir(source.path)
        error("Could not find directory \"$(source.path)\".")
    end

    # Package up this directory and calculate its hash
    tarball_pathv = ""
    tarball_hash  = ""
    mktempdir() do tempdir
        tarball_path = joinpath(tempdir, basename(source.path) * ".tar.gz")
        package(source.path, tarball_path)
        tarball_hash = open(tarball_path, "r") do f
            bytes2hex(sha256(f))
        end

        # Move it to a filename that has the hash as a part of it (to avoid name collisions)
        tarball_pathv = storage_dir("downloads", string(tarball_hash, "-", basename(source.path), ".tar.gz"))
        mv(tarball_path, tarball_pathv; force=true)
    end
    return SetupSource(tarball_pathv, tarball_hash, source.unpack_target)
end

"""
    download_source(source::AbstractSource; verbose::Bool = false)

Download the given `source`.  All downloads are cached within the
BinaryBuilder `downloads` storage directory.
"""
download_source

# Add JSON serialization to sources
JSON.lower(fs::FileSource) = Dict("type" => "file", extract_fields(fs)...)
JSON.lower(gs::GitSource) = Dict("type" => "git", extract_fields(gs)...)
JSON.lower(ds::DirectorySource) = Dict("type" => "directory", extract_fields(ds)...)

# When deserialiasing the JSON file, the sources are in the form of
# dictionaries.  This function converts the dictionary back to the appropriate
# AbstractSource.
function sourcify(d::Dict)
    if d["type"] == "directory"
        return DirectorySource(d["path"])
    elseif d["type"] == "git"
        return GitSource(d["url"], d["hash"])
    elseif d["type"] == "file"
        return FileSource(d["url"], d["hash"])
    else
        error("Cannot convert to source")
    end
end

# XXX: compatibility functions.  These are needed until we support old-style
# Pair/String sources specifications.
coerce_source(source::AbstractSource) = source
function coerce_source(source::AbstractString)
    @warn "Using a string as source is deprecated, use DirectorySource instead"
    return DirectorySource(source)
end
function coerce_source(source::Pair)
    src_url, src_hash = source
    if endswith(src_url, ".git")
        @warn "Using a pair as source is deprecated, use GitSource instead"
        return GitSource(src_url, src_hash)
    else
        @warn "Using a pair as source is deprecated, use FileSource instead"
        return FileSource(src_url, src_hash)
    end
end

ext: jl


text:
using Base.Test
using Clang.wrap_c

headername = "cxx/c_struct_typedef.h"
wc = wrap_c.init(;  headers = [headername], header_wrapped=(a,b)->true)
parsed = wrap_c.parse_c_headers(wc)
obuf = wc.output_bufs[headername]
wrap_c.wrap_header(wc, parsed[headername], headername, obuf)
common_buf = wrap_c.dump_to_buf(wc.common_buf)
io = IOBuffer()
wrap_c.print_buffer(io, common_buf)
str = takebuf_string(io)
@test str == "\ntype Str\n    x::Cint\nend\n"


text:
function test_tf_idf()
    println("-- Testing TF-IDF content-based recommender")

    m = [1  0  1  0  1  1  0  0  0  1
         0  1  1  1  0  0  0  1  0  0
         0  0  0  1  1  1  0  0  0  0
         0  0  1  1  0  0  1  1  0  0
         0  1  0  0  0  0  0  0  1  1
         1  0  0  1  0  0  0  0  0  0
         0  0  0  0  0  0  0  1  0  1
         0  0  1  1  0  0  1  0  0  1
         0  0  0  0  0  1  0  0  1  0
         0  1  0  0  1  0  1  0  0  0
         0  0  1  0  1  0  0  0  1  0
         1  0  0  0  0  1  1  0  0  0
         0  0  1  1  1  0  0  1  0  0
         0  1  1  1  0  0  0  0  1  0
         0  0  0  1  0  1  1  1  0  0
         1  0  0  0  0  1  0  0  1  0
         0  1  1  1  0  0  0  1  0  0
         0  0  0  1  0  0  0  0  1  0
         0  1  1  0  1  0  1  0  0  1
         0  0  1  1  0  0  1  0  1  0]

    tf = sum(m, dims=1)
    idf = 1 ./ tf

    users = [ 1  -1
             -1   1
              0   0
              0   1
              0   0
              1   0
              0   0
              0   0
              0   0
              0   0
              0   0
              0  -1
              0   0
              0   0
              0   0
              1   0
              0   1
              0   0
             -1   0
              0   0]
    uv1 = sparse(users[:, 1])
    uv2 = sparse(users[:, 2])

    n_doc = size(m, 1)
    docs = [i for i in 1:20]

    # Case: basic
    data = DataAccessor(m)
    set_user_attribute(data, 1, uv1)
    set_user_attribute(data, 2, uv2)

    recommender = TFIDF(data, tf, ones(size(tf))) # do not use IDF
    rec = recommend(recommender, 1, n_doc, docs)
    @test first(rec[2]) == 1
    @test first(rec[3]) == 12
    @test last(rec[2]) == 4 && last(rec[3]) == 4

    # Case: normalized matrix
    m_normalized = m ./ sqrt.(sum(m.^2, dims=2))
    data = DataAccessor(m_normalized)
    set_user_attribute(data, 1, uv1)
    set_user_attribute(data, 2, uv2)

    recommender = TFIDF(data, tf, ones(size(tf))) # do not use IDF
    rec = recommend(recommender, 1, n_doc, docs)
    @test first(rec[5]) == 1
    @test isapprox(last(rec[5]), 1.0090, atol=1e-4)
    rec = recommend(recommender, 2, n_doc, docs)
    @test first(rec[10]) == 7
    @test isapprox(last(rec[10]), 0.7444, atol=1e-4)
    @test first(rec[11]) == 19
    @test isapprox(last(rec[11]), 0.4834, atol=1e-4)

    # Case: normalized matrix & using IDF
    recommender = TFIDF(data, tf, idf)
    rec = recommend(recommender, 1, n_doc, docs)
    @test first(rec[4]) == 1
    @test isapprox(last(rec[4]), 0.2476, atol=1e-4)
end

test_tf_idf()


text:
using GNSSTools
using Statistics
using PyPlot
pygui(true)

CN0 = 36
prn = 10
t_length = 1
fd_range = 5000
f_s = 5e6
sigtype = define_l1ca_code_type()
# Simulate L1 C/A signal and perform course acquisition
signal = definesignal(sigtype, f_s, t_length; prn=prn, f_d=4910,
                      code_start_idx=2500, include_phase_noise=true, CN0=CN0)
generatesignal!(signal)
replica = definesignal(sigtype, f_s, 1e-3)
Δfd = 1/replica.t_length
fd_est, n0_est, SNR_est, corr_result = courseacquisition(signal, replica, prn;
                                                         return_corrresult=true,
                                                         Δfd=Δfd, fd_range=fd_range,
                                                         M=1)

max_idx = argmax(corr_result)
ps = corr_result[max_idx[1],:]
PS = corr_result[max_idx]
N, M = size(corr_result)
noise_val = sum(corr_result[max_idx[1],:])
PN = (noise_val - (replica.t_length/1e-3)*PS)/(M - 1)
SNR = 10log10(PS/PN)

################################################################################
h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # TCXO
file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/known_oscillator_data/20210327_013800_000005_Aero-North-TCXO_1575.42_25.sc16"
data = loaddata(:sc16, file, 25e6, L1_freq, L1_freq, 1., skip_to=1+0.00031216)

replica = definesignal(define_l1ca_code_type(), 25e6, 2e-3)
for i in 1:32
    fd_est, n0_est, SNR_est = courseacquisition(data, replica, i; M=10)
    println("PRN $i\t$SNR_est")
end

fig = figure()
acquisiton_T = 1e-3
fine_acq_T = 10e-3
tracking_T = 1e-3
for i in 1:32
    ax = fig.add_subplot(8, 4, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $i")
    acqresults, trackresults = process(data, define_l1ca_code_type(), i;
                                       show_plot=false,
                                       fine_acq_method=:fft, state_num=3,
                                       dynamickf=true, h₀=h_parms[3],
                                       h₋₂=h_parms[1], dll_b=1,
                                       acquisition_T=acquisiton_T,
                                       fine_acq_T=fine_acq_T,
                                       tracking_T=tracking_T,
                                       M=10)
    ax.plot(trackresults.t, real.(trackresults.ZP))
    ax.plot(trackresults.t, imag.(trackresults.ZP))
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.93, left=0.08, right=0.93,
                bottom=0.01)
plot_title = string("Course Acq. T=$(floor(Int, acquisiton_T*1000))ms;    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms;    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms")
suptitle(plot_title)


text:
import Base.*

export S¹
export U1

"""
    Represents a point in a one-sphere.
"""
abstract type S¹ end

"""
    Represents a point in U(1).

field: α.
"""
struct U1 <: S¹
    z::Complex
    U1(z::Complex) = begin
        @assert(isapprox(abs(z), 1), "The magnitude must be equal to 1.")
        new(z)
    end
    U1(α::Real) = begin
        @assert(-pi ≤ α ≤ pi, "The phase angle must be in the interval [-π, π].")
        U1(exp(im * α))
    end
end

Base.angle(u::U1) = angle(u.z)
*(u1::U1, u2::U1) = U1(angle(u1.z * u2.z))
Base.isapprox(u1::U1, u2::U1) = isapprox(u1.z, u2.z)


text:
import JSON
import DataFrames
import CategoricalArrays
using Printf

abstract type AbstractDataFrame end

abstract type StatVariate end

abstract type AbstractFactor{T<:Unsigned} <: StatVariate end

abstract type AbstractCovariate{T<:AbstractFloat} <: StatVariate end

abstract type AbstractBoolVariate <: StatVariate end

abstract type AbstractDateTimeVariate <: StatVariate end

abstract type AbstractIntVariate{T<:Signed} <: StatVariate end

struct DataFrame <: AbstractDataFrame
    length::Int64
    factors::AbstractVector{<:AbstractFactor}
    covariates::AbstractVector{<:AbstractCovariate}
    boolvariates::AbstractVector{<:AbstractBoolVariate}
    datetimevariates::AbstractVector{<:AbstractDateTimeVariate}
    intvariates::AbstractVector{<:AbstractIntVariate}
    headerpath::String
end

mutable struct CovariateStats
    obscount::Int64
    nancount::Int64
    nanpcnt::Float64
    sum::Float64
    sum2::Float64
    mean::Float64
    std::Float64
    min::Float64
    max::Float64
end

mutable struct DateTimeVariateStats
    obscount::Int64
    nancount::Int64
    nanpcnt::Float64
    min::Int64
    max::Int64
end

mutable struct LevelStats
    level::String
    freq::Int64
    freqpcnt::Float64
end

mutable struct FactorStats
    obscount::Int64
    missingfreq::Int64
    missingpcnt::Float64
    levelstats::Vector{LevelStats}
end

function Base.length(dataframe::AbstractDataFrame)
    dataframe.length
end

getlevels(factor::AbstractFactor{T}) where {T<:Unsigned} = factor.levels

getname(factor::AbstractFactor{T}) where {T<:Unsigned} = factor.name

getname(covariate::AbstractCovariate{T}) where {T<:AbstractFloat} = covariate.name

getname(intvariate::AbstractIntVariate{T}) where {T<:Signed} = intvariate.name

getname(boolvar::AbstractBoolVariate) = boolvar.name

getname(datetimevar::AbstractDateTimeVariate) = datetimevar.name

Base.length(covariate::AbstractCovariate{T}) where {T<:AbstractFloat} = covariate.length

Base.length(intvariate::AbstractIntVariate{T}) where {T<:Signed} = intvariate.length

Base.length(boolvar::AbstractBoolVariate) = boolvar.length

Base.length(factor::AbstractFactor{T}) where {T<:Unsigned} = factor.length

function DataFrame(path::String; preload::Bool = true)
    path = abspath(path)
    headerpath = isfile(path) ? path : joinpath(path, "header.txt")
    headerjson = open(headerpath) do f
        read(f, String)
    end
    header = JSON.parse(headerjson)
    factors = Vector{AbstractFactor}()
    covariates = Vector{AbstractCovariate}()
    dtvariates = Vector{AbstractDateTimeVariate}()
    intvariates = Vector{AbstractIntVariate}()
    datacols = header["datacolumns"]
    len = 0
    for datacol in datacols
        datatype = datacol["datatype"]
        len = datacol["length"]
        name = datacol["name"]
        datpath = joinpath(dirname(headerpath), datacol["filename"])

        if datatype == "Float32"
            if preload
                push!(covariates, Covariate{Float32}(name, len, datpath))
            else
                push!(covariates, FileCovariate{Float32}(name, len, datpath))
            end
        end

        if datatype == "DateTime"
            if preload
                push!(dtvariates, DateTimeVariate(name, len, datpath))
            else
                push!(dtvariates, FileDateTimeVariate(name, len, datpath))
            end
        end

        if datatype == "Int64"
            if preload
                push!(intvariates, IntVariate{Int64}(name, len, datpath))
            else
                push!(intvariates, FileIntVariate{Int64}(name, len, datpath))
            end
        end

        if datatype == "UInt8"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt8}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt8}(name, len, levels, datpath))
            end
        end

        if datatype == "UInt16"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt16}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt16}(name, len, levels, datpath))
            end
        end

        if datatype == "UInt32"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt32}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt32}(name, len, levels, datpath))
            end
        end
    end
    DataFrame(len, factors, covariates, AbstractBoolVariate[], dtvariates, intvariates, headerpath)
end

function DataFrame(df::DataFrames.DataFrame)
    factors = Vector{AbstractFactor}()
    covariates = Vector{AbstractCovariate}()
    dtvariates = Vector{AbstractDateTimeVariate}()
    intvariates = Vector{AbstractIntVariate}()
    colnames = DataFrames.names(df)
    for i = 1:length(colnames)
        colname = string(colnames[i])
        col = df[i]
        if typeof(col) <: Vector{<:AbstractFloat}
            T = eltype(col)
            push!(covariates, Covariate{T}(colname, col))
        elseif typeof(col) <: CategoricalArrays.CategoricalArray
            levels = [string(level) for level in CategoricalArrays.index(col.pool)]
            T = eltype(col.refs)
            push!(factors, Factor{T}(colname, levels, col.refs))
        elseif typeof(col) <: AbstractVector{<:Integer}
            T = eltype(col)
            col = convert(Vector{Float32}, collect(col))
            push!(covariates, Covariate{Float32}(colname, col))
        end
    end
    DataFrame(size(df)[1], factors, covariates, AbstractBoolVariate[], dtvariates, intvariates, "")
end

function Base.getindex(df::AbstractDataFrame, name::String)
    for factor in df.factors
        if getname(factor) == name
            return factor
        end
    end
    for cov in df.covariates
        if getname(cov) == name
            return cov
        end
    end
    for dt in df.datetimevariates
        if getname(dt) == name
            return dt
        end
    end

    for ivar in df.intvariates
        if getname(ivar) == name
            return ivar
        end
    end
end

function Base.getindex(df::AbstractDataFrame, name::Symbol)
    df[string(name)]
end

function Base.summary(factor::AbstractFactor{T}; selector::AbstractBoolVariate = BoolVariate("", BitArray{1}())) where {T<:Unsigned}
    io = IOBuffer()
    factorstats = getstats(factor; selector = selector)
    println(io, @sprintf("%-16s%15d", "Obs Count", factorstats.obscount))
    println(io, @sprintf("%-16s%15s%15s", "Level", "Frequency", "Frequency(%)"))
    println(io, @sprintf("%-16s%15d%15G", MISSINGLEVEL, factorstats.missingfreq, factorstats.missingpcnt))
    for levelstats in factorstats.levelstats
        println(io, @sprintf("%-16s%15d%15G", levelstats.level, levelstats.freq, levelstats.freqpcnt))
    end
    print(String(take!(io)))
end

function Base.summary(boolvar::AbstractBoolVariate)
    io = IOBuffer()
    len = length(boolvar)
    slices = slice(boolvar, 1, len, SLICELENGTH)
    truefreq = fold(0, slices) do acc, slice
        res = acc
        for v in slice
            if v
                res += 1
            end
        end
        res
    end
    println(io, @sprintf("%-15s%12d", "Obs Count", len))
    println(io, @sprintf("%-15s%12d", "True Freq", truefreq))
    println(io, @sprintf("%-15s%12G", "True Freq (%)", 100.0 * truefreq / len))
    print(String(take!(io)))
end

function getstats(factor::AbstractFactor{T}; selector::AbstractBoolVariate = BoolVariate("", BitArray{1}())) where {T<:Unsigned}
    len = length(factor)
    levels = getlevels(factor)
    levelcount = length(levels)
    slices = slice(factor, 1, len, SLICELENGTH)
    noselector = length(selector) == 0
    if noselector
        init = zeros(Int64, levelcount + 1)
        freq = fold(init, slices) do frq, slice
            for levelindex in slice
                frq[levelindex + 1] += 1
            end
            frq
        end
    else
        selslices = slice(selector, 1, len, SLICELENGTH)
        zipslices = zip(selslices, slices)
        stats0 = zeros(Int64, levelcount + 1)
        freq, selcount = fold((stats0, 0), zipslices) do acc, zipslice
            frq, k = acc
            selslice, slice = zipslice
            for i in 1:length(slice)
                if selslice[i]
                    k += 1
                    levelindex = slice[i]
                    frq[levelindex + 1] += 1
                end
            end
            frq, k
        end
    end
    missingfreq = freq[1]
    missingpcnt = 100.0 * missingfreq / (noselector ? len : selcount)
    levelstats = [LevelStats(levels[i], freq[i + 1], 100.0 * freq[i + 1] / (noselector ? len : selcount)) for i in 1:levelcount]
    FactorStats(len, missingfreq, missingpcnt, levelstats)
end

function getstats(covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    len = length(covariate)
    init = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
    slices = slice(covariate, 1, len, SLICELENGTH)
    stats = fold(init, slices) do s, slice
        for v in slice
            if isnan(v)
                s.nancount += 1
            else
                if isnan(s.sum)
                    s.sum = v
                    s.sum2 = v * v
                    s.min = v
                    s.max = v
                else
                    s.sum += v
                    s.sum2 += v * v
                    if v < s.min
                        s.min = v
                    end
                    if v > s.max
                        s.max = v
                    end
                end
            end
        end
        s
    end
    stats.obscount = len
    stats.nanpcnt = 100.0 * stats.nancount / len
    stats.mean = stats.sum / (len - stats.nancount)
    stats.std = sqrt(((stats.sum2 - stats.sum * stats.sum / (len - stats.nancount)) / (len - stats.nancount - 1)))
    stats
end

function getstats(factor::AbstractFactor{T}, covariate::AbstractCovariate{S}) where {T<:Unsigned} where {S<:AbstractFloat}
    len = length(covariate)
    levelcount = length(getlevels(factor))
    levelsinit = [CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64) for i in 1:levelcount]
    missinit = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
    covslices = slice(covariate, 1, len, SLICELENGTH)
    factorslices = slice(factor, 1, len, SLICELENGTH)
    zipslices = zip(factorslices, covslices)
    missacc, levelacc = fold((missinit, levelsinit), zipslices) do acc, slice
        missstat, levelsstat = acc
        fslice, cslice = slice
        for (i, levelindex) in enumerate(fslice)
            v = cslice[i]
            if levelindex == zero(T)
                missstat.obscount += 1
                if isnan(v)
                    missstat.nancount += 1
                else
                    if isnan(missstat.sum)
                        missstat.sum = v
                        missstat.sum2 = v * v
                        missstat.min = v
                        missstat.max = v
                    else
                        missstat.sum += v
                        missstat.sum2 += v * v
                        if v < missstat.min
                            missstat.min = v
                        end
                        if v > missstat.max
                            missstat.max = v
                        end
                    end
                end
            else
                levelsstat[levelindex].obscount += 1
                if isnan(v)
                    levelsstat[levelindex].nancount += 1
                else
                    if isnan(levelsstat[levelindex].sum)
                        levelsstat[levelindex].sum = v
                        levelsstat[levelindex].sum2 = v * v
                        levelsstat[levelindex].min = v
                        levelsstat[levelindex].max = v
                    else
                        levelsstat[levelindex].sum += v
                        levelsstat[levelindex].sum2 += v * v
                        if v < levelsstat[levelindex].min
                            levelsstat[levelindex].min = v
                        end
                        if v > levelsstat[levelindex].max
                            levelsstat[levelindex].max = v
                        end
                    end
                end
            end
        end
        acc
    end
    missacc.nanpcnt = 100.0 * missacc.nancount / missacc.obscount
    missacc.mean = missacc.sum / (missacc.obscount - missacc.nancount)
    missacc.std = sqrt(((missacc.sum2 - missacc.sum * missacc.sum / (missacc.obscount - missacc.nancount)) / (missacc.obscount - missacc.nancount - 1)))
    for i in 1:levelcount
        levelacc[i].nanpcnt = 100.0 * levelacc[i].nancount / levelacc[i].obscount
        levelacc[i].mean = levelacc[i].sum / (levelacc[i].obscount - levelacc[i].nancount)
        levelacc[i].std = sqrt(((levelacc[i].sum2 - levelacc[i].sum * levelacc[i].sum / (levelacc[i].obscount - levelacc[i].nancount)) / (levelacc[i].obscount - levelacc[i].nancount - 1)))
    end
    missacc, levelacc
end

function getstats(dtvariate::AbstractDateTimeVariate)
    len = length(dtvariate)
    init = DateTimeVariateStats(0, 0, NaN64, 0, 0)
    slices = slice(dtvariate, 1, len, SLICELENGTH)
    stats = fold(init, slices) do s, slice
        for v in slice
            if v == zero(Int64)
                s.nancount += 1
            else
                if s.min == zero(Int64)
                    s.min = v
                    s.max = v
                else
                    if v < s.min
                        s.min = v
                    end
                    if v > s.max
                        s.max = v
                    end
                end
            end
        end
        s
    end
    stats.obscount = len
    stats.nanpcnt = 100.0 * stats.nancount / len
    stats
end

function Base.summary(covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    io = IOBuffer()
    stats = getstats(covariate)
    println(io, @sprintf("%-15s%12d", "Obs Count", stats.obscount))
    println(io, @sprintf("%-15s%12d", "NaN Freq", stats.nancount))
    println(io, @sprintf("%-15s%12G", "NaN %", stats.nanpcnt))
    println(io, @sprintf("%-15s%12G", "Min", stats.min))
    println(io, @sprintf("%-15s%12G", "Max", stats.max))
    println(io, @sprintf("%-15s%12G", "Mean", stats.mean))
    println(io, @sprintf("%-15s%12G", "Std", stats.std))
    print(String(take!(io)))
end

function Base.summary(dtvariate::AbstractDateTimeVariate)
    io = IOBuffer()
    stats = getstats(dtvariate)
    println(io, @sprintf("%-15s%12d", "Obs Count", stats.obscount))
    println(io, @sprintf("%-15s%12d", "NaN Freq", stats.nancount))
    println(io, @sprintf("%-15s%12G", "NaN %", stats.nanpcnt))
    println(io, @sprintf("%-15s%12s", "Min", Dates.format(Dates.epochms2datetime(stats.min), "yyyy-mm-ddTHH:MM:SS")))
    println(io, @sprintf("%-15s%12s", "Max", Dates.format(Dates.epochms2datetime(stats.max), "yyyy-mm-ddTHH:MM:SS")))
    print(String(take!(io)))
end

function Base.show(io::IO, covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    slices = slice(covariate, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(covariate)
    if slice1 !== nothing
        datahead = join([isnan(v) ? "." : string(v) for v in slice1], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "Covariate $(getname(covariate)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "Covariate $(getname(covariate)) with $(len) obs")
    end
end

function Base.show(io::IO, intvariate::AbstractIntVariate{T}) where {T<:Signed}
    slices = slice(intvariate, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(intvariate)
    if slice1 !== nothing
        datahead = join([v == typemin(T) ? "." : string(v) for v in slice1], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "IntVariate $(getname(intvariate)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "IntVariate $(getname(intvariate)) with $(len) obs")
    end
end

function Base.show(io::IO, boolvar::AbstractBoolVariate)
    slices = slice(boolvar, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(boolvar)
    if slice1 !== nothing
        datahead = join([string(v) for v in slice1], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "BoolVar $(getname(boolvar)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "BoolVar $(getname(boolvar)) with $(len) obs")
    end
end

function Base.show(io::IO, datetimevar::AbstractDateTimeVariate)
    slices = slice(datetimevar, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(datetimevar)
    if slice1 !== nothing
        datahead = join([v == zero(Int64) ? MISSINGLEVEL : string(Dates.epochms2datetime(v)) for v in slice1], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "DateTimeVar $(getname(datetimevar)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "DateTimeVar $(getname(datetimevar)) with $(len) obs")
    end
end

function Base.show(io::IO, factor::AbstractFactor{T}) where {T<:Unsigned}
    slices = slice(factor, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(factor)
    levels = getlevels(factor)
    levelcount = length(levels)
    if slice1 !== nothing
        datahead = join([index == 0 ? MISSINGLEVEL : levels[index] for index in slice1], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "Factor $(getname(factor)) with $(len) obs and $(levelcount) levels: $(datahead)$dataend")
    else
        println(io, "Factor $(getname(factor)) with $(len) obs and $(levelcount) levels")
    end
end

function Base.convert(::Type{Vector{T}}, factor::AbstractFactor{T}) where {T<:Unsigned}
    slices = slice(factor, 1, length(factor), SLICELENGTH)
    data = Vector{T}(undef, length(factor))
    fold(0, slices) do offset, slice
        n = length(slice)
        view(data, (1 + offset):(n + offset)) .= slice
        offset += n
    end
    data
end

function Base.convert(::Type{Vector{T}}, covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    slices = slice(covariate, 1, length(covariate), SLICELENGTH)
    data = Vector{T}(undef, length(covariate))
    fold(0, slices) do offset, slice
        n = length(slice)
        view(data, (1 + offset):(n + offset)) .= slice
        offset += n
    end
    data
end

function Base.convert(::Type{Vector{Int64}}, datetimevariate::AbstractDateTimeVariate)
    slices = slice(datetimevariate, 1, length(datetimevariate), SLICELENGTH)
    data = Vector{Int64}(undef, length(datetimevariate))
    fold(0, slices) do offset, slice
        n = length(slice)
        view(data, (1 + offset):(n + offset)) .= slice
        offset += n
    end
    data
end

function Base.convert(::Type{BitArray{1}}, boolvariate::AbstractBoolVariate)
    slices = slice(boolvariate, 1, length(boolvariate), SLICELENGTH)
    data = BitArray{1}(length(boolvariate))
    fold(0, slices) do offset, slice
        n = length(slice)
        view(data, (1 + offset):(n + offset)) .= slice
        offset += n
    end
    data
end

function Base.convert(::Type{Vector{Bool}}, boolvariate::AbstractBoolVariate)
    slices = slice(boolvariate, 1, length(boolvariate), SLICELENGTH)
    data = Vector{Bool}(undef, length(boolvariate))
    fold(0, slices) do offset, slice
        n = length(slice)
        view(data, (1 + offset):(n + offset)) .= slice
        offset += n
    end
    data
end

function isordinal(factor::AbstractFactor{T}) where {T<:Unsigned}
    false
 end

function slicestring(factor::AbstractFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    levels = getlevels(factor)
    f = (i::T) -> i == zero(T) ? "" : levels[i]
    slices = slice(factor, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, String)
end

function slicestring(covariate::AbstractCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    f = (x::T) -> isnan(x) ? "" : @sprintf("%F", x)
    slices = slice(covariate, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, String)
end

function slicestring(intvariate::AbstractIntVariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Signed}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    f = (x::T) -> x == typemin(T) ? "" : @sprintf("%D", x)
    slices = slice(intvariate, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, String)
end

function tocsv(path::String, dataframe::AbstractDataFrame)
    path = abspath(path)
    iostream = open(path, "w")
    factors = dataframe.factors
    covariates = dataframe.covariates
    intvariates = dataframe.intvariates
    variates = [factors; covariates; intvariates]
    len = length(dataframe)
    strslices = zip(Tuple([slicestring(v, 1, len, SLICELENGTH) for v in variates]))
    ncol = length(variates)
    colnames = [getname(v) for v in variates]
    headersline = join(colnames, ",") * "\r\n"
    write(iostream, headersline)
    foreach(strslices) do x
        slicelen = length(x[1])
        for i in 1:slicelen
            line = join([x[j][i] for j in 1:ncol], ",") * "\r\n"
            write(iostream, line)
        end
    end
    close(iostream)
end

function Base.eltype(x::AbstractFactor{T}) where {T<:Unsigned}
    T
end

function Base.eltype(x::AbstractCovariate{T}) where {T<:AbstractFloat}
    T
end

function Base.eltype(x::AbstractBoolVariate)
    Bool
end

function Base.eltype(x::AbstractDateTimeVariate)
    Int64
end

function Base.eltype(x::AbstractIntVariate{T}) where {T<:Signed}
    T
end

function Base.filter(selector::AbstractBoolVariate, dataframe::DataFrame)
    headerpath = dataframe.headerpath
    selectorname = getname(selector)
    outtempfolder = joinpath(normpath(joinpath(dirname(headerpath), "..")), randstring(10))
    outfolder = joinpath(normpath(joinpath(dirname(headerpath), "..")), selectorname)
    mkpath(outtempfolder)
    headerjson = open(headerpath) do f
        readstring(f)
    end
    header = JSON.parse(headerjson)
    datacols = header["datacolumns"]
    selector = convert(BitArray{1}, selector)
    for datacol in datacols
        datatype = datacol["datatype"]
        len = datacol["length"]
        name = datacol["name"]
        datpath = joinpath(dirname(headerpath), datacol["filename"])
        newfilename = "$(randstring(10)).dat"
        outdatpath = joinpath(outtempfolder, newfilename)

        T =
            if datatype == "Float32"
                Float32
            elseif datatype == "Float64"
                    Float64
            elseif datatype == "Int32"
                Int32
            elseif datatype == "DateTime"
                Int64
            elseif datatype == "UInt8"
                UInt8
            elseif datatype == "UInt16"
                UInt16
            else
                UInt32
            end
        data = Vector{T}(undef, len)
        open(datpath) do f
            read!(f, data)
        end
        seldata = data[selector]
        open(outdatpath, "a") do f
            @inbounds for i in 1:(length(seldata))
                write(f, seldata[i])
            end
        end
        datacol["length"] = length(seldata)
        datacol["filename"] = newfilename
    end

    newheader = DataHeader([DataColumnInfo(x["name"], x["length"], x["filename"], x["datatype"], [string(level) for level in x["levels"]]) for x in datacols])
    newheaderjson = JSON.json(newheader)
    newheaderpath = joinpath(outtempfolder, "header.txt")
    open(newheaderpath, "a") do f
        write(f, newheaderjson)
    end
    mv(outtempfolder, outfolder; remove_destination = true)
end

function clearjumltmpdir()
    tmpjuml = joinpath(tempdir(), JUMLDIR)
    rm(tmpjuml; force = true, recursive = true)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Test

@info "Testing TimeSeriesEcon:"
using TimeSeriesEcon
cd(joinpath(dirname(pathof(TimeSeriesEcon)), "..", "test")) do
    include(joinpath(pwd(), "runtests.jl"))
end

@info "Testing ModelBaseEcon:"
using ModelBaseEcon
cd(joinpath(dirname(pathof(ModelBaseEcon)), "..", "test")) do
    include(joinpath(pwd(), "runtests.jl"))
end

@info "Testing StateSpaceEcon:"
using StateSpaceEcon
cd(joinpath(dirname(pathof(StateSpaceEcon)), "..", "test")) do
    include(joinpath(pwd(), "runtests.jl"))
end

@info "Testing Tutorials"
@testset "US_SW07" begin
    include(joinpath("..", "src", "Tutorials", "US_SW07", "main.jl"))
end

@info "All done."


text:
using Images

B = rand(1:20,3,5)
cmap = uint8(repmat(linspace(12,255,20),1,3))
img = ImageCmap(copy(B),cmap,["colorspace" => "RGB", "pixelspacing" => [2.0, 3.0], "spatialorder" => Images.yx])
imgd = convert(Image, img)

# basic information
@assert size(img) == (3,5)
@assert size(imgd) == (3,5,3)
@assert ndims(img) == 2
@assert ndims(imgd) == 3

# copy/similar
A = randn(3,5,3)
imgc = copy(img)
@assert imgc.data == img.data
imgc = copy(imgd, A)
@assert imgc.data == A
img2 = similar(img)
@assert isa(img2, ImageCmap)
@assert img2.data != img.data
img2 = similar(imgd)
@assert isa(img2, Image)

# ref/assign
prev = img[4]
@assert prev == B[4]
img[4] = prev+1
@assert img.data[4] == prev+1
@assert img[4] == prev+1

# properties
@assert colorspace(img) == "RGB"
@assert colordim(img) == 0
@assert colordim(imgd) == 3
@assert timedim(img) == 0
@assert limits(img) == (12,255)
@assert limits(imgd) == (0,255)
@assert pixelspacing(img) == [2.0, 3.0]
@assert pixelspacing(imgd) == [2.0, 3.0]

@assert sdims(img) == sdims(imgd)
@assert coords_spatial(img) == coords_spatial(imgd)
@assert size_spatial(img) == size_spatial(imgd)

@assert storageorder(img) == Images.yx
@assert storageorder(imgd) == [Images.yx, "color"]

# sub/slice
s = sub(img, 2, 1:4)
@assert ndims(s) == 2
@assert sdims(s) == 2
@assert size(s) == (1,4)
s = subim(img, 2, 1:4)
@assert ndims(s) == 2
@assert sdims(s) == 2
@assert size(s) == (1,4)
s = sliceim(img, 2, 1:4)
@assert ndims(s) == 1
@assert sdims(s) == 1
@assert size(s) == (4,)
s = sliceim(imgd, 2, 1:4, 1:3)
@assert ndims(s) == 2
@assert sdims(s) == 1
@assert colordim(s) == 2

# reslicing
D = randn(3,5,4)
sd = SliceData(D, 2)
C = slice(D, sd, 2)
@assert C == reshape(D[1:end, 2, 1:end], size(C))
reslice!(C, sd, 3)
@assert C == reshape(D[1:end, 3, 1:end], size(C))
sd = SliceData(D, 3)
C = slice(D, sd, 2)
@assert C == reshape(D[1:end, 1:end, 2], size(C))

sd = SliceData(imgd, 2)
s = sliceim(imgd, sd, 2)
@assert colordim(s) == 2
@assert colorspace(s) == "RGB"
@assert spatialorder(s) == ["y"]
@assert s.data == reshape(imgd[:,2,:], size(s))
sd = SliceData(imgd, 3)
s = sliceim(imgd, sd, 2)
@assert colordim(s) == 0
@assert colorspace(s) == "Gray"
@assert spatialorder(s) == Images.yx
@assert s.data == imgd[:,:,2]
reslice!(s, sd, 3)
@assert s.data == imgd[:,:,3]

# named indexing
@assert dimindex(imgd, "color") == 3
@assert dimindex(imgd, "y") == 1
@assert dimindex(imgd, "z") == 0
imgdp = permutedims(imgd, [3,1,2])
@assert dimindex(imgdp, "y") == 2
@assert img["y", 2, "x", 4] == B[2,4]
@assert img["x", 4, "y", 2] == B[2,4]
chan = imgd["color", 2]
Blookup = reshape(cmap[B[:],1], size(B))
@assert chan == Blookup

sd = SliceData(imgd, "x")
s = sliceim(imgd, sd, 2)
@assert spatialorder(s) == ["y"]
@assert s.data == reshape(imgd[:,2,:], size(s))
sd = SliceData(imgd, "y")
s = sliceim(imgd, sd, 2)
@assert spatialorder(s) == ["x"]
@assert s.data == reshape(imgd[2,:,:], size(s))
sd = SliceData(imgd, "x", "y")
s = sliceim(imgd, sd, 2, 1)
@assert s.data == reshape(imgd[1,2,:], 3)

# spatial order, width/height, and permutations
@assert spatialpermutation(Images.yx, imgd) == [1,2]
@assert widthheight(imgd) == (5,3)
C = convert(Array, imgd)
@assert C == imgd.data
imgd.properties["spatialorder"] = ["x", "y"]
@assert spatialpermutation(Images.xy, imgd) == [1,2]
@assert widthheight(imgd) == (3,5)
C = convert(Array, imgd)
@assert C == permutedims(imgd.data, [2,1,3])
imgd.properties["spatialorder"] = ["y", "x"]
@assert spatialpermutation(Images.xy, imgd) == [2,1]
imgd.properties["spatialorder"] = ["x", "L"]
@assert spatialpermutation(Images.xy, imgd) == [1,2]
imgd.properties["spatialorder"] = ["L", "x"]
@assert spatialpermutation(Images.xy, imgd) == [2,1]
@assert spatialpermutation(Images.xy, A) == [2,1]
@assert spatialpermutation(Images.yx, A) == [1,2]

imgd.properties["spatialorder"] = Images.yx
imgp = permutedims(imgd, ["x", "y", "color"])
@assert imgp.data == permutedims(imgd.data, [2,1,3])
imgp = permutedims(imgd, ("color", "x", "y"))
@assert imgp.data == permutedims(imgd.data, [3,2,1])


text:
module TestSmoke
include("preamble_plots.jl")

using Bifurcations.Codim1: resolved_points, SpecialPoint
using Bifurcations: examples

@testset "smoke $name" for (name, ex) in examples()
    prob = ex.prob
    solver = init(prob)
    solve!(solver)
    points = resolved_points(solver)
    @test all(isa.(points, SpecialPoint))

    @testset "show" begin
        smoke_test_solver_show(solver)
    end

    @testset "plot" begin
        sol = solve(prob)
        @test_nothrow nullshow(plot(sol))
        @test_nothrow nullshow(plot(sol; include_points=true))
        @test_nothrow nullshow(plot(sol; bif_style=Dict()))

        @test_nothrow nullshow(plot(solver))
        @test_nothrow nullshow(plot(solver; include_points=false))

        for p in points
            @test_nothrow nullshow(plot(p))
        end
    end
end

@testset "solving!" begin
    (_, ex), = examples()
    prob = ex.prob
    solver = init(prob)
    points = []
    solving!(solver) do x
        push!(points, (; map(n -> n => getproperty(x, n), propertynames(x))...))
    end
    @test length(points) > 0
    @test keys(points[1]) == (:i_sweep, :i_point, :u, :direction, :simple_bifurcation)
    @test all(p.u isa AbstractVector for p in points)
end

@testset "include(gpu.jl)" begin
    m = Module()
    @test begin
        Base.include(m, "../examples/gpu.jl")
        m.solver isa Bifurcations.Continuations.ContinuationSolver
    end
end

end  # module


text:
module Day4

using Advent: eachlinegroup, save_input
using Base.Iterators: map # shadows Base.map

const fields = Set(("byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"))
isvalid(passport) = issubset(fields, keys(passport))

function parse_passports(filename)
    passports = Set(
        Dict{String,String}(
            k => v
            for (k, v) in map(Base.Fix2(split, ":"), Iterators.flatten(map(split, group)))
        )
        for group in eachlinegroup(filename)
    )

    filter!(isvalid, passports)
end

solve_a(filename) = filename |> parse_passports |> length

function validate_passport(current)
    (parse(Int, current["byr"]) in 1920:2002) &&
    (parse(Int, current["iyr"]) in 2010:2020) &&
    (parse(Int, current["eyr"]) in 2020:2030) &&
    (
        (current["hgt"][end - 1:end] == "cm" && parse(Int, current["hgt"][1:end-2]) in 150:193) ||
        (current["hgt"][end - 1:end] == "in" && parse(Int, current["hgt"][1:end-2]) in 59:76)

    ) &&
    match(r"^#[0-9a-f]{6}$", current["hcl"]) !== nothing &&
    current["ecl"] in ("amb", "blu", "brn", "gry", "grn", "hzl", "oth") &&
    match(r"^[0-9]{9}$", current["pid"]) !== nothing
end

function solve_b(filename)
    count(parse_passports(filename)) do current
        try
            validate_passport(current)
        catch
            false
        end
    end
end

solve(filename = save_input(2020, 4)) = solve_a(filename), solve_b(filename)
solve_big() = solve("day4_big.in")

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module CoordinateDescentTest

const  NUMBER_REPEAT = 1

using CoordinateDescent
using Test
using ProximalBase
using Random
using LinearAlgebra
using SparseArrays

# test |Y - β|_2^2 + λ⋅|β|_1
@testset "Small ProxL1" begin
  Y = [1., 1.5]

  λ = 1.2
  g = ProximalBase.ProxL1(λ)
  f = CDQuadraticLoss(Matrix(1.0I, 2, 2), -Y)

  opt = CDOptions(;maxIter=100, optTol=1e-8, warmStart=true, randomize=false)
  x = SparseIterate(2)
  coordinateDescent!(x, f, g, opt)

  @test Vector(x) ≈ Vector([0., 0.3])
end

# check that the warm start and non warm start produce the same result
@testset "ProxL1" begin

    for i=1:NUMBER_REPEAT
      n = 500
      p = 50
      s = 5

      X = randn(n, p)
      β = randn(s)
      Y = X[:,1:s] * β + randn(n)

      λ = 0.01
      g = ProximalBase.ProxL1(λ)
      f = CDLeastSquaresLoss(Y, X)

      opt1 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=false)
      opt2 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=true)
      opt3 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=false)
      opt4 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=true)

      x1 = SparseIterate(sprand(p, 0.6))
      x2 = SparseIterate(sprand(p, 0.6))
      x3 = SparseIterate(sprand(p, 0.6))
      x4 = SparseIterate(sprand(p, 0.6))

      coordinateDescent!(x1, f, g, opt1)
      coordinateDescent!(x2, f, g, opt2)
      coordinateDescent!(x3, f, g, opt3)
      coordinateDescent!(x4, f, g, opt4)

      @test Vector(x1) ≈ Vector(x2) atol=1e-5
      @test Vector(x3) ≈ Vector(x2) atol=1e-5
      @test Vector(x4) ≈ Vector(x2) atol=1e-5
    end
end

@testset "AProxL1" begin

    for i=1:NUMBER_REPEAT
      n = 500
      p = 50
      s = 10

      X = randn(n, p)
      β = randn(s)
      Y = X[:,1:s] * β + randn(n)

      λ = 0.01
      g = ProximalBase.ProxL1(λ, rand(p))
      f = CDLeastSquaresLoss(Y, X)

      opt1 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=false)
      opt2 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=true)
      opt3 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=false)
      opt4 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=true)

      x1 = SparseIterate(sprand(p, 0.6))
      x2 = SparseIterate(sprand(p, 0.6))
      x3 = SparseIterate(sprand(p, 0.6))
      x4 = SparseIterate(sprand(p, 0.6))

      coordinateDescent!(x1, f, g, opt1)
      coordinateDescent!(x2, f, g, opt2)
      coordinateDescent!(x3, f, g, opt3)
      coordinateDescent!(x4, f, g, opt4)

      @test Vector(x1) ≈ Vector(x2) atol=1e-5
      @test Vector(x3) ≈ Vector(x2) atol=1e-5
      @test Vector(x4) ≈ Vector(x2) atol=1e-5
    end
end

# @testset "weighted least squares loss" begin
#
#   n, s, p = 500, 2, 50
#   gk = GaussianKernel(0.3)
#   w = zeros(n)
#   sw = zeros(n)
#
#   for i=1:NUMBER_REPEAT
#     for degree=0:2
#       Y, X, Z, betaMat = genData(n, s)
#       X = [X zeros(n, p-s)]
#       z0 = rand()
#
#       cp = p*(degree+1)
#       eX = zeros(n, cp)
#       _expand_X!(eX, X, Z, z0, degree)
#
#       @. w = evaluate(gk, Z, z0)
#       @. sw = sqrt(w)
#
#       λ = 0.001
#       g = ProximalBase.ProxL1(λ)
#       f1 = CDLeastSquaresLoss(diagm(sw) * Y, diagm(sw)*X)
#       f2 = CDWeightedLSLoss(Y, X, w)
#
#       opt1 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=false)
#       opt2 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=true, randomize=true)
#       opt3 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=false)
#       opt4 = CDOptions(;maxIter=5000, optTol=1e-8, warmStart=false, randomize=true)
#
#       x1 = SparseIterate(sprand(p, 0.6))
#       x2 = SparseIterate(sprand(p, 0.6))
#       x3 = SparseIterate(sprand(p, 0.6))
#       x4 = SparseIterate(sprand(p, 0.6))
#       x5 = SparseIterate(sprand(p, 0.6))
#       x6 = SparseIterate(sprand(p, 0.6))
#       x7 = SparseIterate(sprand(p, 0.6))
#       x8 = SparseIterate(sprand(p, 0.6))
#
#       coordinateDescent!(x1, f1, g, opt1)
#       coordinateDescent!(x2, f1, g, opt2)
#       coordinateDescent!(x3, f1, g, opt3)
#       coordinateDescent!(x4, f1, g, opt4)
#
#       coordinateDescent!(x5, f2, g, opt1)
#       coordinateDescent!(x6, f2, g, opt2)
#       coordinateDescent!(x7, f2, g, opt3)
#       coordinateDescent!(x8, f2, g, opt4)
#
#
#       @test Vector(x1) ≈ Vector(x2) atol=1e-5
#       @test Vector(x3) ≈ Vector(x2) atol=1e-5
#       @test Vector(x4) ≈ Vector(x2) atol=1e-5
#       @test Vector(x5) ≈ Vector(x2) atol=1e-5
#       @test Vector(x6) ≈ Vector(x2) atol=1e-5
#       @test Vector(x7) ≈ Vector(x2) atol=1e-5
#       @test Vector(x8) ≈ Vector(x2) atol=1e-5
#     end
#   end
#
# end

end


text:
# IndPolyhedral: OSQP implementation

using OSQP

struct IndPolyhedralOSQP{R} <: IndPolyhedral
    l::AbstractVector{R}
    A::AbstractMatrix{R}
    u::AbstractVector{R}
    mod::OSQP.Model
    function IndPolyhedralOSQP{R}(
        l::AbstractVector{R}, A::AbstractMatrix{R}, u::AbstractVector{R}
    ) where R
        m, n = size(A)
        mod = OSQP.Model()
        if !all(l .<= u)
            error("function is improper (are some bounds inverted?)")
        end
        OSQP.setup!(mod; P=SparseMatrixCSC{R}(I, n, n), l=l, A=sparse(A), u=u, verbose=false,
            eps_abs=eps(R), eps_rel=eps(R),
            eps_prim_inf=eps(R), eps_dual_inf=eps(R))
        new(l, A, u, mod)
    end
end

# properties

is_prox_accurate(::IndPolyhedralOSQP) = false

# constructors

IndPolyhedralOSQP(
    l::AbstractVector{R}, A::AbstractMatrix{R}, u::AbstractVector{R}
) where R =
    IndPolyhedralOSQP{R}(l, A, u)

IndPolyhedralOSQP(
    l::AbstractVector{R}, A::AbstractMatrix{R}, u::AbstractVector{R},
    xmin::AbstractVector{R}, xmax::AbstractVector{R}
) where R =
    IndPolyhedralOSQP([l; xmin], [A; I], [u; xmax])

IndPolyhedralOSQP(
    l::AbstractVector{R}, A::AbstractMatrix{R}, args...
) where R =
    IndPolyhedralOSQP(
        l, SparseMatrixCSC(A), R(Inf).*ones(R, size(A, 1)), args...
    )

IndPolyhedralOSQP(
    A::AbstractMatrix{R}, u::AbstractVector{R}, args...
) where R =
    IndPolyhedralOSQP(
        R(-Inf).*ones(R, size(A, 1)), SparseMatrixCSC(A), u, args...
    )

# function evaluation

function (f::IndPolyhedralOSQP{R})(x::AbstractVector{R}) where R
    Ax = f.A * x
    return all(f.l .<= Ax .<= f.u) ? zero(R) : Inf
end

# prox

function prox!(y::AbstractVector{R}, f::IndPolyhedralOSQP{R}, x::AbstractVector{R}, gamma::R=one(R)) where R
    OSQP.update!(f.mod; q=-x)
    results = OSQP.solve!(f.mod)
    y .= results.x
    return zero(R)
end

# naive prox

# we want to compute the projection p of a point x
#
# primal problem is: minimize_p (1/2)||p-x||^2 + g(Ap)
# where g is the indicator of the box [l, u]
#
# dual problem is: minimize_y (1/2)||-A'y||^2 - x'A'y + g*(y)
# can solve with (fast) dual proximal gradient method

function prox_naive(f::IndPolyhedralOSQP{R}, x::AbstractVector{R}, gamma::R=one(R)) where R
    y = zeros(R, size(f.A, 1)) # dual vector
    y1 = y
    g = IndBox(f.l, f.u)
    gstar = Conjugate(g)
    gstar_y = zero(R)
    stepsize = one(R)/opnorm(Matrix(f.A*f.A'))
    for it = 1:1e6
        w = y + (it-1)/(it+2)*(y - y1)
        y1 = y
        z = w - stepsize * (f.A * (f.A'*w - x))
        y, = prox(gstar, z, stepsize)
        if norm(y-w)/(1+norm(w)) <= 1e-12 break end
    end
    p = -f.A'*y + x
    return p, zero(R)
end


text:
module ModelExploration
export Explore, Dim, Free

using Catlab.CategoricalAlgebra

"""
A dimension implicitly specifies a metric space of structures, (or, think of it
as a preorder if we don't care about distances) ordered by "complexity" in some
context-dependent sense. We generally want to explore forward in the dimension
only to the extent necessary to model a phenomenon.

The data of a Dim allows us to lazily explore this preorder (the set of
structures implied may be infinite, such as an "N-city model" in epidemiology).
"""
abstract type Dim end

"""

"""
struct Finite <: Dim
  vals::Vector{StructACSetTransformation}
  ordering::Vector{Vector{Int}}
  function Finite(vs::Vector{StructACSetTransformation})
    new(vs, [Int[] for _ in vs])
  end
end

"""
A sequence of models indexed by ℕ
"""
struct Linear <: Dim
  fun::Function # ℕ → CSetTransformation
end

"""
All possible models up to ACSet isomorphism
"""
struct Free <: Dim end # use ModelEnumeration.jl

struct Rewrite <: Dim
  rules::Vector{Pair{ACSetTransformation}}
end

"""
Reductions of model space from all possible ACSets of a schema
"""
abstract type Constraint end

struct Chase <: Constraint
  rules::Vector{EqConstraint}
end

"""
Last resort for a constraint, because model is hidden implicit in code
"""
struct Filter <: Constraint
  f::Function # StructACsetTransformation -> Bool
end

"""
We can only take pullbacks between slices, but we can default to a trivial
slice if none is explicitly provided.
"""
struct Explore
  csettype::Type
  slice::StructACSet
  dims::Vector{Dim}
  constr::Vector{Constraint}
end

struct ModelSpace
  dims::Vector{Dim}
end

struct ModelPoint
  coords::Vector{UInt128} # model hashes in each dimension
end

struct ExploreState
  seen::Set{UInt128}
  space::ModelSpace
end

# Magnitude and direction of greatest improvement?
struct Loss
end

eval(mp::ModelPoint)::Loss = Loss()

stopping_criteria(es::ExploreState)::Bool = false

next_mp(es::ExploreState)::Vector{ModelPoint} = []

end # module


text:
using Printf

_drop_moi(s) = replace(string(s), "MathOptInterface." => "")

struct _PrintOptions{T<:MIME}
    simplify_coefficients::Bool
    default_name::String
    print_types::Bool

    """
        _PrintOptions(
            mime::MIME;
            simplify_coefficients::Bool = false,
            default_name::String = "v",
            print_types::Bool = true,
        )

    A struct to control options for printing.

    ## Arguments

     * `simplify_coefficients` : Simplify coefficients if possible by omitting
       them or removing trailing zeros.
     * `default_name` : The name given to variables with an empty name.
     * `print_types` : Print the MOI type of each function and set for clarity.
    """
    function _PrintOptions(
        mime::MIME;
        simplify_coefficients::Bool = false,
        default_name::String = "v",
        print_types::Bool = true,
    )
        return new{typeof(mime)}(
            simplify_coefficients,
            default_name,
            print_types,
        )
    end
end

function _to_string(mime::MIME, args...; kwargs...)
    return _to_string(_PrintOptions(mime), args...; kwargs...)
end

#------------------------------------------------------------------------
# Math Symbols
#------------------------------------------------------------------------

# REPL-specific symbols
# Anything here: https://en.wikipedia.org/wiki/Windows-1252
# should probably work fine on Windows

_to_string(::_PrintOptions, ::typeof(*)) = "*"
_to_string(::_PrintOptions{MIME"text/latex"}, ::typeof(*)) = "\\times "

function _to_string(::_PrintOptions, ::typeof(^), n::Int)
    return n == 2 ? "²" : string('^', n)
end

function _to_string(::_PrintOptions{MIME"text/latex"}, ::typeof(^), n::Int)
    return string('^', n)
end

_to_string(::_PrintOptions, ::typeof(in)) = @static Sys.iswindows() ? "in" : "∈"

#------------------------------------------------------------------------
# Functions
#------------------------------------------------------------------------

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    v::MOI.VariableIndex,
)
    var_name = MOI.get(model, MOI.VariableName(), v)
    if isempty(var_name)
        return string(options.default_name, "[", v.value, "]")
    else
        return var_name
    end
end

function _to_string(
    options::_PrintOptions{MIME"text/latex"},
    model::MOI.ModelLike,
    v::MOI.VariableIndex,
)
    var_name = MOI.get(model, MOI.VariableName(), v)
    if isempty(var_name)
        return string(options.default_name, "_{", v.value, "}")
    end
    # We need to escape latex math characters that appear in the name.
    # However, it's probably impractical to catch everything, so let's just
    # escape the common ones:
    # Escape underscores to prevent them being treated as subscript markers.
    var_name = replace(var_name, "_" => "\\_")
    # Escape carets to prevent them being treated as superscript markers.
    var_name = replace(var_name, "^" => "\\^")
    # Convert any x[args] to x_{args} so that indices on x print as subscripts.
    m = match(r"^(.*)\[(.+)\]$", var_name)
    if m !== nothing
        var_name = m[1] * "_{" * m[2] * "}"
    end
    return var_name
end

function _shorten(options::_PrintOptions, x::Float64)
    if options.simplify_coefficients && isinteger(x)
        return string(round(Int, x))
    end
    return string(x)
end

_shorten(::_PrintOptions, x) = string("(", x, ")")

"""
    _to_string(options::_PrintOptions, c::Real, x::String)

Write a coefficient-name pair to string. There are a few cases to handle.

          | is_first | !is_first
    -----------------------------
    +2.1x | "2.1 x"  | " + 2.1 x"
    -2.1x | "-2.1 x" | " - 2.1 x"
    -----------------------------
    +2.0x | "2 x"    | " + 2 x"
    -2.0x | "-2 x"   | " - 2 x"
    +1.0x | "x"      | " + x"
    -1.0x | "-x"     | " - x"
"""
function _to_string(options::_PrintOptions, c::Real, x::String; is_first::Bool)
    prefix = if is_first
        c < 0 ? "-" : ""
    else
        c < 0 ? " - " : " + "
    end
    s = _shorten(options, abs(c))
    if options.simplify_coefficients && s == "1"
        return string(prefix, x)
    else
        return string(prefix, s, " ", x)
    end
end

function _to_string(
    options::_PrintOptions,
    c::Number,
    x::String;
    is_first::Bool,
)
    return string(is_first ? "" : " + ", _shorten(options, c), " ", x)
end

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    term::MOI.ScalarAffineTerm;
    is_first::Bool,
)
    name = _to_string(options, model, term.variable)
    return _to_string(options, term.coefficient, name; is_first = is_first)
end

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    f::MOI.ScalarAffineFunction,
)
    s = _shorten(options, f.constant)
    if options.simplify_coefficients && iszero(f.constant)
        s = ""
    end
    is_first = isempty(s)
    for term in f.terms
        s *= _to_string(options, model, term; is_first = is_first)
        is_first = false
    end
    return s
end

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    term::MOI.ScalarQuadraticTerm;
    is_first::Bool,
)
    name_1 = _to_string(options, model, term.variable_1)
    name_2 = _to_string(options, model, term.variable_2)
    # Be careful here when printing the coefficient. ScalarQuadraticFunction
    # assumes an additional 0.5 factor!
    coef = term.coefficient
    name = if term.variable_1 == term.variable_2
        coef /= 2
        string(name_1, _to_string(options, ^, 2))
    else
        string(name_1, _to_string(options, *), name_2)
    end
    return _to_string(options, coef, name; is_first = is_first)
end

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    f::MOI.ScalarQuadraticFunction,
)
    s = _shorten(options, f.constant)
    if options.simplify_coefficients && iszero(f.constant)
        s = ""
    end
    is_first = isempty(s)
    for term in f.affine_terms
        s *= _to_string(options, model, term; is_first = is_first)
        is_first = false
    end
    for term in f.quadratic_terms
        s *= _to_string(options, model, term; is_first = is_first)
        is_first = false
    end
    return s
end

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction,
)
    rows = map(fi -> _to_string(options, model, fi), eachscalar(f))
    max_length = maximum(length.(rows))
    s = join(map(r -> string("│", rpad(r, max_length), "│"), rows), '\n')
    return string(
        "┌",
        rpad("", max_length),
        "┐\n",
        s,
        "\n└",
        rpad("", max_length),
        "┘",
    )
end

function _to_string(
    options::_PrintOptions{MIME"text/latex"},
    model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction,
)
    return string(
        "\\begin{bmatrix}\n",
        join(
            map(fi -> _to_string(options, model, fi), eachscalar(f)),
            "\\\\\n",
        ),
        "\\end{bmatrix}",
    )
end

#------------------------------------------------------------------------
# Sets
#------------------------------------------------------------------------

function _to_string(options::_PrintOptions, set::MOI.LessThan)
    return string("<= ", _shorten(options, set.upper))
end

function _to_string(options::_PrintOptions{MIME"text/latex"}, set::MOI.LessThan)
    return string("\\le ", _shorten(options, set.upper))
end

function _to_string(options::_PrintOptions, set::MOI.GreaterThan)
    return string(">= ", _shorten(options, set.lower))
end

function _to_string(
    options::_PrintOptions{MIME"text/latex"},
    set::MOI.GreaterThan,
)
    return string("\\ge ", _shorten(options, set.lower))
end

function _to_string(options::_PrintOptions, set::MOI.EqualTo)
    return string("== ", _shorten(options, set.value))
end

function _to_string(options::_PrintOptions{MIME"text/latex"}, set::MOI.EqualTo)
    return string("= ", _shorten(options, set.value))
end

function _to_string(options::_PrintOptions, set::MOI.Interval)
    return string(
        _to_string(options, in),
        " [",
        _shorten(options, set.lower),
        ", ",
        _shorten(options, set.upper),
        "]",
    )
end

function _to_string(options::_PrintOptions{MIME"text/latex"}, set::MOI.Interval)
    return string(
        "\\in \\[",
        _shorten(options, set.lower),
        ", ",
        _shorten(options, set.upper),
        "\\]",
    )
end

function _to_string(options::_PrintOptions, ::MOI.ZeroOne)
    return string(_to_string(options, in), " {0, 1}")
end

_to_string(::_PrintOptions{MIME"text/latex"}, ::MOI.ZeroOne) = "\\in \\{0, 1\\}"

function _to_string(options::_PrintOptions, ::MOI.Integer)
    return string(_to_string(options, in), " ℤ")
end

function _to_string(::_PrintOptions{MIME"text/latex"}, ::MOI.Integer)
    return "\\in \\mathbb{Z}"
end

function _to_string(options::_PrintOptions, set::MOI.AbstractSet)
    return string(_to_string(options, in), " ", _drop_moi(set))
end

function _to_string(::_PrintOptions{MIME"text/latex"}, set::MOI.AbstractSet)
    set_str = replace(replace(_drop_moi(set), "{" => "\\{"), "}" => "\\}")
    return string("\\in \\text{", set_str, "}")
end

#------------------------------------------------------------------------
# Constraints
#------------------------------------------------------------------------

function _to_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    cref::MOI.ConstraintIndex,
)
    f = MOI.get(model, MOI.ConstraintFunction(), cref)
    s = MOI.get(model, MOI.ConstraintSet(), cref)
    return string(_to_string(options, model, f), " ", _to_string(options, s))
end

#------------------------------------------------------------------------
# Nonlinear constraints
#------------------------------------------------------------------------

"""
    _VariableNode

A type used to work-around the default printing of Julia expressions.

Without this type, if we subsititued the variable names into the expression
and then converted to a string, each variable would be printed with enclosing
`"`.

To work-around this, create a new type and overload `show`.
"""
struct _VariableNode
    x::String
end
Base.show(io::IO, x::_VariableNode) = print(io, x.x)

_replace_names(::_PrintOptions, ::MOI.ModelLike, x, ::Any) = x
function _replace_names(
    options::_PrintOptions,
    model::MOI.ModelLike,
    x::Expr,
    lookup,
)
    if x.head == :ref
        return get!(lookup, x.args[2]) do
            return _VariableNode(_to_string(options, model, x.args[2]))
        end
    else
        for (i, arg) in enumerate(x.args)
            x.args[i] = _replace_names(options, model, arg, lookup)
        end
    end
    return x
end

function _replace_nonlinear_latex(::_PrintOptions{MIME"text/latex"}, s::String)
    s = replace(s, " * " => " \\times ")
    s = replace(s, " >= " => " \\ge ")
    s = replace(s, " <= " => " \\le ")
    s = replace(s, " == " => " = ")
    return s
end
_replace_nonlinear_latex(::_PrintOptions, s::String) = s

function _print_nonlinear_constraints(
    io::IO,
    options::_PrintOptions{MIME"text/plain"},
    model::MOI.ModelLike,
    block::MOI.NLPBlockData,
)
    lookup = Dict{MOI.VariableIndex,_VariableNode}()
    if options.print_types
        println(io, "\nNonlinear")
    end
    has_expr = :ExprGraph in MOI.features_available(block.evaluator)
    for (i, bound) in enumerate(block.constraint_bounds)
        if has_expr
            ex = MOI.constraint_expr(block.evaluator, i)
            println(io, " ", _replace_names(options, model, ex, lookup))
        else
            println(io, " ", bound.lower, " <= g_$(i)(x) <= ", bound.upper)
        end
    end
end

function _print_nonlinear_constraints(
    io::IO,
    options::_PrintOptions{MIME"text/latex"},
    model::MOI.ModelLike,
    block::MOI.NLPBlockData,
)
    lookup = Dict{MOI.VariableIndex,_VariableNode}()
    if options.print_types
        println(io, " & \\text{Nonlinear} \\\\")
    end
    has_expr = :ExprGraph in MOI.features_available(block.evaluator)
    for (i, bound) in enumerate(block.constraint_bounds)
        if has_expr
            ex = MOI.constraint_expr(block.evaluator, i)
            nl_c = string(_replace_names(options, model, ex, lookup))
            println(io, " & ", _replace_nonlinear_latex(options, nl_c), " \\\\")
        else
            println(
                io,
                "& ",
                bound.lower,
                " \\le g_$(i)(x) \\le ",
                bound.upper,
                " \\\\",
            )
        end
    end
end

function _print_nonlinear_constraints(
    ::IO,
    ::_PrintOptions,
    ::MOI.ModelLike,
    ::Nothing,
)
    return
end

#------------------------------------------------------------------------
# ObjectiveFunction
#------------------------------------------------------------------------

function _objective_function_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    ::Nothing,
)
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    f = MOI.get(model, MOI.ObjectiveFunction{F}())
    return _drop_moi(F), _to_string(options, model, f)
end

function _objective_function_string(
    options::_PrintOptions,
    model::MOI.ModelLike,
    block::MOI.NLPBlockData,
)
    lookup = Dict{MOI.VariableIndex,_VariableNode}()
    if block.has_objective
        f = "f(x)"
        if :ExprGraph in MOI.features_available(block.evaluator)
            ex = MOI.objective_expr(block.evaluator)
            f = string(_replace_names(options, model, ex, lookup))
            f = _replace_nonlinear_latex(options, f)
        end
        return "Nonlinear", f
    else
        return _objective_function_string(options, model, nothing)
    end
end

#------------------------------------------------------------------------
# MOI.ModelLike
#------------------------------------------------------------------------

function _nlp_block(model::MOI.ModelLike)
    try
        block = MOI.get(model, MOI.NLPBlock())
        if block === nothing
            return
        end
        if :ExprGraph in MOI.features_available(block.evaluator)
            MOI.initialize(block.evaluator, [:ExprGraph])
        end
        return block
    catch
        return nothing
    end
end

"""
    _print_model(
        io::IO,
        options::_PrintOptions{MIME"text/plain"},
        model::MOI.ModelLike,
    )

Print a plain-text formulation of `model` to `io`.
"""
function _print_model(
    io::IO,
    options::_PrintOptions{MIME"text/plain"},
    model::MOI.ModelLike,
)
    nlp_block = _nlp_block(model)
    sense = MOI.get(model, MOI.ObjectiveSense())
    if sense == MOI.FEASIBILITY_SENSE
        println(io, "Feasibility")
    else
        F, f = _objective_function_string(options, model, nlp_block)
        sense_s = sense == MOI.MIN_SENSE ? "Minimize" : "Maximize"
        if options.print_types
            println(io, sense_s, " ", F, ":\n ", f)
        else
            println(io, sense_s, ": ", f)
        end
    end
    println(io, "\nSubject to:")
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        if options.print_types
            println(io, "\n$(_drop_moi(F))-in-$(_drop_moi(S))")
        end
        for cref in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            s = _to_string(options, model, cref)
            println(io, " ", replace(s, '\n' => "\n "))
        end
    end
    _print_nonlinear_constraints(io, options, model, nlp_block)
    return
end

"""
    _print_model(
        io::IO,
        options::_PrintOptions{MIME"text/latex"},
        model::MOI.ModelLike,
    )

Print a LaTeX formulation of `model` to `io`.
"""
function _print_model(
    io::IO,
    options::_PrintOptions{MIME"text/latex"},
    model::MOI.ModelLike,
)
    nlp_block = _nlp_block(model)
    println(io, "\$\$ \\begin{aligned}")
    sense = MOI.get(model, MOI.ObjectiveSense())
    if sense == MOI.FEASIBILITY_SENSE
        println(io, "\\text{feasibility}\\\\")
    else
        F, f = _objective_function_string(options, model, nlp_block)
        sense_s = sense == MOI.MIN_SENSE ? "min" : "max"
        println(io, "\\", sense_s, "\\quad & ", f, " \\\\")
    end
    println(io, "\\text{Subject to}\\\\")
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        if options.print_types
            println(io, " & \\text{$(_drop_moi(F))-in-$(_drop_moi(S))} \\\\")
        end
        for cref in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            println(io, " & ", _to_string(options, model, cref), " \\\\")
        end
    end
    _print_nonlinear_constraints(io, options, model, nlp_block)
    return print(io, "\\end{aligned} \$\$")
end

#------------------------------------------------------------------------
# Latex
#------------------------------------------------------------------------

struct _LatexModel{T<:MOI.ModelLike}
    model::T
    kwargs::Any
end

"""
    latex_formulation(model::MOI.ModelLike; kwargs...)

Wrap `model` in a type so that it can be pretty-printed as `text/latex` in a
notebook like IJulia, or in Documenter.

To render the model, end the cell with `latex_formulation(model)`, or call
`display(latex_formulation(model))` in to force the display of the model from
inside a function.

Possible keyword arguments are:

 * `simplify_coefficients` : Simplify coefficients if possible by omitting
   them or removing trailing zeros.
 * `default_name` : The name given to variables with an empty name.
 * `print_types` : Print the MOI type of each function and set for clarity.
"""
latex_formulation(model::MOI.ModelLike; kwargs...) = _LatexModel(model, kwargs)

function Base.show(io::IO, model::_LatexModel)
    return _print_model(
        io,
        _PrintOptions(MIME("text/latex"); model.kwargs...),
        model.model,
    )
end

Base.show(io::IO, ::MIME"text/latex", model::_LatexModel) = show(io, model)

function Base.print(model::MOI.ModelLike; kwargs...)
    for d in Base.Multimedia.displays
        if Base.Multimedia.displayable(d, "text/latex") &&
           startswith("$(typeof(d))", "IJulia.")
            return display(d, "text/latex", latex_formulation(model; kwargs...))
        end
    end
    return print(stdout, model; kwargs...)
end

function Base.print(io::IO, model::MOI.ModelLike; kwargs...)
    return _print_model(io, _PrintOptions(MIME("text/plain"); kwargs...), model)
end


text:
module UnitlessFlatten

using Reexport, Unitful

@reexport using Flatten

import Flatten: flatten, reconstruct, retype, update!

flatten(x::Unitful.Quantity) = (x.val,)
reconstruct(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1
retype(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1
update!(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1

end # module


text:
using LinearAlgebra
using Distributed

############
# Utilities
############
"""
split(N, P, i) -> from:to

Find the ith range when 1:N is split into P consecutive parts of roughly equal size
"""
function split(N, P, i)
    base, rem = divrem(N, P)
    from = (i - 1) * base + min(rem, i - 1) + 1
    to = from + base - 1 + (i ≤ rem ? 1 : 0)
    from : to
end

#########################
# Modeling (demigration)
#########################

# Auxiliary functions
function kirchmod_kernel!(data, model, tt, ot, dt, nt)
    its = @. round(Int, (tt-ot)/dt)
    its = @. min(max(its, 0), nt) + 1
    for ixyz = 1:length(its)
        data[its[ixyz]] += model[ixyz]
    end
end

function kirchmod_get_axes(model::AbstractArray{T, N},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,N,M}
    Nzxy = size(model)

    nr = size(trav_r)[end]
    ns = size(trav_s)[end]
    nzxy = size(trav_r)[1:M-1]
    NZXY = size(trav_s)[1:M-1]

    ot, dt, nt = t[1], t[2]-t[1], length(t)

    if nzxy != Nzxy
        throw(DimensionMismatch("1st $(M-1) dimensions of trav_r must be the same as that of model"))
    end
    if nzxy != NZXY
        throw(DimensionMismatch("1st $(M-1) dimensions of trav_r and trav_s must be the same"))
    end
    return ot, dt, nt, nr, ns, nzxy
end

# Threaded
kirchmod_thread(model::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmod_thread(model, t, trav_r, trav_r)

function kirchmod_thread(model::AbstractArray{T, N},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,N,M}
    ot, dt, nt, nr, ns, nzxy = kirchmod_get_axes(model, t, trav_r, trav_s)

    colons = [Colon() for i=1:M-1]

    nthreads = Threads.nthreads()
    v   = Vector{Vector{Matrix{T}}}(undef, nthreads)
    ind = Vector{Vector{Int}}(undef, nthreads)
    Threads.@threads for i in 1:nthreads
        v[i] = Matrix{T}[]
        ind[i] = Int[]
    end
    Threads.@threads for i in 1:nthreads
        P = Threads.threadid()
        range = split(nr, nthreads, P)
        @fastmath @inbounds for ir in range
            rec = zeros(T, ns, nt+1)
            for is=1:ns
                kirchmod_kernel!(view(rec, is, :), model,
                                 trav_r[colons...,ir] .+ trav_s[colons...,is],
                                 ot, dt, nt)
            end
            push!(v[P], rec[:,1:nt])
            push!(ind[P],ir)
        end
    end
    v = vcat(v...)
    ind = vcat(ind...)
    data = zeros(T, nr, ns, nt)
    @inbounds @simd for i=1:nr
        data[ind[i], :, :] = v[i]
    end
    return data
end

# Parallel
kirchmod_par(model::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmod_par(model, t, trav_r, trav_r)
function kirchmod_par(model::AbstractArray{T, N},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,N,M}
    ot, dt, nt, nr, ns, nzxy = kirchmod_get_axes(model, t, trav_r, trav_s)

    colons = [Colon() for i=1:M-1]
    data = @distributed vcat for ir=1:nr
        rec = zeros(T, ns, nt+1)
        @fastmath @inbounds @simd for is=1:ns
            kirchmod_kernel!(view(rec, is, :), model,
                             trav_r[colons...,ir] .+ trav_s[colons...,is],
                             ot, dt, nt)
        end
        rec
    end
    return permutedims(reshape(data[:,1:nt], ns, nr, nt), (2, 1, 3))
end

# Serial
kirchmod(model::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmod(model, t, trav_r, trav_r)
function kirchmod(model::AbstractArray{T, N},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,N,M}
    ot, dt, nt, nr, ns, nzxy = kirchmod_get_axes(model, t, trav_r, trav_s)

    colons = [Colon() for i=1:M-1]
    data = zeros(T, nr, ns, nt+1)
    @fastmath @inbounds for ir=1:nr
        @simd for is=1:ns
            kirchmod_kernel!(view(data, ir, is, :), model,
                             trav_r[colons...,ir] .+ trav_s[colons...,is],
                             ot, dt, nt)
        end
    end
    return data[:,:,1:nt]
end

##############
# Migration
##############

# Auxiliary functions
function kirchmig_kernel!(model, data, tt, ot, dt, nt)
    its = @. round(Int, (tt-ot)/dt)
    its = @. min(max(its, 0), nt) + 1
    model .+= data[its[:]]
end

function kirchmig_get_axes(data::AbstractArray{T, 3},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,M}
    nR = size(trav_r)[end]
    nS = size(trav_s)[end]
    nzxy = size(trav_r)[1:M-1]
    NZXY = size(trav_s)[1:M-1]

    nr, ns, nt = size(data)
    ot, dt, nT = t[1], t[2]-t[1], length(t)

    if nr != nR
        throw(DimensionMismatch("1st dimension of data (receivers) must be the same as last dimension of trav_r"))
    end
    if ns != nS
        throw(DimensionMismatch("2nd dimension of data (sources) must be the same as last dimension of trav_s"))
    end
    if nzxy != NZXY
        throw(DimensionMismatch("1st $(M-1) dimensions of trav_r and trav_s must be the same"))
    end

    return ot, dt, nt, nr, ns, nzxy
end

# Threaded
kirchmig_thread(model::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmig_thread(model, t, trav_r, trav_r)

function kirchmig_thread(data::AbstractArray{T, 3},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,M}

    ot, dt, nt, nr, ns, nzxy = kirchmig_get_axes(data, t, trav_r, trav_s)

    data_ = zeros(T, nr, ns, nt+1)
    data_[:,:,1:nt] = data
    colons = [Colon() for i=1:M-1]

    nthreads = Threads.nthreads()
    v = Vector{Vector{Vector{T}}}(undef, nthreads)
    Threads.@threads for i in 1:nthreads
        v[i] = Vector{T}[]
    end
    Threads.@threads for i in 1:nthreads
        mod = zeros(T, prod(nzxy))
        P = Threads.threadid()
        range = split(nr, nthreads, P)
        @fastmath @inbounds for ir in range
            for is=1:ns
                kirchmig_kernel!(mod, data_[ir, is, :],
                                 trav_r[colons..., ir] .+ trav_s[colons..., is],
                                 ot, dt, nt)
            end
        end
        push!(v[P], mod)
    end
    return reshape(sum(vcat(v...)), nzxy...)
end

# Parallel
kirchmig_par(data::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmig_par(data, t, trav_r, trav_r)

function kirchmig_par(data::AbstractArray{T, 3},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,M}
    ot, dt, nt, nr, ns, nzxy = kirchmig_get_axes(data, t, trav_r, trav_s)

    data_ = zeros(T, nr, ns, nt+1)
    data_[:,:,1:nt] = data
    colons = [Colon() for i=1:M-1]
    model = @distributed (+) for ir=1:nr
        mod = zeros(T, prod(nzxy))
        @fastmath @inbounds @simd for is=1:ns
            kirchmig_kernel!(mod, data_[ir, is, :],
                             trav_r[colons..., ir] .+ trav_s[colons..., is],
                             ot, dt, nt)
        end
        mod
    end
    return reshape(model, nzxy...)
end

# Serial
kirchmig(model::AbstractArray, t::AbstractVector,
         trav_r::AbstractArray) = kirchmig(model, t, trav_r, trav_r)

function kirchmig(data::AbstractArray{T, 3},
                  t::AbstractVector,
                  trav_r::AbstractArray{<:Real, M},
                  trav_s::AbstractArray{<:Real, M}) where {T,M}
    ot, dt, nt, nr, ns, nzxy = kirchmig_get_axes(data, t, trav_r, trav_s)

    data_ = zeros(T, nr, ns, nt+1)
    data_[:,:,1:nt] = data
    colons = [Colon() for i=1:M-1]
    model = zeros(T, prod(nzxy))
    @fastmath @inbounds for ir=1:nr
        @simd for is=1:ns
            kirchmig_kernel!(model, data_[ir, is, :],
                             trav_r[colons..., ir] .+ trav_s[colons..., is],
                             ot, dt, nt)
        end
    end
    return reshape(model, nzxy...)
end


text:
#using Random
#using RandomExtensions
#
#const rng = MersenneTwister()
#const rand_seed = rand(UInt128)

@testset "Elliptic curves over finite fields" begin

  R1 = GF(23)
  R2, a2 = FlintFiniteField(23, 1, "a")
  R3, a3 = FlintFiniteField(fmpz(23), 1, "a")
  R4, a4 = FlintFiniteField(23, 2, "a")

  E1 = EllipticCurve(R1, [2, 3])
  E2 = EllipticCurve(R2, [2, 3])
  E3 = EllipticCurve(R3, [2, 3])
  E4 = EllipticCurve(R4, [2, 3])
  E5 = EllipticCurve(R4, [1, 2, 3, 0, 6])

  @testset "Random point construction" begin
    @inferred rand(E1)
    @inferred rand(E2)
    @inferred rand(E3)
    @inferred rand(E4)
    @inferred rand(E5)

    T = EllCrvPt{gfp_elem}
    @test rand(rng, E1) isa T
    @test rand(rng, E1, 3) isa Vector{T}

    Random.seed!(rng, rand_seed)
    a = rand(rng, E1)
    Random.seed!(rng, rand_seed)
    @test a == rand(rng, E1)
  end

  @testset "Order computation (Exhaustive_search)" begin
    @test 24 == @inferred order_via_exhaustive_search(E1)
  end

  @testset "Order computation (Legendre)" begin
    @test 24 == @inferred order_via_legendre(E1)
  end

  @testset "Order computation (BSGS)" begin
    @test 24 in @inferred order_via_bsgs(E1)
    @test 24 in @inferred order_via_bsgs(E2)
    @test 24 in @inferred order_via_bsgs(E3)
    @test 576 in @inferred order_via_bsgs(E4)
  end

  @testset "Hasse interval" begin
    l = @inferred hasse_interval(E1)
    @test l[1] <= 24 && 24 <= l[2]
    l = @inferred hasse_interval(E2)
    @test l[1] <= 24 && 24 <= l[2]
    l = @inferred hasse_interval(E3)
    @test l[1] <= 24 && 24 <= l[2]
    l = @inferred hasse_interval(E4)
    @test l[1] <= 576 && 576 <= l[2]
  end

  @testset "Trace of Frobenius" begin
   E = EllipticCurve(GF(7,2), [1, 2, 3, 4, 5])
   @test -13 == @inferred trace_of_frobenius(E)
  end

  @testset "Schoofs algorithm" begin
    @test 24 == @inferred order_via_schoof(E1)
    @test 24 == @inferred order_via_schoof(E2)
    @test 24 == @inferred order_via_schoof(E3)
    @test 576 == @inferred order_via_schoof(E4)
  end

  @testset "Point counting" begin
    RR = GF(2)
    E = EllipticCurve(RR,[1, 1, 0, 0, 1])
    @test 2 == @inferred order(E)
    RR = GF(3)
    E = EllipticCurve(RR, [1, 1])
    @test 4 == @inferred order(E)

    RR = GF(3,6)
    E = EllipticCurve(RR, [1,2,0,1,1])
    @test 784 == @inferred order(E)

    @test 24 == @inferred order(E1)
    @test 24 == @inferred order(E2)
    @test 24 == @inferred order(E3)
    @test 576 == @inferred order(E4)
  end
end


text:
# # Wiring diagrams in Graphviz
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/graphics/graphviz_wiring_diagrams.ipynb)
#
# Catlab can draw wiring diagrams using the `dot` program in
# [Graphviz](https://www.graphviz.org/). This feature requires that Graphviz be
# installed, but does not require any additional Julia packages.

using Catlab.WiringDiagrams, Catlab.Graphics

# ### Symmetric monoidal category

using Catlab.Doctrines

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A)
h = Hom(:h, otimes(A,B), otimes(A,B));

# To start, here are a few very simple examples.

to_graphviz(f)
#-
to_graphviz(compose(f,g))
#-
to_graphviz(otimes(f,g))

# In the next example, notice how Graphviz automatically "untwists" the double
# braiding to minimize edge crossings.

to_graphviz(compose(braid(A,A), otimes(f,f), braid(B,B)))

# Here is a larger composite morphism.

composite = compose(otimes(g,f), h, otimes(f,g))
to_graphviz(composite)

# We can tell Graphviz to lay out the diagrams in the horizontal (left-to-right)
# direction.

to_graphviz(composite, direction=:horizontal)

# When working with very large diagrams (larger than the ones shown here), it is
# sometimes convenient to omit the ports of the outer box and any wires attached
# to them.

to_graphviz(composite, outer_ports=false)

# ### Biproduct category

A, B = Ob(FreeBiproductCategory, :A, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A);

# By default, copies and merges are drawn the way they are represented
# internally, as multiple wires.

f1 = compose(mcopy(A), otimes(f,f))
to_graphviz(f1)
#-
f2 = compose(mcopy(A), otimes(f,f), mmerge(B))
to_graphviz(f2)

# To draw nodes for copies and merges, we need to add junctions to the wiring
# diagram.

to_graphviz(add_junctions!(to_wiring_diagram(f1)))
#-
to_graphviz(add_junctions!(to_wiring_diagram(f2)))


text:
function obimPravougaonika(a,b)
    return 2*a+2*b
end

function matrice(A,B)
    return A + B, A - B, A * B
end

a = 10
if a < 5
    print("oops")
    print("ppps")
else
    print("aahahahah")
end


text:
# Rock 🗿, Paper 📃, Scissors ✂ Game in Julia

function play_rock_paper_scissors()
    moves = ["🗿", "📃", "✂"]
    computer_move = moves[rand(1:3)]

    # Base.prompt is similar to readline which we used before
    human_move = Base.prompt("Please enter 🗿, 📃, or ✂")
    # Appends a ": " to the end of the line by default

    print("Rock...")
    sleep(0.8)

    print("Paper...")
    sleep(0.8)

    print("Scissors...")
    sleep(0.8)

    print("Shoot!\n")

    if computer_move == human_move
        print("You tied, please try again")
    elseif computer_move == "🗿" && human_move == "✂"
        print("You lose, the computer won with 🗿, please try again")
    elseif computer_move == "📃" && human_move == "🗿"
        print("You lose, the computer won with 📃, please try again")
    elseif computer_move == "✂" && human_move == "📃"
        print("You lose, the computer won with ✂, please try again")
    else
        print("You won, the computer lost with $computer_move, nice work!")
    end

end