text:
################################################################################
### For generic methods
################################################################################

# Heuristic for selecting the level with the (presumably) most reliable
# standard error estimate:
# Take the highest lvl with at least 32 bins.
# (Chose 32 based on https://doi.org/10.1119/1.3247985)
function _reliable_level(B::ErrorPropagator{T,N})::Int64 where {T,N}
    isempty(B) && (return 1)                # results in NaN in std_error
    i = findlast(x -> x >= 32, B.count)
    return something(i, 1)
end

_eachlevel(B::ErrorPropagator) = 1:findlast(x -> x > 1, B.count)

################################################################################
### Statistics
################################################################################

"""
    mean(B::ErrorPropagator, i[, lvl = 1])

Calculates the mean for the i-th argument of the error propagator at an optional
binning level.
"""
function mean(B::ErrorPropagator, i::Integer, lvl = 1)
    return B.sums1D[lvl][i] / B.count[lvl]
end

"""
    var(B::ErrorPropagator, i[, lvl])

Calculates the variance for the i-th argument of the error propagator at a given
binning level.
"""
function var(B::ErrorPropagator{<: Real}, i::Integer, lvl = _reliable_level(B))
    n = B.count[lvl]
    X = B.sums1D[lvl][i]
    X2 = B.sums2D[lvl][i, i]

    # lvl = 1 <=> original values
    # correct variance:
    # (∑ xᵢ^2) / (N-1) - (∑ xᵢ)(∑ xᵢ) / (N(N-1))
    return X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(B::ErrorPropagator{<: Complex}, i::Integer, lvl = _reliable_level(B))
    n = B.count[lvl]
    X = B.sums1D[lvl][i]
    X2 = B.sums2D[lvl][i, i]

    return (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

function var(B::ErrorPropagator{<: AbstractArray{<: Real}}, i::Integer, lvl = _reliable_level(B))
    n = B.count[lvl]
    X = B.sums1D[lvl][i]
    X2 = B.sums2D[lvl][i, i]

    return @. X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(B::ErrorPropagator{<: AbstractArray{<: Complex}}, i::Integer, lvl = _reliable_level(B))
    n = B.count[lvl]
    X = B.sums1D[lvl][i]
    X2 = B.sums2D[lvl][i, i]

    return @. (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

"""
    varN(B::ErrorPropagator, i[, lvl])

Calculates the variance/N for the i-th argument of the error propagator at a
given binning level.
"""
function varN(B::ErrorPropagator, i::Integer, lvl = _reliable_level(B))
    return var(B, i, lvl) / B.count[lvl]
end

"""
    tau(B::ErrorPropagator, i[, lvl])

Calculates the autocorrelation time tau for the i-th argument of the error
propagator at a given binning level.
"""
function tau(B::ErrorPropagator{<: Number}, i::Integer, lvl = _reliable_level(B))
    return 0.5 * (varN(B, i, lvl) / varN(B, i, 1) - 1)
end
function tau(B::ErrorPropagator{<: AbstractArray}, i::Integer, lvl = _reliable_level(B))
    return 0.5 * (varN(B, i, lvl) ./ varN(B, i, 1) .- 1)
end
autocorrelation(B::ErrorPropagator, i::Integer, lvl = _reliable_level(B)) = tau(B, i, lvl)
autocorrelation_time(B::ErrorPropagator, i::Integer, lvl = _reliable_level(B)) = tau(B, i, lvl)

"""
    std_error(B::ErrorPropagator, i[, lvl])

Calculates the standard error of the mean for the i-th argument of the error
propagator at a given binning level.
"""
function std_error(B::ErrorPropagator{<: Number}, i::Integer, lvl=_reliable_level(B))
    return sqrt(varN(B, i, lvl))
end
function std_error(B::ErrorPropagator{<: AbstractArray}, i::Integer, lvl=_reliable_level(B))
    return sqrt.(varN(B, i, lvl))
end

################################################################################
### Additional methods
################################################################################

for name in [:mean, :varN, :var, :tau, :std_error, :autocorrelation, :autocorrelation_time]
    # generates functions fs(B[, lvl]) = [f(B, 1, lvl), .., f(B, N_args, lvl)]
    @eval begin
        function $(Symbol(name, :s))(B::ErrorPropagator, lvl=_reliable_level(B))
            return [$name(B, i, lvl) for i in eachindex(B.sums1D[1])]
        end
    end

    # all_fs(ep) = [[f(ep, 1, lvl), .., f(ep, N_args, lvl)] for lvl in eachlvl]
    @eval begin
        function $(Symbol(:all_, name, :s))(B::ErrorPropagator)
            return [$(Symbol(name, :s))(B, lvl) for lvl in _eachlevel(B)]
        end
    end
end

# Docs
@doc """
    varNs(B::ErrorPropagator[, lvl])

Calculates the variance/N for each argument of the error propagator at a given
binning level.
""" varNs
@doc """
    vars(B::ErrorPropagator[, lvl])

Calculates the variance for each argument of the error propagator at a given
binning level.
""" vars
@doc """
    means(B::ErrorPropagator[, lvl])

Calculates the mean for each argument of the error propagator at a given
binning level.
""" means
@doc """
    std_errors(B::ErrorPropagator[, lvl])

Calculates the standard error of the mean for each argument of the error
propagator at a given binning level.
""" std_errors

"""
    covmat(B::ErrorPropagator[, lvl])

Returns the covariance matrix for a given level of the error propgator.
"""
function covmat(
        B::ErrorPropagator{T, N}, lvl = _reliable_level(B)
    ) where {N , T <: Number}

    invN = 1.0 / B.count[lvl]
    invN1 = 1.0 / (B.count[lvl] - 1)
    return [
        (
            B.sums2D[lvl][i, j] -
            B.sums1D[lvl][i] * conj(B.sums1D[lvl][j]) * invN
        ) * invN1
        for i in eachindex(B.sums1D[lvl]), j in eachindex(B.sums1D[lvl])
    ]
end
function covmat(
        B::ErrorPropagator{T, N}, lvl = _reliable_level(B)
    ) where {N , T <: AbstractArray}

    invN = 1.0 / B.count[lvl]
    invN1 = 1.0 / (B.count[lvl] - 1)
    return [
        (
            B.sums2D[lvl][i, j] -
            B.sums1D[lvl][i] .* conj(B.sums1D[lvl][j]) * invN
        ) * invN1
        for i in eachindex(B.sums1D[lvl]), j in eachindex(B.sums1D[lvl])
    ]
end

"""
    var(B::ErrorPropagator, gradient[, lvl])

Gives the first-order variance estimate of a function `f` acting on the
arguments of the error propagator. `gradient` is either the gradient of `f` (a
function) or a vector `∇f(means(B))`. To get an estimate mean value of `f`,
`mean(B, f)` can be used.
"""
function var(
        B::ErrorPropagator{T, N},
        gradient::Vector,
        lvl = _reliable_level(B)
    ) where {T <: Real, N}

    result = 0.0
    invN = 1.0 / B.count[lvl]
    # invNN1 = 1.0 / (B.count[lvl] * (B.count[lvl] - 1))
    invN1 = 1.0 / (B.count[lvl] - 1)
    for i in eachindex(B.sums1D[lvl])
        for j in eachindex(B.sums1D[lvl])
            result += gradient[i] * gradient[j] * (
                B.sums2D[lvl][i, j] -
                B.sums1D[lvl][i] * B.sums1D[lvl][j] * invN
            ) * invN1
        end
    end

    return result
end

function var(
        B::ErrorPropagator{T, N},
        gradient::Vector,
        lvl = _reliable_level(B)
    ) where {T <: Complex, N}

    result = 0.0 * 0.0im
    invN = 1.0 / B.count[lvl]
    # invNN1 = 1.0 / (B.count[lvl] * (B.count[lvl] - 1))
    invN1 = 1.0 / (B.count[lvl] - 1)
    for i in eachindex(B.sums1D[lvl])
        for j in eachindex(B.sums1D[lvl])
            result += gradient[i] * conj(gradient[j]) * (
                B.sums2D[lvl][i, j] -
                B.sums1D[lvl][i] * conj(B.sums1D[lvl][j]) * invN
            ) * invN1
        end
    end

    return abs(result)
end

# Wrappers
function var(B::ErrorPropagator, gradient::Function, lvl = _reliable_level(B))
    grad = gradient(means(B))
    if typeof(grad) <: AbstractArray
        return var(B, grad, lvl)
    else
        return var(B, [grad], lvl)
    end
end

"""
    varN(B::ErrorPropagator, gradient[, lvl])

    Gives the first-order variance/N estimate of a function `f` acting on the
    arguments of the error propagator. `gradient` is either the gradient of `f` (a
    function) or a vector `∇f(means(B))`. To get an estimate mean value of `f`,
    `mean(B, f)` can be used.
"""
function varN(B::ErrorPropagator, gradient::Function, lvl = _reliable_level(B))
    return var(B, gradient, lvl) / B.count[lvl]
end

"""
    std_error(B::ErrorPropagator, gradient[, lvl])

Gives the first-order standard error estimate of a function `f` acting on the
arguments of the error propagator. `gradient` is either the gradient of `f` (a
function) or a vector `∇f(means(B))`. To get an estimate mean value of `f`,
`mean(B, f)` can be used.
"""
function std_error(B::ErrorPropagator, gradient::Function, lvl = _reliable_level(B))
    return sqrt(var(B, gradient, lvl) / B.count[lvl])
end

"""
    mean(B, f[, lvl=1])

Returns an estimate for the mean value of `f`, where `f` is a function acting
on the sample pushed the error porpagator. `f` must be of the form `f(v)`, where
`v = [mean_arg1, mean_arg2, ..., mean_argN]` is a vector containing the averages
of each argument pushed to the error propagator.
"""
function mean(B::ErrorPropagator, f::Function, lvl = 1)
    return f(means(B, lvl))
end

ext: jl

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
export cross, backcross

function sample_genotype(geno::T, n_alleles::Int) where T<:Genotype
    sample(SVector(geno), n_alleles, replace = false)
end

function sample_genotype(geno::Missing, n_alleles::Int)
    return [missing]
end

function haploid_cross!(data::DataFrame, p1::T, p2::T; n::Int) where T <: GenoArray
    iter_df = DataFrames.groupby(data, :name)
    for simulation in iter_df
        all_alleles = getindex.(collect.(zip.(p1,p2)), 1)
        offspring_geno = Tuple(Base.Iterators.flatten(rand.(all_alleles, 1)) |> collect)
        miss_idx = reduce(union, findall.(i -> i == (missing,), offspring_geno))
        simulation.genotype[Not(miss_idx)] .= offspring_geno[Not(miss_idx)]
    end
    return data
end

function polyploid_cross!(data::DataFrame, p1::T, p2::T; n::Int, ploidy::Int) where T <: GenoArray
    iter_df = DataFrames.groupby(data, :name)
    n_alleles = ploidy ÷ 2
    for simulation in iter_df
        p1_contrib = Base.Iterators.flatten(sample_genotype.(p1, n_alleles)) |> collect
        p2_contrib = Base.Iterators.flatten(sample_genotype.(p2, n_alleles)) |> collect
        miss_idx = reduce(union, findall.(ismissing, [p1_contrib, p2_contrib]))
        offspring_geno = zip(p1_contrib, p2_contrib) |> collect
        simulation.genotype[Not(miss_idx)] .= offspring_geno[Not(miss_idx)]
    end
    return data
end

"""
    cross(data::PopData, parent1::String, parent2::String; n::Int = 100, generation::String = "F1")

Simulate a breeding cross between individuals `parent1` and `parent2` from a `PopData` object.
Returns PopData consisting of `n` offspring resulting from the cross.
#### Keyword Arguments
- `n` : Integer of number of offspring to generate (default: `100`)
- `generation` : A string to assign `population` identity to the offspring (default: `"F1"`)
"""
function cross(data::PopData, parent1::String, parent2::String; n::Int = 100, generation::String = "F1")
    # check for presence of parents
    err = ""
    err *= parent1 ∉ (@view data.meta[!, :name]) ? "$parent1 " : ""
    err *= parent2 ∉ (@view data.meta[!, :name]) ? "$parent2" : ""
    err != "" && error("One or more parents not found in PopData: $err")

    # get parental genotypes
    p1 = get_genotypes(data, parent1)
    p2 = get_genotypes(data, parent2)

    # check for parents not having mixed ploidy
    length(unique(length.(skipmissing(p1)))) != 1 && error("Parent $parent1 has mixed ploidy, which is unsupported")
    length(unique(length.(skipmissing(p2)))) != 1 && error("Parent $parent2 has mixed ploidy, which is unsupported")

    # Get the ploidy value & check for equal ploidy
    p1_ploidy = length.(skipmissing(p1)) |> first
    p2_ploidy = length.(skipmissing(p2)) |> first
    p1_ploidy != p2_ploidy && error("Parents must have identical ploidy. Parent1 = $p1_ploidy | Parent2 = $p2_ploidy")

    loci = unique(data.loci.locus)

    # pre-allocate all output information
    out_loci_names = fill(loci, n) |> Base.Iterators.flatten |> collect
    #parents = Vector{Tuple{String, String}}(undef, n)

    out_offspring = fill.(["$generation" * "_offspring_$i" for i in 001:n], length(loci)) |> Base.Iterators.flatten |> collect
    out_population = fill(generation, n * length(loci))
    out_geno = similar(p1, n * length(loci))
    out_loci = DataFrame(:name => out_offspring, :population => out_population, :locus => out_loci_names, :genotype => out_geno)
    #categorical!(out_loci, [:name, :population, :locus], compress = true)
    out_loci.name = PooledArray(out_loci.name)
    out_loci.population = PooledArray(out_loci.population)
    out_loci.locus = PooledArray(out_loci.locus)

    # perform the cross
    if p1_ploidy == 1
        haploid_cross!(data, parent1, parent2, n = n)
    elseif p1_ploidy ∈  [2, 4, 6, 8]
        polyploid_cross!(out_loci, p1, p2, n = n, ploidy = p1_ploidy)
    else
        error("Currently supported ploidy: 1, 2, 4, 6, 8")
    end
    out_meta = DataFrame(
        :name => unique(out_loci.name),
        :ploidy => fill(p1_ploidy, n),
        :population => fill(generation, n),
        :latitude => Vector{Union{Missing, Float32}}(undef, n),
        :longitude => Vector{Union{Missing, Float32}}(undef, n),
        :parents => fill((parent1,parent2), n)
    )
    PopData(out_meta, out_loci)
end

"""
    cross(parent_1::Pair, parent_2::Pair, n::Int = 100, generation::String = "F1")

Simulate a breeding cross between individuals `parent` and `parent2` from two different `PopData` objects.
Returns PopData consisting of `n` offspring resulting from the cross. `parent_1_data` and `parent_2_data`
are positional arguments, therefore they must be written without keywords and in the order of parents 1, parent 2.
#### Keyword Arguments
- `parent_1` : Pair of `PopData => "Parent1Name"`
- `parent_2` : Pair of `PopData => "Parent1Name"`
- `n` : Integer of number of offspring to generate (default: `100`)
- `generation` : A string to assign `population` identity to the offspring (default: `"F1"`)
"""
function cross(parent_1::Pair, parent_2::Pair; n::Int = 100, generation::String = "F1")
    parent_1_data = parent_1.first
    parent_2_data = parent_2.first

    parent1 = parent_1.second
    parent2 = parent_2.second

    # check for presence of parents
    parent1 ∉ (@view parent_1_data.meta[!, :name]) && error("$parent1 not found in PopData")
    parent2 ∉ (@view parent_2_data.meta[!, :name]) && error("$parent2 not found in PopData")

    # get parental genotypes
    p1 = get_genotypes(parent_1_data, parent1)
    p2 = get_genotypes(parent_2_data, parent2)

    # check for parents not having mixed ploidy
    length(unique(length.(skipmissing(p1)))) != 1 && error("Parent $parent1 has mixed ploidy, which is unsupported")
    length(unique(length.(skipmissing(p2)))) != 1 && error("Parent $parent2 has mixed ploidy, which is unsupported")

    # Get the ploidy value & check for equal ploidy
    p1_ploidy = length.(skipmissing(p1)) |> first
    p2_ploidy = length.(skipmissing(p2)) |> first
    p1_ploidy != p2_ploidy && error("Parents must have identical ploidy. Parent1 = $p1_ploidy | Parent2 = $p2_ploidy")

    # verify identical loci
    loci = unique(parent_1_data.loci.locus)
    loci_p2 = unique(parent_2_data.loci.locus)
    length(loci) != length(loci_p2) && error("Both parents must have the same number of loci. $parent1 : $length(loci) | $parent2 : $length(loci_p2")
    loci_check = loci .!= loci_p2
    culprits_p1 = loci[loci_check]
    culprits_p2 = loci_p2[loci_check]
    culp_print = "Parent 1\tParent 2" * "\n---------\t---------\n" * join("$i\t$j\n" for (i,j) in zip(culprits_p1, culprits_p2))
    length(culprits_p1) > 0 && error("Both datasets must have loci in the same order. Loci causing this error:\n" * culp_print)

    # pre-allocate all output information
    out_loci_names = fill(loci, n) |> Base.Iterators.flatten |> collect
    #parents = Vector{Tuple{String, String}}(undef, n)

    out_offspring = fill.(["$generation" * "_offspring_$i" for i in 001:n], length(loci)) |> Base.Iterators.flatten |> collect
    out_population = fill(generation, n * length(loci))
    out_geno = similar(p1, n * length(loci))
    out_loci = DataFrame(:name => out_offspring, :population => out_population, :locus => out_loci_names, :genotype => out_geno)
    #categorical!(out_loci, [:name, :population, :locus], compress = true)
    out_loci.name = PooledArray(out_loci.name)
    out_loci.population = PooledArray(out_loci.population)
    out_loci.locus = PooledArray(out_loci.locus)

    # perform the cross
    if p1_ploidy == 1
        haploid_cross!(data, parent1, parent2, n = n)
    elseif p1_ploidy ∈  [2, 4, 6, 8]
        polyploid_cross!(out_loci, p1, p2, n = n, ploidy = p1_ploidy)
    else
        error("Currently supported ploidy: 1, 2, 4, 6, 8")
    end
    out_meta = DataFrame(
        :name => unique(out_loci.name),
        :ploidy => fill(p1_ploidy, n),
        :population => fill(generation, n),
        :latitude => Vector{Union{Missing, Float32}}(undef, n),
        :longitude => Vector{Union{Missing, Float32}}(undef, n),
        :parents => fill((parent1,parent2), n)
    )
    PopData(out_meta, out_loci)
end


text:
using Code
using Documenter

DocMeta.setdocmeta!(Code, :DocTestSetup, :(using Code); recursive=true)

makedocs(;
    modules=[Code],
    authors="Alfonso Landeros <alanderos@ucla.edu> and contributors",
    repo="https://github.com/alanderos91/Code.jl/blob/{commit}{path}#{line}",
    sitename="Code.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
@testset "Path Sampler" begin
    @test_throws AssertionError StructDualDynProg.SDDP._samplepaths(4, [0.8, 0.3], false, false)
    @test_throws AssertionError StructDualDynProg.SDDP._samplepaths(4, [0.5, 0.3], true, false)
    pmf = [0.5, 0.5]
    @test [2, 2] == StructDualDynProg.SDDP._samplepaths(4, pmf, true, false)
    @test [0.5, 0.5] == pmf # It shouldn't have been modified
    @test [2, 3] == sort(StructDualDynProg.SDDP._samplepaths(5, pmf, true, false))
end


text:
using SQLite, LightXML, JSON, CSV, DataFrames, Dates

# load the SQLite DB. eventually the name will be configurable
db = SQLite.DB("transcode.sqlite")
# defines the tables in SQLite. See readme for details.
function make_tables()
	tables = """
	CREATE TABLE IF NOT EXISTS Transcodes (
	id INTEGER PRIMARY KEY,
	ingest_id INTEGER NOT NULL,
	output_id INTEGER NOT NULL,
	transcode_date INTEGER,
	UNIQUE(ingest_id, output_id)
	);"""
	ingests = """
	CREATE TABLE IF NOT EXISTS Ingests (
	id INTEGER PRIMARY KEY,
	transcode_id INTEGER,
	title TEXT NOT NULL,
	authors TEXT,
	event TEXT NOT NULL,
	video_file TEXT NOT NULL,
	sub_file TEXT NOT NULL,
	metadata_file TEXT NOT NULL,
	UNIQUE(video_file)
	);"""
	outputs = """
	CREATE TABLE IF NOT EXISTS Outputs (
	id INTEGER PRIMARY KEY,
	transcode_id INTEGER,
	title TEXT NOT NULL,
	authors TEXT NOT NULL,
	track TEXT NOT NULL,
	room TEXT NOT NULL,
	event_id TEXT NOT NULL,
	base_path TEXT NOT NULL,
	UNIQUE(event_id)
	);
	"""
	metadata = """
	CREATE TABLE IF NOT EXISTS Metadata (
	id INTEGER PRIMARY KEY,
	ingest_id INTEGER,
	description TEXT,
	social_tags TEXT,
	social_handles TEXT,
	summary TEXT,
	preprint TEXT,
	doi TEXT
	);
	"""
	DBInterface.execute(db, tables)
	DBInterface.execute(db, ingests)
	DBInterface.execute(db, outputs)
	DBInterface.execute(db, metadata)
end
make_tables()

function load_conf_pub_dir(db, dir, ingest_dir)
	modified = (0,0)
	for file in readdir(dir)
		extension = last(split(file, "."))
		if extension == "xml"
			modified = modified .+ load_conf_pub(db, joinpath(dir, file), ingest_dir)
		end
	end
	return modified
end

function load_manual_csv(db, csv)
	function check_file_exists(file, ftyp)
		if !Base.Filesystem.ispath(file)
			println("$ftyp file does not exist at $file")
			return false
		end
		return true
	end
	data = CSV.File(csv) |> DataFrame
	insert_ingest = SQLite.Stmt(db, "INSERT OR REPLACE INTO Ingests (title, authors, event, video_file, sub_file, metadata_file) VALUES (:title, :authors, :event, :video, :subs, :metadata)")
	lookup_inserted = SQLite.Stmt(db, "SELECT id FROM Ingests WHERE title = :title")
	transcode_exists = SQLite.Stmt(db, "SELECT id FROM Transcodes WHERE ingest_id=:iid AND output_id = :oid")
	insert_transcode = SQLite.Stmt(db, "INSERT OR IGNORE INTO Transcodes (ingest_id, output_id) VALUES (:iid, :oid)")
	lookup_binding = SQLite.Stmt(db, "SELECT id FROM Outputs WHERE event_id = :eid")
	for (eid, title, track, authors, video_file, sub_file, metadata_file) in eachrow(data)
		if ismissing(video_file) || isempty(strip(video_file))
			println("event $title does not have a video file, continuing")
			continue
		end
		if (!check_file_exists(sub_file, "subtitle") || !check_file_exists(metadata_file, "metadata"))
			continue
		end
		res = DataFrame(DBInterface.execute(lookup_binding, (eid=eid, )))
		if isempty(res)
			println("Invalid event id $eid, continuing")
			continue
		end
		output_id = first(res["id"])
		if isempty(DBInterface.execute(lookup_inserted, (title=title, )))
			DBInterface.execute(insert_ingest, (title=title, authors=authors, event=track, video=video_file, subs=sub_file, metadata=metadata_file))
		end
		inserted = DataFrame(DBInterface.execute(lookup_inserted, (title=title, )))
		ingest_id = first(inserted["id"])
		if isempty(DBInterface.execute(transcode_exists, (iid=ingest_id, oid=output_id)))
			DBInterface.execute(insert_transcode, (iid=ingest_id, oid=output_id))
		end
	end
end

function load_conf_pub(db, pub_xml, ingest_dir)
	existsq = SQLite.Stmt(db, "SELECT 1 FROM Ingests WHERE video_file = :video AND event = :event")
	insert = SQLite.Stmt(db, "INSERT OR IGNORE INTO Ingests (title, authors, event, video_file, sub_file, metadata_file) VALUES (:title, :authors, :event, :video, :subs, :metadata)")
	xdoc = nothing
	try
		xdoc = parse_file(pub_xml)
	catch e
		error("Could not parse conference publishing XML file")
		rethrow(e)
	end
	DBInterface.execute(db, "BEGIN TRANSACTION;")
	inserted = 0
	skipped = 0
	try
		function add_definition(article, event_name)
			confpub_id = attribute(article, "id")
			title = content(first(article["title"]))
			authors = join(map(author_xml -> "$(content(first(author_xml["givennames"]))) $(content(first(author_xml["surname"])))", first(article["authors"])["author"]), ", ")
			video_file = "$ingest_dir/$confpub_id-Video.mp4"
			sub_file = "$ingest_dir/$confpub_id-Video.srt"
			metadata_file = "$ingest_dir/$confpub_id-Video.json"
			if isempty(DBInterface.execute(existsq, (video=video_file, event=event_name)))
				DBInterface.execute(insert, (title=title, authors = authors, event=event_name, video=video_file, subs=sub_file, metadata=metadata_file))
				inserted += 1
			else
				println(article)
				skipped += 1
			end
		end
		xroot = root(xdoc)
		event_name = content(first(xroot["eventmain"]))
		for eventsub in xroot["eventsub"]
			for session in eventsub["session"]
				for article in session["article"]
					add_definition(article, event_name)
				end
			end
			for track in eventsub["track"]
				for article in track["article"]
					add_definition(article, event_name)
				end
			end
		end
		DBInterface.execute(db, "COMMIT;")
	catch e
		DBInterface.execute(db, "ROLLBACK;")
		rethrow(e)
	end
	return inserted, skipped
end

function process_scheduler(db, required, base_path)
	data = JSON.parsefile(required)
	existsq = SQLite.Stmt(db, "SELECT 1 FROM Outputs WHERE event_id = :eid")
	insert = SQLite.Stmt(db, "INSERT OR IGNORE INTO Outputs (title, track, authors, room, event_id, base_path) VALUES (:title, :track, :authors, :room, :eventid, :basepath)")
	DBInterface.execute(db, "BEGIN TRANSACTION;")
	inserted = 0
	skipped = 0
	try
		for req in data
			if isempty(DBInterface.execute(existsq, (eid=req["id"], )))
				DBInterface.execute(insert, (title=req["title"], track=req["track"], room=req["room"], authors=req["authors"], eventid=req["id"], basepath=base_path))
				inserted += 1
			else
				skipped += 1
			end
		end
		DBInterface.execute(db, "COMMIT;")
	catch e
		DBInterface.execute(db, "ROLLBACK;")
		rethrow(e)
	end
	return inserted, skipped
end

function load_metadata(db, pub_xml, ingest_dir)
	check = SQLite.Stmt(db, "SELECT id FROM Metadata WHERE ingest_id = :iid")
	lookup_id = SQLite.Stmt(db, "SELECT id FROM Ingests WHERE title = :title")
	insert = SQLite.Stmt(db, """
		INSERT INTO Metadata (ingest_id, title, abstract, url, urlinfo, event, paper_keywords, description, social_tags, social_handles, summary, doi)
		VALUES (:iid, :title, :abst, :url, :urlinfo, :event, :paper_keywords, :description, :social_tags, :social_handles, :summary, :doi)""")
	update = SQLite.Stmt(db, """
		UPDATE Metadata SET
			abstract=:abst,
			url=:url, urlinfo=:urlinfo, event=:event,
			paper_keywords=:paper_keywords, description=:description, social_tags=:social_tags,
			social_handles=:social_handles, summary=:summary, doi=:doi
		WHERE ingest_id = :iid""")
	xdoc = nothing
	try
		xdoc = parse_file(pub_xml)
	catch e
		error("Could not parse conference publishing XML file")
		rethrow(e)
	end
	function add_definition(article, event_name)
		confpub_id = attribute(article, "id")
		if !Base.Filesystem.ispath(joinpath(ingest_dir, "$confpub_id-Video.json"))
			return
		end
		title = content(first(article["title"]))
		abstr = content(first(article["abstract"]))
		url = content(first(article["url"]))
		urlinfo = content(first(article["urlinfo"]))
		doi = content(first(article["doi"]))
		paper_keywords = content(first(article["keywords"]))

		metadata_file = joinpath(ingest_dir, "$confpub_id-Video.json")
		metadata_json = JSON.parsefile(metadata_file)
		desc = metadata_json["submissionmaterialdesc"]
		tags = if haskey(metadata_json, "submissionmaterialsocialmediatags") metadata_json["submissionmaterialsocialmediatags"] else "" end
		handles = if haskey(metadata_json, "submissionmaterialsocialmediahandles") metadata_json["submissionmaterialsocialmediahandles"] else "" end
		summary = if haskey(metadata_json, "submissionmaterialtwittersummary") metadata_json["submissionmaterialtwittersummary"] else "" end

		(igid, ) = first(DBInterface.execute(lookup_id, (title=title, )))
		if isempty(DBInterface.execute(check, (iid=igid, )))
			DBInterface.execute(insert, (iid=igid,title=title,abst=abstr,url=url,urlinfo=urlinfo,event=event_name,paper_keywords=paper_keywords,description=desc,social_tags=tags,social_handles=handles,summary=summary, doi=doi))
		else
			DBInterface.execute(update, (iid=igid,title=title,abst=abstr,url=url,urlinfo=urlinfo,event=event_name,paper_keywords=paper_keywords,description=desc,social_tags=tags,social_handles=handles,summary=summary, doi=doi))
		end
	end
	xroot = root(xdoc)
	event_name = content(first(xroot["eventmain"]))
	for eventsub in xroot["eventsub"]
		for session in eventsub["session"]
			for article in session["article"]
				add_definition(article, event_name)
			end
		end
		for track in eventsub["track"]
			for article in track["article"]
				add_definition(article, event_name)
			end
		end
	end
end

function csv_ingest(db, csv)
	df = CSV.File(csv) |> DataFrame
	return df
end

function title_bind(db)
	add_mappingsl = SQLite.Stmt(db, """
	   INSERT INTO Transcodes (output_id, ingest_id) SELECT Outputs.id, Ingests.id FROM Outputs
       LEFT JOIN Ingests ON lower(Outputs.title) LIKE "%" || lower(Ingests.title) || "%"
       LEFT JOIN Transcodes ON Outputs.id = Transcodes.output_id
       WHERE Transcodes.id IS NULL AND Ingests.id NOT NULL""")
	DBInterface.execute(add_mappingsl)
	add_mappingsr = SQLite.Stmt(db, """
	   INSERT INTO Transcodes (output_id, ingest_id) SELECT Outputs.id, Ingests.id FROM Outputs
       LEFT JOIN Ingests ON lower(Ingests.title) LIKE "%" || lower(Outputs.title) || "%"
       LEFT JOIN Transcodes ON Outputs.id = Transcodes.output_id
       WHERE Transcodes.id IS NULL AND Ingests.id NOT NULL""")
	DBInterface.execute(add_mappingsr)
	authors_match = SQLite.Stmt(db, """
	INSERT INTO Transcodes (output_id, ingest_id) SELECT oid, iid FROM (SELECT Outputs.id as oid, Ingests.id as iid, COUNT(DISTINCT(Ingests.id)) AS count FROM Outputs
      LEFT JOIN Transcodes ON Transcodes.output_id = Outputs.id
      LEFT JOIN Ingests ON Ingests.authors LIKE Outputs.authors
      WHERE Ingests.id IS NOT NULL
      GROUP BY Ingests.id) LEFT JOIN Transcodes ON Transcodes.output_id = oid
      WHERE Transcodes.id IS NULL AND count = 1
	""")
	DBInterface.execute(authors_match)
end

function bind_events(db, title_output, title_ingest)
	find_output = SQLite.Stmt(db, "SELECT id FROM Outputs where Outputs.title = :title")
	find_ingest = SQLite.Stmt(db, "SELECT id FROM Ingests where Ingests.title = :title")
	insert = SQLite.Stmt(db, "INSERT INTO Transcodes (output_id, ingest_id) VALUES (:output_id, :ingest_id)")
	out = DataFrame(DBInterface.execute(find_output, (title=title_output, )))
	if isempty(out)
		throw("Could not find output $(title_output)")
	end

	ing = DataFrame(DBInterface.execute(find_ingest, (title=title_ingest, )))
	if isempty(ing)
		throw("Could not find input $(title_ingest)")
	end

	DBInterface.execute(insert, (output_id=first(out["id"]), ingest_id=first(ing["id"])))
end

function get_fps(file)
	ffdata = read(`ffprobe -loglevel 8 -print_format json -show_streams $file`, String)
	data = JSON.parse(ffdata)
	for stream in data["streams"]
		if !haskey(stream, "codec_type") || stream["codec_type"] != "video"
			continue
		end
		rfr = stream["avg_frame_rate"]
		num,denom = parse.(Int, split(rfr, "/"))
		return num/denom
	end
end

function get_duration(file)
	ffdata = read(`ffprobe -loglevel 8 -print_format json -show_streams $file`, String)
	data = JSON.parse(ffdata)
	for stream in data["streams"]
		if !haskey(stream, "codec_type") || stream["codec_type"] != "video"
			continue
		end
		rfr = stream["duration"]
		return parse.(Float64, rfr)
	end
end

function normalization_firstpass(video)
	normalize_task = `ffmpeg -y -i $video -max_muxing_queue_size 400 -af loudnorm=I=-15:LRA=9:tp=-1:print_format=json -f null -`
	out = Pipe()
	err = Pipe()
	process = run(pipeline(normalize_task, stdout=out, stderr=err))
	close(out.in)
	close(err.in)
	read_data = String(read(err))
	jsondata = join(split(read_data, "\n")[end-12:end], "\n")
	loudness_data = JSON.parse(jsondata)
	return loudness_data
end

function process_outstanding_videos(db)
	# build the workqueue
	transcodes = DataFrame(DBInterface.execute(db, """
		SELECT video_file, sub_file, metadata_file, event_id, base_path, transcode_date, ingest_id, output_id from "Transcodes"
		INNER JOIN Ingests ON Ingests.id = Transcodes.ingest_id
		INNER JOIN Outputs ON Outputs.id = Transcodes.output_id"""))
	done_stmt = SQLite.Stmt(db, "UPDATE Transcodes SET transcode_date=:date WHERE output_id=:oid AND ingest_id=:iid")
	for (inp_vid, inp_sub, inp_metadata, eventid, base_path, transcode_date, iid, oid) in eachrow(transcodes)
		output_video = joinpath(base_path, "$eventid.mp4")
		if !Base.Filesystem.ispath(base_path)
			Base.Filesystem.mkpath(base_path)
		end
		if !Base.Filesystem.ispath(inp_vid)
			println("Missing input video $inp_vid, continuing.")
			continue
		end
		inp_mod = Base.Filesystem.mtime(inp_vid)
		if !ismissing(transcode_date) && inp_mod <= transcode_date
			continue
		end
		computed_fps = get_fps(inp_vid)
		fps_cmd = ``
		if computed_fps < 23 || computed_fps > 35
			computed_fps = 30
			fps_cmd = `-r 30`
		end
		gop_size = floor(Int, computed_fps*2)
		loudness_data = normalization_firstpass(inp_vid)

		transcode_task = `ffmpeg -y -i $inp_vid $fps_cmd -vf scale=1920:1080 -pix_fmt yuv420p -threads 0 -vcodec libx264 -g $gop_size -sc_threshold 0 -b:v 3000k -max_muxing_queue_size 400
								-bufsize 1216k -maxrate 6000k -preset medium -profile:v high -tune film
								-acodec aac -b:a 128k -ac 2 -ar 44100
								-af "loudnorm=I=-15:LRA=9:tp=-1:measured_I=$(loudness_data["input_i"]):measured_LRA=$(loudness_data["input_lra"]):measured_tp=$(loudness_data["input_tp"]):offset=$(loudness_data["target_offset"]),aresample=async=1:min_hard_comp=0.100000:first_pts=0" $output_video`
		proc = run(transcode_task)
		wait(proc)
		if proc.exitcode != 0
			println("Errored while processing video $inp_vid with ID $eventid")
			continue
		end
		DBInterface.execute(done_stmt, (date=floor(Int, datetime2unix(Dates.now())), oid=oid, iid=iid))
	end
end

function reset_transcodes(db)
	DBInterface.execute(db, "UPDATE Transcodes SET transcode_date=NULL")
end

function copy_metadata(db)
	# build the workqueue
	transcodes = DataFrame(DBInterface.execute(db, """
		SELECT video_file, sub_file, metadata_file, event_id, base_path, transcode_date, ingest_id, output_id from "Transcodes"
		INNER JOIN Ingests ON Ingests.id = Transcodes.ingest_id
		INNER JOIN Outputs ON Outputs.id = Transcodes.output_id"""))
	for (inp_vid, inp_sub, inp_metadata, eventid, base_path, transcode_date, iid, oid) in eachrow(transcodes)
		output_subs = joinpath(base_path, "$eventid.srt")
		output_meta = joinpath(base_path, "$eventid.json")
		if !Base.Filesystem.ispath(inp_sub)
			println("Missing $inp_sub")
			continue
		end
		if !Base.Filesystem.ispath(inp_metadata)
			println("Missing $inp_metadata")
			continue
		end
		Base.Filesystem.cp(inp_sub, output_subs, force=true)
		Base.Filesystem.cp(inp_metadata, output_meta, force=true)
	end
end

function check_conference(db)
	#first, identify outputs that haven't been matched to inputs
	unmatched = DataFrame(DBInterface.execute(db, """
		SELECT Outputs.title, Outputs.authors, Outputs.track, Outputs.room from Outputs
		LEFT JOIN Transcodes ON Transcodes.output_id = Outputs.id
		WHERE Transcodes.id IS NULL"""))
	if !isempty(unmatched)
		println("Unmatched required outputs!")
		for (title, authors, track, room) in eachrow(unmatched)
			println("Unmatched: $title by $authors at $track in $room")
		end
	end
	# next, identify inputs which are mapped to outputs that don't correspond to a video
	required_inp_files = DataFrame(DBInterface.execute(db, """
		SELECT Ingests.title, Ingests.event, Ingests.video_file, Ingests.sub_file, Ingests.metadata_file from Ingests
		INNER JOIN Transcodes ON Transcodes.ingest_id = Ingests.id"""))
	for (title, event, vid, sub, meta) in eachrow(required_inp_files)
		missing_els = []
		hasmissing = false
		if !Base.Filesystem.ispath(vid)
			push!(missing_els, "video")
			hasmissing = true
		end
		if !Base.Filesystem.ispath(sub)
			push!(missing_els, "subs")
			hasmissing = true
		end
		if !Base.Filesystem.ispath(meta)
			push!(missing_els, "meta")
			hasmissing = true
		end
		if hasmissing
			println("Missing $(join(missing_els, ", ")) for $title at $event w/ name $vid")
		end
	end
	# check to see if there are any outstanding transcodes
	transcodes_to_do = DataFrame(DBInterface.execute(db, """
		SELECT Transcodes.transcode_date, Ingests.video_file from Transcodes
		INNER JOIN Ingests ON Transcodes.ingest_id = Ingests.id"""))
	for (tcd, ingest_file) in eachrow(transcodes_to_do)
		if Base.Filesystem.ispath(ingest_file) && (ismissing(tcd) || tcd < Base.Filesystem.mtime(ingest_file))
			println("Need to do transcode $ingest_file")
		end
	end
end

using DataFramesMeta

function schedule_status(db)
	events = DataFrame(DBInterface.execute(db, """
		SELECT Outputs.title, Outputs.event_id, Outputs.base_path, Transcodes.id FROM Outputs
		LEFT JOIN Transcodes ON Transcodes.output_id = Outputs.id"""))
	file_exists(base_path, id) = Base.Filesystem.ispath(joinpath(base_path, "$id.mp4"))
	get_duration_l(base_path, id) = if file_exists(base_path, id) get_duration(joinpath(base_path, "$id.mp4")) else -1 end
	CSV.write("status.csv", @linq events |> transform(has_file=file_exists.(:base_path, :event_id), (duration=get_duration_l.(:base_path, :event_id))))
end

function write_sv_csv(db)
	tovalidate = DataFrame(DBInterface.execute(db, """
		SELECT Outputs.title, Outputs.authors, Outputs.track, Outputs.event_id, Outputs.base_path, Transcodes.id FROM Outputs
		LEFT JOIN Transcodes ON Transcodes.output_id = Outputs.id"""))
	make_path(base_path, id) = return "rtmp://3.18.146.46/vod/mp4:store/$base_path/$id.mp4"
	file_exists(base_path, id) = Base.Filesystem.ispath(joinpath(base_path, "$id.mp4"))
	towrite = @linq tovalidate |> transform(url=make_path.(:base_path, :event_id)) |> where(file_exists.(:base_path, :event_id))
	CSV.write("sv_csv.csv", towrite)
end

#=
if length(ARGS) == 0
	action = "default"
else
	action = ARGS[1]
end

if action == "metadata"
	process_metadata(ARGS[2], ARGS[3], ARGS[4], ARGS[5])
elseif action == "sidechain"
	add_sidechain(ARGS[2], ARGS[3], ARGS[4])
else
	error("Invalid action $action; valid values are metadata and sidechain")
end
=#

max_stars_repo_name: SPLASH-infrastructure/batch-transcoder

max_issues_repo_name: SPLASH-infrastructure/batch-transcoder

max_forks_repo_name: SPLASH-infrastructure/batch-transcoder


text:
# # Building the BIOCLIM model

# **Justification for this use case:** `SImpleSDMLayers` can be used as a
# platform to build your own species distribution models. In this example, which
# assumes that you have read the vignettes on GBIF integration and variable
# selection through VIF, we will build our own version of the BIOCLIM model, and
# apply it to the distribution of *Hypomyces lactifluorum* in North America.

using SimpleSDMLayers
using GBIF
using Plots
using GLM
using StatsBase
using Statistics
using GeometryBasics

# BIOCLIM is a very simple model, which only requires presence information. The
# first step is therefore to get occurrences of *Hypomyces lactifluorum* in
# North America. Because the data in GBIF is only as good as the original data
# source, sometimes searching by `continent` gives fewer results than searching
# by country.

observations = occurrences(
    GBIF.taxon("Hypomyces lactifluorum"; strict=true),
    "hasCoordinate" => "true",
    "country" => "CA",
    "country" => "US",
    "limit" => 300,
)

# We will now page through additional results (300 at a time).

while length(observations) < size(observations)
    occurrences!(observations)
end

# At this point, we could read the whole predictor variables directly, and then
# clip them. This would be fairly wasteful, as we need a small area. For this
# reason, we will calculate the bounding box first, and then use it to only read
# the section we want.

left, right = extrema(longitudes(observations)) .+ (-5, 5)
bottom, top = extrema(latitudes(observations)) .+ (-5, 5)
boundaries = (left=left, right=right, bottom=bottom, top=top)

# With this information in hand, we can start getting our variables. In this
# example, we will take all of the BioClim data from WorldClim, at the 10 arc
# minute resolution, and add the elevation layer. Note that using the bounding
# box coordinates when calling the layers is *much* faster than clipping after
# the fact (assuming that you already have the files downloaded).

predictors =
    convert.(
        Float32, SimpleSDMPredictor(WorldClim, BioClim, 1:19; resolution=10.0, boundaries...)
    );

# We will add the elevation to the stack of variables we use -- we need to
# convert everything to `Float32` layers, because elevation is originally an
# `Int16` one and a number of operations we will make will require floating
# points

push!(
    predictors,
    convert(
        Float32, SimpleSDMPredictor(WorldClim, Elevation; resolution=10.0, boundaries...)
    ),
);

# It is not a bad idea to plot all of the predictors:

plot(plot.(predictors, grid=:none, axes=false, frame=:none, leg=false, c=:imola)...)

# Clearly, some of them show strong autocorrelation; we will therefore re-use
# our VIF code to select a subset that has uncorrelated variables.

function vif(model)
    R² = r2(model)
    return 1 / (1-R²)
end

function stepwisevif(
    layers::Vector{T}, selection=collect(1:length(layers)), threshold::Float64=5.0
) where {T<:SimpleSDMLayer}
    x = hcat([layer[keys(layer)] for layer in layers[selection]]...)
    X = (x .- mean(x; dims=1)) ./ std(x; dims=1)
    vifs = zeros(Float64, length(selection))
    for i in eachindex(selection)
        vifs[i] = vif(lm(X[:, setdiff(eachindex(selection), i)], X[:, i]))
    end
    all(vifs .<= threshold) && return selection
    drop = last(findmax(vifs))
    popat!(selection, drop)
    @info "Variables remaining: $(selection)"
    return stepwisevif(layers, selection, threshold)
end

# We will apply this function with the default parameters:

layers_to_keep = stepwisevif(predictors)

# When this is done, we can plot  the layers again to check that they are all
# more or less unique:

plot(
    plot.(
        predictors[layers_to_keep], grid=:none, axes=false, frame=:none, leg=false, c=:imola
    )...,
)

# The point of BIOCLIM (the model, not the dataset) is that the score assigned
# to a pixel is maximal if this pixel is the *median* value for a given
# variable. Therefore, we need to measure the cumulative density function for
# every pixel in every variable, and transform it with:

_pixel_score(x) = 2.0(x > 0.5 ? 1.0 - x : x);

# The *actual* model generation is fairly straightforward, as we will need to
# get the values of the layers in the cells occupied by an observation. Because
# sampling bias is very real, we will grid the observations by transforming them
# into a boolean layer:

presences = mask(predictors[1], observations, Bool)
plot(convert(Float32, presences); c=cgrad([:lightgrey, :black]), leg=false)

# This step is very important so as not to bias the estimation of quantiles,
# which overcounting observations within the same cell would do. We can now
# define the model:

function SDM(predictor::T1, observations::T2) where {T1<:SimpleSDMLayer,T2<:SimpleSDMLayer}
    _tmp = mask(observations, predictor)
    qf = ecdf(convert(Vector{Float32}, _tmp[keys(_tmp)])) # We only want the observed values
    return (_pixel_score ∘ qf)
end

# Note that we use the ∘ (`\circ`) operator to chain the quantile estimation and
# the pixel scoring, which requires Julia 1.4. This function returns a *model*,
# *i.e.* a function that we can broadcast to a given layer, which might not be
# the same one we used for the training.

# The next step in BIOCLIM is to get the *minimum* suitability across all layers
# for every pixel. Because we have a `min` method defined for a pair of layers,
# we can call `minimum` on an array of layers:

function SDM(predictors::Vector{T}, models) where {T<:SimpleSDMLayer}
    @assert length(models) == length(predictors)
    return minimum([broadcast(models[i], predictors[i]) for i in 1:length(predictors)])
end

# The advantage of this approach is that we can call the `SDM` method for
# prediction on a smaller layer, or a different layer. This can allow us to do
# thing like stitching layers together with `hcat` and `vcat` to use
# multi-threading, or use a different resolution for the prediction than we did
# for the training.

models = [SDM(predictor, presences) for predictor in predictors];

# We now get the prediction:

prediction = SDM(predictors, models)

# It's not a bad idea to look at this prediction, to get a sense of where the
# hotspots of presence would be:

plot(prediction; c=:bamako, frame=:box)
xaxis!("Longitude")
yaxis!("Latitude")

# Just because we may want to visualize this result in a transformed way, *i.e.*
# by looking at the quantiles of suitability, we can call the `rescale`
# function:

prediction_quantile = rescale(prediction, collect(0.0:0.005:1.0))

# As this map now represents the quantiles of suitability, we may want to remove
# the lower 5%. For this, we need to create a boolean mask, which we can do by
# broadcasting a conditional:

cutoff = broadcast(x -> x > 0.05, prediction_quantile)

# The raw prediction, minus the 5% bottom quantiles, can the be plotted:

plot(prediction; frame=:box, c=:lightgrey)
plot!(mask(cutoff, prediction); c=:bamako)
xaxis!("Longitude")
yaxis!("Latitude")

# And there it is! A simple way to write the BIOCLIM model by building on the
# integration between SimpleSDMLayers and GBIF. BIOCLIM has a tendency to
# underfit the distribution quite a bit - in fact, the range returned here is
# larger than what other methods (like BRT) would return. The next vignettes in
# this section will be focused on using additional functionalities of
# `SimpleSDMLayers` until we are able to make a better model for this species.


text:
#! format: off
import QuantumControl
import Documenter

"""Return a list of symbols for the names directly defined in `pkg`.

This filters out re-exported names and sub-modules. By default, for `all=true`,
both public (exported) and private names are included. With `all=false`, the
list is filtered to include only exported names.
"""
function get_local_members(pkg; all=true)
    return [
        m for m in names(pkg, all=all) if !(
            (startswith(String(m), "#")) ||       # compiler-generated names
            (getfield(pkg, m) isa Union{Dict,Array,Set}) ||  # global variable
            (m == Symbol(pkg)) ||                 # the package itself
            (m == :eval) ||                       # compiler-injected "eval"
            (m == :include) ||                    # compiler-injected "include"
            ((getfield(pkg, m)) isa Module) ||    # sub-modules
            (_parentmodule(getfield(pkg, m), pkg)) ≠ pkg  # re-exported
        )
    ]
end

_parentmodule(m, pkg) = parentmodule(m)
_parentmodule(m::Number, pkg) = pkg

"""Return a list of symbols for all the sub-modules of `pkg`.
"""
function get_submodules(pkg)
    return [
        m for m in names(pkg, all=true)
        if (getfield(pkg, m) isa Module) && !(m == Symbol(pkg))
    ]
end

"""Return the canonical fully qualified name of an object (function, type).

Given e.g. a function object, it returns a string containing the canonical
fully qualified name of the original function definition.
"""
function canonical_name(obj)
    mod = parentmodule(obj)
    modpath = fullname(mod)
    modname = join((String(sym) for sym in modpath), ".")
    objname = String(nameof(obj))
    return "$modname.$objname"
end

quantum_control_members = [
    m for m in names(QuantumControl)
    if m ≠ :QuantumControl
]

quantum_control_local_members = get_local_members(QuantumControl, all=false)

quantum_control_reexported_members = [
    m for m in quantum_control_members
    if m ∉ quantum_control_local_members
]

quantum_control_sub_modules = get_submodules(QuantumControl)

subpackages = [
    (:QuantumPropagators, "quantum_propagators.md"),
    (:QuantumControlBase, "quantum_control_base.md"),
    (:Krotov, "krotov.md"),
    (:GRAPE, "grape.md"),
]

outfile = joinpath(@__DIR__, "src", "api", "quantum_control.md")
println("Generating API for QuantumControl in $outfile")
open(outfile, "w") do out
    write(out, "```@meta\n")
    write(out, "EditURL = \"../../generate_api.jl\"\n")
    write(out, "```\n\n")
    write(out, "# [QuantumControl](@id QuantumControlAPI)\n\n")
    if length(quantum_control_local_members) > 0
        println(out, "```@docs")
        for name ∈ quantum_control_local_members
            println(out, name)
        end
        println(out, "```")
    end
    write(out, """

    QuantumControl re-exports the following members:

    """)
    for name ∈ quantum_control_reexported_members
        obj = getfield(QuantumControl, name)
        ref = canonical_name(obj)
        println(out, "* [`$name`](@ref $ref)")
    end

    for submod in quantum_control_sub_modules
        write(out, "\n\n### [`QuantumControl.$submod`](@id QuantumControl$(submod)API)\n\n")
        for name in names(getfield(QuantumControl, submod))
            if name ≠ submod
                obj = getfield(getfield(QuantumControl, submod), name)
                ref = canonical_name(obj)
                println(out, "* [`QuantumControl.$submod.$name`](@ref $ref)")
            end
        end
    end
end

for (pkgname::Symbol, outfilename) in subpackages

    local outfile = joinpath(@__DIR__, "src", "api", outfilename)
    println("Generating API for $pkgname in $outfile")
    open(outfile, "w") do out

        pkg = getfield(QuantumControl, pkgname)
        all_local_members = get_local_members(pkg)
        public_members = get_local_members(pkg, all=false)
        documented_members = [
            k.var for k in keys(Documenter.DocSystem.getmeta(pkg))
        ]
        documented_private_members = [
            name for name in documented_members
            if (name ∉ public_members) && (name ∈ all_local_members)
        ]
        write(out, "```@meta\n")
        write(out, "EditURL = \"../../generate_api.jl\"\n")
        write(out, "```\n\n")
        write(out, "\n\n# $pkgname Package\n\n")
        write(out, "## Index\n\n")
        write(out, raw"""
        ``\gdef\tgt{\text{tgt}}``
        ``\gdef\tr{\operatorname{tr}}``
        ``\gdef\Re{\operatorname{Re}}``
        ``\gdef\Im{\operatorname{Im}}``
        """)
        write(out, """

        ```@index
        Pages   = ["$outfilename"]
        ```

        """)
        write(out, "\n\n## [$pkgname](@id $(pkgname)API)\n\n")
        if length(public_members) > 0
            write(out, "\n### Public\n\n")
            write(out, "```@docs\n")
            for name in public_members
                write(out, "$pkgname.$name\n")
            end
            write(out, "```\n\n")
        end
        if length(documented_private_members) > 0
            write(out, "\n### Private\n\n")
            write(out, "```@docs\n")
            for name in documented_private_members
                write(out, "$pkgname.$name\n")
            end
            write(out, "```\n\n")
        end

        sub_modules = get_submodules(pkg)
        for submodname in sub_modules
            submod = getfield(pkg, submodname)
            all_local_members = get_local_members(submod)
            public_members = get_local_members(submod, all=false)
            documented_members = [
                k.var for k in keys(Documenter.DocSystem.getmeta(submod))
            ]
            documented_private_members = [
                name for name in documented_members
                if (name ∉ public_members) && (name ∈ all_local_members)
            ]
            if length(public_members) + length(documented_private_members) > 0
                write(out, "\n## [$pkgname.$submodname](@id $(pkgname)$(submodname)API)\n\n")
            end
            if length(public_members) > 0
                write(out, "\n### Public\n\n")
                write(out, "```@docs\n")
                for name in public_members
                    write(out, "$pkgname.$submodname.$name\n")
                end
                write(out, "```\n\n")
            end
            if length(documented_private_members) > 0
                write(out, "\n### Private\n\n")
                write(out, "```@docs\n")
                for name in documented_private_members
                    write(out, "$pkgname.$submodname.$name\n")
                end
                write(out, "```\n\n")
            end
        end

    end

end


text:
using OptimalTransport
using Distances
using PythonOT: PythonOT

using LinearAlgebra
using Random
using Test

const POT = PythonOT

Random.seed!(100)

@testset "sinkhorn_unbalanced.jl" begin
    M = 250
    N = 200

    @testset "example" begin
        μ = fill(1 / N, M)
        ν = fill(1 / N, N)

        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

        # compute optimal transport plan
        eps = 0.01
        lambda = 1
        γ = sinkhorn_unbalanced(μ, ν, C, lambda, lambda, eps)

        # compare with POT
        γ_pot = POT.sinkhorn_unbalanced(μ, ν, C, eps, lambda; stopThr=1e-9)
        @test γ_pot ≈ γ

        # compute optimal transport cost
        c = sinkhorn_unbalanced2(μ, ν, C, lambda, lambda, eps; maxiter=5_000)

        # compare with POT
        c_pot = POT.sinkhorn_unbalanced2(
            μ, ν, C, eps, lambda; numItermax=5_000, stopThr=1e-9
        )[1]
        @test c_pot ≈ c

        # ensure that provided plan is used
        c2 = sinkhorn_unbalanced2(similar(μ), similar(ν), C, rand(), rand(), rand(); plan=γ)
        @test c2 ≈ c
    end

    @testset "proxdiv operators" begin
        μ = fill(1 / N, M)
        ν = fill(1 / N, N)
        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

        # compute optimal transport plan and cost with real-valued parameters
        eps = 0.01
        lambda1 = 0.4
        lambda2 = 0.5
        γ = sinkhorn_unbalanced(μ, ν, C, lambda1, lambda2, eps)
        c = sinkhorn_unbalanced2(μ, ν, C, lambda1, lambda2, eps)
        @test c ≈ dot(γ, C)

        # compute optimal transport plan and cost with manual "proxdiv" functions
        function proxdivF1!(s, p, ε)
            return (s .= s .^ (ε / (ε + 0.4)) .* p .^ (0.4 / (ε + 0.4)) ./ s)
        end
        function proxdivF2!(s, p, ε)
            return (s .= s .^ (ε / (ε + 0.5)) .* p .^ (0.5 / (ε + 0.5)) ./ s)
        end
        γ_proxdiv = sinkhorn_unbalanced(μ, ν, C, proxdivF1!, proxdivF2!, eps)
        c_proxdiv = sinkhorn_unbalanced2(μ, ν, C, proxdivF1!, proxdivF2!, eps)
        @test γ_proxdiv ≈ γ
        @test c_proxdiv ≈ c
    end

    @testset "consistency with balanced OT" begin
        μ = fill(1 / M, M)
        ν = fill(1 / N, N)
        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

        # compute optimal transport plan and cost with manual "proxdiv" functions for
        # balanced OT
        ε = 0.01
        proxdivF!(s, p, ε) = (s .= p ./ s)
        γ = sinkhorn_unbalanced(μ, ν, C, proxdivF!, proxdivF!, ε)
        c = sinkhorn_unbalanced2(μ, ν, C, proxdivF!, proxdivF!, ε)
        @test c ≈ dot(γ, C)

        # compute optimal transport plan and cost for balanced OT
        γ_balanced = sinkhorn(μ, ν, C, ε)
        c_balanced = sinkhorn2(μ, ν, C, ε)
        @test γ_balanced ≈ γ rtol = 1e-4
        @test c_balanced ≈ c rtol = 1e-4
    end

    @testset "deprecations" begin
        μ = fill(1 / N, M)
        ν = fill(1 / N, N)
        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

        # compute optimal transport plan and cost with real-valued parameters
        ε = 0.01
        λ1 = 0.4
        λ2 = 0.5
        γ = sinkhorn_unbalanced(μ, ν, C, λ1, λ2, ε)
        c = sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε)

        # compute optimal transport plan and cost with manual "proxdiv" functions
        # as keyword arguments
        function proxdivF1(s, p)
            return s .^ (0.01 / (0.01 + 0.4)) .* p .^ (0.4 / (0.01 + 0.4)) ./ s
        end
        function proxdivF2(s, p)
            return s .^ (0.01 / (0.01 + 0.5)) .* p .^ (0.5 / (0.01 + 0.5)) ./ s
        end
        γ_proxdiv = @test_deprecated sinkhorn_unbalanced(
            μ, ν, C, rand(), rand(), ε; proxdiv_F1=proxdivF1, proxdiv_F2=proxdivF2
        )
        c_proxdiv = @test_deprecated sinkhorn_unbalanced2(
            μ, ν, C, rand(), rand(), ε; proxdiv_F1=proxdivF1, proxdiv_F2=proxdivF2
        )
        @test γ_proxdiv ≈ γ
        @test c_proxdiv ≈ c

        # deprecated `tol` keyword argument
        γ = sinkhorn_unbalanced(μ, ν, C, λ1, λ2, ε; atol=1e-7)
        c = sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε; atol=1e-7)
        γ_tol = @test_deprecated sinkhorn_unbalanced(μ, ν, C, λ1, λ2, ε; tol=1e-7)
        c_tol = @test_deprecated sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε; tol=1e-7)
        @test γ_tol == γ
        @test c_tol == c

        # deprecated `max_iter` keyword argument
        γ = sinkhorn_unbalanced(μ, ν, C, λ1, λ2, ε; maxiter=50)
        c = sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε; maxiter=50)
        γ_max_iter = @test_deprecated sinkhorn_unbalanced(μ, ν, C, λ1, λ2, ε; max_iter=50)
        c_max_iter = @test_deprecated sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε; max_iter=50)
        @test γ_max_iter == γ
        @test c_max_iter == c
    end
