text:
# Protein Data Bank (PDB) molecular structure file format

module PDB

using Printf
using ..Dorothy
using ..Dorothy.Geometry
using ..Dorothy.PBC
using ..Dorothy.Properties
using ..Dorothy.Utils
using Formats

export register_pdb, readpdb!, writepdb

function register_pdb()
	Formats.addformat("structure/x-pdb")
	Formats.addextension("structure/x-pdb", ".pdb")
	Formats.addsignature("structure/x-pdb", "HEADER")
	Formats.addsignature("structure/x-pdb", "TITLE ")
	Formats.addsignature("structure/x-pdb", "CRYST1")
	Formats.addsignature("structure/x-pdb", "ATOM  ")
	Formats.addsignature("structure/x-pdb", "HETATM")
	Formats.addreader("structure/x-pdb", DorothyIO())
	Formats.addwriter("structure/x-pdb", DorothyIO())
end

Base.read(::DorothyIO, ::MIME"structure/x-pdb", io::IO, args...; kwargs...) =
		readpdb!(io, MolecularModel(), args...; kwargs...)

Base.read!(::DorothyIO, ::MIME"structure/x-pdb", io::IO, model::MolecularModel,
		args...; kwargs...) = readpdb!(io, model, args...; kwargs...)

Base.write(::DorothyIO, ::MIME"structure/x-pdb", io::IO,
		model::ParticleCollection, args...; kwargs...) =
		writepdb(io, model, args...; kwargs...)

pdb_particle_name(resname::AbstractString) =
		length(resname) > 3 ? resname : @sprintf(" %-3s", resname)

function readpdb!(io::IO, model::MolecularModel;
		lnoffset::Ref{<:Integer} = Ref(0))
	ids = get!(model, :ids, undef)
	names = get!(model, :names, undef)
	resids = get!(model, :resids, undef)
	resnames = get!(model, :resnames, undef)
	chainids = get!(model, :chainids, undef)
	R = get!(model, :R, undef)
	occupancies = get!(model, :occupancies, undef)
	bfactors = get!(model, :bfactors, undef)
	elements = get!(model, :elements, undef)
	readpdb!(io, model, ids, names, resids, resnames, chainids, R, occupancies,
			bfactors, elements, lnoffset)
end

function readpdb!(io::IO, model::MolecularModel, ids::AbstractVector{<:Integer},
		names::AbstractVector{<:AbstractString},
		resids::AbstractVector{<:Integer},
		resnames::AbstractVector{<:AbstractString},
		chainids::AbstractVector{<:AbstractString},
		R::AbstractVector{Vector3D}, occupancies::AbstractVector{<:Real},
		bfactors::AbstractVector{<:Real},
		elements::AbstractVector{<:AbstractString},
		lnoffset::Ref{<:Integer})

	@debug "reading PDB molecular structure"
	titlelines = String[]
	cell = nothing
	atomrecords = Tuple{Int,String}[]
	ln = lnoffset[]
	while !eof(io)
		ln += 1
		line = readline(io)
		if startswith(line, "ATOM") || startswith(line, "HETATM")
			@debug "reading particle record" ln
			push!(atomrecords, (ln, line))
		elseif startswith(line, "TITLE")
			@debug "reading title" ln
			push!(titlelines, strip(line[11:end]))
		elseif startswith(line, "CRYST1")
			@debug "reading cell" ln
			if cell != nothing
				error("multiple cell definitions")
			else
				tokens = split(line[7:end])
				a, b, c = [parse(Float64, t) for t in tokens[1:3]]
				α, β, γ = [deg2rad(parse(Float64, t)) for t in tokens[4:6]]
				cell = TriclinicCell((a,b,c), (α,β,γ))
			end
		elseif startswith(line, "END")
			@debug "stopping" ln
			break
		else
			@debug "ignoring unknown record" ln
		end
	end
	lnoffset[] = ln

	@debug "preparing molecular model"
	n = length(atomrecords)
	n > 0 || error("no particle record found")
	resize!(model, n)
	if length(titlelines) > 0
		model.header.title = join(titlelines, " ")
	end
	if cell != nothing
		model.header.cell = cell
	end

	@debug "parsing particle records" n
	prevline = repeat("\n" , 78)
	for (i, (ln, line)) in enumerate(atomrecords)
		@debug "parsing particle record" i ln
		@debug "parsing id"
		ids[i] = parse(Int, SubString(line, 7, 11))
		@debug "parsing name"
		names[i] = substrip(line, 13, 16)
		@debug "parsing resname"
		if subsetsequal(line, prevline, 18, 21)
			resnames[i] = resnames[i-1]
		else
			resnames[i] = substrip(line, 18, 21)
		end
		@debug "parsing chainid"
		if subsetsequal(line, prevline, 22, 22)
			chainids[i] = chainids[i-1]
		else
			chainids[i] = substrip(line, 22, 22)
		end
		@debug "parsing resid"
		resids[i] = parse(Int, SubString(line, 23, 26))
		@debug "parsing Rx"
		x = parse(Float64, SubString(line, 31, 38))
		@debug "parsing Ry"
		y = parse(Float64, SubString(line, 39, 46))
		@debug "parsing Rz"
		z = parse(Float64, SubString(line, 47, 54))
		R[i] = Vector3D(x,y,z)
		@debug "parsing occupancy"
		occupancies[i] = parse(Float64, SubString(line, 55, 60))
		@debug "parsing B"
		bfactors[i] = parse(Float64, SubString(line, 61, 66))
		@debug "parsing element"
		if subsetsequal(line, prevline, 77, 78)
			elements[i] = elements[i-1]
		else
			elements[i] = substrip(line, 77, 78)
		end
		prevline = line
	end
	unwrapids!(ids, 99999)
	unwrapnames!(names)
	unwrapids!(resids, 9999)

	model
end

function writepdb(io::IO, model::ParticleCollection;
		modelid::Union{Integer,Nothing} = nothing, endrecord::Bool = true,
		lnoffset::Ref{<:Integer} = Ref(0))
	n = length(model)
	writepdb(io, get(model.header, :title, nothing),
			get(model.header, :cell, nothing), get(model, :ids, 1:n),
			model.names, model.resids, model.resnames,
			get(model, :chainids, Repeated("", n)), model.R,
			get(model, :occupancies, Repeated(1.0, n)),
			get(model, :bfactors, Repeated(0.0, n)),
			get(model, :elements, Repeated("", n)), modelid,
			endrecord, lnoffset)
end

function writepdb(io::IO, title::Union{AbstractString,Nothing},
		cell::Union{TriclinicPBC,Nothing}, ids::AbstractVector{<:Integer},
		names::AbstractVector{<:AbstractString},
		resids::AbstractVector{<:Integer},
		resnames::AbstractVector{<:AbstractString},
		chainids::AbstractVector{<:AbstractString}, R::AbstractVector{Vector3D},
		occupancies::AbstractVector{<:Real}, bfactors::AbstractVector{<:Real},
		elements::AbstractVector{<:AbstractString},
		modelid::Union{Integer,Nothing}, endrecord::Bool,
		lnoffset::Ref{<:Integer})
	@debug "writing PDB molecular structure"
	n = length(ids)

	@debug "checking model properties"
	@boundscheck begin
		n > 0 || error("number of particles must be strictly positive")
		length(names) == length(resids) == length(resnames) ==
				length(chainids) == length(R) == length(occupancies) ==
				length(bfactors) == length(elements) ||
				error("size mismatch between property arrays")
		i = findfirst(x -> length(x) > 4, names)
		i == nothing || error("particle $(i): name exceeds 4 chars")
		i = findfirst(x -> length(x) > 4, resnames)
		i == nothing || error("particle $(i): resname exceeds 4 chars")
		i = findfirst(x -> length(x) > 1, chainids)
		i == nothing || error("particle $(i): chainid exceeds 1 char")
		i = findfirst(x -> length(x) > 2, elements)
		i == nothing || error("particle $(i): element exceeds 2 chars")
	end

	ln = lnoffset[]

	if title != nothing
		ln += 1
		@debug "writing title" ln
		for (i, line) in enumerate(wraptext(title, 70))
			if i == 1
				@printf(io, "TITLE     %s\n", line)
			else
				ln += 1
				@debug "continuing title" ln
				@printf(io, "TITLE %4d %s\n", i, line)
			end
		end
	end

	if cell != nothing
		ln += 1
		@debug "writing cell" ln
		edges, angles = pbcgeometry(cell)
		@printf(io, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n",
				edges..., rad2deg.(angles)...)
	end

	if modelid != nothing
		ln += 1
		@debug "start of model" ln
		@printf(io, "MODEL     %4d\n", modelid)
	end

	@debug "writing particle records" n
	for i = 1:n
		ln += 1
		@debug "writing particle record" i ln
		id = wrapid(ids[i], 99999)
		name = pdb_particle_name(names[i])
		resid = wrapid(resids[i], 9999)
		x, y, z = R[i]
		element = (elements[i] == "H" ? "" : elements[i])
		@printf(io, "ATOM  %5d %4s %-4s%1s%4d    %8.3f%8.3f%8.3f%6.2f%6.2f          %2s  \n",
				id, name, resnames[i], chainids[i], resid, x, y, z,
				occupancies[i], bfactors[i], element)
	end

	if modelid != nothing
		ln += 1
		@debug "end of model" ln
		write(io, "ENDMDL\n")
	end

	if endrecord
		ln += 1
		@debug "end of PDB" ln
		write(io, "END\n")
	end

	lnoffset[] = ln
end

end # module


text:
const _emptyismissing = Ref(false)

"""
    Mill.emptyismissing()

Get the current value of the `emptyismissing` parameter.

See also: [`Mill.emptyismissing!`](@ref).
"""
emptyismissing() = _emptyismissing[]

"""
    Mill.emptyismissing!(::Bool)

Set the new value to the `emptyismissing` parameter.

See also: [`Mill.emptyismissing`](@ref).
"""
emptyismissing!(a) = _emptyismissing[] = a

const _string_start_code = Ref(0x02) # STX in ascii

"""
    Mill.string_start_code()

Get the current value of the `string_start_code` parameter used as a code point of the abstract string-start character.
The default value of the parameter is `0x02`, which corresponds to the `STX` character in ASCII encoding.

See also: [`Mill.string_start_code!`](@ref), [`Mill.string_end_code`](@ref), [`Mill.string_end_code!`](@ref).
"""
string_start_code() = _string_start_code[]

"""
    Mill.string_start_code!(c::UInt8)

Set the new value to the `string_start_code` parameter used as a code point of the abstract string-start character to `c`.
The default value of the parameter is `0x02`, which corresponds to the `STX` character in ASCII encoding.

See also: [`Mill.string_start_code`](@ref), [`Mill.string_end_code`](@ref), [`Mill.string_end_code!`](@ref).
"""
string_start_code!(c) = _string_start_code[] = c

const _string_end_code = Ref(0x03) # ETX in ascii

"""
    Mill.string_end_code()

Get the current value of the `string_end_code` parameter used as a code point of the abstract string-end character.
The default value of the parameter is `0x03`, which corresponds to the `ETX` character in ASCII encoding.

See also: [`Mill.string_end_code!`](@ref), [`Mill.string_start_code`](@ref), [`Mill.string_start_code!`](@ref).
"""
string_end_code() = _string_end_code[]

"""
    Mill.string_end_code!(c::UInt8)

Set the new value to the `string_end_code` parameter used as a code point of the abstract string-end character to `c`.
The default value of the parameter is `0x03`, which corresponds to the `ETX` character in ASCII encoding.

See also: [`Mill.string_end_code`](@ref), [`Mill.string_start_code`](@ref), [`Mill.string_start_code!`](@ref).
"""
string_end_code!(c) = _string_end_code[] = c

const _wildcard_code = Ref(0x00) # NUL in ascii

wildcard_code() = _wildcard_code[]
wildcard_code!(c) = _wildcard_code[] = c

ext: jl