end


text:
fermi(x) = x <= 0 ? (1 / (1 + exp(x))) : 1 - fermi(-x)
dfermi(x) = - fermi(x) * fermi(-x)

fermi(ϵ, β) = fermi(ϵ*β)
dfermi(ϵ, β) = β * dfermi(ϵ*β)


text:
## GCMGridSpec function with default gridName argument:

GCMGridSpec() = GCMGridSpec("LLC90")

## GCMGridSpec function with gridName argument:

function GCMGridSpec(gridName)

global grDir, nFaces, grTopo, ioSize, facesSize, ioPrec;

if gridName=="LLC90";
    grDir="GRID_LLC90/";
    nFaces=5;
    grTopo="llc";
    ioSize=[90 1170];
    facesSize=[[90 270];[90 270];[90 90];[270 90];[270 90]];
    ioPrec=Float64;
elseif gridName=="CS32";
    grDir="GRID_CS32/";
    nFaces=6;
    grTopo="cs";
    ioSize=[32 192];
    facesSize=[[32 32];[32 32];[32 32];[32 32];[32 32];[32 32]];
    ioPrec=Float32;
elseif gridName=="LL360";
    grDir="GRID_LL360/";
    nFaces=1;
    grTopo="ll";
    ioSize=[360 160];
    facesSize=[360 160];
    ioPrec=Float32;