text:
# ```@setup covmat
# # the example command line input for this script
# ARGS = ["example.toml", "143", "143", "143", "143", "T", "T", "T", "T", "1", "2", "1", "2", "--plot"]
# ```
#
# # [Covariance Matrix Estimation (covmat.jl)](@id covmat)
# This part of the pipeline dominates the total computational cost. Many of the steps like
# spectra estimation could have been run on a laptop, but estimating the covariance for the
# full Planck data vector involves computing ``\sim 10^3`` covariance matrix blocks. The
# covariance matrix involves four fields, and those fields each have a frequency, channel,
# and map split. We thus
# need every combination of frequencies ``f_1``, ``f_2``, ``f_3``, ``f_4`` for
# ``f_i \in \{ 100, 143, 217 \}`` with channels ``A``, ``B``, ``C``, ``D``
# for ``X_i \in \{T, E\}``, and splits
# ``s_1``, ``s_2``, ``s_3``, ``s_4`` for ``s_i \in \{ \mathrm{hm1}, \mathrm{hm2} \}``.
#
# These are the inputs for the command line interface for this script, `covmat.jl`.
#
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia covmat.jl global.toml f1 f2 f3 f3 X1 X2 X3 X4 s1 s2 s3 s4 --plot</code></pre>
# ```
# As usual, the `--plot` argument is optional.
#
#

# ## Setup
# As usual, we need to load up our dependencies and then unpack
# the command line arguments into variables. We'll keep the
# `freqs`, `specs`, and `splits` as just arrays of strings for now.

freqs = [ARGS[2], ARGS[3], ARGS[4], ARGS[5]]
specs = [ARGS[6], ARGS[7], ARGS[8], ARGS[9]]
splits = [ARGS[10], ARGS[11], ARGS[12], ARGS[13]]

using TOML
using PowerSpectra
using Healpix
using JLD2
using Plots
include("util.jl")

# As usual, we get some general information from the configuration file
# specified as the first command line argument.

configfile = ARGS[1]
config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
lmax = nside2lmax(nside)
lmax_planck = min(2508, lmax)
run_name = config["general"]["name"]

# ## Change to canonical ordering
# The covariance matrix is over fields ``F_1``, ``F_2``, ``F_3``, ``F_4`` and represents
# the covariance between the ``\langle F_1, F_2 \rangle`` and
# ``\langle F_3 \times F_4 \rangle`` cross-spectra.
# However, the covariance matrix is invariant when swapping ``F_1 \leftrightarrow F_2``
# as well as ``F_3 \leftrightarrow F_4``. Under the transformation that swaps the
# first pair for the second pair (i.e. ``F_1 \leftrightarrow F_3`` and
# ``F_2 \leftrightarrow F_4`` simultaneously), the covariance matrix is transposed.
#
# Given these convenient properties, we define a canonical ordering of the fields
# in a covariance matrix. We insist that no temperature field can ever be after a
# polarization field in an ordering. Thus, we don't need to implement an expression
# for a ``TETT`` covariance, because it can always be transformed into the canonical
# ``TTTE`` (this example requires a transposition pf the matrix at the end).
# We'll write a few utility functions to enable these transformations here.
#
## This utility tests if it's not a temperature channel
ispol(spec_str) = (spec_str == "T") ? 0 : 1

function swap!(specs, freqs, splits, i, j)
    specs[i], specs[j] = specs[j], specs[i]
    freqs[i], freqs[j] = freqs[j], freqs[i]
    splits[i], splits[j] = splits[j], splits[i]
end

## canonical form: T as far left as possible

"""
    canonical!(specs, freqs, splits)

Given four fields, described by length-4 arrays of specs (T, E), freqs, and splits,
modify these input arrays so they are in the canonical ordering, i.e. with
temperature fields before polarization.

### Arguments:
- `specs`: array of four strings, containing "T" and "E" elements
- `freqs`: array of four frequencies
- `splits`: array of four choices of data split

### Returns:
- `Boolean`: some canonical orderings are a matrix transpose of the inputs. If this
    function returns true, then one needs to perform this transposition.
"""
function canonical!(specs, freqs, splits)
    transposed = false

    if sum(ispol.(specs[1:2])) > sum(ispol.(specs[3:4]))
        swap!(specs, freqs, splits, 1, 3)
        swap!(specs, freqs, splits, 2, 4)
        transposed = true
    end
    if specs[1] == "E" && specs[2] == "T"
        swap!(specs, freqs, splits, 1, 2)
    end
    if specs[3] == "E" && specs[4] == "T"
        swap!(specs, freqs, splits, 3, 4)
    end

    return transposed
end

# We re-arrange the `specs`, `freqs`, and `splits` spectra so they are in the
# canonical ordering. We save whether we need to transpose the result.

transposed = canonical!(specs, freqs, splits)
@show specs freqs splits transposed

## convenience names since we looking for covariance between AB and CD spectra
specAB = Symbol(specs[1]*specs[2])
specCD = Symbol(specs[3]*specs[4])
@show specAB specCD

# Next, we need to create dictionaries that contain the fruits of our labors from
# earlier parts of the pipeline. These include the theory spectra, the ratio of the
# noise power spectrum to the white noise power spectrum.

## store info needed for covariance calculations
spectra = Dict{SpectrumName, SpectralVector{Float64, Vector{Float64}}}()
noiseratios = Dict{SpectrumName, SpectralVector{Float64, Vector{Float64}}}()
identity_spectrum = SpectralVector(ones(lmax+1));

# The point source holes in the Planck likelihood maps are not sufficiently apodized,
# which causes problems in the covariance matrix estimation. Recall that we ran
# signal-only simulations in `signalsim.jl` so that we could correct for them here.
# We read in the correction curve estimated in `corrections.jl` for those simuulations,
# for use in the covariance matrix.

function extend_signal(s::Vector{T}, nside) where T
    lmax = nside2lmax(nside)
    es = SpectralVector(zeros(T, lmax+1))
    for ℓ ∈ 0:min(length(s) - 1, lmax)
        es[ℓ] = s[ℓ+1]
    end
    return es
end

function get_correction(freq1, freq2, spec)
    spec_str = string(spec)
    if parse(Float64, freq1) > parse(Float64, freq2)
        freq1, freq2 = freq2, freq1
        spec_str = reverse(spec_str)
    end
    correction_path = joinpath(config["scratch"], "point_source_corrections")
    correction_file = joinpath(correction_path, "$(freq1)_$(freq2)_$(spec_str)_corr.dat")
    correction_gp = readdlm(correction_file)[:,1]
    correction_gp[1:3] .= 0.0
    return extend_signal(correction_gp, nside)
end

# We want to supply the ratios of the noise power spectrum and the white noise power spectrum.
# To do that, we need to read in the estimates of the white noise power spectrum that
# were obtained from the maps.
function whitenoiselevel(config, freq::String, split::String)
    whitenoisefile = joinpath(config["scratch"], "whitenoise.dat")
    df = CSV.read(whitenoisefile, DataFrame)
    for row in 1:nrow(df)
        if (string(df[row,:freq]) == freq) && (string(df[row,:split]) == split)
            return df[row, :noiseT], df[row, :noiseP]
        end
    end
    throw(ArgumentError("freq and split combo missing from white noise file"))
end

# We save the noise power spectrum models in terms of parameters of the smooth
# model fit. We use the expression from CamSpec.
@. camspec_model(ℓ, α) =  α[1] * (100. / ℓ)^α[2] + α[3] * (ℓ / 1000.)^α[4] / ( 1 + (ℓ / α[5])^α[6] )^α[7]

# Now to put it all together. We need to put all of this information into our
# various dictionaries.

ells = collect(0:(lmax))
coefficientpath = joinpath(config["scratch"], "noise_model_coeffs")
beampath = joinpath(config["scratch"], "beams")

## loop over combinations and put stuff into the dictionaries
for (f1, s1) in zip(freqs, splits)
    for (f2, s2) in zip(freqs, splits)
        f1_name = "P$(f1)hm$(s1)"
        f2_name = "P$(f2)hm$(s2)"
        signal, theory_dict = signal_and_theory(f1, f2, config)

        WlTT = util_planck_beam_Wl(f1, "hm"*s1, f2, "hm"*s2, :TT, :TT;
            lmax=lmax, beamdir=beampath)
        WlTE = util_planck_beam_Wl(f1, "hm"*s1, f2, "hm"*s2, :TE, :TE;
            lmax=lmax, beamdir=beampath)
        WlEE = util_planck_beam_Wl(f1, "hm"*s1, f2, "hm"*s2, :EE, :EE;
            lmax=lmax, beamdir=beampath)

        spectra[(:TT, f1_name, f2_name)] = extend_signal(signal["TT"], nside) .* WlTT .*
            sqrt.(1 .+ get_correction(f1, f2, "TT"))
        spectra[(:TE, f1_name, f2_name)] = extend_signal(signal["TE"], nside) .* WlTE .*
            sqrt.(1 .+ get_correction(f1, f2, "TE"))
        spectra[(:EE, f1_name, f2_name)] = extend_signal(signal["EE"], nside) .* WlEE .*
            sqrt.(1 .+ get_correction(f1, f2, "EE"))

        if f1_name == f2_name
            whitenoiseT, whitenoiseP = whitenoiselevel(config, f1, s1)
            coeffTT = readdlm(joinpath(coefficientpath, "$(run_name)_$(f1)_TT_hm$(s1).dat"))[:,1]
            coeffEE = readdlm(joinpath(coefficientpath, "$(run_name)_$(f1)_EE_hm$(s1).dat"))[:,1]
            nlTT = camspec_model(ells, coeffTT)
            nlEE = camspec_model(ells, coeffEE)
            ratioTT = SpectralVector(nlTT ./ whitenoiseT)
            ratioEE = SpectralVector(nlEE ./ whitenoiseP)
            ratioTT[0:1] .= 1.0
            ratioEE[0:1] .= 1.0
            noiseratios[(:TT, f1_name, f2_name)] = sqrt.(ratioTT)
            noiseratios[(:EE, f1_name, f2_name)] = sqrt.(ratioEE)
        else
            noiseratios[(:TT, f1_name, f2_name)] = identity_spectrum
            noiseratios[(:EE, f1_name, f2_name)] = identity_spectrum
        end
    end
end

fnames = ["P$(freq)hm$(split)" for (freq, split) in zip(freqs, splits)]

if "--plot" ∈ ARGS
    plot(noiseratios[(:TT, fnames[1], fnames[1])].parent, label="$(freqs[1]) hm1 TT")
    plot!(noiseratios[(:EE, fnames[1], fnames[1])].parent, label="$(freqs[1]) hm1 EE",
        ylim=(0.0,2), xlim=(0,200), size=(600,300))
end

# Each of the four fields involved with the covariance matrix involves a polarized map,
# with associated mask and variance map. These are packed toegether in the `CovField` struct
# in `PowerSpectra.jl`.

## pack all the information that the covmat needs into a structure
function loadcovfield(freq, split)
    ## read maskT, maskP, covariances
    mapid = "P$(freq)hm$(split)"
    maskfileT = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid)_maskT.fits")
    maskfileP = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid)_maskP.fits")
    mapfile = joinpath(config["scratch"], "maps", config["map"][mapid])
    maskT = readMapFromFITS(maskfileT, 1, Float64)
    maskP = readMapFromFITS(maskfileP, 1, Float64)
    covII = nest2ring(readMapFromFITS(mapfile, 5, Float64)) * 1e12
    covQQ = nest2ring(readMapFromFITS(mapfile, 8, Float64)) * 1e12
    covUU = nest2ring(readMapFromFITS(mapfile, 10, Float64)) * 1e12
    return CovField(mapid, maskT, maskP, PolarizedHealpixMap(covII, covQQ, covUU))
end

fields = [loadcovfield(freqs[i], splits[i]) for i in 1:4];