else;
    error("unknown gridName case");
end;

return "GCMGridSpec: passed"

end

## GCMGridLoad function

function GCMGridLoad()

    global XC, XG, YC, YG, RAC, RAZ, DXC, DXG, DYC, DYG, hFacC, hFacS, hFacW, Depth;

    list0=("XC","XG","YC","YG","RAC","RAZ","DXC","DXG","DYC","DYG","hFacC","hFacS","hFacW","Depth");
    for ii=1:length(list0);
        tmp1=read_bin(grDir*list0[ii]*".data",ioPrec);
        tmp2=Symbol(list0[ii]);
        @eval (($tmp2) = ($tmp1))
    end

    return "GCMGridLoad: passed"

end

function GCMGridOnes(grTp,nF,nP)

    global grDir, nFaces, grTopo, ioSize, facesSize, ioPrec;
    global XC, XG, YC, YG, RAC, RAZ, DXC, DXG, DYC, DYG, hFacC, hFacS, hFacW, Depth;

    grDir=""
    grTopo=grTp
    nFaces=nF
    ioSize=[nP nP*nF];
    facesSize=fill(nP,(nF,2));
    ioPrec=Float32;

    list0=("XC","XG","YC","YG","RAC","RAZ","DXC","DXG","DYC","DYG","hFacC","hFacS","hFacW","Depth");
    for ii=1:length(list0);
        tmp1=fill(1.,nP,nP*nF);
        tmp1=convert2gcmfaces(tmp1);
        tmp2=Symbol(list0[ii]);
        @eval (($tmp2) = ($tmp1))
    end

    return "GCMGridOnes: passed"