# Next, we need to compute the mode-coupling matrices. These are multiplied with
# the covariance matrices in order to correct for the effects of the mask, just
# like in the case of the spectra.
## compute mode-coupling matrices
mcm_type_AB = (specAB == :EE) ? :M⁺⁺ : specAB
mcm_type_CD = (specCD == :EE) ? :M⁺⁺ : specCD

@time M₁₂ = mcm(mcm_type_AB, fields[1], fields[2])
@time M₃₄ = mcm(mcm_type_CD, fields[3], fields[4]);

## compute mode-coupled covariance matrix
workspace = CovarianceWorkspace(fields...);
@time C = coupledcov(Symbol(specs[1]*specs[2]), Symbol(specs[3]*specs[4]),
    workspace, spectra, noiseratios)

## this is the product of A and B pixel windows
function spec2pixwin(XY::Symbol, nside)
    pixwinT, pixwinP = pixwin(nside; pol=true)
    pixwinP[1:2] .= 1.0
    pX = (first(string(XY)) == 'T') ? pixwinT : pixwinP
    pY =  (last(string(XY)) == 'T') ? pixwinT : pixwinP
    return SpectralVector(pX .* pY)
end

pixwinAB = spec2pixwin(specAB, nside)
pixwinCD = spec2pixwin(specCD, nside)

# Next, we need decouple the modes from the effect of the mask, correct for the
# instrumental beam and pixel window function, and then finally bin the covariance matrix.
#

## load binnings
binfile = joinpath(@__DIR__, "../", "input", "binused.dat")
P, lb = util_planck_binning(binfile; lmax=lmax)
beampath = joinpath(config["scratch"], "beams")

C₀ = decouple_covmat(C, M₁₂, M₃₄)
WlAB = util_planck_beam_Wl(freqs[1], "hm"*splits[1], freqs[2], "hm"*splits[2], specAB, specAB;
    lmax=lmax, beamdir=beampath)
WlCD = util_planck_beam_Wl(freqs[3], "hm"*splits[3], freqs[4], "hm"*splits[4], specCD, specCD;
    lmax=lmax, beamdir=beampath)

for l1 in 0:lmax
    for l2 in 0:lmax
        C₀[l1, l2] /= WlAB[l1] * WlCD[l2] * pixwinAB[l1] * pixwinCD[l2]
    end
end

Cbb = P * parent(C₀) * (P')

## transpose if needed, from turning into the canonical form
if transposed
    Cbb = Array(transpose(Cbb))
end

spectrapath = joinpath(config["scratch"], "rawspectra")
covpath = joinpath(config["scratch"], "covmatblocks")
mkpath(covpath)

covmat_filename = joinpath(covpath, join(specs) * "_" *
    join(freqs) * "_" * join(splits) * ".jld2")
@save covmat_filename Cbb

max_forks_repo_name: sgiardie/PSpipe


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiSchedulingV1PriorityClass <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    globalDefault::Any # spec type: Union{ Nothing, Bool } # spec name: globalDefault
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    preemptionPolicy::Any # spec type: Union{ Nothing, String } # spec name: preemptionPolicy
    value::Any # spec type: Union{ Nothing, Int32 } # spec name: value

    function IoK8sApiSchedulingV1PriorityClass(;apiVersion=nothing, description=nothing, globalDefault=nothing, kind=nothing, metadata=nothing, preemptionPolicy=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("globalDefault"), globalDefault)
        setfield!(o, Symbol("globalDefault"), globalDefault)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("preemptionPolicy"), preemptionPolicy)
        setfield!(o, Symbol("preemptionPolicy"), preemptionPolicy)
        validate_property(IoK8sApiSchedulingV1PriorityClass, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiSchedulingV1PriorityClass

const _property_map_IoK8sApiSchedulingV1PriorityClass = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("description")=>Symbol("description"), Symbol("globalDefault")=>Symbol("globalDefault"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("preemptionPolicy")=>Symbol("preemptionPolicy"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiSchedulingV1PriorityClass = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("description")=>"String", Symbol("globalDefault")=>"Bool", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("preemptionPolicy")=>"String", Symbol("value")=>"Int32")
Base.propertynames(::Type{ IoK8sApiSchedulingV1PriorityClass }) = collect(keys(_property_map_IoK8sApiSchedulingV1PriorityClass))
Swagger.property_type(::Type{ IoK8sApiSchedulingV1PriorityClass }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiSchedulingV1PriorityClass[name]))}
Swagger.field_name(::Type{ IoK8sApiSchedulingV1PriorityClass }, property_name::Symbol) =  _property_map_IoK8sApiSchedulingV1PriorityClass[property_name]

function check_required(o::IoK8sApiSchedulingV1PriorityClass)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiSchedulingV1PriorityClass }, name::Symbol, val)
end


text:
# ---
# title: 398. Random Pick Index
# id: problem398
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Reservoir Sampling
# link: <https://leetcode.com/problems/random-pick-index/description/>
# hidden: true
# ---
#
# Given an array of integers with possible duplicates, randomly output the index
# of a given target number. You can assume that the given target number must
# exist in the array.
#
# **Note:**
# The array size can be very large. Solution that uses too much extra space will
# not pass the judge.
#
# **Example:**
#
#
#
#     int[] nums = new int[] {1,2,3,3,3};
#     Solution solution = new Solution(nums);
#
#     // pick(3) should return either index 2, 3, or 4 randomly. Each index should have equal probability of returning.
#     solution.pick(3);
#
#     // pick(1) should return 0. Since in the array only nums[0] is equal to 1.
#     solution.pick(1);
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "wind.jl" begin

        @test Vground(110.0, 40.0, 315.0, 245.0) ≈ 89.7 atol = 0.1
        @test Vground(120.0, 40.0, 315.0, 245.0) ≈ 100.3 atol = 0.1

        @test head_wind(20.0, 60.0, 30.0) ≈ 17.32 atol = 0.01
        @test cross_wind(20.0, 60.0, 30.0) ≈ 10.00 atol = 0.01

        @test isnan(Vground(110.0, 120.0, 0.0, 90.0))
end


text:
module IO

	using JLD2
	using Flux
	using TinyML

	"""
		load(::String)

	Loads a jld2 model from the informed current path
	"""
	load(filename::String) = JLD2.load(filename)

	"""
		save(model::Chain, training_set::Genetic.TrainingSet)

	Saves a jld2 model to the current path
	"""
	function save(filename::String, model::Chain, trainingset::Genetic.TrainingSet)
		tosave = Dict("model" => model,
					"trainingset" => trainingset)
		JLD2.save(filename, tosave)
	end

	export load, save

end