end


text:
module Forms

using ComputedFieldTypes
using LinearAlgebra
using StaticArrays

using ..Defs

"""
    Form{D,R,T}

A differential form living in `D` dimensions (`D>=0`) with Rank `R`
(`0<=D<=R`) holding elements of type `T`.
"""
Form

export Form
@computed struct Form{D,R,T}
    elts::SVector{binomial(Val(D), Val(R)),T}
end

# Constructor without explicit type
Form{D,R}(x::Form{D,R}) where {D,R} = x
Form{D,R}(elts::SVector{N,T}) where {D,R,N,T} = Form{D,R,T}(elts)

# Constructor with added computed type (which must match)
Form{D,R,T,X}(args...) where {D,R,T,X} = Form{D,R,T}(args...)::Form{D,R,T,X}

const IteratorTypes = Union{Base.Generator,Iterators.Flatten}
function Form{D,R,T}(gen::IteratorTypes) where {D,R,T}
    N = length(Form{D,R})
    return Form{D,R,T}(SVector{N,T}(gen))
end
function Form{D,R}(gen::IteratorTypes) where {D,R}
    @assert IteratorEltype(typeof(gen)) == HasEltype()
    N = length(Form{D,R})
    T = eltype(gen)
    return Form{D,R,T}(SVector{N,T}(gen))
end

function Form{D,R,T}(tup::Tuple) where {D,R,T}
    N = length(Form{D,R})
    return Form{D,R,T}(SVector{N,T}(tup))
end
function Form{D,R}(tup::Tuple) where {D,R}
    N = length(Form{D,R})
    return Form{D,R}(SVector{N}(tup))
end
function Form{D,R}(tup::Tuple{}) where {D,R}
    return throw(TypeError("Cannot create Form from emtpy tuple without specifying type"))
end

# Conversions
function Form{D,R,T}(f::Form{D,R}) where {D,R,T}
    return Form{D,R,T}(SVector{length(Form{D,R}),T}(f.elts))
end

Base.convert(::Type{T}, x::Form{D,0,T}) where {D,T} = x.elts[1]
Base.convert(::Type{SVector}, x::Form{D,1,T}) where {D,T} = x.elts
Base.convert(::Type{SVector{D}}, x::Form{D,1,T}) where {D,T} = x.elts
Base.convert(::Type{SVector{D,T}}, x::Form{D,1,T}) where {D,T} = x.elts

export fdim
fdim(::Type{<:Form{D}}) where {D} = D
fdim(x::Form) = fdim(typeof(x))
export frank
frank(::Type{<:Form{D,R}}) where {D,R} = R
frank(x::Form) = frank(typeof(x))

################################################################################

# I/O

function Base.show(io::IO, x::Form{D,R,T}) where {D,R,T}
    print(io, "$T{$D,$R}⟦")
    for n in 1:length(x)
        n > 1 && print(io, ", ")
        inds = lin2lst(Val(D), Val(R), n)
        print(io, "[")
        for ind in inds
            print(io, ind)
        end
        print(io, "]:", x[n])
    end
    print(io, "⟧")
    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", x::Form{D,R,T}) where {D,R,T}
    skiptype = get(io, :typeinfo, Any) <: Form{D,R,T}
    if !skiptype
        print(io, "$T{$D,$R}")
    end
    print(io, "⟦")
    for n in 1:length(x)
        n > 1 && print(io, ", ")
        inds = lin2lst(Val(D), Val(R), n)
        if !get(io, :compact, false)
            print(io, "[")
            for ind in inds
                print(io, ind)
            end
            print(io, "]:")
        end
        show(IOContext(io, :compact => true, :typeinfo => T), mime, x[n])
    end
    print(io, "⟧")
    return nothing
end

################################################################################

# Comparisons

Base.:(==)(x1::Form{D,R}, x2::Form{D,R}) where {D,R} = x1.elts == x2.elts
Base.:(<)(x1::Form{D,R}, x2::Form{D,R}) where {D,R} = x1.elts < x2.elts
Base.isequal(x1::Form, x2::Form) = isequal(x1.elts, x2.elts)
Base.isless(x1::Form, x2::Form) = isless(x1.elts, x2.elts)
Base.hash(x1::Form, h::UInt) = hash(hash(x1.elts, h), UInt(0xc060e76f))
function Base.isapprox(x1::Form{D,R}, x2::Form{D,R}; kw...) where {D,R}
    scale = max(norm(x1), norm(2))
    return isapprox(scale + norm(x1 - x2), scale; kw...)
end

################################################################################

# Special constructors

export fscalar
"""
    fscalar(Val(D), x)
    fscalar(D, x)

Create a `D`-dimensional scalar (a `0`-form)
"""
fscalar(::Val{D}, x) where {D} = Form{D,0}((x,))
fscalar(D::Integer, x) = fscalar(Val(Int(D)), x)

export fvector
"""
    fvector(x...)

Create a `D`-dimensional vector (a `1`-form) from its `D` elements
"""
fvector(x...) = Form{length(x),1}(SVector(x))

export fpseudoscalar
"""
    fpseudoscalar(Val(D), x)
    fpseudoscalar(D, x)

Create a `D`-dimensional pseudo-scalar (a `D`-form)
"""
fpseudoscalar(::Val{D}, x) where {D} = Form{D,D}((x,))
fpseudoscalar(D::Integer, x) = fpseudoscalar(Val(Int(D)), x)

# Constructors from collections

export form

# TODO: Use `@generated` only to generate the index lists; access the
# array in a regular funcion
@generated function form(::Val{D}, ::Val{R}, ::Type{T}, arr::AbstractArray{T,R}) where {D,R,T}
    N = binomial(D, R)
    quote
        @assert all(==(D), size(arr))
        elts = SVector{$N,T}($([:(arr[$(lin2lst(Val(D), Val(R), n)...)]) for n in 1:N]...))
        return Form{D,R,T}(elts)
    end
end
function form(::Val{D}, ::Val{R}, ::Type{T}, arr::AbstractArray{U,R}) where {D,R,T,U}
    return Form{D,R,T}(Form{D,R,U}(arr))
end
function form(::Val{D}, ::Val{R}, arr::AbstractArray) where {D,R}
    T = eltype(arr)
    return Form{D,R,T}(arr)
end