text:
function simSimpleWithDifferentAgents()
    srand(1234)
    for acc in 0:10:100
        fixed = mean([simSimple(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        div10 = mean([simSimple(1000, setupRandomAgents(100, acc, 10), [benchmarkRegistryMean]) for i in 1:10])[1]
        div20 = mean([simSimple(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:10])[1]
        mixed = mean([simSimple(1000, setupMixedAgents(90, acc, 10, 100-acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        println("$acc, $fixed, $div10, $div20, $mixed")
    end
end

function simChallengeWithDifferentAgents()
    srand(1234)
    for acc in 0:10:100
        fixed = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        div10 = mean([simChallenge(1000, setupRandomAgents(100, acc, 10), [benchmarkRegistryMean]) for i in 1:10])[1]
        div20 = mean([simChallenge(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:10])[1]
        mixed = mean([simChallenge(1000, setupMixedAgents(100-acc, 0, acc, 100), [benchmarkRegistryMean]) for i in 1:10])[1]

        println("$acc, $fixed, $div10, $div20, $mixed")
    end
end

function simTokenWithDifferentAgents()
    srand(1234)
    for acc in 0:10:100
        fixed = mean([simToken(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        div10 = mean([simToken(1000, setupRandomAgents(100, acc, 10), [benchmarkRegistryMean]) for i in 1:10])[1]
        div20 = mean([simToken(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:10])[1]
        mixed = mean([simToken(1000, setupMixedAgents(100-acc, 10, acc, 90), [benchmarkRegistryMean]) for i in 1:10])[1]

        println("$acc, $fixed, $div10, $div20, $mixed")
    end
end

function compareRegistryQualityOnDifferentSetups()
    for acc in 10:5:100
        noChallenge = mean([simSimple(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        challenge = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:10])[1]
        divAgents = mean([simChallenge(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:10])[1]
        token = mean([simToken(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:10])[1]

        println("$acc, $noChallenge, $challenge, $divAgents, $token")
    end
end

"""
Accuracy is fixed to 50, Number of steps to 1000
"""
function manipulatingNoAgents()
    for noAgents in 10:10:200
        noToken = [simChallenge(1000, setupRandomAgents(noAgents, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
        token = [simToken(1000, setupRandomAgents(noAgents, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
        boost = mean(token) - mean(noToken)
        test = pvalue(UnequalVarianceTTest(token, noToken))
        significant = test < 0.05
        println("$noAgents, $boost, $test, $significant")
    end
end

"""
Accuracy is fixed to 50, Number of agents to 100
"""
function manipulatingNoSteps()
    for noSteps in 100:100:2000
        noToken = [simChallenge(noSteps, setupRandomAgents(100, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
        token = [simToken(noSteps, setupRandomAgents(100, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
        boost = mean(token) - mean(noToken)
        test = pvalue(UnequalVarianceTTest(token, noToken))
        significant = test < 0.05
        println("$noSteps, $boost, $test, $significant")
    end
end

function efficiencyArea()
    for noAgents in 10:10:200
        for noSteps in 100:100:2000
            noToken = [simChallenge(noSteps, setupRandomAgents(noAgents, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
            token = [simToken(noSteps, setupRandomAgents(noAgents, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
            boost = mean(token) - mean(noToken)
            println("$noSteps, $noAgents, $boost")
        end
    end
end

function accuracyAndBalanceCorrelation()
    for acc in 0:5:100
        results = mean([simToken(1500, setupRandomAgents(100, acc, 20), [benchmarkWeightedAccuracyBoost, benchmarkAccuracyBoost, benchmarkAccuracyBalanceCorrelation]) for i in 1:10])
        weighted = results[1]
        binary = results[2]
        correlation = results[3]
        println("$acc, $weighted, $binary, $correlation")
    end
end

function curationProcess()
    for noSteps in 50:50:2000
        results = mean([simToken(noSteps, setupMixedAgents(80, 0, 20, 100), [benchmarkWeightedAccuracyBoost, benchmarkAccuracyBalanceCorrelation, benchmarkRegistryMean]) for i in 1:30])
        boost = results[1]
        correlation = results[2]
        quality = results[3]
        println("$noSteps, $correlation, $quality")
    end
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
using AWSClusterManagers: AWSBatchManager
using Distributed

const NUM_WORKERS = 4

# Must start workers before any @everywhere for them to be available on the workers
addprocs(AWSBatchManager(NUM_WORKERS))

@everywhere function throw_darts()
    MAX_THROW = 250_000
    dart_count = rand(1:MAX_THROW)
    hits = 0

    for _ in 1:dart_count
        x = rand()
        y = rand()

        if ((x^2 + y^2) < 1)
            hits += 1
        end
    end

    return hits, dart_count
end

calc_approx_pi(hits::Integer, thrown::Integer) = 4 * hits / thrown

function main()
    total_hits = 0
    total_throws = 0

    futures = [remotecall(throw_darts, worker) for worker in workers()]

    for (worker, future) in enumerate(futures)
        hits, throws = fetch(future)
        percentile = (hits / throws)

        println("Worker $(worker) landed $(hits) of $(throws), ($(percentile)%)!")
        total_hits += hits
        total_throws += throws
    end

    approx_pi = calc_approx_pi(total_hits, total_throws)
    println("$(NUM_WORKERS) workers threw a total of $(total_throws) darts.")
    println("They calculated PI to be approximately,\n$(approx_pi)")
end


text:
module twister.parse
using pest
using pest.parse

def parse
    println ""
end

export parse
export twister.parse

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
## Presence absence distribution functions

"""
    presence_absence(species::AbstractDataFrame, layer::SimpleSDMLayer)

Convert the occurrences the `species` DataFrame to a presence-absence layer
based on the dimensions and coordinates from `layer`.
"""
function presence_absence(species::AbstractDataFrame, layer::SimpleSDMLayer)
    @assert all(in(names(species)), ["longitude", "latitude"]) "Coordinates must be stored in columns called longitude and latitude"

    # Create empty grid for presence-absence data (with NaN)
    distribution_grid = copy(layer.grid)
    replace!(x -> !isnothing(x) ? 0.0 : NaN, distribution_grid)

    # Get unique sites/coordinates
    unique_sites = unique(species, [:longitude, :latitude])

    # Loop for all sites
    lon_all = collect(longitudes(layer))
    lat_all = collect(latitudes(layer))
    for (lon, lat) in zip(unique_sites.longitude, unique_sites.latitude)
        # Get grid position for each site
        i_lon = findmin(abs.(lon .- lon_all))[2]
        j_lat = findmin(abs.(lat .- lat_all))[2]
        # Add 1 per species presence
        distribution_grid[j_lat, i_lon] += 1.0
    end

    # Reduce sites with more than 1 presence to binary value (default)
    replace!(x -> x > 1.0 ? 1.0 : x, distribution_grid)
    # Replace zeros (absences) by nothing
    replace!(distribution_grid, 0.0 => nothing)
    replace!(distribution_grid, NaN => nothing)
    # Create SimpleSDMLayer
    distribution_layer = SimpleSDMResponse(distribution_grid, layer)
    return distribution_layer
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"""
Filter out cells which are not satisfy specific criteria.
"""
function filter_cells!(p::Profile; min_counts::Real=0, max_counts::Real=maximum(sum(p.data,dims=1)),
                       min_genes::Real=0, max_genes::Real=nrow(p))
      check_cell_args(min_counts, max_counts, min_genes, max_genes)
      return _filter_cells!(p, min_counts, max_counts, min_genes, max_genes)
end

function filter_cells(p::Profile; min_counts::Real=0, max_counts::Real=maximum(sum(p.data,dims=1)),
                      min_genes::Real=0, max_genes::Real=nrow(p))
      check_cell_args(min_counts, max_counts, min_genes, max_genes)
      return _filter_cells!(copy(p), min_counts, max_counts, min_genes, max_genes)
end

function _filter_cells!(p::Profile, min_counts::Real, max_counts::Real, min_genes::Real, max_genes::Real)
      expressed_genes = sum(p.data .!= 0, dims=1)
      bool_idx = min_genes .<= expressed_genes .<= max_genes
      library_size = sum(p.data, dims=1)
      bool_idx .&= min_counts .<= library_size .<= max_counts
      filter_layers!(p, obs_idx=vec(bool_idx))
      idx = findall(vec(bool_idx))
      p.data = p.data[:, idx]
      p.obs = p.obs[idx, :]
      return p
end

function check_cell_args(min_counts::Real, max_counts::Real, min_genes::Real, max_genes::Real)
      (min_counts < 0) && throw(ArgumentError("`min_counts` must not lower than zero but got $min_counts."))
      (min_genes < 0) && throw(ArgumentError("`min_genes` must not lower than zero but got $min_genes."))
      (max_genes < min_genes) &&
            throw(ArgumentError("`max_genes` must not lower than minimum number of genes ($min_genes) but got $max_genes."))
      (max_counts < min_counts) &&
            throw(ArgumentError("`max_counts` must not lower than minimum count ($min_counts) but got $max_counts."))
end

"""
Filter out genes which are not satisfy specific criteria.
"""
function filter_genes!(p::Profile; min_counts::Real=0, max_counts::Real=maximum(sum(p.data,dims=2)),
                       min_cells::Real=0, max_cells::Real=ncol(p))
      check_genes_args(min_counts, max_counts, min_cells, max_cells)
      return _filter_genes!(p, min_counts, max_counts, min_cells, max_cells)
end

function filter_genes(p::Profile; min_counts::Real=0, max_counts::Real=maximum(sum(p.data,dims=2)),
                      min_cells::Real=0, max_cells::Real=ncol(p))
      check_genes_args(min_counts, max_counts, min_cells, max_cells)
      return _filter_genes!(copy(p), min_counts, max_counts, min_cells, max_cells)
end

function _filter_genes!(p::Profile, min_counts::Real, max_counts::Real, min_cells::Real, max_cells::Real)
      expressed_cells = sum(p.data .!= 0, dims=2)
      bool_idx = min_cells .<= expressed_cells .<= max_cells
      expressed_counts = sum(p.data, dims=2)
      bool_idx .&= min_counts .<= expressed_counts .<= max_counts
      filter_layers!(p, var_idx=vec(bool_idx))
      idx = findall(vec(bool_idx))
      p.data = p.data[idx, :]
      p.var = p.var[idx, :]
      return p
end

function check_genes_args(min_counts::Real, max_counts::Real, min_cells::Real, max_cells::Real)
      (min_counts < 0) && throw(ArgumentError("`min_counts` must not lower than zero but got $min_counts."))
      (min_cells < 0) && throw(ArgumentError("`min_cells` must not lower than zero but got $min_cells."))
      (max_cells < min_cells) &&
            throw(ArgumentError("`max_cells` must not lower than minimum number of cells ($min_cells) but got $max_cells."))
      (max_counts < min_counts) &&
            throw(ArgumentError("`max_counts` must not lower than minimum count ($min_counts) but got $max_counts."))
end


text:
module gtpmip

using PDDL
using JuMP
using Gurobi
using SymbolicPlanners
using LinearAlgebra
using mutex
using HybridFunnelGraphs

include("types.jl")
include("solvers/hpd_solver.jl")
include("solvers/legacies/pure_symbolic_gtpmip.jl")

export solve_hpd,
       pure_symbolic_solve

end


text:
using Convex: AbstractExpr, ConicObj
using LinearAlgebra
using SparseArrays

# It's not super easy to capture the output
# I ended up using this pattern from Suppressor:
# https://github.com/JuliaIO/Suppressor.jl/blob/b4ff08f0fe795a2ce9e592734a758c9e6d8e2bc4/src/Suppressor.jl#L124-L152
function solve_and_return_output(problem, solver; kwargs...)
        original_stdout = stdout
        rd, wr = redirect_stdout()
        out_task = @async read(rd, String)
        try
            solve!(problem, solver; kwargs...)
        finally
            Base.Libc.flush_cstdio() #  https://github.com/JuliaLang/julia/issues/31236
            redirect_stdout(original_stdout)
            close(wr)
        end
        return fetch(out_task)
end

@testset "Utilities" begin

    @testset "`solve!` does not return anything" begin
        x = Variable()
        p = satisfy(x >= 0)
        output = solve!(p, () -> SCS.Optimizer(verbose=0, eps=1e-6))
        @test output === nothing
    end

    @testset "`silent_solver` works" begin
        x = Variable()
        p = satisfy(x >= 0)
        output_non_silent = solve_and_return_output(p, () -> SCS.Optimizer(eps=1e-6))
        @test output_non_silent != ""
        output_silent = solve_and_return_output(p, () -> SCS.Optimizer(eps=1e-6), silent_solver=true)
        @test output_silent == ""
    end

    # This might get deprecated later.
    @testset "`solve!` can take an optimizer directly" begin
        x = Variable()
        p = satisfy(x >= 0)
        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
        @test output === nothing
    end

    @testset "`solve!` with MPB solver errors" begin
        x = Variable()
        p = satisfy(x >= 0)
        @test_throws ErrorException solve!(p, SCSSolver())
    end

    @testset "Complex objective function errors" begin
        x = Variable()
        @test_throws ErrorException minimize(x + im*x)
    end

    @testset "`optval` is nothing before `solve!`" begin
        x = Variable()
        p = minimize(x, x >= 0)
        @test p.optval === nothing
        solve!(p, () -> SCS.Optimizer(verbose=false))
        @test p.optval ≈ 0.0 atol=1e-3
        @test Convex.termination_status(p) == MOI.OPTIMAL
        @test Convex.objective_value(p) ≈ 0.0 atol=1e-3
    end

    @testset "Default problem type is `Float64`" begin
        x = Variable()
        p = minimize(x, x >= 0)
        @test p isa Convex.Problem{Float64}
    end

    @testset "`set_value!` doesn't convert to `Float64`" begin
        x = Variable()
        set_value!(x, big"1.0")
        @test evaluate(x) isa BigFloat

        x = Variable(2)
        set_value!(x, big.([1.0, 2.0]))
        @test evaluate(x) isa Vector{BigFloat}

        x = Variable(2, 2)
        set_value!(x, big.([1.0 2.0; 3.0 4.0]))
        @test evaluate(x) isa Matrix{BigFloat}
    end

    @testset "Show" begin
        x = Variable()
        @test sprint(show, x) == """
        Variable
        size: (1, 1)
        sign: real
        vexity: affine
        $(Convex.show_id(x))"""
        fix!(x, 1.0)
        @test sprint(show, x) == """
        Variable
        size: (1, 1)
        sign: real
        vexity: constant
        $(Convex.show_id(x))
        value: 1.0"""

        @test sprint(show, 2*x) == """
        * (constant; real)
        ├─ 2
        └─ real variable (fixed) ($(Convex.show_id(x)))"""

        free!(x)
        p = maximize( log(x), x >= 1, x <= 3 )

        @test sprint(show, p) == """
        maximize
        └─ log (concave; real)
           └─ real variable ($(Convex.show_id(x)))
        subject to
        ├─ >= constraint (affine)
        │  ├─ real variable ($(Convex.show_id(x)))
        │  └─ 1
        └─ <= constraint (affine)
           ├─ real variable ($(Convex.show_id(x)))
           └─ 3

        status: `solve!` not called yet"""

        x = ComplexVariable(2,3)
        @test sprint(show, x) == """
        Variable
        size: (2, 3)
        sign: complex
        vexity: affine
        $(Convex.show_id(x))"""

        # test `MAXDEPTH`
        # We construct a binary tree of depth >= 3
        # to make sure it gets truncated appropriately.
        x = Variable(2)
        y = Variable(2)
        level3 = hcat(x,y)
        level2 = hcat(level3, level3)
        root = hcat(level2, level2)
        p = minimize(sum(x), root == root)
        @test sprint(show, p) == """
        minimize
        └─ sum (affine; real)
           └─ 2-element real variable ($(Convex.show_id(x)))
        subject to
        └─ == constraint (affine)
           ├─ hcat (affine; real)
           │  ├─ hcat (affine; real)
           │  │  ├─ …
           │  │  └─ …
           │  └─ hcat (affine; real)
           │     ├─ …
           │     └─ …
           └─ hcat (affine; real)
              ├─ hcat (affine; real)
              │  ├─ …
              │  └─ …
              └─ hcat (affine; real)
                 ├─ …
                 └─ …

        status: `solve!` not called yet"""

        # test `MAXWIDTH`
        x = Variable()
        p = satisfy([ x == i for i = 1:100])
        old_maxwidth = Convex.MAXWIDTH[]
        Convex.MAXWIDTH[] = 2
        @test sprint(show, p) == """
            minimize
            └─ 0
            subject to
            ├─ == constraint (affine)
            │  ├─ real variable ($(Convex.show_id(x)))
            │  └─ 1
            ├─ == constraint (affine)
            │  ├─ real variable ($(Convex.show_id(x)))
            │  └─ 2
            ⋮

            status: `solve!` not called yet"""
        Convex.MAXWIDTH[] = old_maxwidth

        # solved problem
        x = Variable()
        p = satisfy(x >= 0)
        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
        @test sprint(show, p) == """
                minimize
                └─ 0
                subject to
                └─ >= constraint (affine)
                   ├─ real variable ($(Convex.show_id(x)))
                   └─ 0

                termination status: OPTIMAL
                primal status: FEASIBLE_POINT
                dual status: FEASIBLE_POINT"""

        # test small `MAXDIGITS`
        x = Variable()
        old_maxdigits = Convex.MAXDIGITS[]
        Convex.MAXDIGITS[] = 2
        @test length(Convex.show_id(x)) == length("id: ") + 5
        Convex.MAXDIGITS[] = old_maxdigits

        # test large `MAXDIGITS`
        x = Variable()
        old_maxdigits = Convex.MAXDIGITS[]
        Convex.MAXDIGITS[] = 100
        @test length(Convex.show_id(x)) == length("id: ") + length(string(x.id_hash))
        Convex.MAXDIGITS[] = old_maxdigits

    end

    @testset "vartype and set_vartype" begin
        for x in (Variable(), Variable(1), ComplexVariable(2, 2))
            @test vartype(x) == ContVar

            vartype!(x, BinVar)
            @test vartype(x) == BinVar
            @test x.vartype == BinVar

            vartype!(x, IntVar)
            @test vartype(x) == IntVar
            @test x.vartype == IntVar

            vartype!(x, ContVar)
            @test vartype(x) == ContVar
            @test x.vartype == ContVar

        end
    end

    @testset "Constructors" begin

        # Constructors with sign
        for sgn in (Positive(), NoSign())
            for x in    [   # tuple size
                            Variable((2, 2), sgn),
                            Variable((2, 2), sgn, BinVar),
                            Variable((2, 2), sgn, :Bin),
                            # individual size
                            Variable(2, 2, sgn),
                            Variable(2, 2, sgn, BinVar),
                            Variable(2, 2, sgn, :Bin),
                            # single dimension
                            Variable(2, sgn),
                            Variable(2, sgn, BinVar),
                            Variable(2, sgn, :Bin),
                            # no dimension
                            Variable(sgn),
                            Variable(sgn, BinVar),
                            Variable(sgn, :Bin),  ]
                @test x isa Variable
                @test x isa Convex.AbstractVariable
                @test sign(x) == sgn
                @test x.sign == sgn
            end
        end

        # constructors without sign
        for x in    [   # tuple size
                        Variable((2, 2)),
                        Variable((2, 2), BinVar),
                        Variable((2, 2), :Bin),
                        # individual size
                        Variable(2, 2),
                        Variable(2, 2, BinVar),
                        Variable(2, 2, :Bin),
                        # single dimension
                        Variable(2),
                        Variable(2, BinVar),
                        Variable(2, :Bin),
                        # no dimension
                        Variable(),
                        Variable(BinVar),
                        Variable(:Bin),  ]
            @test x isa Variable
            @test x isa Convex.AbstractVariable
            @test sign(x) == NoSign()
            @test x.sign == NoSign()

            Convex.sign!(x, Positive())
            @test sign(x) == Positive()
            @test x.sign == Positive()
        end

        # ComplexVariable
        for x in    [   # tuple size
                    ComplexVariable((2, 2)),
                    Variable((2, 2), ComplexSign()),
                    ComplexVariable((2, 2), :Semidefinite),
                    # individual size
                    ComplexVariable(2, 2),
                    Variable(2, 2, ComplexSign()),
                    ComplexVariable(2, 2, :Semidefinite),
                    # single dimension
                    ComplexVariable(2),
                    Variable(2, ComplexSign()),
                    # no dimension
                    ComplexVariable(),
                    Variable(ComplexSign()),  ]
            @test x isa Variable
            @test x isa Convex.AbstractVariable
            @test sign(x) == ComplexSign()
            @test x.sign == ComplexSign()
        end

        for vt in (BinVar, IntVar), V in (ComplexVariable, Semidefinite, HermitianSemidefinite)
            @test_throws Any V(2; vartype=vt)
        end

        for vt in (:Bin, :Int), V in (Semidefinite, HermitianSemidefinite, ComplexVariable)
            @test_throws Any V(2, vt)
        end

        # Semidefinite
        for x in [
                Variable((2,2), :Semidefinite),
                Variable(2,2, :Semidefinite),
                ComplexVariable((2,2), :Semidefinite),
                ComplexVariable(2,2, :Semidefinite),
                HermitianSemidefinite((2,2)),
                HermitianSemidefinite(2, 2),
                HermitianSemidefinite(2),
                Semidefinite((2,2)),
                Semidefinite(2, 2),
                Semidefinite(2),
            ]
            @test length(constraints(x)) == 1
            @test constraints(x)[] isa Convex.SDPConstraint
        end

        @test_throws ErrorException HermitianSemidefinite(2,3)
        @test_throws ErrorException Semidefinite(2,3)
    end

    @testset "ConicObj" for T = [UInt32, UInt64]
        c = ConicObj()
        z = zero(T)
        @test !haskey(c, z)
        c[z] = (1, 1)
        @test c[z] == (1, 1)
        x = T[]
        for (k, v) in c
            push!(x, k)
        end
        @test x == collect(keys(c))
        d = copy(c)
        @test d !== c
    end

    @testset "length and size" begin
        x = Variable(2, 3)
        @test length(x) == 6
        @test size(x) == (2, 3)
        @test size(x, 1) == 2
        @test size(x, 2) == 3

        x = Variable(3)
        @test length(x) == 3
        @test size(x) == (3, 1)

        x = Variable()
        @test length(x) == 1
        @test size(x) == (1, 1)
    end

    @testset "lastindex and axes" begin
        x = Variable(2, 3)
        @test axes(x) == (Base.OneTo(2), Base.OneTo(3))
        @test axes(x, 1) == Base.OneTo(2)
        @test lastindex(x) == 6
        @test lastindex(x, 2) == 3
        y = x[:,end]
        @test y isa AbstractExpr
        @test size(y) == (2, 1)
    end

    @testset "Parametric constants" begin
        z = Constant([1.0 0.0im; 0.0 1.0])
        @test z isa Constant{Matrix{Complex{Float64}}}

        # Helper functions
        @test Convex.ispos(1)
        @test Convex.ispos(0)
        @test !Convex.ispos(-1)
        @test Convex.ispos([0,1,0])
        @test !Convex.ispos([0,-1,0])
        @test Convex.isneg(-1)
        @test Convex.isneg(0)
        @test !Convex.isneg(1)
        @test Convex.isneg([0,-1,0])
        @test !Convex.isneg([0,1,0])
        @test Convex._size(3) == (1, 1)
        @test Convex._sign(3) == Positive()
        @test Convex._size([-1,1,1]) == (3, 1)
        @test Convex._sign([-1,1,1]) == NoSign()
        @test Convex._sign([-1,-1,-1]) == Negative()
        @test Convex._size([0 0; 0 0]) == (2, 2)
        @test Convex._sign([0 0; 0 0]) == Positive()
        @test Convex._size(0 + 1im) == (1, 1)
        @test Convex._sign(0 + 1im) == ComplexSign()

        @test Convex.imag_conic_form(Constant(1.0)) == [0.0]
        @test Convex.imag_conic_form(Constant([1.0, 2.0])) == [0.0, 0.0]
    end

    @testset "#341: Evaluate for constants" begin
        A = rand(4,4)
        @test evaluate(Constant(A)) ≈ copy(A)
        @test Constant(A).size == (4,4)
        b = rand(4)
        @test evaluate(Constant(b)) ≈ copy(b)
        @test Constant(b).size == (4,1)
        c = 1.0
        @test evaluate(Constant(c)) ≈ c
        @test Constant(c).size == (1,1)

        @test evaluate(sumlargesteigs(Variable(4, 4), 0)) == 0
        @test evaluate(sumlargest(Variable(4), 0)) == 0
        @test evaluate(sumsmallest(Variable(4), 0)) == 0
    end

    @testset "Base.vect" begin
    # Issue #223: ensure we can make vectors of variables
        @test size([Variable(2), Variable(3, 4)]) == (2,)
    end

    @testset "Iteration" begin
        x = Variable(2, 3)
        s = sum([xi for xi in x])
        set_value!(x, [1 2 3; 4 5 6])
        # evaluate(s) == [21] (which might be wrong? expected 21)
        # but [21][1] === 21[1] === 21
        # so this should pass even after "fixing" that
        @test evaluate(s)[1] == 21

        x = Variable(4)
        @test [xi.inds for xi in x] == [1:1, 2:2, 3:3, 4:4]

        x = Variable(0)
        @test [xi for xi in x] == []
        @test iterate(x) == nothing
    end
    # returns [21]; not sure why
    # context("iteration") do
    #     x = Variable(2,3)
    #     s = sum([xi for xi in x])
    #     x.value = [1 2 3; 4 5 6]
    #     @fact evaluate(s) --> 21
    # end

    @testset "DCP warnings" begin
        # default is to log
        @test_logs (:warn, r"not DCP compliant") Convex.NotDcp()

        @eval Convex.emit_dcp_warnings() = false
        @test_logs  Convex.NotDcp()
        @eval Convex.emit_dcp_warnings() = true
        @test_logs (:warn, r"not DCP compliant") Convex.NotDcp()

    end

    @testset "`add_constraints!` (#380)" begin
        x = Variable(3, 3)
        p = minimize(norm_1(x))
        y = randn(3, 3)
        c = (norm2(x-y) < 1)
        @test length(p.constraints) == 0
        add_constraint!(p, c)
        @test length(p.constraints) == 1
        empty!(p.constraints)
        add_constraints!(p, c)
        @test length(p.constraints) == 1
        empty!(p.constraints)
        c2 = (norm2(x-rand(3,3)) < 3)
        add_constraints!(p, [c, c2])
        @test length(p.constraints) == 2
    end

    @testset "`diagm` (#401)" begin
        x = Variable(3)
        @test diagm(x) isa AbstractExpr
    end

    @testset "`is_psd` with type $T" for T in (Float64, Float32, Int, ComplexF64, BigFloat, Rational{BigInt}, Complex{Rational{BigInt}})
        A = zeros(T, 3, 3)
        A[1,1] = one(T)
        @test Convex.is_psd(A)
        @test Convex.is_psd(sparse(A))
        B = A .- one(T) / T(5000)
        @test !Convex.is_psd(B)
        @test !Convex.is_psd(sparse(B))

        # See https://github.com/jump-dev/Convex.jl/issues/452 for details
        C = [70.12718378756115 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -70.12718378756115 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 277.05985725905305 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 70.12718378756115 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 347.1870410466142 -103.46633673574595 0.0 -70.12718378756115 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 103.46633673574595 -103.46633673574595 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -70.12718378756115 0.0 -103.46633673574595 277.05985725905305 -103.46633673574595; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -103.46633673574595 103.46633673574595]
        @test Convex.is_psd(C)
        @test Convex.is_psd(sparse(C))
    end

    @testset "`assume_psd` option for `quadform`" begin
        A = [-1 0; 0 1] # neither PSD nor negative semidefinite
        x = Variable(2)

        @test_throws ErrorException quadform(x, A) # default
        @test_throws ErrorException quadform(x, A; assume_psd=false)
        @test quadform(x, A; assume_psd=true) isa Convex.AbstractExpr
    end

    @testset "`logsumexp` stability" begin
        v = Convex.Constant([1000, 1000, 1000])
        @test Convex.evaluate(Convex.logsumexp(v)) ≈ 1001.098612
    end
end


text:
# Example for modeling simultaneous sources and reverse-time migration
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using JUDI, PyPlot

## Set up model structure
n = (120, 100)	# (x, y, z) or (x, z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32, n) .+ 0.4f0
v0 = ones(Float32, n) .+ 0.4f0
v[:, Int(round(end/2)):end] .= 4.0f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup model structure
nsrc = 1    # one simultaneous source experiment
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

## Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(100f0, stop=100f0, length=nxrec)
timeR = 2000f0
dtR = 4f0

# Setup receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt = dtR, t = timeR)

## Set up source geometry
xsrc = [250f0, 500f0, 750f0, 1000f0]	# four simultaneous sources
ysrc = 0f0
zsrc = [50f0, 50f0, 50f0, 50f0]

# Source sampling and number of time steps
timeS = 2000f0
dtS = 4f0

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt = dtS, t = timeS)

# setup wavelets
f0 = 0.01	# source peak frequencies
q = ricker_wavelet(500f0, dtS, f0)  # 500 ms wavelet

# Create four source functions with different time shifts of q
wavelet1 = zeros(srcGeometry.nt[1])
wavelet1[1:1+length(q)-1] = q
wavelet2 = zeros(srcGeometry.nt[1])
wavelet2[41:41+length(q)-1] = q
wavelet3 = zeros(srcGeometry.nt[1])
wavelet3[121:121+length(q)-1] = q
wavelet4 = zeros(srcGeometry.nt[1])
wavelet4[201:201+length(q)-1] = q
wavelet = [wavelet1 wavelet2 wavelet3 wavelet4]	# Collect all wavelets
###################################################################################################

# Setup operators
Pr = judiProjection(recGeometry)
A_inv = judiModeling(model)
A0_inv = judiModeling(model0)
Ps = judiProjection(srcGeometry)
q = judiVector(srcGeometry, wavelet)

# Linearized modeling
d = Pr*A_inv*adjoint(Ps)*q  # simultaneous shot
J = judiJacobian(Pr*A0_inv*adjoint(Ps), q)
dD = J*dm   # simultaneous reflection data
rtm = adjoint(J)*dD # RTM

# Plot results
figure(); imshow(d.data[1], vmin = -1e1, vmax = 1e1, cmap = "seismic"); title("Simultaneous shot w/ four sources")
figure(); imshow(copy(adjoint(reshape(rtm,model.n))), ColorMap("gray")); title("Migrated simultaneous shot")


text:
pomdp = BabyPOMDP()

# testing constructor
b0 = DiscreteBelief(pomdp, [0.5,0.5])
@test pdf(b0,true) == 0.5
@test pdf(b0,false) == 0.5

println("There should be a warning below:")
DiscreteBelief(pomdp, [0.6, 0.5])
println("There should be a warning below:")
DiscreteBelief(pomdp, [-0.1, 1.1])

println("There should NOT be a warning below:")
DiscreteBelief(pomdp, [-0.1, 1.1], check=false)

# testing uniform belief
b1 = uniform_belief(pomdp)
@test pdf(b1,true) == 0.5
@test pdf(b1,false) == 0.5

# testing iterator
@test support(b1) == ordered_states(pomdp)

# testing equality (== function)
b2 = uniform_belief(pomdp)
b3 = DiscreteBelief(pomdp, [0.0,1.0])
@test b2 == b1
@test b2 == b0
@test b2 != b3

# testing hashing
@test hash(b0) == hash(b0)
@test hash(b0) == hash(b1)
@test hash(b1) == hash(b2)
@test hash(b2) != hash(b3)

# testing updater initialization
up = DiscreteUpdater(pomdp)
isd = initialstate_distribution(pomdp) # can be replaced with initialstate when pomdps 0.8 support is dropped
b4 = initialize_belief(up, isd)
@test pdf(b4,true) == pdf(isd,true)
@test pdf(b4,false) == pdf(isd,false)

# testing update function; if we feed baby, it won't be hungry
a = true
o = true
b4p = update(up, b4, a, o)
@test pdf(b4p,true) == 0.0
@test pdf(b4p,false) == 1.0

# if we don't feed the baby and observe crying
a = false
o = true
b4p = update(up, b4, false, true)
@test isapprox(pdf(b4p,true), 0.470588, atol=1e-4)
@test isapprox(pdf(b4p,false), 0.52941, atol=1e-4)

# testing that it works in a solve/simulation loop
# I'm not sure I need this test (could eliminate FIB dependency if not)
# r = test_solver(FIBSolver(), pomdp, max_steps=100)
# @test isapprox(r, -20.414855)

# Some more tests with tiger problem (old tests, but still work)
pomdp = TigerPOMDP()
up = DiscreteUpdater(pomdp)
bold = initialize_belief(up, initialstate_distribution(pomdp))

a = 0
o = true
bnew = update(up, bold, a, o)
@test isapprox(bnew.b, [0.15, 0.85])
@test isapprox(pdf(bnew, false), 0.15)
@test isapprox(pdf(bnew, true), 0.85)

# test mean and mode
b5 = DiscreteBelief(pomdp, [0.4, 0.6])
@test @inferred(mean(b5)) == 0.6
@test @inferred(mode(b5)) == true


text:
@testset "IntSetDomain.jl" begin
    @testset "IntSetDomain()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 20)

    end

    @testset "exclude!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 20)

        SeaPearl.saveState!(trailer)

        SeaPearl.exclude!(dom, 11)

        @test !SeaPearl.is_possible(dom, 11)

        SeaPearl.restoreState!(trailer)

        @test SeaPearl.is_possible(dom, 11)
    end

    @testset "require!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 20)

        SeaPearl.saveState!(trailer)

        @test SeaPearl.is_possible(dom, 10)
        @test !SeaPearl.is_required(dom, 10)

        SeaPearl.require!(dom, 10)

        @test SeaPearl.is_possible(dom, 10)
        @test SeaPearl.is_required(dom, 10)

        SeaPearl.restoreState!(trailer)

        @test SeaPearl.is_possible(dom, 10)
        @test !SeaPearl.is_required(dom, 10)
    end

    @testset "is_possible()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 20)

        SeaPearl.require!(dom, 11)
        SeaPearl.exclude!(dom, 10)

        @test SeaPearl.is_possible(dom, 11)
        @test SeaPearl.is_possible(dom, 12)
        @test !SeaPearl.is_possible(dom, 10)
        @test !SeaPearl.is_possible(dom, 9)
    end

    @testset "is_required()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 20)

        SeaPearl.require!(dom, 11)
        SeaPearl.exclude!(dom, 10)

        @test SeaPearl.is_required(dom, 11)
        @test !SeaPearl.is_required(dom, 12)
        @test !SeaPearl.is_possible(dom, 10)
        @test !SeaPearl.is_possible(dom, 9)
    end

    @testset "required_values()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 13)

        @test SeaPearl.required_values(dom) == Set{Int}()

        SeaPearl.require!(dom, 11)
        SeaPearl.require!(dom, 12)

        @test SeaPearl.required_values(dom) == Set{Int}([11, 12])
    end

    @testset "possible_not_required_values()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 13)

        @test SeaPearl.possible_not_required_values(dom) == Set{Int}([10, 11, 12, 13])

        SeaPearl.require!(dom, 11)
        SeaPearl.require!(dom, 12)

        @test SeaPearl.possible_not_required_values(dom) == Set{Int}([10, 13])
    end

    @testset "exchangePositions!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 5, 10)

        @test dom.values == [1, 2, 3, 4, 5, 6]
        @test dom.indexes == [1, 2, 3, 4, 5, 6]

        SeaPearl.exchangePositions!(dom, 6, 9)

        @test dom.values == [1, 5, 3, 4, 2, 6]
        @test dom.indexes == [1, 5, 3, 4, 2, 6]

        SeaPearl.exchangePositions!(dom, 9, 5)

        @test dom.values == [5, 1, 3, 4, 2, 6]
        @test dom.indexes == [2, 5, 3, 4, 1, 6]
    end

    @testset "require_all!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 13)
        SeaPearl.exclude!(dom, 12)
        SeaPearl.require!(dom, 11)

        SeaPearl.saveState!(trailer)

        SeaPearl.require_all!(dom)

        @test SeaPearl.required_values(dom) == Set{Int}([10, 11, 13])
        @test SeaPearl.possible_not_required_values(dom) == Set{Int}()

        SeaPearl.restoreState!(trailer)

        @test SeaPearl.required_values(dom) == Set{Int}([11])
        @test SeaPearl.possible_not_required_values(dom) == Set{Int}([10, 13])
    end

    @testset "exclude_all!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntSetDomain(trailer, 10, 13)
        SeaPearl.exclude!(dom, 12)
        SeaPearl.require!(dom, 11)

        SeaPearl.saveState!(trailer)

        SeaPearl.exclude_all!(dom)

        @test SeaPearl.required_values(dom) == Set{Int}([11])
        @test SeaPearl.possible_not_required_values(dom) == Set{Int}()

        SeaPearl.restoreState!(trailer)

        @test SeaPearl.required_values(dom) == Set{Int}([11])
        @test SeaPearl.possible_not_required_values(dom) == Set{Int}([10, 13])
    end

    @testset "reset_domain!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntSetVar(5, 10, "x", trailer)

        SeaPearl.exclude!(x.domain, 6)
        SeaPearl.require!(x.domain, 8)

        @test x.domain.values == [4, 6, 3, 1, 5, 2]
        @test x.domain.indexes == [4, 6, 3, 1, 5, 2]
        @test SeaPearl.length(SeaPearl.required_values(x.domain)) == 1
        @test SeaPearl.length(SeaPearl.possible_not_required_values(x.domain)) == 4
        SeaPearl.reset_domain!(x.domain)
        @test x.domain.values == [1, 2, 3, 4, 5, 6]
        @test x.domain.indexes == [1, 2, 3, 4, 5, 6]
        @test SeaPearl.length(SeaPearl.required_values(x.domain)) == 0
        @test SeaPearl.length(SeaPearl.possible_not_required_values(x.domain)) == 6
    end
end


text:
using LinearAlgebra
using SparseArrays
using PartitionedArrays
using Test
using IterativeSolvers
using Distances

function test_interfaces(parts)

  nparts = num_parts(parts)
  @assert nparts == 4

  parts2 = get_part_ids(parts)
  map_parts(parts,parts2) do part1, part2
    @test part1 == part2
  end

  parts_rcv = map_parts(parts) do part
    if part == 1
      [2,3]
    elseif part == 2
      [4,]
    elseif part == 3
      [1,2]
    else
      [1,3]
    end
  end

  parts_snd = map_parts(parts) do part
    if part == 1
      [3,4]
    elseif part == 2
      [1,3]
    elseif part == 3
      [1,4]
    else
      [2]
    end
  end

  data_snd = map_parts(i->10*i,parts_snd)
  data_rcv = map_parts(similar,parts_rcv)

  exchange!(
    data_rcv,
    data_snd,
    parts_rcv,
    parts_snd)

  map_parts(parts,data_rcv) do part, data_rcv
    if part == 1
      r = [10,10]
    elseif part == 2
      r = [20]
    elseif part == 3
      r = [30,30]
    else
      r= [40,40]
    end
    @test r == data_rcv
  end

  a = reduce_main(+,parts,init=0)
  @test get_main_part(a) == 1+2+3+4
  b = reduce_all(+,parts,init=0)
  map_parts(b) do b
    @test b == 1+2+3+4
  end
  @test reduce(+,parts,init=0) == 1+2+3+4
  @test sum(parts) == 1+2+3+4

  a = map_parts(parts) do part
    if part == 1
      4
    elseif part == 2
      2
    elseif part == 3
      6
    else
      3
    end
  end
  b = iscan(+,a,init=0)
  map_parts(parts,b) do part,b
    if part == 1
      @test b == 4
    elseif part == 2
      @test b == 6
    elseif part == 3
      @test b == 12
    else
      @test b == 15
    end
  end
  b, n = iscan(+,reduce,a,init=0)
  @test n == 15
  b, n = iscan_all(+,reduce,a,init=0)
  @test n == 15
  map_parts(b) do b
    @test b == [4,6,12,15]
  end

  b = xscan(+,a,init=1)
  map_parts(parts,b) do part,b
    if part == 1
      @test b == 1
    elseif part == 2
      @test b == 5
    elseif part == 3
      @test b == 7
    else
      @test b == 13
    end
  end

  b, n = xscan(+,reduce,a,init=1)
  @test n == 15+1
  b, n = xscan_all(+,reduce,a,init=1)
  @test n == 15+1
  map_parts(b) do b
    @test b == [1,5,7,13]
  end

  data_rcv, t = async_exchange(
    data_snd,
    parts_rcv,
    parts_snd)

  map_parts(schedule,t)
  map_parts(wait,t)

  map_parts(parts,data_rcv) do part, data_rcv
    if part == 1
      r = [10,10]
    elseif part == 2
      r = [20]
    elseif part == 3
      r = [30,30]
    else
      r= [40,40]
    end
    @test r == data_rcv
  end

  data_rcv = exchange(
    data_snd,
    parts_rcv,
    parts_snd)

  map_parts(parts,data_rcv) do part, data_rcv
    if part == 1
      r = [10,10]
    elseif part == 2
      r = [20]
    elseif part == 3
      r = [30,30]
    else
      r= [40,40]
    end
    @test r == data_rcv
  end

  parts_snd_2 = discover_parts_snd(parts_rcv)
  map_parts(parts_snd,parts_snd_2) do parts_snd, parts_snd_2
    @test parts_snd == parts_snd_2
  end

  PartitionedArrays.ERROR_DISCOVER_PARTS_SND[] = true
  @test_throws ErrorException parts_snd_2 = discover_parts_snd(parts_rcv)
  PartitionedArrays.ERROR_DISCOVER_PARTS_SND[] = false

  n = 10

  partition = map_parts(parts) do part
    if part == 1
      IndexSet(part,[1,2,3,5,7,8],Int32[1,1,1,2,3,3])
    elseif part == 2
      IndexSet(part,[2,4,5,10],Int32[1,2,2,4])
    elseif part == 3
      IndexSet(part,[6,7,8,5,4,10],Int32[3,3,3,2,2,4])
    else
      IndexSet(part,[1,3,7,9,10],Int32[1,1,3,4,4])
    end
  end

  exchanger = Exchanger(partition)

  map_parts(parts,exchanger.parts_snd,exchanger.lids_snd) do part, parts_snd, lids_snd
    if part == 1
      _parts = [2,4]
      ids = [[2],[1,3]]
    elseif part == 2
      _parts = [1,3]
      ids = [[3],[3,2]]
    elseif part == 3
      _parts = [1,4]
      ids = [[2,3],[2]]
    else
      _parts = [2,3]
      ids = [[5],[5]]
    end
    @test _parts == parts_snd
    @test ids == lids_snd
  end

  values = map_parts(parts,partition) do part, partition
    values = fill(0.0,num_lids(partition))
    for lid in 1:length(partition.lid_to_part)
      owner = partition.lid_to_part[lid]
      if owner == part
        values[lid] = 10*part
      end
    end
    values
  end

  exchange!(values,exchanger)

  map_parts(values,partition) do values, partition
    for lid in 1:length(partition.lid_to_part)
      owner = partition.lid_to_part[lid]
      @test values[lid] == 10*owner
    end
  end

  values_rcv = map_parts(values) do values
    fill(10.0,length(values))
  end

  values_snd = map_parts(values) do values
    fill(20.0,length(values))
  end

  exchange!(values_rcv,values_snd,exchanger)

  map_parts(values_rcv,partition) do values_rcv, partition
    for lid in 1:length(partition.lid_to_part)
      if partition.lid_to_part[lid] == partition.part
        @test values_rcv[lid] == 10.0
      else
        @test values_rcv[lid] == 20.0
      end
    end
  end

  map_parts(values_snd) do values_snd
    @test all(values_snd .== 20.0)
  end

  values = map_parts(parts,partition) do part, partition
    values = [ zeros(Int,3) for lid in 1:num_lids(partition)]
    for lid in partition.oid_to_lid
      gid = partition.lid_to_gid[lid]
      for i in 1:3
        values[lid][i] = 100*part + 10*gid + i
      end
    end
    Table(values)
  end

  exchange!(values,exchanger)

  map_parts(parts,partition,values) do part, partition, values
    for lid in 1:num_lids(partition)
      gid = partition.lid_to_gid[lid]
      owner = partition.lid_to_part[lid]
      for i in 1:3
        @test values[lid][i] == 100*owner + 10*gid + i
      end
    end
  end

  exchanger_rcv = exchanger # receives data at ghost ids from remote parts
  exchanger_snd = reverse(exchanger_rcv) # sends data at ghost ids to remote parts

  values = map_parts(parts,partition) do part, partition
    values = fill(0.0,num_lids(partition))
    for lid in 1:length(partition.lid_to_part)
      values[lid] = 10*part
    end
    values
  end
  exchange!(+,values,exchanger_snd)
  exchange!(values,exchanger_rcv)

  ids = PRange(n,partition)
  ids2 = copy(ids)

  @test ids !== ids2
  @test ids.partition !== ids2.partition
  @test ids.exchanger !== ids2.exchanger

  map_parts(ids.partition) do partition
    @test get_part(partition) == partition.part
    @test get_lid_to_gid(partition) == partition.lid_to_gid
    @test get_lid_to_part(partition) == partition.lid_to_part
    @test get_oid_to_lid(partition) == partition.oid_to_lid
    @test get_hid_to_lid(partition) == partition.hid_to_lid
    @test get_lid_to_ohid(partition) == partition.lid_to_ohid
    @test get_gid_to_lid(partition) == partition.gid_to_lid
  end
  @test num_parts(ids) == nparts
  @test num_gids(ids) == n

  ids2 = PRange(parts,n)
  @test ids2.ghost == false

  gids = map_parts(parts) do part
    if part == 1
      [1,4,6]
    elseif part == 2
      [3,1,2,8]
    elseif part == 3
      [1,9,6]
    else
      [3,2,8,10]
    end
  end

  i_to_part = map_parts(getindex,ids2.gid_to_part,gids)
  ids3 = add_gids(ids2,gids,i_to_part)
  @test ids3.ghost == true

  gids2 = map_parts(parts) do part
    if part == 1
      [4,6]
    elseif part == 2
      [1,2]
    elseif part == 3
      [1,9]
    else
      [3]
    end
  end

  hids = touched_hids(ids3,gids2)
  map_parts(hids,gids2,ids3.partition) do hids, gids, ids3
    @test gids == ids3.lid_to_gid[ids3.hid_to_lid[hids]]
  end

  ids3 = add_gids(ids2,gids)
  @test ids3.ghost == true
  to_lids!(gids,ids3)
  to_gids!(gids,ids3)

  a = map_parts(parts) do part
    if part == 1
      4
    elseif part == 2
      2
    elseif part == 3
      6
    else
      3
    end
  end
  ids5 = PRange(parts,a)
  map_parts(parts,ids5.partition,ids5.gid_to_part) do part, ids5, gid_to_part
    if part == 1
      @test ids5.lid_to_gid == [1, 2, 3, 4]
    elseif part == 2
      @test ids5.lid_to_gid == [5, 6]
    elseif part == 3
      @test ids5.lid_to_gid == [7, 8, 9, 10, 11, 12]
    else
      @test ids5.lid_to_gid == [13, 14, 15]
    end
    @test gid_to_part == [1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4]
  end
  ids5 = PRange(parts,reduce(+,a,init=0),a)
  ids5 = PRange(parts,reduce(+,a,init=0),a,xscan(+,a,init=1))

  hid_to_gid, hid_to_part = map_parts(parts) do part
    Int[], Int32[]
  end
  ids5 = PRange(parts,reduce(+,a,init=0),a,xscan(+,a,init=1),hid_to_gid,hid_to_part)

  if ndims(parts) > 1

    ids4 = PRange(parts,(5,4))
    @test ids4.ghost == false
    @test num_gids(ids4) == 4*5
    map_parts(parts,ids4.partition,ids4.gid_to_part) do part, ids4, gid_to_part
      if part == 1
        @test ids4.lid_to_gid == [1, 2, 6, 7]
      elseif part == 2
        @test ids4.lid_to_gid == [3, 4, 5, 8, 9, 10]
      elseif part == 3
        @test ids4.lid_to_gid == [11, 12, 16, 17]
      else
        @test ids4.lid_to_gid == [13, 14, 15, 18, 19, 20]
      end
      @test gid_to_part == [1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 3, 3, 4, 4, 4]
    end

    pcis = PCartesianIndices(parts,(5,4))
    map_parts(parts,pcis) do part, cis
      if part == 1
        @test cis == CartesianIndices((1:2,1:2))
      elseif part == 2
        @test cis == CartesianIndices((3:5,1:2))
      elseif part == 3
        @test cis == CartesianIndices((1:2,3:4))
      else
        @test cis == CartesianIndices((3:5,3:4))
      end
    end

    pcis = PCartesianIndices(parts,(5,4),no_ghost)
    pcis = PCartesianIndices(parts,(5,4),with_ghost)
    map_parts(parts,pcis) do part, cis
      if part == 1
        @test cis == CartesianIndices((1:3,1:3))
      elseif part == 2
        @test cis == CartesianIndices((2:5,1:3))
      elseif part == 3
        @test cis == CartesianIndices((1:3,2:4))
      else
        @test cis == CartesianIndices((2:5,2:4))
      end
    end

    ids4 = PRange(parts,(5,4),no_ghost)
    ids4 = PRange(parts,(5,4),with_ghost)
    map_parts(parts,ids4.partition,ids4.gid_to_part) do part, ids4,gid_to_part
      if part == 1
        @test ids4.lid_to_gid == [1, 2, 3, 6, 7, 8, 11, 12, 13]
      elseif part == 2
        @test ids4.lid_to_gid == [2, 3, 4, 5, 7, 8, 9, 10, 12, 13, 14, 15]
      elseif part == 3
        @test ids4.lid_to_gid == [6, 7, 8, 11, 12, 13, 16, 17, 18]
      else
        @test ids4.lid_to_gid == [7, 8, 9, 10, 12, 13, 14, 15, 17, 18, 19, 20]
      end
      @test gid_to_part == [1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 3, 3, 4, 4, 4]
    end
    @test ids4.ghost = true

  end

  v = PVector(gids,map_parts(copy,gids),ids3;ids=:global)
  v = PVector(gids,map_parts(copy,gids),ids3;ids=:local)
  v = PVector(gids,map_parts(copy,gids),n;ids=:global)
  u = 2*v
  map_parts(u.values,v.values) do u,v
    @test u == 2*v
  end
  u = v + u
  map_parts(u.values,v.values) do u,v
    @test u == 3*v
  end

  @test any(i->i>4,v) == true
  @test any(i->i>10,v) == false
  @test all(i->i<10,v) == true
  @test all(i->i<4,v) == false
  @test maximum(v) == 5
  @test minimum(v) == 0
  @test maximum(i->i-1,v) == 4
  @test minimum(i->i-1,v) == -1

  w = copy(v)
  rmul!(w,-1)
  @test all(i->i==0,v+w)

  @test w == w
  @test w != v

  @test sqeuclidean(w,v) ≈ (norm(w-v))^2
  @test euclidean(w,v) ≈ norm(w-v)

  w = similar(v)
  w = similar(v,Float64)
  w = similar(v,Float64,ids3)
  w = similar(typeof(v),ids3)
  w = zero(v)
  @test norm(w) == 0
  @test sum(w) == 0

  w = v .- u
  @test isa(w,PVector)
  w =  1 .+ v
  @test isa(w,PVector)
  w =  v .+ 1
  @test isa(w,PVector)
  w =  v .+ w .- u
  @test isa(w,PVector)
  w =  v .+ 1 .- u
  @test isa(w,PVector)

  w .= v .- u
  w .= v .- 1 .- u
  w .= u
  map_parts(w.values,u.values) do w,u
    @test w == u
  end

  u = PVector(1.0,ids2)
  w = PVector(3.0,ids3)

  @test oids_are_equal(u.rows,u.rows)
  @test hids_are_equal(u.rows,u.rows)
  @test lids_are_equal(u.rows,u.rows)
  @test oids_are_equal(u.rows,w.rows)
  @test !hids_are_equal(u.rows,w.rows)
  @test !lids_are_equal(u.rows,w.rows)

  w = v .- u
  @test isa(w,PVector)
  w =  v .+ w .- u
  @test isa(w,PVector)
  w =  v .+ 1 .- u
  @test isa(w,PVector)

  w .= v .- u
  w .= v .- 1 .- u
  w .= u
  map_parts(w.owned_values,u.owned_values) do w,u
    @test w == u
  end

  map_parts(parts,local_view(v,v.rows)) do part,v
    if part == 3
      v[1] = 6
    end
  end

  map_parts(parts,local_view(v,v.rows)) do part,v
    if part == 3
      @test v[1] == 6
    end
  end

  subrows_partition = map_parts(v.rows.partition) do rows
    perm = [3,1,2]
    IndexSet(rows.part,rows.lid_to_gid[perm],rows.lid_to_part[perm])
  end
  subrows = PRange(length(v.rows),subrows_partition)
  subv = PVector(1.0,subrows)
  map_parts(local_view(subv,v.rows)) do v
    v[[1,2,3]] = [6,7,8]
  end
  map_parts(subv.values) do v
    perm = [2,3,1]
    @test v[perm] == [6,7,8]
  end
  map_parts(local_view(subv,v.rows)) do v
    @test v[[1,2,3]] == [6,7,8]
  end

  map_parts(parts,global_view(v,v.rows)) do part,v
    if part == 4
      v[9] = 6
    end
  end

  map_parts(parts,global_view(v,v.rows)) do part,v
    if part == 4
      @test v[9] == 6
    end
  end

  v = PVector{Float64}(undef,ids)
  fill!(v,1.0)

  v = PVector{Float64}(undef,ids)
  map_parts(parts,v.values,v.rows.partition) do part, values, partition
    for lid in 1:length(partition.lid_to_part)
      owner = partition.lid_to_part[lid]
      if owner == part
        values[lid] = 10*part
      end
    end
  end
  exchange!(v)
  map_parts(parts,v.values,v.rows.partition) do part, values, partition
    for lid in 1:length(partition.lid_to_part)
      owner = partition.lid_to_part[lid]
      @test values[lid] == 10*owner
    end
  end

  assemble!(v)

  cols = ids
  rows = cols

  values = map_parts(rows.partition,cols.partition) do rows, cols
    i = collect(1:num_lids(rows))
    j = i
    v = fill(2.0,length(i))
    a=sparse(i,j,v,num_lids(rows),num_lids(cols))
    a
  end

  x = PVector{Float64}(undef,cols)
  fill!(x,3)
  b = PVector{Float64}(undef,rows)

  A = PSparseMatrix(values,rows,cols)
  A = PSparseMatrix(values,rows,cols,A.exchanger)
  mul!(b,A,x)

  map_parts(b.owned_values) do values
    @test all( values .== 6 )
  end

  exchange!(b)
  map_parts(b.values) do values
    @test all( values .== 6 )
  end

  LinearAlgebra.fillstored!(A,1.0)
  fill!(x,3.0)
  mul!(b,A,x)
  exchange!(b)
  map_parts(b.values) do values
    @test all( values .== 3 )
  end

  exchange!(A)
  assemble!(A)

  I,J,V = map_parts(parts) do part
    if part == 1
      [1,2,1,2], [2,6,1,2], [1.0,2.0,30.0,10.0]
    elseif part == 2
      [3,3,4], [3,8,4], [10.0,2.0,30.0]
    elseif part == 3
      [5,5,6,7], [5,6,6,7], [10.0,2.0,30.0,1.0]
    else
      [9,9,8,10], [9,2,8,10], [10.0,2.0,30.0,50.0]
    end
  end
  A = PSparseMatrix(I,J,V,n,n;ids=:global)
  local_view(A,A.rows,A.cols)
  local_view(A,copy(A.rows),copy(A.cols))
  global_view(A,A.rows,A.cols)

  x = PVector{Float64}(undef,A.cols)
  fill!(x,1.0)
  y = A*x
  dy = y - y

  y = PVector(1.0,A.rows)
  x = IterativeSolvers.cg(A,y)

  x = PVector(0.0,A.cols)
  IterativeSolvers.cg!(x,A,y)
  fill!(x,0.0)

  x = A\y
  @test isa(x,PVector)
  r = A*x-y
  @test norm(r) < 1.0e-9

  to_gids!(I,A.rows)
  to_gids!(J,A.cols)
  exchange!(I,J,V,A.cols)

  factors = lu(A)
  x .= 0
  ldiv!(x,factors,y)
  r = A*x-y
  @test norm(r) < 1.0e-9

  lu!(factors,A)
  x .= 0
  ldiv!(x,factors,y)
  r = A*x-y
  @test norm(r) < 1.0e-9

end


text:
module MetaBalls

include("metaball.jl")
include("visualize.jl")

end


text:
using GeometricFlux
using Flux
using Flux: onehotbatch, onecold, crossentropy, throttle
using JLD2
using StatsBase
using SparseArrays
using LightGraphs.SimpleGraphs

@load "data/cora_features.jld2" features
@load "data/cora_labels.jld2" labels
@load "data/cora_graph.jld2" g

num_nodes = 2708
num_features = 1433
target_catg = 7

## Preprocessing data
train_X = features  # dim: num_features * num_nodes
train_y = labels  # dim: target_catg * num_nodes

## Model
model = Chain(GCNConv(g, num_features=>1000, relu),
              GCNConv(g, 1000=>500, relu),
              GCNConv(g, 500=>7),
              softmax)
# test model
# model(train_X)

## Loss
loss(x, y) = crossentropy(model(x), y)
accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

## Training
ps = Flux.params(model)
train_data = [(train_X, train_y)]
opt = ADAM(0.01)
evalcb() = @show(accuracy(train_X, train_y))

Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))


text:
"""
BFS
---
    compute breadth first search distances and returns a distance(d),
    the discover time(dt), predecessor array(pred) in the tuple (d,dt,pred)
    pred[i] = 0 if vertex i is in a component not reachable from u and i != u.
    Search stops when it hits the vertex target.

Functions
---------
* (d,dt,pred) = bfs(A::MatrixNetwork,u::Int64,target::Int64)
* (d,dt,pred) = bfs{T}(A::SparseMatrixCSC{T,Int64}),u::Int64,target::Int64)
* (d,dt,pred) = bfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64,target::Int64)\n
If target is not specified, it is assigned to 0

Example
-------
~~~
A = load_matrix_network("bfs_example")
(d,dt,pred) = bfs(A,1)
~~~
"""

function bfs(A::MatrixNetwork,u::Int64,target::Int64)
    (rp,ci) = (A.rp,A.ci)
    n=length(rp)-1
    d=-1*ones(Int64,n)
    dt=-1*ones(Int64,n)
    pred=zeros(Int64,n)
    sq=zeros(Int64,n)
    sqt=0
    sqh=0 # search queue and search queue tail/head

    # start bfs at u
    sqt=sqt+1
    sq[sqt]=u
    t=0
    d[u]=0
    dt[u]=t
    t=t+1
    pred[u]=u
    while sqt-sqh>0
        sqh=sqh+1
        v=sq[sqh] # pop v off the head of the queue
        for ri=rp[v]:rp[v+1]-1
            w=ci[ri]
            if d[w]<0
                sqt=sqt+1
                sq[sqt]=w
                d[w]=d[v]+1
                dt[w]=t
                t=t+1
                pred[w]=v
                if w==target
                    return (d,dt,pred)
                end
            end
        end
    end
    return (d,dt,pred)
end

############################
### Additional functions ###
############################

function bfs(A::MatrixNetwork,u::Int64)
    return bfs(A,u,0)
end

## CSC sparse matrices:
function bfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64)
    return bfs(MatrixNetwork(A),u)
end

function bfs{T}(A::SparseMatrixCSC{T,Int64},u::Int64,target::Int64)
    return bfs(MatrixNetwork(A),u,target)
end

## Triplet Format:
function bfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64)
    return bfs(MatrixNetwork(ei,ej),u)
end

function bfs(ei::Vector{Int64},ej::Vector{Int64},u::Int64,target::Int64)
    return bfs(MatrixNetwork(ei,ej),u,target)
end

## CSR sparse matrices:
function bfs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64,u::Int64)
    return bfs(MatrixNetwork(n,rp,ci,vals),u)
end

function bfs{T}(rp::Vector{Int64},ci::Vector{Int64},vals::Vector{T},n::Int64,u::Int64,target::Int64)
    return bfs(MatrixNetwork(n,rp,ci,vals),u,target)
end


text:
#Assign a value of 13 to the variable x
x = 13;
println("Variable x has the value $x");

#Perform math with x's value
sum = x + 5;
println("The sum of x + 5 is $sum");

#Assing a new type to existing variable
x = "x has been changed to this string instead of a number";
println("Variables can be assinged to with newe types: ['$x']");

# Flexiable naming of Variables
# Variable names are case-sensitive, and have no sematic meaning
# (that is, the language will not streat variables differently based on their names)

z = "This is the variable z";
Z = "This is the variable Z";
println("z:['$z'] and Z:['$Z']");

number_five = 5;
println("The variable number_five has value: [$number_five]");

alsoNumberFive = 5;
println("The variable alsoNumberFive has value: [$alsoNumberFive]");

#Universal Declaration Of Human Rights Start
UniversalLanguageSupport = "人人生而自由，在尊严和权利上一律平等。";
println("Universal Language Support: [$UniversalLanguageSupport]");

#Unicode names (UTF-8 encoding) are allowed:
δ = 0.00001;
println("δ = $δ");

안녕하세요 = "Hello";
println("안녕하세요 variable holds the value: ['$안녕하세요']");

#Redefine built-in constants and functions if needed (although this is not recommended to avoid potential confusions):

#redefine built-in constant pi
println("pi constant $pi");

#However, if you try to redefine a built-in constant or function already in use, Julia will give you an error:
#ERROR: cannot assign a value to variable MathConstants.pi from module Main
#pi = 3
#println("pi constant redefined $pi");

#Stylistic Conventions
#While Julia imposes few restrictions on valid names, it has become useful to adopt the following conventions:

#Names of variables are in lower case.
#Word separation can be indicated by underscores ('_'), but use of underscores is discouraged unless the name would be hard to read otherwise.
#Names of Types and Modules begin with a capital letter and word separation is shown with upper camel case instead of underscores.
#Names of functions and macros are in lower case, without underscores.
#Functions that write to their arguments have names that end in !. These are sometimes called "mutating" or "in-place" functions because they are intended to produce changes in their arguments after the function is called, not just return a value.

max_stars_repo_name: CubicleJockey/julia-playground

max_issues_repo_name: CubicleJockey/julia-playground

max_forks_repo_name: CubicleJockey/julia-playground


text:
"""
Save the Workout at the end of every epoch.

# Usage
```julia
fit!(workout, data, cb=EpochSave())
```

"""
struct EpochSave
end

function (c::EpochSave)(workout::Workout, phase::Symbol)
    if phase == :valid
        saveWorkout(workout)
    end
end

"""
Save the Workout if a certain metric has improved

# Usage
```julia

# save as long as the validation loss is declining
fit!(workout, data, cb=AutoSave(:val_loss))
```

"""
mutable struct AutoSave
    value::Float64
    metric::Symbol
    AutoSave(metric::Symbol) = new(Inf, metric)
end

function (c::AutoSave)(workout::Workout, phase::Symbol)
    if phase == :valid
        getmetricvalue(workout, c.metric) do x
            if x < c.value
                saveWorkout(workout)
                c.value = x
            end
        end
    end
end

"""
Stop the training if a certain metric didn't improve
"""
mutable struct EarlyStop
    value::Float64
    metric::Symbol
    EarlyStop(metric::Symbol) = new(Inf, metric)
end

function (c::EarlyStop)(workout::Workout, phase::Symbol)
    if phase == :valid
        getmetricvalue(workout, c.metric) do x
            if x > c.value
                stop(workout,"$(c.metric) didn't improve anymore")
            else
                c.value = x
            end
        end
    end
end