function form(::Val{D}, ::Val{R}, ::Type{T}, tup::Tuple) where {D,R,T}
    N = binomial(Val(D), Val(R))
    return Form{D,R,T}(SVector{N,T}(tup))
end
function form(::Val{D}, ::Val{R}, tup::Tuple{}) where {D,R}
    return error("Cannot create Form from emtpy tuple without specifying type")
end
function form(::Val{D}, ::Val{R}, tup::Tuple) where {D,R}
    N = binomial(Val(D), Val(R))
    return Form{D,R}(SVector{N}(tup))
end

function form(::Val{D}, ::Val{R}, ::Type{T}, gen::IteratorTypes) where {D,R,T}
    N = binomial(Val(D), Val(R))
    elts = SVector{N,T}(gen)
    return Form{D,R,T}(elts)
end

# TODO: Use `@generated` only to generate the index lists; call the
# function in a regular funcion
@generated function form(::Val{D}, ::Val{R}, ::Type{T}, fun::Function) where {D,R,T}
    N = binomial(D, R)
    quote
        elts = SVector{$N,T}($([:(fun($(lin2lst(Val(D), Val(R), n)...))::T) for n in 1:N]...))
        return Form{D,R,T}(elts)
    end
end

export MakeForm
struct MakeForm{D,R,T} end
MakeForm{D,R,T}(elts) where {D,R,T} = form(Val(D), Val(R), T, elts)
MakeForm{D,R}(elts) where {D,R} = form(Val(D), Val(R), elts)

################################################################################

"""
We use three different representations for the "tensor indices" for
differential forms. Let us consider a rank `R` form in `D` dimensions.

- The first representation is a list of indices `[i,j,k]`. There are
  `R` indices, and each index is in the range `1:D`. Since the form is
  antisymmetric, we can choose a canonical representations where the
  indices are strictly increasing (`i < j < k`). There can be zero
  indices (for `D = 0`), and there are at most `D` indices (since `R
  <= D`). This representation is easily understood and preferred by
  humans.

- The second representation is a bit mask. Since all indices are
  different and they are ordered, each bit represents whether the
  corresponding index is present or not. For example, for `D = 4`, the
  index list `[1,2,3]` is represented by the bit mask `0111`, and
  `[2,4]` is represented by `1010`. The empty index list `[]` is
  represented by `0000`. This representation makes index manipulations
  (e.g. determining how the elements in a wedge product multiply)
  straightforward.

  This bitwise representation also clarifies that a `D`-dimensional
  form can have `2^D` elements (when summed over all ranks), and that
  a rank-`R` `D`-dimensional form has `R`-choose-`D` components.

- The third representation is a linear index. We number all bit masks
  for a particular rank, and this defines the conversion between a bit
  mask and a linear index. This representation is used to store the
  element of a form in an array.

This is how the elements of a `D = 2` form look in the various
representations:

    rank      linear    bits    indices
    ------------------------------------
    R = 0:    1         0,0     []
    R = 0:    1         0,0     []
    R = 1:    1         0,1     [1]
              2         1,0     [2]
    R = 2:    1         1,1     [1,2]

And here for `D = 3`:

    rank      linear    bits     indices
    ------------------------------------
    R = 0:    1         0,0,0    []
    R = 1:    1         0,0,1    [1]
              2         0,1,0    [2]
              3         1,0,0    [3]
    R = 2:    1         0,1,1    [1,2]
              2         1,0,1    [1,3]
              3         1,1,0    [2,3]
    R = 3:    1         1,1,1    [1,2,3]
"""
indices

# A fast type-level version of `binomial`
@generated function Base.binomial(::Val{D}, ::Val{R}) where {D,R}
    @assert 0 <= R <= D
    return binomial(D, R)
end

function bits2uint(bits::SVector{D,Bool}) where {D}
    ubits = Unsigned(0)
    for d in 1:D
        ubits += Unsigned(bits[d]) << (d - 1)
    end
    return ubits::Unsigned
end

function uint2bits(::Val{D}, ubits::Unsigned) where {D}
    bits = zeros(SVector{D,Bool})
    for d in 1:D
        bits = Base.setindex(bits, (ubits & (Unsigned(1) << (d - 1))) != 0, d)
    end
    return bits::SVector{D,Bool}
end
function uint2bits(::Val{D}, ubits::Integer) where {D}
    return uint2bits(Val(D), Unsigned(ubits))
end

"""
all bit indices for dimension `D` with rank `R`
"""
function dim_bitindices_slow(D::Int, R::Int)
    @assert 0 <= R <= D
    bitindices = Vector{SVector{D,Bool}}()
    for ibits in CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> 1, D))
        bits = SVector{D,Bool}(Tuple(ibits))
        if count(bits) == R
            push!(bitindices, bits)
        end
    end
    @assert length(bitindices) == binomial(D, R)
    return bitindices
end
@generated function dim_bitindices(::Val{D}, ::Val{R}) where {D,R}
    return dim_bitindices_slow(D, R)
end

"""
all linear indices for dimension `D` with rank `R`
"""
function dim_linindices_slow(D::Int, R::Int)
    @assert 0 <= R <= D
    bitindices = dim_bitindices(Val(D), Val(R))
    linindices = zeros(Int, 1 << D)
    for (lin, bits) in enumerate(bitindices)
        ubits = bits2uint(bits)
        linindices[ubits + 1] = lin
    end
    @assert count(!=(0), linindices) == binomial(D, R)
    return linindices
end
@generated function dim_linindices(::Val{D}, ::Val{R}) where {D,R}
    return dim_linindices_slow(D, R)
end

"""
convert a linear index to a bit index
"""
function lin2bit(::Val{D}, ::Val{R}, lin::Int) where {D,R}
    return dim_bitindices(Val(D), Val(R))[lin]::SVector{D,Bool}
end

"""
Convert a bit index to a linear index. Return 0 for failure.
"""
function bit2lin(::Val{D}, ::Val{R}, bits::SVector{D,Bool}) where {D,R}
    return dim_linindices(Val(D), Val(R))[bits2uint(bits) + 1]::Int
end

"""
Convert a bit index to an index list
"""
function bit2lst(::Val{D}, ::Val{R}, bits::SVector{D,Bool}) where {D,R}
    list = zero(SVector{R,Int})
    r = 0
    for d in 1:D
        if bits[d]
            r += 1
            list = Base.setindex(list, d, r)
        end
    end
    @assert r == R
    return list::SVector{R,Int}
end

"""
Convert an index list to a bit index
"""
function lst2bit(::Val{D}, ::Val{R}, list::SVector{R,Int}) where {D,R}
    bits = zero(SVector{D,Bool})
    for r in 1:(R - 1)
        @assert list[r] < list[r + 1]
    end
    for r in 1:R
        bits = Base.setindex(bits, true, list[r])
    end
    return bits::SVector{D,Bool}
end

function lin2lst(::Val{D}, ::Val{R}, lin::Int) where {D,R}
    return bit2lst(Val(D), Val(R), lin2bit(Val(D), Val(R), lin))
end

function lst2lin(::Val{D}, ::Val{R}, list::SVector{R,Int}) where {D,R}
    return bit2lin(Val(D), Val(R), lst2bit(Val(D), Val(R), list))
end

################################################################################

# Forms are collections

Base.eltype(::Type{<:Form{D,R,T}}) where {D,R,T} = T
Base.firstindex(::Type{<:Form}) = 1
Base.firstindex(x::Form) = firstindex(typeof(x))
Base.iterate(x::Form, state...) = iterate(x.elts, state...)
Base.ndims(::Type{<:Form{D,R}}) where {D,R} = R
Base.ndims(x::Form) = ndims(typeof(x))
Base.lastindex(::Type{<:Form{D,R}}) where {D,R} = length(Form{D,R})
Base.lastindex(x::Form) = lastindex(typeof(x))
Base.length(::Type{<:Form{D,R}}) where {D,R} = binomial(Val(D), Val(R))
Base.length(x::Form) = length(typeof(x))
Base.size(::Type{<:Form{D,R}}) where {D,R} = ntuple(d -> D, R)
Base.size(x::Form) = size(typeof(x))
Base.size(::Type{<:Form{D}}, r) where {D} = D
Base.size(x::Form, r) = size(typeof(x))

Base.getindex(x::Form, ind::Integer) = x.elts[ind]
function Base.getindex(x::Form{D,R}, inds::SVector{R}) where {D,R}
    return x.elts[lst2lin(Val(D), Val(R), inds)]
end
Base.getindex(x::Form, inds::Tuple{}) = x[SVector{0,Int}()]
Base.getindex(x::Form, inds::Tuple) = x[SVector(inds)]
Base.getindex(x::Form, inds::Integer...) = x[inds]

function Base.setindex(x::Form{D,R}, val, ind::Integer) where {D,R}
    return Form{D,R}(Base.setindex(x.elts, val, ind))
end
function Base.setindex(x::Form{D,R}, val, inds::SVector{R}) where {D,R}
    return Base.setindex(x, val, lst2lin(Val(D), Val(R), inds))
end
function Base.setindex(x::Form, val, inds::Tuple{})
    return Base.setindex(x, val, SVector{0,Int}())
end
Base.setindex(x::Form, val, inds::Tuple) = Base.setindex(x, val, SVector(inds))
Base.setindex(x::Form, val, inds::Integer...) = Base.setindex(x, val, inds)

function Base.map(f, x::Form{D,R}, ys::Form{D,R}...) where {D,R}
    return Form{D,R}(map(f, x.elts, map(y -> y.elts, ys)...))
end
Base.reduce(op, x::Form) = reduce(op, x.elts)
Base.mapreduce(f, op, x::Form{D,R}, ys::Form{D,R}...) where {D,R} = mapreduce(op, f, x.elts, map(y -> y.elts, ys)...)

################################################################################

# Forms form a vector space

function Base.rand(::Type{<:Form{D,R,T}}) where {D,R,T}
    N = binomial(Val(D), Val(R))
    return Form{D,R}(rand(SVector{N,T}))
end
function Base.zero(::Type{<:Form{D,R,T}}) where {D,R,T}
    N = binomial(Val(D), Val(R))
    return Form{D,R}(zero(SVector{N,T}))
end
Base.zero(::Type{<:Form{D,R}}) where {D,R} = zero(Form{D,R,Float64})
Base.zero(x::Form) = zero(typeof(x))
Base.iszero(x::Form) = all(iszero, x.elts)

# Deprecated
@inline Base.zeros(::Type{<:Form{D,R,T}}) where {D,R,T} = zero(Form{D,R,T})
@inline Base.zeros(::Type{<:Form{D,R}}) where {D,R} = zeros(Form{D,R,Float64})
@inline Base.zeros(x::Form) = zeros(typeof(x))

function Defs.unit(::Type{<:Form{D,R,T}}, ind::Integer) where {D,R,T}
    return setindex(zero(Form{D,R,T}), one(T), ind)
end
function Defs.unit(::Type{<:Form{D,R}}, ind::Integer) where {D,R}
    return unit(Form{D,R,Float64}, ind)
end
function Defs.unit(::Type{<:Form{D,R,T}}, inds::SVector{R}) where {D,R,T}
    return setindex(zero(Form{D,R,T}), one(T), inds)
end
function Defs.unit(::Type{<:Form{D,R}}, inds::SVector{R}) where {D,R}
    return unit(Form{D,R,Float64}, inds)
end
Defs.unit(F::Type{<:Form}, inds::Tuple{}) = unit(F, SVector{0,Int}())
Defs.unit(F::Type{<:Form}, inds::Tuple) = unit(F, SVector(inds))
Defs.unit(F::Type{<:Form}, inds::Integer...) = unit(F, inds)

Base.:+(x::Form{D,R}) where {D,R} = Form{D,R}(+x.elts)
Base.:-(x::Form{D,R}) where {D,R} = Form{D,R}(-x.elts)
Base.conj(x::Form{D,R}) where {D,R} = Form{D,R}(conj(x.elts))
Base.:+(x::Form{D,R}, y::Form{D,R}) where {D,R} = Form{D,R}(x.elts + y.elts)
Base.:-(x::Form{D,R}, y::Form{D,R}) where {D,R} = Form{D,R}(x.elts - y.elts)
Base.:*(x::Form{D,R}, a) where {D,R} = Form{D,R}(x.elts * a)
Base.:/(x::Form{D,R}, a) where {D,R} = Form{D,R}(x.elts / a)
Base.div(x::Form{D,R}, a) where {D,R} = Form{D,R}(map(b -> div(b, a), x.elts))
Base.mod(x::Form{D,R}, a) where {D,R} = Form{D,R}(map(b -> mod(b, a), x.elts))
Base.:*(a, x::Form{D,R}) where {D,R} = Form{D,R}(a * x.elts)
Base.:\(a, x::Form{D,R}) where {D,R} = Form{D,R}(a \ x.elts)

################################################################################

# Forms form an algebra

@inline Base.one(::Type{<:Form{D,R,T}}) where {D,R,T} = one(Form{D,0,T})
@inline Base.one(::Type{<:Form{D,0,T}}) where {D,T} = Form{D,0,T}((one(T),))
@inline Base.one(::Type{<:Form{D,R}}) where {D,R} = one(Form{D,0})
@inline Base.one(::Type{<:Form{D,0}}) where {D} = one(Form{D,0,Float64})
@inline Base.one(::Type{<:Form{D}}) where {D} = one(Form{D,0})
@inline Base.one(x::Form) = one(typeof(x))
@inline Base.isone(x::Form{D,0}) where {D} = isone(x[])

"""
    reverse(x)
    ~x

Reverse
"""
reverse
@inline Base.reverse(x::Form{D,R}) where {D,R} = bitsign((R - 1) * R ÷ 2) * x
@inline Base.:~(x::Form) = reverse(x)
@inline Base.inv(::typeof(~)) = ~

"""
    cycle_basis(x)

Cycle basis: `e_i => e_{i+1}`
"""
cycle_basis
@generated function cycle_basis(x1::Form{D,R}) where {D,R}
    @assert 0 <= R <= D
    U = typeof(one(eltype(x1)))
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n1 in 1:length(x1)
        lst1 = lin2lst(Val(D), Val(R), n1)
        lstr0 = mod1.((lst1 .+ 1), D)
        lstr, parity = sort_perm(lstr0)
        s = bitsign(parity)
        op = s > 0 ? :+ : :-
        ind = lst2lin(Val(D), Val(R), lstr)
        elts[ind] = :($op(x1[$n1]))
    end
    @assert !any(==(nothing), elts)
    quote
        Form{D,$R,$U}(SVector{$N,$U}($(elts...)))
    end
end
export cycle_basis

"""
    reverse_basis(x)

Reverse basis: `e_i => e_{D+1-i}`
"""
reverse_basis
@generated function reverse_basis(x1::Form{D,R}) where {D,R}
    @assert 0 <= R <= D
    U = typeof(one(eltype(x1)))
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n1 in 1:length(x1)
        lst1 = lin2lst(Val(D), Val(R), n1)
        lstr0 = (D + 1) .- lst1
        lstr, parity = sort_perm(lstr0)
        s = bitsign(parity)
        op = s > 0 ? :+ : :-
        ind = lst2lin(Val(D), Val(R), lstr)
        elts[ind] = :($op(x1[$n1]))
    end
    @assert !any(==(nothing), elts)
    quote
        Form{D,$R,$U}(SVector{$N,$U}($(elts...)))
    end
end
export reverse_basis

"""
    hodge(x)
    ⋆x   (typed: \\star<tab>)

Hodge dual
"""
hodge
@generated function hodge_algorithm(::Val{D}, ::Val{R1}) where {D,R1}
    @assert 0 <= R1 <= D
    R = D - R1
    @assert 0 <= R <= D
    N1 = binomial(D, R1)
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(R1), n1)
        bitsr = .~bits1
        _, parity = sort_perm(SVector{R1 + R,Int}(bit2lst(Val(D), Val(R1), bits1)..., bit2lst(Val(D), Val(R), bitsr)...))
        s = isodd(parity)
        ind = bit2lin(Val(D), Val(R), bitsr)
        elts[ind] = (s, n1)
    end
    @assert !any(==(nothing), elts)
    return SVector{N,Tuple{Bool,Int}}(elts)
end
@inline @inbounds function eval_hodge_term(term::Tuple{Bool,Int}, x1)
    s, i = term
    return bitsign(s) * x1[i]
end
function hodge(x1::Form{D,R1}) where {D,R1}
    @assert 0 <= R1 <= D
    R = D - R1
    @assert 0 <= R <= D
    algorithm = hodge_algorithm(Val(D), Val(R1))::SVector
    return Form{D,R}(map(term -> eval_hodge_term(term, x1), algorithm))
end
export hodge, ⋆
const ⋆ = hodge

"""
    invhodge(x)
    inv(⋆)x   (typed: \\star<tab>)

Inverse of Hodge dual: `inv(⋆)⋆x = x`
"""
invhodge
@generated function invhodge_algorithm(::Val{D}, ::Val{R1}) where {D,R1}
    @assert 0 <= R1 <= D
    R = D - R1
    @assert 0 <= R <= D
    N1 = binomial(D, R1)
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n1 in 1:N1
        bits1 = lin2bit(Val(D), Val(R1), n1)
        bitsr = .~bits1
        # The "opposite" parity as `hodge`
        _, parity = sort_perm(SVector{R + R1,Int}(bit2lst(Val(D), Val(R), bitsr)..., bit2lst(Val(D), Val(R1), bits1)...))
        s = isodd(parity)
        ind = bit2lin(Val(D), Val(R), bitsr)
        elts[ind] = (s, n1)
    end
    @assert !any(==(nothing), elts)
    return SVector{N,Tuple{Bool,Int}}(elts)
end
@inline function eval_invhodge_term(term::Tuple{Bool,Int}, x1)
    s, i = term
    @inbounds bitsign(s) * x1[i]
end
function invhodge(x1::Form{D,R1}) where {D,R1}
    @assert 0 <= R1 <= D
    R = D - R1
    @assert 0 <= R <= D
    algorithm = invhodge_algorithm(Val(D), Val(R1))::SVector
    return Form{D,R}(map(term -> eval_invhodge_term(term, x1), algorithm))
end
export invhodge
@inline Base.inv(::typeof(hodge)) = invhodge

@inline Base.conj(x::Form{D,R}) where {D,R} = Form{D,R}(conj.(x.elts))

"""
    wedge(x, y)
    x ∧ y   (typed: \\wedge<tab>)

Outer producxt
"""
wedge
@generated function wedge_algorithm(::Val{D}, ::Val{R1}, ::Val{R2}) where {D,R1,R2}
    @assert 0 <= R1 <= D
    @assert 0 <= R2 <= D
    R = R1 + R2
    @assert 0 <= R <= D
    N1 = binomial(D, R1)
    N2 = binomial(D, R2)
    N = binomial(D, R)
    elts = [Any[] for n in 1:N]
    for n1 in 1:N1, n2 in 1:N2
        bits1 = lin2bit(Val(D), Val(R1), n1)
        bits2 = lin2bit(Val(D), Val(R2), n2)
        if !any(bits1 .& bits2)
            bitsr = bits1 .| bits2
            _, parity = sort_perm(SVector{R1 + R2,Int}(bit2lst(Val(D), Val(R1), bits1)..., bit2lst(Val(D), Val(R2), bits2)...))
            s = isodd(parity)
            ind = bit2lin(Val(D), Val(R), bitsr)
            push!(elts[ind], (s, n1, n2))
        end
    end
    M = length(elts[1])
    @assert all(elt -> length(elt) == M, elts)
    elts = Tuple.(elts)::Vector{NTuple{M,Tuple{Bool,Int,Int}}}
    return SVector{N,NTuple{M,Tuple{Bool,Int,Int}}}(elts)
end
@inline @inbounds function eval_wedge_term(term::NTuple{M,Tuple{Bool,Int,Int}}, x1, x2) where {M}
    U = typeof(one(eltype(x1)) * one(eltype(x2)))
    M == 0 && return zero(U)
    s, i, j = term[1]
    r = bitsign(s) * x1[i] * x2[j]
    for m in 2:M
        s, i, j = term[m]
        r += bitsign(s) * x1[i] * x2[j]
    end
    return r
end
function wedge(x1::Form{D,R1}, x2::Form{D,R2}) where {D,R1,R2}
    @assert 0 <= R1 <= D
    @assert 0 <= R2 <= D
    R = R1 + R2
    @assert 0 <= R <= D
    algorithm = wedge_algorithm(Val(D), Val(R1), Val(R2))::SVector
    return Form{D,R}(map(term -> eval_wedge_term(term, x1, x2), algorithm))
end
@inline wedge(x::Form) = x
@inline wedge(x1::Form, x2::Form, x3s::Form...) = wedge(wedge(x1, x2), x3s...)
@inline wedge(xs::SVector{0,<:Form{D,R,T}}) where {D,R,T} = one(Form{D,0,T})
@inline wedge(xs::SVector{1,<:Form}) = xs[1]
@inline wedge(xs::SVector{N,<:Form}) where {N} = ∧(pop(xs)) ∧ last(xs)
export wedge, ∧
const ∧ = wedge

"""
    vee(x, y)
    x ∨ y   (typed: \\vee<tab>)

Regressive product: `⋆(x ∨ y) = ⋆x ∧ ⋆y`
(Inspired by Grassmann.jl)
"""
vee
@inline vee(x1::Form, x2::Form) = inv(⋆)(⋆x1 ∧ ⋆x2)
@inline vee(x::Form) = x
@inline vee(x1::Form, x2::Form, x3s::Form...) = vee(vee(x1, x2), x3s...)
@inline vee(xs::SVector{0,<:Form{D,R,T}}) where {D,R,T} = ⋆one(Form{D,0,T})
@inline vee(xs::SVector{1,<:Form}) = xs[1]
@inline vee(xs::SVector{N,<:Form}) where {N} = ∨(pop(xs)) ∨ last(xs)
export vee, ∨
const ∨ = vee

"""
    dot(x, y)
    x ⋅ y   (typed: \\cdot<tab>)

Dot product: `x ⋅ y = x ∨ ⋆y`
(Inspired by Grassmann.jl)
"""
@inline LinearAlgebra.dot(x1::Form, x2::Form) = x1 ∨ ⋆x2
export dot, ⋅

# Base.abs2(x::Form) = (x ⋅ x)[]
# Base.abs(x::Form) = sqrt(abs2(x))
@inline norm2(x::Missing) = missing
@inline norm2(x::Number) = abs2(x)
@inline norm2(x::AbstractArray) = sum(norm2.(x))
# @inline norm2(x::Form) = (x ⋅ x)[]
@inline norm2(x::Form) = norm2(x.elts)
export norm2
@inline LinearAlgebra.norm(x::Form) = sqrt(norm2(x))
export norm

"""
    cross(x, y)
    x × y   (typed: \\times<tab>)

Dot product: `x × y = ⋆(x ∧ y)`
(Inspired by Grassmann.jl)
"""
@inline LinearAlgebra.cross(x1::Form, x2::Form) = ⋆(x1 ∧ x2)
export cross, ×

################################################################################

"""
    tensorsum(x, y)
    x ⊕ y   (typed: \\oplus<tab>)

Tensor sum
"""
tensorsum
@generated function tensorsum_algorithm(::Val{D1}, ::Val{D2}, ::Val{R}) where {D1,D2,R}
    @assert 0 < R <= D1
    @assert 0 < R <= D2
    D = D1 + D2
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n in 1:N
        bitsr = lin2bit(Val(D), Val(R), n)
        bits1 = SVector{D1,Bool}(bitsr[1:D1])
        bits2 = SVector{D2,Bool}(bitsr[(D1 + 1):end])
        n1 = bit2lin(Val(D1), Val(R), bits1)
        n2 = bit2lin(Val(D2), Val(R), bits2)
        if n1 > 0
            @assert n2 == 0
            elts[n] = (Val(1), n1)
        elseif n2 > 0
            @assert n1 == 0
            elts[n] = (Val(2), n2)
        else
            elts[n] = nothing
        end
    end
    return Tuple(elts)
end
@inline function eval_tensorsum_term(term::Nothing, x1, x2)
    U = typeof(one(eltype(x1)) * one(eltype(x2)))
    return zero(U)
end
@inline function eval_tensorsum_term(term::Tuple{Val{I},Int}, x1, x2) where {I}
    _, i = term
    I == 1 && return (@inbounds x1[i])
    return I == 2 && return (@inbounds x2[i])
end
function tensorsum(x1::Form{D1,R}, x2::Form{D2,R}) where {D1,D2,R}
    @assert 0 < R <= D1
    @assert 0 < R <= D2
    D = D1 + D2
    algorithm = tensorsum_algorithm(Val(D1), Val(D2), Val(R))
    return Form{D,R}(map(term -> eval_tensorsum_term(term, x1, x2), algorithm))
end
tensorsum(x::Form) = x
function tensorsum(x1::Form, x2::Form, x3s::Form...)
    return tensorsum(tensorsum(x1, x2), x3s...)
end
const ⊕ = tensorsum
export tensorsum, ⊕

"""
    tensorproduct(x, y)
    x ⊗ y   (typed: \\otimes<tab>)

Tensor product
"""
tensorproduct
@generated function tensorproduct_algorithm(::Val{D1}, ::Val{R1}, ::Val{D2}, ::Val{R2}) where {D1,R1,D2,R2}
    @assert 0 <= R1 <= D1
    @assert 0 <= R2 <= D2
    D = D1 + D2
    R = R1 + R2
    N = binomial(D, R)
    elts = Any[nothing for n in 1:N]
    for n in 1:N
        bitsr = lin2bit(Val(D), Val(R), n)
        bits1 = SVector{D1,Bool}(bitsr[1:D1])
        bits2 = SVector{D2,Bool}(bitsr[(D1 + 1):end])
        n1 = bit2lin(Val(D1), Val(R1), bits1)
        n2 = bit2lin(Val(D2), Val(R2), bits2)
        if n1 > 0 && n2 > 0
            elts[n] = (n1, n2)
        else
            elts[n] = nothing
        end
    end
    return Tuple(elts)
end
@inline function eval_tensorproduct_term(term::Nothing, x1, x2)
    U = typeof(one(eltype(x1)) * one(eltype(x2)))
    return zero(U)
end
@inline function eval_tensorproduct_term(term::Tuple{Int,Int}, x1, x2) where {I}
    i, j = term
    return (@inbounds x1[i] * x2[j])
end
function tensorproduct(x1::Form{D1,R1}, x2::Form{D2,R2}) where {D1,R1,D2,R2}
    @assert 0 <= R1 <= D1
    @assert 0 <= R2 <= D2
    D = D1 + D2
    R = R1 + R2
    algorithm = tensorproduct_algorithm(Val(D1), Val(R1), Val(D2), Val(R2))
    return Form{D,R}(map(term -> eval_tensorproduct_term(term, x1, x2), algorithm))
end
tensorproduct(x::Form) = x
function tensorproduct(x1::Form, x2::Form, x3s::Form...)
    return tensorproduct(tensorproduct(x1, x2), x3s...)
end
const ⊗ = tensorproduct
export tensorproduct, ⊗

end


text:
using Documenter

using STMOZOO
using STMOZOO.LocalSearch

makedocs(sitename="STMO ZOO",
    format = Documenter.HTML(),
    modules=[LocalSearch], # add your module
    pages=Any[
        "Sudoku_local" => "man/Sudoku_local.md"
    ])

#=
deploydocs(
            repo = "github.com/michielstock/STMOZOO.jl.git",
        )
=#


text:
"""
    Intermodes <: AbstractImageBinarizationAlgorithm
    Intermodes()

    binarize([T,] img, f::AdaptiveThreshold; [window_size])
    binarize!([out,] img, f::AdaptiveThreshold; [window_size])

Under the assumption that the image histogram is bimodal the image histogram is
smoothed using a length-3 mean filter until two modes remain. The binarization
threshold is then set to the average value of the two modes.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, Intermodes())
```

## Reference

1. C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” CVGIP: Graphical Models and Image Processing, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
"""
struct Intermodes <: AbstractImageBinarizationAlgorithm end

function (f::Intermodes)(out::GenericGrayImage, img::GenericGrayImage)
    edges, counts = build_histogram(img,  256)
    t = find_threshold(HistogramThresholding.Intermodes(), counts[1:end], edges)
    @simd for i in CartesianIndices(img)
      out[i] = img[i] < t ? 0 : 1
    end
    out
end

(f::Intermodes)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
function IterationEvent(arg0::Object, arg1::jint)
    return IterationEvent((Object, jint), arg0, arg1)
end

function get_iterations(obj::IterationEvent)
    return jcall(obj, "getIterations", jint, ())
end

function get_source(obj::EventObject)
    return jcall(obj, "getSource", Object, ())
end

function to_string(obj::EventObject)
    return jcall(obj, "toString", JString, ())
end


text:
mutable struct ProjMPS
  lpos::Int
  rpos::Int
  nsite::Int
  M::MPS
  LR::Vector{ITensor}
  ProjMPS(M::MPS) = new(0,
                        length(M)+1,
                        2,
                        M,
                        Vector{ITensor}(undef, length(M)))
  ProjMPS(nsite::Int, M::MPS) = new(0,
                        length(M)+1,
                        nsite,
                        M,
                        Vector{ITensor}(undef, length(M)))
end

nsite(P::ProjMPS) = P.nsite

Base.length(P::ProjMPS) = length(P.M)

function lproj(P::ProjMPS)
  (P.lpos <= 0) && return nothing
  return P.LR[P.lpos]
end

function rproj(P::ProjMPS)
  (P.rpos >= length(P)+1) && return nothing
  return P.LR[P.rpos]
end

function product(P::ProjMPS,
                 v::ITensor)::ITensor
  if P.rpos - P.lpos != nsite(P) + 1
    error("P.lpos and P.rpos values inconsistent with nsite (must satisfy P.rpos - P.lpos = nsite + 1")
  end

  tensor_list = ITensor[]
  for j = P.lpos+1:P.rpos-1
    if j == P.lpos+1
      Lpm = dag(prime(P.M[j],"Link"))
      !isnothing(lproj(P)) && (Lpm *= lproj(P))
      push!(tensor_list, copy(Lpm))
    elseif j == P.rpos-1
      Rpm = dag(prime(P.M[j],"Link"))
      !isnothing(rproj(P)) && (Rpm *= rproj(P))
      push!(tensor_list, copy(Rpm))
    else
      Mpm = dag(prime(P.M[j],"Link"))
      push!(tensor_list, copy(Mpm))
    end
  end
  pm = prod(tensor_list)

  pv = scalar(pm*v)

  Mv = pv*dag(pm)

  return noprime(Mv)
end

#function Base.eltype(P::ProjMPS)
#  elT = eltype(P.M[P.lpos+1])
#  for j = P.lpos+2:P.rpos-1
#    elT = promote_type(elT,eltype(P.M[j]))
#  end
#  if !isnull(lproj(P))
#    elT = promote_type(elT,eltype(lproj(P)))
#  end
#  if !isnull(rproj(P))
#    elT = promote_type(elT,eltype(rproj(P)))
#  end
#  return elT
#end

(P::ProjMPS)(v::ITensor) = product(P,v)

#function Base.size(P::ProjMPS)::Tuple{Int,Int}
#  d = 1
#  if P.lpos > 0
#    d *= dim(linkind(M,P.lpos))
#  end
#  for j=P.lpos+1:P.rpos-1
#    d *= dim(siteind(P.M,j))
#  end
#  if P.rpos-1 < N
#    d *= dim(linkind(M,P.rpos-1))
#  end
#  return (d,d)
#end

function makeL!(P::ProjMPS,
                psi::MPS,
                k::Int)
  while P.lpos < k
    ll = P.lpos
    if ll <= 0
      P.LR[1] = psi[1]*dag(prime(P.M[1],"Link"))
      P.lpos = 1
    else
      P.LR[ll+1] = P.LR[ll]*psi[ll+1]*dag(prime(P.M[ll+1],"Link"))
      P.lpos += 1
    end
  end
end

function makeR!(P::ProjMPS,
                psi::MPS,
                k::Int)
  N = length(P.M)
  while P.rpos > k
    rl = P.rpos
    if rl >= N+1
      P.LR[N] = psi[N]*dag(prime(P.M[N],"Link"))
      P.rpos = N
    else
      P.LR[rl-1] = P.LR[rl]*psi[rl-1]*dag(prime(P.M[rl-1],"Link"))
      P.rpos -= 1
    end
  end
end

function position!(P::ProjMPS,
                   psi::MPS,
                   pos::Int)
  makeL!(P,psi,pos-1)
  makeR!(P,psi,pos+nsite(P))

  #These next two lines are needed
  #when moving lproj and rproj backward
  P.lpos = pos-1
  P.rpos = pos+nsite(P)
end


text:
module TestDataArray
	using Base.Test
	using DataArrays

	v = [1, 2, 3, 4]
	dv = DataArray(v, falses(size(v)))

	m = [1 2; 3 4]
	dm = DataArray(m, falses(size(m)))

	t = Array(Int, 2, 2, 2)
	t[1:2, 1:2, 1:2] = 1
	dt = DataArray(t, falses(size(t)))

	dv = DataArray(v)
	dv = DataArray(v, [false, false, false, false])

	dv = DataArray(Int, 2)
	dm = DataArray(Int, 2, 2)
	dt = DataArray(Int, 2, 2, 2)

	similar(dv)
	similar(dm)
	similar(dt)

	similar(dv, 2)
	similar(dm, 2, 2)
	similar(dt, 2, 2, 2)

    x = DataArray([9, 9, 8])
    y = DataArray([1, 9, 3, 2, 2])
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = DataArray([9, 9, 8])
    y = [1, 9, 3, 2, 2]
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = @data [1, 2, NA]
    y = @data [3, NA, 5]
    @test isequal(copy(x), x)
    @test isequal(copy!(y, x), x)

	# Test vecbind
	# a = [1:4]
	# d = DataArray(a)
	# @assert isequal(vecbind(a,a), [a,a])
	# @assert isequal(vecbind(a,1.0 * a), 1.0 * [a,a])
	# @assert isequal(vecbind(d,a), [d,d])
	# @assert isequal(vecbind(a,d), [d,d])
	# @assert isequal(vecbind(1.0 * a,d), 1.0 * [d,d])
	# @assert isequal(vecbind(d,d), [d,d])
	# @assert isequal(vecbind(a,IndexedVector(a)), [a,a])
	# @assert isequal(vecbind(a,IndexedVector(d)), [d,d])
	# @assert isequal(vecbind(PooledDataArray(a),IndexedVector(d)), [d,d])
	# @assert isequal(vecbind(PooledDataArray(a),IndexedVector(a)), [a,a])
	# @assert isequal(vecbind(a,RepeatedVector(a,2)), [a,a,a])
	# @assert (vecbind(a,StackedVector({a,1.0*a})) == [a,a,a])
end


text:
## Constructor when not calibrated
function Param(name :: Symbol, description :: T2, symbol :: T2, defaultValue :: T1) where {T1 <: Any,  T2 <: AbstractString}

    return Param(name, description, symbol, defaultValue, defaultValue,
        defaultValue .- 0.001, defaultValue .+ 0.001, false)
end

function validate(p :: Param{F1}; silent = true) where F1
    sizeV = size(default_value(p));
    isValid = true;
    if !Base.isempty(value(p))
        (size(value(p)) == sizeV)  ||  (isValid = false);
    end
    if !Base.isempty(p.lb)
        (size(p.lb) == sizeV)  ||  (isValid = false);
        (size(p.ub) == sizeV)  ||  (isValid = false);
    end
    if !isValid  &&  !silent
        @warn """
            Invalid Param $p
            default value:  $(default_value(p))
            value:          $(value(p))
            lb:             $(p.lb)
            ub:             $(p.ub)
        """
    end
    return isValid
end

## -----------  Access

Base.size(p :: AbstractParam) = size(default_value(p));
Base.length(p :: AbstractParam) = length(default_value(p));

name(p :: AbstractParam) = p.name;
lsymbol(p :: AbstractParam) = p.symbol;

"""
	$(SIGNATURES)

Is this parameter calibrated?
"""
is_calibrated(p :: AbstractParam) = p.isCalibrated;

"""
	$(SIGNATURES)

Retrieve value of a `Param`.
"""
value(p :: Param{F1}) where F1 = p.value;

default_value(p :: AbstractParam) = p.defaultValue;
lb(p :: AbstractParam) = p.lb;
ub(p :: AbstractParam) = p.ub;

# This is what the numerical optimizer sees (only the calibrated entries).
calibrated_value(p :: Param{F1}) where F1 =
    is_calibrated(p) ? value(p) : nothing;
calibrated_lb(p :: AbstractParam) = lb(p);
calibrated_ub(p :: AbstractParam) = ub(p);

# Is a parameter value close to lower or upper bounds?
close_to_lb(p :: AbstractParam; rtol = 0.01)  =
    any((value(p) .- lb(p)) ./ (ub(p) .- lb(p)) .< rtol);

close_to_ub(p :: AbstractParam; rtol = 0.01)  =
    any((ub(p) .- value(p)) ./ (ub(p) .- lb(p)) .< rtol);

close_to_bounds(p :: AbstractParam; rtol = 0.01) where F1 =
    close_to_lb(p; rtol = rtol) || close_to_ub(p; rtol = rtol);

## ------------  Show

# Short string that summarizes the `Param`
# For reporting parameters in a table
function show_string(p :: AbstractParam)
    if is_calibrated(p)
        calStr = "calibrated";
    else
        calStr = "fixed";
    end
    return string(p.name) * ": " * value_string(p) * "  ($calStr)"
end

# Summary of the value
function value_string(p :: AbstractParam)
    pType = eltype(value(p));
    if isa(value(p), Real)
        outStr = string(round(value(p), digits = 3));
    elseif isa(value(p), Array)
        outStr = "Array{$pType} of size $(size(value(p)))";
    else
        outStr = "of type $pType";
    end
    return outStr
end

Base.show(io :: IO,  p :: Param{F1}) where F1 =
    print(io, "Param:  " * show_string(p));

function short_string(p :: AbstractParam)
    vStr = formatted_value(value(p));
    return "$(name(p)): $vStr"
end

"""
    report_param

Short summary of parameter and its value.
Can be used to generate a simple table of calibrated parameters.
"""
function report_param(p :: AbstractParam)
    vStr = formatted_value(value(p));
    println("\t$(p.description):\t$(name(p)) = $vStr")
end

## ----------  Change / update

function set_calibration_status!(p :: AbstractParam, isCalibrated :: Bool)
    if isCalibrated
        calibrate!(p);
    else
        fix!(p);
    end
end

"""
    $(SIGNATURES)

Change calibration status to `true`
"""
function calibrate!(p :: Param{F1}) where F1
    p.isCalibrated = true
    return nothing
end

"""
    $(SIGNATURES)

Change calibration status to `false`
"""
function fix!(p :: Param{F1}; pValue = nothing) where F1
    p.isCalibrated = false;
    if !isnothing(pValue)
        set_value!(p, pValue);
        set_default_value!(p, pValue);
    end
    return nothing
end

"""
	$(SIGNATURES)

Set bounds for a `Param`.
"""
function set_bounds!(p :: AbstractParam; lb = nothing, ub = nothing)
    isnothing(lb)  ||  (p.lb = lb);
    isnothing(ub)  ||  (p.ub = ub);
    @assert size(ModelParams.lb(p)) == size(ModelParams.ub(p)) == size(p);
end

"""
    $(SIGNATURES)

Set parameter value. Not used during calibration.
"""
function set_value!(p :: AbstractParam, vIn;
    skipInvalidSize = false) where F1

    oldValue = value(p);
    if size(default_value(p)) == size(vIn)
        p.value = deepcopy(vIn);
    else
        @warn("""
            Wrong size for $p
            Given: $(size(vIn))
            Expected: $(size(default_value(p)))
            """);
        if !skipInvalidSize
            error("Stopped");
        end
    end
    return oldValue
end

# This is used during the calibration
set_calibrated_value!(p :: AbstractParam, vIn; skipInvalidSize = false) =
    set_value!(p, vIn; skipInvalidSize);

function set_default_value!(p :: AbstractParam, vIn)
    @assert size(default_value(p)) == size(vIn)  "Size invalid for $(p.name): $(size(vIn)). Expected $(size(default_value(p)))"
    p.defaultValue = vIn
end

"""
    $(SIGNATURES)

Update a parameter with optional arguments.
"""
function update!(p :: AbstractParam; value = nothing, defaultValue = nothing,
    lb = nothing, ub = nothing, isCalibrated = nothing)
    if !isnothing(value)
        set_value!(p, value)
    end
    if (!isnothing(lb))  ||  (!isnothing(ub))
        set_bounds!(p; lb, ub);
    end
    if !isnothing(defaultValue)
        set_default_value!(p, defaultValue);
    end
    if !isnothing(isCalibrated)
        set_calibration_status!(p, isCalibrated);
    end
    return nothing
end

# """
# 	$(SIGNATURES)

# Compare two Params. Report differences as Vector{Symbol}
# """
# function param_diffs(p1 :: Param{F1}, p2 :: Param{F2};
#     ignoreCalibrationStatus :: Bool = true) where {F1, F2}

#     diffs = Vector{Symbol}();
#     if (F1 == F2)
#         if !isequal(value(p1), value(p2))
#             push!(diffs, :value);
#         end
#     else
#         push!(diffs, :type);
#     end
#     return diffs
# end

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