text:
function loc2cs3D(loc1::Union{Int64,Array{Int64}},loc2::Union{Int64,Array{Int64}},loc3::Union{Int64,Array{Int64}},n::Array{Int64,1})
@inbounds cs = loc1 .+ (loc2.-1)*n[1] .+ (loc3.-1)*n[1]*n[2];
return cs;
end

function generateProjectionOperator(param::SfSParam)
Mesh = param.ObjMesh;
CamLoc = param.CamLoc;
XScreenLoc = param.XScreenLoc;
ScreenMesh = param.ScreenMeshX2X3;
theta_phi = param.theta_phi_rad;

X = getCellCenteredGrid(Mesh);
mid_domain = [(Mesh.domain[1]+Mesh.domain[2])/2,(Mesh.domain[3]+Mesh.domain[4])/2,(Mesh.domain[5]+Mesh.domain[6])/2];
x_mesh = (Mesh.domain[1] + 0.5*Mesh.h[1]):Mesh.h[1]:(Mesh.domain[2] - 0.5*Mesh.h[1]);

X[:,1] .= X[:,1] .- mid_domain[1];
X[:,2] .= X[:,2] .- mid_domain[2];
X[:,3] .= X[:,3] .- mid_domain[3];

X_screen = getCellCenteredGrid(ScreenMesh);
midScreen = [(ScreenMesh.domain[1]+ScreenMesh.domain[2])/2,(ScreenMesh.domain[3]+ScreenMesh.domain[4])/2];
X_screen[:,1] .-= midScreen[1];
X_screen[:,2] .-= midScreen[2];

half_screen_domain_length = [(ScreenMesh.domain[2]-ScreenMesh.domain[1])/2,(ScreenMesh.domain[4]-ScreenMesh.domain[3])/2];

offset = 0;
n =  Mesh.n;
invh = 1.0 ./Mesh.h;

volume = prod(Mesh.n);
n_tup = tuple(Mesh.n...);
n_screen_tup = tuple(ScreenMesh.n...);
n_screen = ScreenMesh.n;
numelScreen = prod(ScreenMesh.n);

I = ones(Int32,(n[1]*numelScreen));
J = ones(Int32,(n[1]*numelScreen));
V = zeros(Float32,(n[1]*numelScreen));
offset = 0;
g_idx = 1;
for k=1:numelScreen
	## definition of the line between screen point and camera point
	X_screen_k = X_screen[k,:];
	## projecting point Z to the screen:
	# ProjectionLoc(t) = CamLoc + t*(X_screen - CamLoc)
	# We need all the mesh points which are on this line;
	# For a given x location in the mesh, we search for the YZ locations on the line

	# For the first coordinate: x_mesh_jj = CamLoc[1] + t*(X_screen - CamLoc[1]);
	for ii = 1:n[1]
		x_mesh_ii = x_mesh[ii];
		t = (x_mesh_ii - CamLoc[1])/(XScreenLoc - CamLoc[1]);
		y_mesh = CamLoc[2] + t*(X_screen[k,1] - CamLoc[2]);
		z_mesh = CamLoc[3] + t*(X_screen[k,2] - CamLoc[3]);
		y_mesh += mid_domain[2];
		z_mesh += mid_domain[3];
		jj = round(Int64,y_mesh*invh[2]+0.5);
		kk = round(Int64,z_mesh*invh[3]+0.5);
		if jj <= n[2] && jj > 0 && kk <= n[3] && kk > 0
			J[g_idx] = loc2cs3D(ii,jj,kk,Mesh.n);
			I[g_idx] = k;
			V[g_idx] = 1.0;
			g_idx += 1;
		end
	end
end
RT = sparse(J,I,V,volume,numelScreen);

return RT;
end

# function generateProjectionOperatorWithRotation(param::SfSParam)
# Mesh = param.ObjMesh;
# CamLoc = param.CamLoc;
# XScreenLoc = param.XScreenLoc;
# ScreenMesh = param.ScreenMeshX2X3;
# theta_phi = param.theta_phi_rad;

# X = convert(Array{Float32,2},getCellCenteredGrid(Mesh));
# mid_domain = [(Mesh.domain[1]+Mesh.domain[2])/2,(Mesh.domain[3]+Mesh.domain[4])/2,(Mesh.domain[5]+Mesh.domain[6])/2];
# X[:,1] = X[:,1] - mid_domain[1];
# X[:,2] = X[:,2] - mid_domain[2];
# X[:,3] = X[:,3] - mid_domain[3];

# midScreen = [(ScreenMesh.domain[1]+ScreenMesh.domain[2])/2,(ScreenMesh.domain[3]+ScreenMesh.domain[4])/2];
# half_screen_domain_length = [(ScreenMesh.domain[2]-ScreenMesh.domain[1])/2,(ScreenMesh.domain[4]-ScreenMesh.domain[3])/2];

# num_rotations = size(theta_phi,1);
# offset = 0;
# n =  Mesh.n;
# volume = prod(Mesh.n);

# n_screen_tup = tuple(ScreenMesh.n...);

# numelScreen = prod(ScreenMesh.n);

# globalRowIdx = zeros(Int32,(volume*num_rotations+5));
# globalColIdx = zeros(Int32,(volume*num_rotations+5));

# calcSR = 1;

# R = 0;
# offset = 0;
# for k=1:num_rotations

	# ## Rotating the original locations X.
	# Rk = convert(Array{Float32,2},getRotate3D(deg2rad.(theta_phi[k,1]),deg2rad.(theta_phi[k,2])));
	# Xk = X*Rk';
	# Zk = Xk; Xk = 0;

	# ## projecting point Z to the screen:
	# # ProjectionLoc(t) = CamLoc + t*(X - CamLoc)
	# for k=1:3  ## Zk = Xk - CamLoc
		# Zk[:,k] = Zk[:,k] - CamLoc[k];
	# end;
	# # ProjectionLoc(t) at screelLock will have first coordinate equals XScreenLoc
	# tk = (XScreenLoc - CamLoc[1])./Zk[:,1];
	# Sk = zeros(eltype(Zk),size(Zk,1),2);
	# Sk[:,1] = CamLoc[2] + tk.*Zk[:,2];
	# Sk[:,2] = CamLoc[3] + tk.*Zk[:,3];

	# Sk[:,1] = (Sk[:,1] + half_screen_domain_length[1])*(1 ./ScreenMesh.h[1]) + 0.5;
	# Sk[:,2] = (Sk[:,2] + half_screen_domain_length[2])*(1 ./ScreenMesh.h[2]) + 0.5;
	# Sk = round.(Int32,Sk);
	# II = find((Sk[:,1].>=1) .& (Sk[:,1].<=ScreenMesh.n[1]) .& (Sk[:,2].>=1) .& (Sk[:,2].<=ScreenMesh.n[2]) );
	# Sk = Sk[II,:];
	# RowIdx_k = Base.sub2ind(n_screen_tup,Sk[:,1],Sk[:,2]);
	# globalColIdx[(offset+1):(offset + length(II))] = II;
	# globalRowIdx[(offset+1):(offset + length(II))] = RowIdx_k;
	# offset = offset + length(II);
# end
# X = 0;
# II = globalRowIdx .> 0;
# globalRowIdx = globalRowIdx[II];
# globalColIdx = globalColIdx[II];
# II = 0;
# V = ones(Float32,size(globalColIdx));
# RT = sparse(globalColIdx,globalRowIdx,V,volume,num_rotations*numelScreen);
# return RT;
# end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Variables #
# ========================== #
function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, set::DecisionSet)
    return VariableConstrainedOnCreation(variable, decision_set(variable, set))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, set::DecisionSet)
    return VariablesConstrainedOnCreation(variables, decision_set(variables, set))
end

function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, set::KnownSet)
    return VariableConstrainedOnCreation(variable, decision_set(variable, set))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, set::KnownSet)
    return VariablesConstrainedOnCreation(variables, decision_set(variables, set))
end

function JuMP.add_variable(model::Model, variable::VariableConstrainedOnCreation{<:SingleDecisionSet}, name::String)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Create a regular JuMP variable if decisions are not handled
        return JuMP.add_variable(model, variable.scalar_variable, name)
    end
    var_index, con_index = MOI.add_constrained_variable(backend(model), variable.set)
    set_stage!(decisions, var_index, variable.set.stage)
    # Map to model decisions after indices are known
    if !has_decision(decisions, var_index)
        # Store decision if is seen for the first time
        set_decision!(decisions, var_index, variable.set.decision)
    else
        # Reuse if decision has been created already
        MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, decision(decisions, var_index)))
    end
    # Add any given decision constraints
    if state(variable.set.decision) != Known
        _moi_constrain_decision(backend(model), decisions, var_index, variable.scalar_variable.info, variable.set)
    end
    # Finally, set any given name
    if !isempty(name)
        MOI.set(backend(model), MOI.VariableName(), var_index, name)
    end
    # Return created decision as DecisionRef
    return DecisionRef(model, var_index)
end

function JuMP.add_variable(model::Model, variable::VariablesConstrainedOnCreation{<:MultipleDecisionSet}, names)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Create regular JuMP variables if decisions are not handled
        var_refs = map(zip(variable.scalar_variables, names)) do (scalar_variable, name)
            JuMP.add_variable(model, scalar_variable, name)
        end
        return reshape_vector(var_refs, variable.shape)
    end
    if variable.set.constraint isa NoSpecifiedConstraint
        # Sufficient to use single decision bridges if decisions
        # are not initially constrained.
        refs = map(zip(variable.scalar_variables, variable.set.decisions, names)) do (scalar_variable, decision, name)
            set = SingleDecisionSet(variable.set.stage, decision, variable.set.constraint, variable.set.is_recourse)
            add_variable(model, VariableConstrainedOnCreation(scalar_variable, set), name)
        end
        return reshape_vector(refs, variable.shape)
    else
        var_indices, con_index = MOI.add_constrained_variables(backend(model), variable.set)
        # Map to model decisions after indices are known
        seen_decisions = Vector{Decision{Float64}}()
        for (i, var_index) in enumerate(var_indices)
            set_stage!(decisions, var_index, variable.set.stage)
            if !has_decision(decisions, var_index)
                # Store decision if is seen for the first time
                set_decision!(decisions, var_index, variable.set.decisions[i])
            else
                # Reuse if decision has been created already
                push!(seen_decisions, decision(decisions, var_index))
            end
        end
        if !isempty(seen_decisions)
            # Sanity check
            length(seen_decisions) == length(variable.scalar_variables) || error("Inconsistency in number of seen decisions and created variables.")
            # Update decision set for reuse
            MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, seen_decisions))
        end
        # Add any given decision constraints
        for (index, scalar_variable, decision) in zip(var_indices, variable.scalar_variables, variable.set.decisions)
            if decision.state != Known
                _moi_constrain_decision(backend(model), decisions, index, scalar_variable.info, variable.set)
            end
        end
        # Finally, set any given names
        for (var_index, name) in zip(var_indices, JuMP.vectorize(names, variable.shape))
            if !isempty(name)
                MOI.set(backend(model), MOI.VariableName(), var_index, name)
            end
        end
        # Return created decisions as DecisionRef
        refs = map(var_indices) do index
            DecisionRef(model, index)
        end
        return reshape_vector(refs, variable.shape)
    end
end

# Containers #
# ========================== #
struct DecisionArray{A <: AbstractArray, S <: Union{DecisionSet, KnownSet}}
    array::A
    set::S

    function DecisionArray(array::AbstractArray, set::Union{DecisionSet, KnownSet})
        A = typeof(array)
        S = typeof(set)
        return new{A, S}(array, set)
    end
end

struct DecisionDenseAxisArray{A <: DenseAxisArray, S <: Union{DecisionSet, KnownSet}}
    array::A
    set::S

    function DecisionDenseAxisArray(array::DenseAxisArray, set::Union{DecisionSet, KnownSet})
        A = typeof(array)
        S = typeof(set)
        return new{A, S}(array, set)
    end
end

struct DecisionSparseAxisArray{A <: SparseAxisArray, S <: Union{DecisionSet, KnownSet}}
    array::A
    set::S

    function DecisionSparseAxisArray(array::SparseAxisArray, set::Union{DecisionSet, KnownSet})
        A = typeof(array)
        S = typeof(set)
        return new{A, S}(array, set)
    end
end

function JuMP.build_variable(_error::Function, variables::AbstractArray{<:JuMP.ScalarVariable}, set::Union{DecisionSet,KnownSet})
    return DecisionArray(variables, set)
end

function JuMP.build_variable(_error::Function, variables::DenseAxisArray{<:JuMP.ScalarVariable}, set::Union{DecisionSet,KnownSet})
    return DecisionDenseAxisArray(variables, set)
end

function JuMP.build_variable(_error::Function, variables::SparseAxisArray{<:JuMP.ScalarVariable}, set::Union{DecisionSet,KnownSet})
    return DecisionSparseAxisArray(variables, set)
end

function JuMP.add_variable(model::Model, variable::DecisionArray, names::AbstractArray{String})
    array = variable.array
    refs = Array{DecisionRef}(undef, size(array))
    for idx in eachindex(array)
        var = array[idx]
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, decision_set(var, variable.set)),
                                      names[idx])
    end
    return refs
end

function JuMP.add_variable(model::Model, variable::DecisionDenseAxisArray, names::DenseAxisArray{String})
    array = variable.array
    refs = DenseAxisArray{DecisionRef}(undef, array.axes...)
    for idx in eachindex(array)
        var = array[idx]
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, decision_set(var, variable.set)),
                                      names[idx])
    end
    return refs
end

function JuMP.add_variable(model::Model,
                           variable::DecisionSparseAxisArray{SparseAxisArray{T,N,K}},
                           names::SparseAxisArray{String}) where {T, N, K}
    refs = SparseAxisArray(Dict{K,DecisionRef}())
    for (idx, var) in variable.array.data
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, decision_set(var, variable.set)),
                                      names[idx])
    end
    return refs
end

# Constraints #
# ========================== #
function JuMP._functionize(dref::DecisionRef)
    return convert(DecisionAffExpr{Float64}, dref)
end
function JuMP._functionize(drefs::AbstractArray{DecisionRef})
    return JuMP._functionize.(drefs)
end

function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, set::MOI.AbstractScalarSet)
    offset = constant(aff.variables)
    JuMP.add_to_expression!(aff.variables, -offset)
    shifted_set = MOIU.shift_constant(set, -offset)
    return JuMP.ScalarConstraint(aff, shifted_set)
end

function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, lb::Real, ub::Real)
    JuMP.build_constraint(_error, aff, MOI.Interval(lb, ub))
end

function JuMP.add_constraint(model::Model,
                             constraint::ScalarConstraint{DecisionRef, S},
                             name::String = "") where S <: MOI.AbstractScalarSet
    decisions = get_decisions(model)::Decisions
    check_belongs_to_model(constraint, model)
    ci = CI{SingleDecision, S}(moi_function(constraint).decision.value)
    inner = moi_add_constraint(backend(model), moi_function(constraint), moi_set(constraint))
    map_constraint!(decisions, ci, inner)
    con_ref = ConstraintRef(model, inner, shape(constraint))
    if !isempty(name)
        set_name(con_ref, name)
    end
    return con_ref
end

# Helper function #
# ========================== #
function _moi_constrain_decision(backend::MOI.ModelLike,
                                 decisions::Decisions,
                                 index::MOI.VariableIndex,
                                 info::VariableInfo,
                                 set::Union{SingleDecisionSet, MultipleDecisionSet})
    nothing_added = true
    if info.has_lb
        ci = CI{SingleDecision, MOI.GreaterThan{Float64}}(index.value)
        inner = MOI.add_constraint(backend, SingleDecision(index),
                                   MOI.GreaterThan{Float64}(info.lower_bound))
        map_constraint!(decisions, ci, inner)
        nothing_added &= false
    end
    if info.has_ub
        ci = CI{SingleDecision, MOI.LessThan{Float64}}(index.value)
        inner = MOI.add_constraint(backend, SingleDecision(index),
                                   MOI.LessThan{Float64}(info.upper_bound))
        map_constraint!(decisions, ci, inner)
        nothing_added &= false
    end
    if info.has_fix
        ci = CI{SingleDecision, MOI.EqualTo{Float64}}(index.value)
        inner = MOI.add_constraint(backend, SingleDecision(index),
                                   MOI.EqualTo{Float64}(info.fixed_value))
        map_constraint!(decisions, ci, inner)
        nothing_added &= false
    end
    if info.binary
        ci = CI{SingleDecision, MOI.ZeroOne}(index.value)
        inner = MOI.add_constraint(backend, SingleDecision(index),
                                   MOI.ZeroOne())
        map_constraint!(decisions, ci, inner)
        nothing_added &= false
    end
    if info.integer
        ci = CI{SingleDecision, MOI.Integer}(index.value)
        inner = MOI.add_constraint(backend, SingleDecision(index), MOI.Integer())
        map_constraint!(decisions, ci, inner)
        nothing_added &= false
    end
    if info.has_start
        MOI.set(backend, MOI.VariablePrimalStart(), index,
                Float64(info.start))
        nothing_added &= false
    end
    return nothing
end

ext: jl


text:
using Registrator.RegEdit: RegEdit,
    DEFAULT_REGISTRY_URL,
    parse_registry,
    showsafe,
    registration_branch,
    get_registry
using LibGit2
using Pkg.TOML
using Pkg.Types: Project

using Test

const TEST_GITCONFIG = Dict(
    "user.name" => "RegistratorTests",
    "user.email" => "ci@juliacomputing.com",
)
const TEST_SIGNATURE = LibGit2.Signature(
    TEST_GITCONFIG["user.name"],
    TEST_GITCONFIG["user.email"],
)

@testset "RegEdit" begin

@testset "Utilities" begin
    @testset "showsafe" begin
        @test string(showsafe(3)) == "3"
        @test string(showsafe(nothing)) == "nothing"
    end

    @testset "registration_branch" begin
        example = Project(Dict(
            "name" => "Example", "version" => "1.10.2",
            "uuid" => "698ec630-83b2-4a6d-81d4-a10176273030"
        ))
        @test registration_branch(example) == "registrator/example/698ec630/v1.10.2"
    end
end

@testset "RegistryCache" begin
    @testset "get_registry" begin
        mktempdir(@__DIR__) do temp_cache_dir
            # test when registry is not in the cache and not downloaded
            cache = RegEdit.RegistryCache(temp_cache_dir)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL

            # test when registry is in the cache but not downloaded
            registry_path = RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            rm(registry_path, recursive=true, force=true)
            @test !ispath(registry_path)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL

            # test when registry is in the cache, downloaded, but mutated
            orig_hash = LibGit2.GitHash()
            LibGit2.branch!(repo, "newbranch", force=true)
            LibGit2.remove!(repo, "Registry.toml")
            LibGit2.commit(
                repo,
                "Removing Registry.toml in Registrator tests";
                author=TEST_SIGNATURE,
                committer=TEST_SIGNATURE,
            )
            @test LibGit2.GitObject(repo, "HEAD") != LibGit2.GitObject(repo, "master")
            @test ispath(registry_path)
            repo = get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)
            @test LibGit2.path(repo) == RegEdit.path(cache, DEFAULT_REGISTRY_URL)
            @test LibGit2.branch(repo) == "master"
            @test !LibGit2.isdirty(repo)
            @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL
        end
    end
end

@testset "RegistryData" begin
    blank = RegEdit.RegistryData("BlankRegistry", "d4e2f5cd-0f48-4704-9988-f1754e755b45")

    example = Project(Dict(
        "name" => "Example", "uuid" => "7876af07-990d-54b4-ab0e-23690620f79a"
    ))

    @testset "I/O" begin
        registry = """
            name = "General"
            uuid = "23338594-aafe-5451-b93e-139f81909106"
            repo = "https://github.com/JuliaRegistries/General.git"

            description = \"\"\"
            Official general Julia package registry where people can
            register any package they want without too much debate about
            naming and without enforced standards on documentation or
            testing. We nevertheless encourage documentation, testing and
            some amount of consideration when choosing package names.
            \"\"\"

            [packages]
            00701ae9-d1dc-5365-b64a-a3a3ebf5695e = { name = "BioAlignments", path = "B/BioAlignments" }
            00718b61-6157-5045-8849-3d4c4093d022 = { name = "Convertible", path = "C/Convertible" }
            0087ddc6-3964-5e57-817f-9937aefb0357 = { name = "MathOptInterfaceMosek", path = "M/MathOptInterfaceMosek" }
            """

        registry_data = parse_registry(IOBuffer(registry))
        @test registry_data isa RegEdit.RegistryData
        written_registry = sprint(TOML.print, registry_data)
        written_registry_data = parse_registry(IOBuffer(written_registry))

        @test written_registry_data == registry_data
        @test written_registry == registry
    end

    @testset "Package Operations" begin
        registry_data = copy(blank)

        @test isempty(registry_data.packages)
        @test push!(registry_data, example) == registry_data
        @test length(registry_data.packages) == 1
        @test haskey(registry_data.packages, string(example.uuid))
        @test registry_data.packages[string(example.uuid)]["name"] == "Example"
        @test registry_data.packages[string(example.uuid)]["path"] == joinpath("E", "Example")
    end
end

@testset "check_version!" begin
    import Registrator.RegEdit: RegBranch, check_version!
    import Pkg.Types: Project

    for ver in ["0.0.2", "0.3.2", "4.3.2"]
        pkg = Project(Dict("name" => "TestTools", "version" => ver))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, VersionNumber[])
        @test haskey(regbr.metadata, "warning")
        @test length(regbr.metadata["warning"]) != 0
        @test !haskey(regbr.metadata, "error")
    end

    for ver in ["0.0.1", "0.1.0", "1.0.0"]
        pkg = Project(Dict("name" => "TestTools", "version" => ver))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, VersionNumber[])
        @test !haskey(regbr.metadata, "warning")
        @test !haskey(regbr.metadata, "error")
    end

    versions_list = [v"0.0.5", v"0.1.0", v"0.1.5", v"1.0.0"]
    let    # Less than least existing version
        pkg = Project(Dict("name" => "TestTools", "version" => "0.0.4"))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, versions_list)
        @test haskey(regbr.metadata, "error")
        @test length(regbr.metadata["error"]) != 0
        @test !haskey(regbr.metadata, "warning")
    end

    let    # Existing version
        pkg = Project(Dict("name" => "TestTools", "version" => "0.0.5"))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, versions_list)
        @test haskey(regbr.metadata, "error")
        @test length(regbr.metadata["error"]) != 0
        @test !haskey(regbr.metadata, "warning")
    end

    let    # Non-existing version
        pkg = Project(Dict("name" => "TestTools", "version" => "0.0.6"))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, versions_list)
        @test !haskey(regbr.metadata, "error")
        @test !haskey(regbr.metadata, "warning")
    end

    let    # Skip a version
        pkg = Project(Dict("name" => "TestTools", "version" => "0.0.7"))
        regbr = RegBranch(pkg, "test")
        check_version!(regbr, versions_list)
        @test !haskey(regbr.metadata, "error")
        @test haskey(regbr.metadata, "warning")
        @test length(regbr.metadata["warning"]) != 0
    end
end

end


text:
import Base: +, *, isapprox

# numerical type for elements of all quantum objects
const Scalar = Union{Float64, ComplexF64}

abstract type AbstractQuantumObject end
abstract type AbstractQuantumMatrix <: AbstractQuantumObject end
abstract type AbstractQuantumVector <: AbstractQuantumMatrix end
abstract type AbstractPropagator end

data(A::AbstractQuantumObject) = A.data
size(A::AbstractQuantumObject) = size(A.data)
getindex(A::AbstractQuantumObject, i::Integer) = A.data[i]
adjoint(A::AbstractQuantumObject) = error("Adjoint is not defined for abstract type", AbstractQuantumObject)

isapprox(A::AbstractQuantumMatrix, obj2::AbstractQuantumMatrix) = size(obj1) == size(obj2) && obj1.data ≈ obj2.data
adjoint(op::AbstractQuantumMatrix) = adjoint(op.data)
ishermitian(op::AbstractQuantumMatrix) = ishermitian(op.data)
isunitary(op::AbstractQuantumMatrix) = isunitary(op.data)

norm(ψ::AbstractQuantumMatrix) = norm(ψ.data)

+(obj1::AbstractQuantumObject, obj2::AbstractQuantumObject) = error("Addition not defined for abstract type", typeof(obj1))
*(obj1::AbstractQuantumObject, obj2::AbstractQuantumObject) = error("Multiplication not defined for abstract type", typeof(obj1))

function check_dimensions_match(obj1::AbstractQuantumObject, obj2::AbstractQuantumObject)
    size(obj1)==size(obj2) || throw(DimensionMismatch("Dimensions of ", obj1, " and ", obj2, " do not match."))
end


text:
"""
    Project(task_durations, links)

 - `task_durations::NamedTuple`: a namedtuple where the keys are the names of the tasks,
   and the values indicate the timing.
    - The timing values can be a `Real` value (like 5.0) to represent a fixed duration;
      or anything sampleable to represent a duration only probabilistically known.
      i.e. any object `dur` that has `rand(rng, dur)::Real` defined on it. e.g. any
      collection or distribution object. Particularly relevent is `PertBeta(min, mode, max)`
      which is the standard used for PERT estimation.
 - `links::Vector{Pair{Symbol, Symbol}}`: a collection of pairs for all the depencencies of
   tasks. e.g. `[:a => :b, :a=>:c, :c=>:d]`. A convienent way to write this is often using
   broadcasting and matrix notation for vcat: `[:a .=> [:b, :c]; :c => :d]`.

The task names `start` and `finish` are special, in that they represent the begining
and ending milestones of the project.
They must be in the `task_durations` with duration `0`, and they respectively must be the
root and only leaf of the graph decribed by the `links`.
"""
struct Project{D<:NamedTuple}
    task_durations::D
    links::Vector{Pair{Symbol, Symbol}}

    function Project(task_durations::D, links) where D
        proj = new{D}(task_durations, links)
        validate_project(proj)
        return proj
    end
end

function validate_project(proj::Project)
    validate_milestones(proj)
    validate_links(proj)
end

function validate_milestones(proj)
    haskey(proj.task_durations, :start) || throw(MissingMilestoneException(:start))
    haskey(proj.task_durations, :finish) || throw(MissingMilestoneException(:finish))

    after_finish = [b for (a, b) in proj.links if a === :finish]
    isempty(after_finish) || throw(PostFinishTasksException(after_finish))

    before_start = [a for (a, b) in proj.links if b === :start]
    isempty(before_start) || throw(PreStartTasksException(before_start))

    validate_milestone_duration(proj, :start)
    validate_milestone_duration(proj, :finish)
end

function validate_milestone_duration(proj, milestone_name)
    dur = proj.task_durations[milestone_name]
    if length(dur) != 1 || !iszero(only(dur))
        throw(BadMilestoneDurationException(milestone_name, dur))
    end
end

"""
    validate_links(proj)

Checks that:
    - the links in `proj` form a directed acyclic graph
    - that every end-point in the `links` is on `task_durations`
    - that every task in `task_durations` occurs in `links`
"""
function validate_links(proj)
    children_lookup = MultiDict(proj.links)
    all_seen = Set{Symbol}()  # for checking at the end that everything was covered

    function inner(node, path)
        push!(all_seen, node)
        node == :finish && return
        children = get(children_lookup, node, nothing)
        children === nothing && throw(DisconnectedPathException(path))
        for child in children
            push!(path, child)
            child ∈ path[1:end-1] && throw(CyclicPathException(path))
            inner(child, path)
            pop!(path)  # take it off before we do the next
        end
    end

    inner(:start, Symbol[])

    missing_links = setdiff(keys(proj.task_durations), all_seen)
    isempty(missing_links) || throw(UnlinkedTasksException(missing_links))

    missing_durations = setdiff(all_seen, keys(proj.task_durations))
    isempty(missing_durations) || throw(UnspecifiedTasksException(missing_durations))

    return nothing
end

struct MissingMilestoneException <: Exception
    milestone
end

struct BadMilestoneDurationException <: Exception
    milestone_name
    duration
end
function Base.showerror(io::IO, err::BadMilestoneDurationException)
    print(io, "The Milestone $(err.milestone_name) has duration $(err.duration). ")
    print(io, "milestones must have singlton duration, that is zero valued.")
end

struct PostFinishTasksException <: Exception
    task
end
struct PreStartTasksException <: Exception
    task
end

struct CyclicPathException <: Exception
    cycle
end
struct DisconnectedPathException <: Exception
    path
end

struct UnlinkedTasksException <: Exception
    tasks
end

struct UnspecifiedTasksException <: Exception
    tasks
end


text:
## ------------------------------------------------------------------
# Property handling
function Base.getproperty(obj::RBObject, p::Symbol)
    props = _get_properties(obj)
    return haskey(props, p) ? props[p] : getfield(obj, p)
end

function Base.getproperty(obj::RBObject, p::Symbol, dflt)
    props = _get_properties(obj)
    get(props, p, dflt)
end

function getproperty!(obj::RBObject, p::Symbol, dflt)
    props = _get_properties(obj)
    get!(props, p, dflt)
end

getproperty!(f::Function, obj::RBObject, p::Symbol) = getproperty!(obj, p, f())

function Base.setproperty!(obj::RBObject, p::Symbol, value)
    props = _get_properties(obj)
    props[p] = value
    return value
end

function Base.propertynames(obj::RBObject, private::Bool = false)
    return collect(keys(_get_properties(obj)))
end

## ------------------------------------------------------------------
# common metas
get_label(obj::RBObject) = getproperty(obj, :label, "")

srcfile!(obj::RBObject, file::String) = setproperty!(obj, :srcfile, realpath(strip(file)))
srcfile(obj::RBObject)::String = getproperty(obj, :srcfile, "")

srcline!(obj::RBObject, line::Int) = setproperty!(obj, :srcline, line)
srcline(obj::RBObject)::Int = getproperty(obj, :srcline, -1)

get_parent(obj::RBObject) = getproperty(obj, :parent, nothing)
set_parent!(obj::RBObject, parent::RBObject) = setproperty!(obj, :parent, parent)

get_book(d::RBObject) = get_book(get_parent(d))
bookdir(q::RBObject) = bookdir(get_book(q))

get_tags(obj::RBObject) = getproperty!(() -> Set{String}(), obj, :tags)
add_tag!(obj::RBObject,  t::String, ts::String...) = (_push_csv!(get_tags(obj), t, ts...); obj)

get_todos(obj::RBObject) = getproperty!(() -> Set{String}(), obj, :todos)
add_todo!(obj::RBObject,  td::String, tds::String...) = (_push_csv!(get_todos(obj), td, tds...); obj)

function localpath(obj::RBObject; rel = false)
    file = srcfile(obj)
    file = rel ? relpath(file) : file
    line = srcline(obj)
    linestr = (line == -1) ? "" : string(":", line)
    return string(file, linestr)
end

localrelpath(obj::RBObject) = localpath(obj; rel = true)


text:
"""
    SubsetEstimates

BLB linear mixed model estimates from one subset.
"""
struct SubsetEstimates{T <: LinearAlgebra.BlasReal}
    # blb parameter
    n_boots::Int64
    # dimension of parameter estimates, for pre-allocation
    p::Int64
    q::Int64
    # blb results
    βs::Matrix{T}
    Σs::Array{T}
    σ²s::Vector{T}
end

"""
    SubsetEstimates(n_boots, p, q)

Constructor for `SubsetEstimates`.
"""
function SubsetEstimates(n_boots::Int64, p::Int64, q::Int64)
    # vector of vector/matrices approach
    # βs  = [Vector{Float64}(undef,  p) for _ in 1:n_boots] #Vector{Vector{Float64}}()
    # Σs  = [Matrix{Float64}(undef, q, q) for _ in 1:n_boots] #Vector{Matrix{Float64}}()

    # multi-dimensional array approach
    βs  = Matrix{Float64}(undef, n_boots, p)
    Σs  = Array{Float64}(undef, q, q, n_boots) # 3d array
    σ²s = Vector{Float64}(undef, n_boots)
    SubsetEstimates(n_boots, p, q, βs, Σs, σ²s)
end

"""
    blbEstimates

BLB linear mixed model estimates, which contains blb parameters and a vector of `SubsetEstimates`.
"""
struct blbEstimates{T <: LinearAlgebra.BlasReal}
    # blb parameters
    n_subsets::Int64
    subset_size::Int64
    n_boots::Int64
    fenames::Vector{String}
    renames::Vector{String}
    # subset estimates from all subsets
    all_estimates::Vector{SubsetEstimates{T}}
    runtime::Vector{Float64}
end

"""
    save_bootstrap_result!(subset_estimates, β, Σ, σ²)

Save the result from one bootstrap iteration to subset_estimates

## Positional arguments
* `subset_estimates`: an object of the SubsetEstimates type
* `i`: update the ith element of βs, Σs and σ²s
* `β`: parameter estimates for fixed effects
* `Σ`: the covariance matrix of variance components
* `σ²`: estimate of error variance
"""
function save_bootstrap_result!(
    subset_estimates::SubsetEstimates{T},
    i::Int64,
    β::Vector{T},
    Σ::Matrix{T},
    σ²::T
    ) where T <: BlasReal
    # using the vector of vectors/matrices approach
    # copyto!(subset_estimates.βs[i], β)
    # copyto!(subset_estimates.Σs[i], Σ)

    # using the multi-dimensional array approach
    subset_estimates.βs[i, :] = β
    subset_estimates.Σs[:, :, i] = Σ
    subset_estimates.σ²s[i] = σ²
end

"""
    blb_one_subset(rng, m; n_boots, solver, verbose, nonparametric_boot)

Performs Bag of Little Bootstraps on a subset.

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `m`: an object of the blblmmModel type

# Keyword arguments
- `n_boots`: number of bootstrap iterations. Default to 1000
- `solver`: solver for the optimization problem.
- `verbose`: Bool, whether to print bootstrap progress
- `nonparametric_boot`: Bool, whether to use nonparametric bootstrap

# Values
- `subset_estimates`: an object of type `SubsetEstimates`
"""
function blb_one_subset(
    rng::Random.AbstractRNG,
    m::blblmmModel{T};
    n_boots::Int = 1000,
    solver,
    verbose::Bool,
    nonparametric_boot::Bool
    ) where T <: BlasReal

    # Initalize model parameters
    init_ls!(m)
    # Fit LMM on the subset
    fit!(m, solver)
    # Initalize an instance of SubsetEstimates type for storing results
    subset_estimates = SubsetEstimates(n_boots, m.p, m.q)

    # construct the simulator
    if nonparametric_boot
        verbose && print("Using Non-parametric Bootstrap\n")
        simulator = NonparametricBootSimulator(m)
    else
        verbose && print("Using Parametric Bootstrap\n")
        simulator = ParametricBootSimulator(m)
    end

    # Bootstrapping
    @inbounds for k = 1:n_boots
        if verbose
            flush(stdout)
            print("Bootstrap iteration ", k, "\n")
        end

        if !nonparametric_boot
            # Parametric bootstrapping. Updates m.data[i].y for all i
            simulate!(rng, m, simulator)
        end

        # Get weights by drawing N i.i.d. samples from multinomial
        rand!(simulator.mult_dist, simulator.ns)

        # Update weights in blblmmModel
        update_w!(m, simulator.ns)

        # Fit model on the bootstrap sample
        fit!(m, solver)

        # Save estimates
        save_bootstrap_result!(subset_estimates, k, m.β, m.Σ, m.σ²[1])

        # Reset model parameter to subset estimates because
        # using the bootstrap estimates from each iteration may be unstable.
        copyto!(m.β, simulator.β_subset)
        copyto!(m.Σ, simulator.Σ_subset)
        copyto!(m.σ², simulator.σ²_subset)
    end
    return subset_estimates
end
blb_one_subset(m::blblmmModel; n_boots::Int = 1000, solver, verbose::Bool, nonparametric_boot::Bool) =
    blb_one_subset(Random.GLOBAL_RNG, m; n_boots = n_boots, solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)

# function replicate(f::Function, n::Integer; use_threads=false)
#     if use_threads
#         # no macro version yet: https://github.com/timholy/ProgressMeter.jl/issues/143
#         # get the type
#         rr = f()
#         # pre-allocate
#         results = [rr for _ in Base.OneTo(n)]
#         Threads.@threads for idx = 2:n
#             results[idx] = f()
#         end
#     else
#         results = [f() for _ in Base.OneTo(n)]
#     end
#     results
# end

# """
#     blb_one_subset(m; n_boots, solver, verbose)

# Performs Bag of Little Bootstraps on a subset.

# # Positional arguments
# - `m`: an object of the blblmmModel type

# # Keyword arguments
# - `n_boots`: number of bootstrap iterations. Default to 1000
# - `solver`: solver for the optimization problem.
# - `verbose`: Bool, whether to print bootstrap progress (percentage completion)

# # Values
# - `subset_estimates`: an object of the SubsetEstimates type
# """
# function blb_one_subset_2(
#     rng::Random.AbstractRNG,
#     m::blblmmModel{T};
#     n_boots::Int64 = 1000,
#     solver = Ipopt.IpoptSolver(print_level = 0),
#     verbose::Bool = false,
#     use_threads::Bool = false
#     ) where T <: Real

#     print("m.β = ", m.β, "\n")
#     print("m.data[1].y = ", m.data[1].y, "\n")

#     βsc, Σsc, σ²sc, p, q, mc = similar(m.β), similar(m.Σ), similar(m.σ²), m.p, m.q, deepcopy(m)

#     init_ls!(mc, verbose)
#     print("mc.β = ", mc.β, "\n")
#     print("mc.data[1].y = ", mc.data[1].y, "\n")
#     simulatorc = Simulator(mc)

#     mc_threads = [m]
#     simulatorc_threads = [simulatorc]
#     βsc_threads = [βsc]
#     Σsc_threads = [Σsc]
#     σ²sc_threads = [σ²sc]

#     if use_threads
#         Threads.resize_nthreads!(mc_threads)
#         Threads.resize_nthreads!(βsc_threads)
#         Threads.resize_nthreads!(Σsc_threads)
#         Threads.resize_nthreads!(σ²sc_threads)
#     end

#     rnglock = ReentrantLock()
#     samp = replicate(n_boots, use_threads = use_threads) do
#         mod = mc_threads[Threads.threadid()]
#         init_ls!(mod, verbose)
#         print("mod.β = ", mod.β, "\n")
#         print("mod.data[1].y = ", mod.data[1].y, "\n")
#         simulator = simulatorc_threads[Threads.threadid()]
#         local βsc = βsc_threads[Threads.threadid()]
#         local Σsc = Σsc_threads[Threads.threadid()]
#         local σ²sc = σ²sc_threads[Threads.threadid()]
#         lock(rnglock)
#         simulate!(rng, mod, simulator)
#         rand!(simulator.mult_dist, simulator.ns)
#         update_w!(mod, simulator.ns)
#         unlock(rnglock)
#         fit!(mod, solver)
#         (
#             #  β = copyto!(βsc, mod.β),
#             #  Σ = copyto!(Σsc, mod.Σ),
#             #  σ² = copyto!(σ²sc, mod.σ²)
#          β = copyto!(βsc, mod.β),
#          Σ = copyto!(Σsc, mod.Σ),
#          σ² = copyto!(σ²sc, mod.σ²)
#         )
#     end

#     return samp
# end

"""
    blblmmobs(datatable)

Constructor for the type `blblmmObs`

# Positional arguments
- `data_obs`: a table object that is compatible with Tables.jl
- `feformula`: the formula for the fixed effects
- `reformula`: the formula for the fixed effects
"""
function blblmmobs(data_obs, feformula::FormulaTerm, reformula::FormulaTerm)
    y, X = StatsModels.modelcols(feformula, data_obs)
    Z = StatsModels.modelmatrix(reformula, data_obs)
    return blblmmObs(y, X, Z)
end
blblmmobs(feformula::FormulaTerm, reformula::FormulaTerm) = data_obs -> blblmmobs(data_obs, feformula, reformula)

"""
    count_levels(data_columns, cat_names)

Count the number of levels of each categorical variable.

# Positional arguments
- `data_columns`: an object of the AbstractColumns type
- `cat_names`: a character vector of the names of the categorical variables

# Values
- `cat_levels`: a dictionary that contains the number of levels of each categorical variable.
"""
function count_levels(data_columns, cat_names::Vector{String})
    cat_levels = Dict{String, Int64}()
    @inbounds for cat_name in cat_names
        cat_levels[cat_name] = length(countmap(Tables.getcolumn(data_columns, Symbol(cat_name))))
    end
    return cat_levels
end
count_levels(cat_names::Vector{String}) = data_columns -> count_levels(data_columns, cat_names)

"""
    subsetting!(rng, subset_id, data_columns, id_name, unique_id, cat_names, cat_levels)

Draw a subset from the full dataset. The IDs of the subset is stored in `subset_id`

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `subset_id`: a vector for storing the IDs of the subset
- `data_columns`: an object of the AbstractColumns type, or a DataFrame
- `unique_id`: a vector of the unique ID in the full data set
- `cat_names`: a character vector of the names of the categorical variables
- `cat_levels`: a dictionary that contains the number of levels of each categorical variable
"""
function subsetting!(
    rng::Random.AbstractRNG,
    subset_id::Vector,
    data_columns,
    id_name::Symbol,
    unique_id::Vector,
    cat_names::Vector{String},
    cat_levels::Dict
    )
    if length(cat_names) == 0
        sample!(rng, unique_id, subset_id; replace = false)
    else
        # there are categorical variables
        good_subset = false
        while !good_subset
            # Sample from the full dataset
            sample!(rng, unique_id, subset_id; replace = false)
            # subset_indices = LinearIndices(id)[findall(in(blb_id_unique), id)]
            cat_levels_subset = data_columns |>
                TableOperations.filter(x -> Tables.getcolumn(x, Symbol(id_name)) .∈ Ref(Set(subset_id))) |>
                Tables.columns |>
                count_levels(cat_names)
            # If the subset levels do not match the full dataset levels,
            # skip the current iteration and take another subset
            if cat_levels_subset == cat_levels
                good_subset = true
            end
        end
    end
    sort!(subset_id)
end
subsetting!(subset_id::Vector, data_columns, id_name::Symbol, unique_id::Vector, cat_names::Vector{String}, cat_levels::Dict) =
    subsetting!(Random.GLOBAL_RNG, subset_id, data_columns, id_name, unique_id, cat_names, cat_levels)

"""
    x_in_y(x, sorted_y, k)

Test whether x is in a sorted vector y or not.

# Positional arguments
- `x`: a value of type T
- `y`: a sorted vector of type T
"""
function x_in_y(x::T, sorted_y::Vector{T}, k::Int) where T
    # test whether x is in sorted_y
    index = searchsortedfirst(sorted_y, x)
    return index ≤ k && sorted_y[index] == x
end

"""
    blb_full_data(rng, datatable; feformula, reformula, id_name, cat_names, subset_size, n_subsets, n_boots, solver, verbose,  nonparametric_boot)

Performs Bag of Little Bootstraps on the full dataset

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `datatable`: a data table type that is compatible with Tables.jl.

# Keyword arguments
- `feformula`: model formula for the fixed effects.
- `reformula`: model formula for the random effects.
- `id_name`: name of the cluster identifier variable. String.
- `cat_names`: a vector of the names of the categorical variables.
- `subset_size`: number of clusters in the subset.
- `n_subsets`: number of subsets.
- `n_boots`: number of bootstrap iterations. Default to 1000
- `solver`: solver for the optimization problem.
- `verbose`: Bool, whether to print bootstrap progress (percentage completion)
- `nonparametric_boot`: Bool, whether to use Nonparametric bootstrap

# Values
- `result`: an object of the blbEstimates type
"""
function blb_full_data(
    rng::Random.AbstractRNG,
    datatable;
    feformula::FormulaTerm,
    reformula::FormulaTerm,
    id_name::String,
    cat_names::Vector{String} = Vector{String}(),
    subset_size::Int,
    n_subsets::Int = 10,
    n_boots::Int = 1000,
    solver = Ipopt.IpoptSolver(print_level=0, mehrotra_algorithm = "yes", warm_start_init_point = "yes"),
    verbose::Bool = false,
    # use_threads::Bool = false,
    nonparametric_boot::Bool = true
    )
    # Create Tables.Columns type for subsequent processing
    datatable_cols = Tables.columns(datatable) # This step does not allocate memory
    # Get the unique ids, which will be used for subsetting
    typeof(id_name) <: String && (id_name = Symbol(id_name))
    unique_id = unique(Tables.getcolumn(datatable_cols, id_name))
    # Same performance as: unique(Tables.columntable(TableOperations.select(datatable, id_name))[id_name])
    N = length(unique_id) # number of individuals/clusters in the full dataset
    if N < subset_size
        error(string("The subset size should not be bigger than the total number of clusters. \n",
                        "Total number of clusters = ", N, "\n",
                        "Subset size = ", subset_size, "\n"))
    end

    # apply df-wide schema
    feformula = apply_schema(feformula, schema(feformula, datatable))
    reformula = apply_schema(reformula, schema(reformula, datatable))

    #fenames = coefnames(feformula)[2]
    if typeof(coefnames(feformula)[2]) == String
        fenames = [coefnames(feformula)[2]]
    else
        fenames = coefnames(feformula)[2]
    end

    if typeof(coefnames(reformula)[2]) == String
        renames = [coefnames(reformula)[2]]
    else
        renames = coefnames(reformula)[2]
    end

    # f = @formula(y ~ 1+a)
    # f = apply_schema(f, schema(f, df))
    # if typeof(coefnames(f)[2]) == String
    #     renames = [coefnames(f)[2]]
    # else
    #     renames = coefnames(f)[2]
    # end

    # By chance, some factors of a categorical variable may not show up in a subset.
    # To make sure this does not happen, we
    # 1. count the number of levels of each categorical variable in the full data;
    # 2. for each sampled subset, we check whether the number of levels match.
    # If they do not match, the subset is resampled.
    if length(cat_names) > 0
        cat_levels = count_levels(datatable_cols, cat_names)
    else
        cat_levels = Dict()
    end

    # Initialize an array to store the unique IDs for the subset
    subset_id = Vector{eltype(unique_id)}(undef, subset_size)
    # Initialize a vector of SubsetEstimates for storing estimates from subsets
    all_estimates = Vector{SubsetEstimates{Float64}}(undef, n_subsets)
    # Initialize a vector of the blblmmObs objects
    obsvec = Vector{blblmmObs{Float64}}(undef, subset_size)
    # Initialize a vector for storing the runtime
    # Currently, runtime doesn't make sense when there are multiple workers
    runtime = Vector{Float64}(undef, n_subsets)
    # do grouping once, and collect
    datatable_grouped = datatable |> @groupby(_.id) |> collect

    if Distributed.nworkers() > 1
        # Using multi-processing
        futures = Vector{Future}(undef, n_subsets)
        # wks_schedule assigns each subset to a worker
        wks_schedule = Vector{Int}(undef, n_subsets)
        # nwks = Distributed.nworkers()
        wk = 2 # note that the worker number starts from 2
        wk_max = maximum(workers())
        @inbounds for i = 1:n_subsets
            wks_schedule[i] = wk
            wk == wk_max ? wk = 2 : wk += 1
        end
        flush(stdout)
        print("wks_schedule = ", wks_schedule, "\n")

        @inbounds for j = 1:n_subsets
            time0 = time_ns()
            # Take a subset
            subsetting!(subset_id, datatable_cols, id_name, unique_id, cat_names, cat_levels)
            # Construct blblmmObs objects
            # if use_groupby
            obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(blblmmobs(_, feformula, reformula)) |> collect |> Array{blblmmObs{Float64}, 1}
            # else
            #     Threads.@threads for i in 1:subset_size
            #         obsvec[i] = datatable_cols |>
            #             TableOperations.filter(x -> Tables.getcolumn(x, id_name) == subset_id[i]) |>
            #             Tables.columns |>
            #             blblmmobs(feformula, reformula)
            #     end
            # end
            # Construct the blblmmModel type
            m = blblmmModel(obsvec, fenames, renames, N)#, use_threads)
            # Process this subset on worker "wks_schedule[j]"
            futures[j] = remotecall(blb_one_subset, wks_schedule[j], rng, m;
                                    n_boots = n_boots, solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)
            # A remote call returns a Future to its result. Remote calls return immediately;
            # the process that made the call proceeds to its next operation while the remote call happens somewhere else.
            # You can wait for a remote call to finish by calling wait on the returned Future,
            # and you can obtain the full value of the result using fetch.
            runtime[j] = (time_ns() - time0) / 1e9
        end
        @inbounds for j = 1:n_subsets
            # Fetch results from workers
            all_estimates[j] = fetch(futures[j])
        end
    else
        # Not using multi-processing
        @inbounds for j = 1:n_subsets
            time0 = time_ns()
            # Take a subset
            subsetting!(subset_id, datatable_cols, id_name, unique_id, cat_names, cat_levels)
            # Construct blblmmObs objects
            # if use_groupby
            obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(blblmmobs(_, feformula, reformula)) |> collect |> Array{blblmmObs{Float64}, 1}
            # else
            #     Threads.@threads for i in 1:subset_size
            #         obsvec[i] = datatable_cols |>
            #             TableOperations.filter(x -> Tables.getcolumn(x, id_name) == subset_id[i]) |>
            #             Tables.columns |>
            #             blblmmobs(feformula, reformula)
            #     end
            # end
            # Construct the blblmmModel type
            m = blblmmModel(obsvec, fenames, renames, N)#, use_threads)
            all_estimates[j] = blb_one_subset(rng, m; n_boots = n_boots, solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)
            runtime[j] = (time_ns() - time0) / 1e9
        end
    end

    # Create a blbEstimates instance for storing results from all subsets
    result = blbEstimates{Float64}(n_subsets, subset_size, n_boots, fenames, renames, all_estimates, runtime)
    return result
end

blb_full_data(datatable; feformula::FormulaTerm, reformula::FormulaTerm, id_name::String,
                cat_names::Vector{String} = Vector{String}(), subset_size::Int, n_subsets::Int = 10, n_boots::Int = 200,
                solver = Ipopt.IpoptSolver(), verbose::Bool = false, nonparametric_boot::Bool = true) =
    blb_full_data(Random.GLOBAL_RNG, datatable; feformula = feformula, reformula = reformula, id_name = id_name,
                    cat_names = cat_names, subset_size = subset_size, n_subsets = n_subsets, n_boots = n_boots,
                    solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)

"""
    confint(subset_ests, level)

Calculate confidence intervals using estimates from one subset.

# Positional arguments
- `subset_ests`: an object of type `SubsetEstimates`
- `level`: confidence level, usually set to 0.95
"""
function confint(subset_ests::SubsetEstimates, level::Real)
    ci_βs = Matrix{Float64}(undef, subset_ests.p, 2) # p-by-2 matrix
    for i in 1:subset_ests.p
        ci_βs[i, :] = StatsBase.percentile(view(subset_ests.βs, :, i), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
    end
    q = subset_ests.q
    ci_Σs = Matrix{Float64}(undef, ◺(q), 2)
    k = 1
    # For Σ, we get the CI for the diagonals first, then the upper off-diagonals
    @inbounds for i in 1:q
        ci_Σs[k, :] = StatsBase.percentile(view(subset_ests.Σs, i, i, :), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
        k += 1
    end
    @inbounds for i in 1:q, j in (i+1):q
        ci_Σs[k, :] = StatsBase.percentile(view(subset_ests.Σs, i, j, :), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
        k += 1
    end
    ci_σ²s = reshape(StatsBase.percentile(subset_ests.σ²s, 100 * [(1 - level) / 2, 1 - (1-level) / 2]), 1, 2)
    return ci_βs, ci_Σs, ci_σ²s
end

"""
    confint(blb_ests, level)

Calculate confidence intervals using estimates from all subsets.

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
- `level`: confidence level, usually set to 0.95
"""
function confint(blb_ests::blbEstimates, level::Real)
    # initialize arrays for storing the CIs from each subset
    cis_βs = Array{Float64}(undef, blb_ests.all_estimates[1].p, 2, blb_ests.n_subsets)
    cis_Σs = Array{Float64}(undef, ◺(blb_ests.all_estimates[1].q), 2, blb_ests.n_subsets)
    cis_σ²s = Array{Float64}(undef, 1, 2, blb_ests.n_subsets)
    @inbounds for i in 1:blb_ests.n_subsets
        cis_βs[:, :, i], cis_Σs[:, :, i], cis_σ²s[:, :, i] = confint(blb_ests.all_estimates[i], level)
    end
    ci_β  = mean(cis_βs, dims = 3)[:, :, 1]
    ci_Σ  = mean(cis_Σs, dims = 3)[:, :, 1]
    ci_σ² = mean(cis_σ²s, dims = 3)[:, :, 1]
    return ci_β, ci_Σ, ci_σ²
end
confint(blb_ests::blbEstimates) = confint(blb_ests, 0.95)

"""
    fixef(blb_ests)

Calculate BLB fixed effect estimates, which are averages of fixed effect estimates from all subsets

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
"""
# returns fixed effect estimates
function fixef(blb_ests::blbEstimates)
    means_βs = Matrix{Float64}(undef, blb_ests.n_subsets, blb_ests.all_estimates[1].p) # n-by-p matrix
    for i in 1:blb_ests.n_subsets
        means_βs[i, :] = mean(blb_ests.all_estimates[i].βs, dims = 1)
    end
    mean_β = mean(means_βs, dims = 1)
    return mean_β
end

"""
    vc(blb_ests)

Calculate BLB variance components estimates, which are averages of variance component estimates from all subsets

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
"""
# returns variance components estimates
function vc(blb_ests::blbEstimates)
    means_Σs = Array{Float64}(undef, blb_ests.all_estimates[1].q, blb_ests.all_estimates[1].q, blb_ests.n_subsets)
    # print(means_Σs, "\n")
    means_σ²s = Vector{Float64}(undef, blb_ests.n_subsets)
    for i in 1:blb_ests.n_subsets
        # print(mean(blb_ests.all_estimates[i].Σs, dims = 3), "\n")
        means_Σs[:,:,i] = mean(blb_ests.all_estimates[i].Σs, dims = 3)
        # view(means_Σs, :, :, i) .= mean(blb_ests.all_estimates[i].Σs, dims = 3)
        means_σ²s[i] = mean(blb_ests.all_estimates[i].σ²s)
    end
    mean_Σ = mean(means_Σs, dims = 3)[:, :, 1]
    mean_σ² = mean(means_σ²s)
    return mean_Σ, mean_σ²
end

function StatsBase.coeftable(ests::Vector, ci::Matrix, varnames::Vector{String})
    CoefTable(
        hcat(reshape(ests, :, 1), ci[:, 1], ci[:, 2]),
        ["Estimate", "CI Lower", "CI Upper"],
        varnames
        # 4 # pvalcol
    )
end

function vectorize_Σ(Σ::Matrix)
    q  = size(Σ, 1)
    v = Vector{Float64}(undef, ◺(q))
    idx = 1
    # Extract the diagonals first
    @inbounds for i in 1:q
        v[idx] = Σ[i, i]
        idx += 1
    end
    # Extract the upper off-diagonals
    @inbounds for i in 1:q, j in (i+1):q
        v[idx] = Σ[i, j]
        idx += 1
    end
    return v
end

function vc_names(blb_ests::blbEstimates)
    q = length(blb_ests.renames)
    re_names = Vector{String}(undef, ◺(q))
    idx = 1
    # Extract the diagonals' names first
    @inbounds for i in 1:q
        re_names[idx] = blb_ests.renames[i]
        idx += 1
    end
    # Extract the upper off-diagonals' names
    @inbounds for i in 1:q, j in (i+1):q
        re_names[idx] = string(blb_ests.renames[i], " : ", blb_ests.renames[j])
        idx += 1
    end
    return vcat(re_names, "Residual")
end

function Base.show(io::IO, blb_ests::blbEstimates)
    println("Bag of Little Boostrap (BLB) for linear mixed models.")
    println("Number of subsets: ", blb_ests.n_subsets)
    println("Number of grouping factors per subset: ", blb_ests.subset_size)
    println("Number of bootstrap samples per subset: ", blb_ests.n_boots)
    println("Confidence interval level: 95%")
    println(io)

    # calculate all CIs and mean estimates
    ci_β, ci_Σ, ci_σ² = confint(blb_ests)
    mean_β = fixef(blb_ests)
    mean_Σ, mean_σ² = vc(blb_ests)

    println("Variance Components parameters")
    show(io, coeftable(vcat(vectorize_Σ(mean_Σ), mean_σ²), vcat(ci_Σ, ci_σ²), vc_names(blb_ests)))
    println(io)
    println(io)

    println("Fixed-effect parameters")
    show(io, coeftable(vec(mean_β), ci_β, blb_ests.fenames))
end


text:
@doc raw"""
    ℛ(z::Number, tableau::AbstractButcherTableau) :: Number
    ℛ(z::Number, solver::AbstractRungeKuttaSolver) :: Number

returns the stability function of an [`AbstractRungeKuttaSolver`](@ref):
```math
    R(z) = \frac{\det(I - z(A - \mathbb{1}b^\intercal))}{\det(I - zA)}.
```
"""
function ℛ(z::Number, tableau::AbstractButcherTableau)
    @↓ A, b, s = tableau
    e = ones(s)
    return det(I - z * (A - e * b')) / det(I - z * A)
end
ℛ(z::Number, solver::AbstractRungeKuttaSolver) = ℛ(z, solver.tableau)

"""
    ℛ(Z::AbstractMatrix, tableau::AbstractButcherTableau) :: AbstractMatrix
    ℛ(Z::AbstractMatrix, solver::AbstractRungeKuttaSolver) :: AbstractMatrix

returns the stability function of an [`AbstractRungeKuttaSolver`](@ref).
"""
function ℛ(Z::AbstractMatrix, tableau::AbstractButcherTableau)
    @↓ A, b, s = tableau
    e = ones(s)
    tmp = kron(e, Z)
    tmp = (I - kron(A, Z)) \ tmp
    tmp = kron(b', Matrix(1.0I, size(Z)...)) * tmp
    return I + tmp
end
ℛ(Z::AbstractMatrix, solver::AbstractRungeKuttaSolver) = ℛ(Z, solver.tableau)


text:
__precompile__()

module NearestNeighbors

using Distances
import Distances: Metric, result_type, eval_reduce, eval_end, eval_op, eval_start, evaluate

using StaticArrays
import Base.show

export NNTree, BruteTree, KDTree, BallTree, DataFreeTree
export knn, inrange # TODOs? , allpairs, distmat, npairs
export injectdata

export Euclidean,
       Cityblock,
       Minkowski,
       Chebyshev,
       Hamming

# Change this to enable debugging
const DEBUG = false

abstract type NNTree{V <: AbstractVector,P <: Metric} end

const MinkowskiMetric = Union{Euclidean,Chebyshev,Cityblock,Minkowski}

function check_input(::NNTree{V1}, ::AbstractVector{V2}) where {V1, V2 <: AbstractVector}
    if length(V1) != length(V2)
        throw(ArgumentError(
            "dimension of input points:$(length(V2)) and tree data:$(length(V1)) must agree"))
    end
end

function check_input(::NNTree{V1}, point::AbstractVector{T}) where {V1, T <: Number}
    if length(V1) != length(point)
        throw(ArgumentError(
            "dimension of input points:$(length(point)) and tree data:$(length(V1)) must agree"))
    end
end

get_T(::Type{T}) where {T <: AbstractFloat} = T
get_T(::T) where {T} = Float64

include("debugging.jl")
include("evaluation.jl")
include("tree_data.jl")
include("datafreetree.jl")
include("knn.jl")
include("inrange.jl")
include("hyperspheres.jl")
include("hyperrectangles.jl")
include("utilities.jl")
include("brute_tree.jl")
include("kd_tree.jl")
include("ball_tree.jl")
include("tree_ops.jl")

end # module


text:
"""
`module JAC.AngularMomentum`
    ... a submodel of JAC that contains various methods to calculate (standard) symbols, coefficients and functions from
        the theory of angular momentum and Racah's algebra.
"""
module AngularMomentum

    using SpecialFunctions, ..Basics
    using GSL: sf_coupling_3j, sf_coupling_6j, sf_coupling_9j

    """
    `AngularMomentum.allowedKappaSymmetries(syma::LevelSymmetry, symb::LevelSymmetry)`
        ... to determine all allowed single-electron symmetries/partial waves kappa (l,j) that can be coupled to the given
            level symmetries. A list::Array{Int64,1} of kappa-values is returned.
    """
    function  allowedKappaSymmetries(syma::LevelSymmetry, symb::LevelSymmetry)
        kappaList = Int64[];    JList = Basics.oplus(syma.J, symb.J);    kappaMax = syma.J.num + symb.J.num + 2
        for  kappa = -kappaMax:kappaMax
             kappa == 0  &&  continue
             if  syma.parity == Basics.plus   la = 0               else   la = 1       end
             if  symb.parity == Basics.plus   lb = 0               else   lb = 1       end
             if  kappa < 0                 l  = abs(kappa) -1   else   l  = kappa   end
             isodd( la + lb + l )          && continue
             for  j in JList
                 if  j == AngularJ64(abs(kappa) - 1//2)   push!( kappaList, kappa)  end
             end
        end
        return( kappaList )
    end

   """
    `AngularMomentum.allowedMultipoleSymmetries(syma::LevelSymmetry, multipole::EmMultipole)`
        ... to determine all allowed level symmetries for which the given multipole can give rise to a non-zero (transition)
            amplitude; a symList::Array{LevelSymmetry,1} is returned.
    """
    function  allowedMultipoleSymmetries(syma::LevelSymmetry, multipole::EmMultipole)
        symList = LevelSymmetry[];    JList = Basics.oplus(syma.J, AngularJ64(multipole.L) )
        for  J in JList
            if      parityEmMultipolePi(syma.parity, multipole, Basics.plus)     push!( symList, LevelSymmetry(J, Basics.plus) )
            elseif  parityEmMultipolePi(syma.parity, multipole, Basics.minus)    push!( symList, LevelSymmetry(J, Basics.minus) )
            else    error("stop a")
            end
        end
        return( symList )
    end

    """
    `AngularMomentum.allowedTotalSymmetries(syma::LevelSymmetry, kappa::Int64)`
        ... to determine all allowed total symmetries J^P that can be constructed by coupling a partial wave kappa (l,j)
            to the given level symmetry syma. A list::Array{LevelSymmetry,1} of total symmetries is returned.
    """
    function allowedTotalSymmetries(syma::LevelSymmetry, kappa::Int64)
        symtList = LevelSymmetry[]
        if  kappa < 0    l  = abs(kappa) -1   else   l  = kappa   end
        j = AngularJ64( abs(kappa) - 1//2 )
        #
        JList = Basics.oplus(syma.J, j)
        for J in JList
            if    iseven(l)   push!( symtList, LevelSymmetry(J, syma.parity) )
            else              push!( symtList, LevelSymmetry(J, Basics.invertParity(syma.parity)) )
            end
        end
        return( symtList )
    end

    """
    `AngularMomentum.bracket(jList::Array{AngularJ64,1})`
        ... to compute the bracket [a, b, c, ... ] = (2a+1) * (2b+1) * (2b+1) * ... of the given angular momenta.
            A value::Int64 is returned.
    """
    function bracket(jList::Array{AngularJ64,1})
        value = 1;    for  j in jList    value = value * (AngularMomentum.twoJ(j) + 1)    end
        return( value )
    end

    """
    `AngularMomentum.ChengI

    + (kapa::Int64, ma::AngularM64, kapb::Int64, mb::AngularM64, L::AngularJ64, M::AngularM64)`
        ... evaluates the angular I (kappa m, kappa' m', LM) integral as defined by Cheng, NATO summerschool (198x),
            Eq. (A4.5), and including the full magnetic (orientational) dependence. A value::Float64 is returned.
    """
    function ChengI(kapa::Int64, ma::AngularM64, kapb::Int64, mb::AngularM64, L::AngularJ64, M::AngularM64)
        ja = Basics.subshell_j( Subshell(9, kapa) );    jb = Basics.subshell_j( Subshell(9, kapb) )   # Use principal QN n=9 arbitrarely here
        la = Basics.subshell_l( Subshell(9, kapa) );    lb = Basics.subshell_l( Subshell(9, kapb) )
        #
        # Test for parity
        if  L.den != 1   error("stop a")                end
        if  isodd( la + lb + L.num )    return( 0. )    end
        #
        wa = AngularMomentum.phaseFactor([jb, +1, L, -1, ja]) * sqrt( (twoJ(jb)+1)*(twoJ(JL)+1) / (4*pi*(twoJ(ja)+1) ) ) *
             AngularMomentum.ClebschGordan(jb, AngularM64(1//2), L, AngularM64(0), ja, AngularM64(1//2)) *
             AngularMomentum.ClebschGordan(jb, mb, L, M, ja, ma)
        return( wa )
    end

    """
    + (kapa::Int64, kapb::Int64, L::AngularJ64)`
        ... evaluates the same angular I (kappa m, kappa' m', LM) integral but without the magnetic (orientational) dependence.
            A value::Float64 is returned.
    """
    function ChengI(kapa::Int64, kapb::Int64, L::AngularJ64)
        ja = Basics.subshell_j( Subshell(9, kapa) );    jb = Basics.subshell_j( Subshell(9, kapb) )   # Use principal QN n=9 arbitrarely here
        la = Basics.subshell_l( Subshell(9, kapa) );    lb = Basics.subshell_l( Subshell(9, kapb) )
        #
        # Test for parity
        if  L.den != 1   error("stop a")                end
        if  isodd( la + lb + L.num )    return( 0. )    end
        #
        # There occurs here two changes with regard to Cheng formulas:
        # i)  We now divide by sqrt(twoJ(jb)+1) instead of sqrt(twoJ(ja)+1) ... which is likely related to change emission - absorption
        # ii) The phase (-1)^(ja + L - jb) is replaced by (-1)^L  to get a proper phase between 1/2 --> 3/2 ME (compared to 3/2 --> 3/2) ...
        #     but which already comes from the Clebsch-Gordan
        ## wa = AngularMomentum.phaseFactor([ja, +1, L, -1, jb]) *
        wa = (-1)^L.num  * sqrt( (twoJ(jb)+1)*(twoJ(L)+1) / (4*pi) ) / sqrt(twoJ(jb)+1) *
             AngularMomentum.ClebschGordan(jb, AngularM64(1//2), L, AngularM64(0), ja, AngularM64(1//2))
        return( wa )
    end

    """
    `AngularMomentum.ClebschGordan(ja, ma, jb, mb, Jab, Mab)`
        ... calculates the Clebsch-Gordan coefficient  <ja, ma, jb, mb; Jab, Mab> for given quantum numbers by
            a proper call to a Wigner 3-j symbol. A value::Float64 is returned.
    """
    function ClebschGordan(ja, ma, jb, mb, Jab, Mab)
        mab = - Basics.twice(Mab) / 2
        pp  = (Basics.twice(ja) - Basics.twice(jb) + Basics.twice(Jab))/2
        cg  = (-1)^pp * sqrt(Basics.twice(Jab) + 1) * sf_coupling_3j(Basics.twice(ja), Basics.twice(jb), Basics.twice(Jab),
                                                                     Basics.twice(ma), Basics.twice(mb), Basics.twice(mab))
        return( cg )
    end

    """
    `AngularMomentum.ClebschGordan_old(ja::AngularJ64, ma::AngularM64, jb::AngularJ64, mb::AngularM64, Jab::AngularJ64, Mab::AngularM64)`
        ... calculates the Clebsch-Gordan coefficient  <ja, ma, jb, mb; Jab, Mab> for given quantum numbers by
            a proper call to a Wigner 3-j symbol; a value::Float64 is returned.
    """
    function ClebschGordan_old(ja::AngularJ64, ma::AngularM64, jb::AngularJ64, mb::AngularM64, Jab::AngularJ64, Mab::AngularM64)
        cg = 0.   ## ja-jb+Mab must be integer
        mab = AngularM64(-Mab.num, Mab.den)
        cg = AngularMomentum.phaseFactor([ja, -1, jb, +1, Mab]) * sqrt( (twoJ(Jab)+1) ) *
             AngularMomentum.Wigner_3j(ja, jb, Jab, ma, mb, mab)
    end

    """
    `AngularMomentum.CL_reduced_me(suba::Subshell, L::Int64, subb::Subshell)`
        ... calculates the reduced matrix element of the C^L spherical tensor <suba || C^(L) || subb>; a value::Float64
            is returned.
    """
    function  CL_reduced_me(suba::Subshell, L::Int64, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        rem(ja2+1, 2) != 0    &&    error("stop a")

        redme = ((-1)^((ja2+1)/2)) * sqrt( (ja2+1)*(jb2+1) ) *
                Wigner_3j(AngularJ64(ja2//2), AngularJ64(L), AngularJ64(jb2//2),  AngularM64(1//2), AngularM64(0), AngularM64(-1//2) )

        return( redme )
    end

    """
    `AngularMomentum.CL_reduced_me_rb(suba::Subshell, L::Int64, subb::Subshell)`
        ... calculates the reduced matrix element of the C^L spherical tensor <suba || C^(L) || subb>; a value::Float64 is returned.
    """
    function  CL_reduced_me_rb(suba::Subshell, L::Int64, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        if rem(la + lb + L, 2) != 0   return 0.     end

        redme = ((-1)^((jb2-2L-1)/2)) * sqrt( (jb2+1) ) *
                Wigner_3j(AngularJ64(ja2//2), AngularJ64(jb2//2), AngularJ64(L), AngularM64(1//2), AngularM64(-1//2), AngularM64(0) )

        return( redme )
    end

    """
    `AngularMomentum.isAllowedMultipole(syma::LevelSymmetry, multipole::EmMultipole, symb::LevelSymmetry)`
        ... evaluates to true if the given multipole may connect the two level symmetries, and false otherwise.
    """
    function  isAllowedMultipole(syma::LevelSymmetry, multipole::EmMultipole, symb::LevelSymmetry)
        if  isTriangle(syma.J, AngularJ64(multipole.L), symb.J)  &&
            parityEmMultipolePi(syma.parity, multipole, symb.parity)    return( true )
        else                                                            return( false )
        end
    end

    """
    `AngularMomentum.isTriangle(ja::AngularJ64, jb::Int64, jc::AngularJ64)`
        ... evaluates to true if Delta(ja,jb,jc) = 1, ie. if the angular momenta ja, jb and jc can couple to each other,
            and false otherwise.
    """
    function isTriangle(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)
        if  ja.den == 1   ja2 = 2ja.num   else   ja2 = ja.num   end
        if  jb.den == 1   jb2 = 2jb.num   else   jb2 = jb.num   end
        if  jc.den == 1   jc2 = 2jc.num   else   jc2 = jc.num   end
        isodd(ja2 + jb2 + jc2)   &&    error("Angular momenta do no fullfill proper coupling rules; 2ja = $ja2, 2jb = $jb2, 2jc = $jc2")
        if  ja2 + jb2 >= jc2     &&    jc2 + ja2 >= jb2   &&   jb2 + jc2 >= ja2    return( true )
        else                                                                       return( false )
        end
    end

    """
    `AngularMomentum.kappa_j(kappa::Int64)`  ... calculates the j::AngularJ64 value of a given kappa.
    """
    function  kappa_j(kappa::Int64)
        j = AngularJ64( abs(kappa) - 1//2 )
        return( j )
    end

    """
    `AngularMomentum.kappa_l(kappa::Int64)`  ... calculates the l::AngularJ64 value of a given kappa.
    """
    function  kappa_l(kappa::Int64)
        if  kappa < 0    l  = abs(kappa) -1   else   l  = kappa   end
        return( AngularJ64(l) )
    end

    """
    `AngularMomentum.oneJ(ja::AngularJ64)`  ... calculates ja; a (positive) value::Float64 is returned.
    """
    function  oneJ(ja::AngularJ64)
        if  ja.den  == 1    ja1 = 1.0 * ja.num   else   ja1 = ja.num / 2.   end
        return( ja1 )
    end

    """
    `AngularMomentum.parityEmMultipolePi(pa::Parity, multipole::EmMultipole, pb::Parity)`
        ... evaluates to true if the given multipole fullfills the parity selection rule pi(a, multipole, b) = 1,
            and false otherwise. This includes a proper test for both, electric and magnetic multipoles, based on
            multipole.electric.
    """
    function  parityEmMultipolePi(pa::Parity, multipole::EmMultipole, pb::Parity)
        # The the multipolarity into account to 'interprete' electric multipoles
        if       isodd(multipole.L)    Le = multipole.electric    else   Le = !(multipole.electric)   end

        if        pa == Basics.plus    &&   Le   &&   pb == Basics.minus        return( true )
        elseif    pa == Basics.minus   &&   Le   &&   pb == Basics.plus         return( true )
        elseif    pa == pb          &&   !(Le)                                          return( true )
        else                                                                            return( false )
        end
    end

    """
    `AngularMomentum.phaseFactor(list::Array{Any,1})`
        ... checks and calculates the phase factor (-1)^(ja + mb -jc ...) that occur frequently in angular momentum theory;
            a value +1. or -1. is returned. Use phaseFactor([ja::Union{AngularJ64,AngularM64), -1, mb::Union{AngularJ64,AngularM64),
            ..., 1, jc::Union{AngularJ64,AngularM64)]) to specify the phase.
    """
    function phaseFactor(list::Array{Any,1})
        if   iseven( length(list) )                                                   error("Wrong number of arguments.")   end
        for  i = 1:2:length(list)
             if  !(typeof(list[i]) == AngularJ64  || typeof(list[i]) == AngularM64)   error("Wrong type of argument $i")    end
             if  i == length(list)                                                    break                                 end
             if  list[i+1] != 1   &&  list[i+1] != -1                error("Wrong type of argument $(i+1); must be +-1")    end
        end
        #
        ja = list[1];        if  ja.den  == 1    jm2 = 2ja.num                  else   jm2 = ja.num                  end
        for  i = 2:2:length(list)
           ja = list[i+1];   if  ja.den  == 1    jm2 = jm2 + list[i]*2ja.num    else   jm2 = jm2 + list[i]*ja.num    end
        end
        #
        if rem(jm2,2) != 0    error("Improper combination of angular momenta")   end

        return( (-1.)^(jm2/2) )
    end

    """
    `AngularMomentum.sigma_reduced_me(suba::Subshell, subb::Subshell)`
        ... calculates the reduced matrix element of the sigma^(1) spherical tensor <suba || sigma^(1) || subb>;
            a value::Float64 is returned.
    """
    function  sigma_reduced_me(suba::Subshell, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        if rem(la + lb, 2) != 0   return 0.     end

        ## redme = ((-1)^((jb2-2L-1)/2)) * sqrt( (jb2+1) ) *
        ##         Wigner_3j(AngularJ64(ja2//2), AngularJ64(jb2//2), AngularJ64(L), AngularM64(1//2), AngularM64(-1//2), AngularM64(0) )
        redme = 0.

        return( redme )
    end

    """
    `AngularMomentum.triangularDelta(ia2::Int64, ib2::Int64, ic2::Int64)`
        ... calculates the tringular Delta(ja,jb,jc). The arguments in this integer function are i2a = 2*ja+1, ...
            The result is 0 if the triangular condition failes and 1 otherwise.
    """
    function triangularDelta(ia2::Int64, ib2::Int64, ic2::Int64)
        i = ib2 - ic2
        if  ia2 >= abs(i) + 1   &&   ia2 <= ib2 + ic2 - 1    return( 1 )   else    return( 0 )    end
    end

    """
    `AngularMomentum.triangularDelta(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)`
        ... calculates the tringular Delta(ja,jb,jc). The result is 0 if the triangular condition failes and 1 otherwise.
    """
    function triangularDelta(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)
        if  abs(ja.num//ja.den - jb.num//jb.den) <= jc.num//jc.den <= ja.num//ja.den + jb.num//jb.den  return( 1 )   else    return( 0 )    end
    end

    """
    `AngularMomentum.twoJ(ja::AngularJ64)`  ... calculates 2*ja; an (positive) value::Int64 is returned.
    """
    function  twoJ(ja::AngularJ64)
        if  ja.den  == 1    ja2 = 2ja.num   else   ja2 = ja.num   end
        return( ja2 )
    end

    """
    `AngularMomentum.twoM(ma::AngularM64)`  ... calculates 2*ma; a value::Int64 is returned.
    """
    function  twoM(ma::AngularM64)
        if  ma.den  == 1    ma2 = 2ma.num   else   ma2 = ma.num   end
        return( ma2 )
    end

    """
    `AngularMomentum.Wigner_3j(a, b, c, m_a, m_b, m_c)`
        ... calculates the value of a Wigner 3-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_3j(a, b, c, m_a, m_b, m_c)
        sf_coupling_3j(Basics.twice(a),   Basics.twice(b),   Basics.twice(c),
                       Basics.twice(m_a), Basics.twice(m_b), Basics.twice(m_c))
    end

    """
    `AngularMomentum.Wigner_3j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, m_a::AngularM64, m_b::AngularM64, m_c::AngularM64)`
         ... calculates the value of a Wigner 3-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_3j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, m_a::AngularM64, m_b::AngularM64, m_c::AngularM64)
        # Use twice the given arguments and integer arithmetik throughout.
        if  a.den == 1      ja = a.num + a.num    else    ja = a.num      end
        if  b.den == 1      jb = b.num + b.num    else    jb = b.num      end
        if  c.den == 1      jc = c.num + c.num    else    jc = c.num      end
        if  m_a.den == 1    ma = 2m_a.num         else    ma = m_a.num    end
        if  m_b.den == 1    mb = 2m_b.num         else    mb = m_b.num    end
        if  m_c.den == 1    mc = 2m_c.num         else    mc = m_c.num    end

        # Test the triangular condition and that for magnetic quantum numbers
        if      ma+mb+mc != 0  ||   triangularDelta(ja+1,jb+1,jc+1) == 0      return( 0. )
        elseif  abs(ma) > ja   ||   abs(mb) > jb   ||  abs(mc) > jc           return( 0. )
        elseif  rem(ma+ja+ja,2) != rem(ja,2)   ||   rem(mb+jb+jb,2) != rem(jb,2)   ||   rem(mc+jc+jc,2) != rem(jc,2)    error("stop a")
        end

        ik = [ ja + jb - jc,    ja - jb + jc,    -ja + jb + jc,    ja + jb + jc + 2,    ja - ma,         ja + ma,         jb - mb,
               jb + mb,         jc - mc,          jc + mc,         jb - jc - ma,        ja - jc + mb,    jc - jb + ma,    ja - jb - mc ]
        for  i  in  1:14
            rem(ik[i],2) == 1   &&   error("stop b")
            ik[i] = ik[i] / 2
        end
        # Calculate the 3-j delta factor
        delta = (factorial( ik[1] ) * factorial( ik[2] ) * factorial( ik[3]) * factorial( ik[5]) * factorial( ik[6] ) * factorial( ik[7] ) *
                 factorial( ik[8] ) * factorial( ik[9] ) * factorial(ik[10] )  /  factorial( ik[4] ) )
        # Find out the intervall of summation  k  and sum up
        kmin = max(0, ik[11], ik[12])
        kmax = min(ik[1], ik[5], ik[8] )
        qsum = 0.
        for  k in kmin:kmax
            sumk = (factorial(k) * factorial( ik[1]-k )  * factorial( ik[5]-k ) * factorial( ik[8]-k ) *
                    factorial( ik[13]+k ) * factorial( k-ik[12]) )
            sumk = 1.0 / sumk
            if  rem(k,2) == 0   qsum = qsum + sumk    else   qsum = qsum - sumk    end
        end

        if  rem(ik[14],2) != 0    w3j = - sqrt(delta) * qsum    else    w3j = sqrt(delta) * qsum    end
        return( w3j )
    end

    """
    `AngularMomentum.Wigner_6j(a, b, c, d, e, f)`
        ... calculates the value of a Wigner 6-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_6j(a, b, c, d, e, f)
        sf_coupling_6j(Basics.twice(a), Basics.twice(b), Basics.twice(c),
                       Basics.twice(d), Basics.twice(e), Basics.twice(f))
    end

    """
    `AngularMomentum.Wigner_6j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, e::AngularJ64, f::AngularJ64)`
         ... calculates the value of a Wigner 6-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_6j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, ee::AngularJ64, f::AngularJ64)
        # Use twice the given arguments and integer arithmetik throughout.
        if  a.den == 1       ja = a.num + a.num    else    ja = a.num      end
        if  b.den == 1       jb = b.num + b.num    else    jb = b.num      end
        if  c.den == 1       jc = c.num + c.num    else    jc = c.num      end
        if  d.den == 1       jd = d.num + d.num    else    jd = d.num      end
        if  ee.den == 1      je = ee.num + ee.num  else    je = ee.num     end
        if  f.den == 1       jf = f.num + f.num    else    jf = f.num      end
        #
        wa = AngularMomentum.triangularDelta(ja+1, jb+1, jc+1) * AngularMomentum.triangularDelta(ja+1, je+1, jf+1) *
             AngularMomentum.triangularDelta(jd+1, jb+1, jf+1) * AngularMomentum.triangularDelta(jd+1, je+1, jc+1)
        if  wa == 0    return( 0. )    end
        #
        # Evaluate the algebraic formulae
        n1 = Int64( (ja + jb + jc)//2 );    n2 = Int64( (je + jd + jc)//2 );       n3 = Int64( (ja + je + jf)//2 )
        n4 = Int64( (jb + jd + jf)//2 );    n5 = Int64( (ja + jb + je + jd)//2 );  n6 = Int64( (ja + jd + jc + jf)//2 )
        n7 = Int64( (jb + je + jc + jf)//2 )
        kmin = max(n1,n2,n3,n4) + 1;   kmax = min(n5,n6,n7) + 1
        #
        w6j = 1.0;   icount = 0
        if  kmin != kmax
            for k = kmin+1:kmax
                ki = kmax - icount
	            w6j = 1.0 - (w6j * ki*(n5 - ki + 2.) * (n6 - ki + 2.) * (n7 - ki + 2.)) /
                            ((ki - 1. - n1) * (ki - 1. - n2) * (ki - 1. - n3) * (ki - 1. - n4))
                icount = icount + 1
            end
        end
        #
        w6j = w6j * exp( (lgamma(kmin+1) - lgamma(kmin-n1) - lgamma(kmin-n2) - lgamma(kmin-n3) - lgamma(kmin-n4) - lgamma(n5+2-kmin) -
	                      lgamma(n6+2-kmin) - lgamma(n7+2-kmin))  +
	                   ( (lgamma(n1+1-ja) + lgamma(n1+1-jb) + lgamma(n1+1-jc) - lgamma(n1+2) + lgamma(n2+1-je) + lgamma(n2+1-jd) +
	                      lgamma(n2+1-jc) - lgamma(n2+2) + lgamma(n3+1-ja) + lgamma(n3+1-je) + lgamma(n3+1-jf) - lgamma(n3+2) +
                          lgamma(n4+1-jb) + lgamma(n4+1-jd) + lgamma(n4+1-jf) - lgamma(n4+2)) / 2. ) )
        #
        if   rem(n5+kmin,2) == 0       w6j = -w6j   end
        if   rem(ja+jb+je+jd,4) != 0   w6j = -w6j   end

        return( w6j )
    end

    """
    `AngularMomentum.Wigner_9j(a, b, c, d, e, f, g, h, i)`
        ... calculates the value of a Wigner 3-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_9j(a, b, c, d, e, f, g, h, i)
        sf_coupling_9j(Basics.twice(a), Basics.twice(b), Basics.twice(c),
                       Basics.twice(d), Basics.twice(e), Basics.twice(f),
                       Basics.twice(g), Basics.twice(h), Basics.twice(i))
    end

    """
    `AngularMomentum.Wigner_9j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, e::AngularJ64, f::AngularJ64,
                                   g::AngularJ64, h::AngularJ64, i::AngularJ64)`
         ... calculates the value of a Wigner 6-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_9j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, ee::AngularJ64, f::AngularJ64,
                       g::AngularJ64, h::AngularJ64, i::AngularJ64)
        Defaults.warn(AddWarning, "AngularMomentum.Wigner_9j(): Here, we wish to call an external Julia package/interface; now set to 2.0")
        w9j = 2.0
        return( w9j )
    end

    #========================  Currently not in use

    """
        clebschgordan(j₁, m₁, j₂, m₂, J, M)

    Calculate the Clebsch–Gordan coefficient ``\\langle j_1 m_1 j_2 m_2 | J M \\rangle``. The
    calculation uses the Wigner 3-j symbol from the GNU Scientific Library (cf.
    [`wigner3j`](@ref)). The Condon–Shortley phase convention is employed.
    """
    function clebschgordan(j₁, m₁, j₂, m₂, J, M)
        J₁ = HalfInt(j₁) # explicit conversion because AngularJ64 does not support arithmetic
        J₂ = HalfInt(j₂) # explicit conversion because AngularJ64 does not support arithmetic
        mM = -HalfInt(M) # explicit conversion because AngularM64 does not support arithmetic
        exp = mM-J₁+J₂
        isinteger(exp) || return 0.0
        phase = ifelse(iseven(Integer(exp)), 1, -1)
        phase * sqrt(Basics.twice(J)+1) * wigner3j(J₁, J₂, J, m₁, m₂, mM)
    end

    """
        wigner3j(j₁, j₂, j₃, m₁, m₂, m₃)

    Calculate the Wigner 3-j symbol
    ```math
    \\begin{pmatrix}
        j_1 & j_2 & j_3 \\\\
        m_1 & m_2 & m_3
    \\end{pmatrix}
    ```
    by calling the corresponding function from the GNU Scientific Library.
    """
    function wigner3j(j₁, j₂, j₃, m₁, m₂, m₃)
        sf_coupling_3j(Basics.twice(j₁), Basics.twice(j₂), Basics.twice(j₃),
                       Basics.twice(m₁), Basics.twice(m₂), Basics.twice(m₃))
    end

    """
        wigner6j(j₁, j₂, j₃, j₄, j₅, j₆)

    Calculate the Wigner 6-j symbol
    ```math
    \\begin{Bmatrix}
        j_1 & j_2 & j_3 \\\\
        j_4 & j_5 & j_6
    \\end{Bmatrix}
    ```
    by calling the corresponding function from the GNU Scientific Library.
    """
    function wigner6j(j₁, j₂, j₃, j₄, j₅, j₆)
        sf_coupling_6j(Basics.twice(j₁), Basics.twice(j₂), Basics.twice(j₃),
                       Basics.twice(j₄), Basics.twice(j₅), Basics.twice(j₆))
    end

    """
        wigner9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)

    Calculate the Wigner 9-j symbol
    ```math
    \\begin{Bmatrix}
        j_1 & j_2 & j_3 \\\\
        j_4 & j_5 & j_6 \\\\
        j_7 & j_8 & j_9
    \\end{Bmatrix}
    ```
    by calling the corresponding function from the GNU Scientific Library.
    """
    function wigner9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
        sf_coupling_9j(Basics.twice(j₁), Basics.twice(j₂), Basics.twice(j₃),
                       Basics.twice(j₄), Basics.twice(j₅), Basics.twice(j₆),
                       Basics.twice(j₇), Basics.twice(j₈), Basics.twice(j₉))
    end
    ===========================#

end # module


text:
const SReal{V} = Union{SInteger{V},SFloat{V},SRational{V}}

SReal(x::BaseInteger) = SInteger(x)
SReal(x::BaseFloat) = SFloat(x)
SReal(x::Rational) = SRational(x)

Base.promote_eltype(x::SReal, y::BaseNumber) = promote_type(eltype(x), eltype(y))
Base.promote_eltype(x::BaseNumber, y::SReal) = promote_type(eltype(x), eltype(y))
Base.promote_eltype(x::SReal, y::SReal) = promote_type(eltype(x), eltype(y))

function SReal(::Val{V}) where V
    if V isa Integer
        if V isa Unsigned
            if V isa UInt8
                SUInt8{V}()
            elseif V isa UInt16
                SUInt16{V}()
            elseif V isa UInt32
                SUInt32{V}()
            elseif V isa UInt128
                SUInt128{V}()
            else
                SUInt64{V}()
            end
        elseif V isa Bool
            SBool{V}()
        elseif V isa Signed
            if V isa Int8
                SInt8{V}()
            elseif V isa Int16
                SInt16{V}()
            elseif V isa Int32
                SInt32{V}()
            elseif V isa Int128
                SInt128{V}()
            else
                SInt64{V}()
            end
        end
    elseif V isa AbstractFloat
        if V isa Float16
            SFloat16{V}()
        elseif V isa Float32
            SFloat32{V}()
        else
            SFloat64{V}()
        end
    elseif V isa Rational
        SRational(val)
    end
end

S2B = Dict(SUInt128 => UInt128,
           SUInt16 => UInt16,
           SUInt32 => UInt32,
           SUInt64 => UInt64,
           SUInt8 => UInt8,

           SInt128 => Int128,
           SInt16 => Int16,
           SInt32 => Int32,
           SInt64 => Int64,
           SInt8 => Int8,

           SFloat64 => Float64,
           SFloat32 => Float32,
           SFloat16 => Float16)

B2S = Dict(UInt128 => SUInt128,
           UInt16 => SUInt16,
           UInt32 => SUInt32,
           UInt64 => SUInt64,
           UInt8 => SUInt8,

           Int128 => SInt128,
           Int16 => SInt16,
           Int32 => SInt32,
           Int64 => SInt64,
           Int8 => SInt8,

           Float64 => SFloat64,
           Float32 => SFloat32,
           Float16 => SFloat16)


text:
module TestGraphviz
using Test

using Catlab.Graphics.Graphviz

# Pretty-print
##############

spprint(expr::Expression) = sprint(pprint, expr)

# Node statement
@test spprint(Node("n")) == "n;"
@test spprint(Node("n"; label="foo")) == "n [label=\"foo\"];"
@test spprint(Node("n"; shape="box", style="filled")) ==
  "n [shape=\"box\",style=\"filled\"];"

# Edge statement
@test spprint(Edge("n1","n2")) == "n1 -- n2;"
@test spprint(Edge("n1","n2"; label="foo")) == "n1 -- n2 [label=\"foo\"];"
@test spprint(Edge("n1","n2"; style="dotted", weight="10")) ==
  "n1 -- n2 [style=\"dotted\",weight=\"10\"];"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"))) == "n1:p1 -- n2:p2;"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"); label="bar")) ==
  "n1:p1 -- n2:p2 [label=\"bar\"];"
@test spprint(Edge(NodeID("n1","p1","w"), NodeID("n2", "p2", "e"))) ==
  "n1:p1:w -- n2:p2:e;"
@test spprint(Edge("n1","n2","n3")) == "n1 -- n2 -- n3;"
@test spprint(Edge(NodeID("n1"), NodeID("n2"), NodeID("n3"))) ==
  "n1 -- n2 -- n3;"

# Graph statement
graph = Graph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
graph G {
  n1;
  n2;
  n1 -- n2;
}
"""

graph = Digraph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
digraph G {
  n1;
  n2;
  n1 -> n2;
}
"""

graph = Digraph("G",
  Node("n1"; label="foo"),
  Node("n2"; label="bar"),
  Edge("n1","n2");
  graph_attrs = Attributes(:rankdir => "LR"),
  node_attrs = Attributes(:shape => "box", :style => "filled"),
  edge_attrs = Attributes(:style => "dotted")
)
@test spprint(graph) == """
digraph G {
  graph [rankdir="LR"];
  node [shape="box",style="filled"];
  edge [style="dotted"];
  n1 [label="foo"];
  n2 [label="bar"];
  n1 -> n2;
}
"""

# Subgraph statement
subgraph = Subgraph("sub",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(subgraph) == """
subgraph sub {
  n1;
  n2;
  n1 -- n2;
}"""

subgraph = Subgraph(
  Node("n1"),
  Node("n2"),
  graph_attrs = Attributes(:rank => "same")
)
@test spprint(subgraph) == """
{
  graph [rank="same"];
  n1;
  n2;
}"""

end


text:
module TailRec
export @tailrec

macro tailrec(func)
    fargs=map(e->if isa(e,Expr) e.args[1] else e end,func.args[1].args)
    fbody=func.args[2]
    fbody=rewrite(fbody,fargs)
    func.args[2]=Expr(:block,:(@label retry),fbody)
    esc(func)
end

function rewrite(expr,args,callflag=false)
    if !isa(expr,Expr)
        expr
    elseif expr.head == :call && expr.args[1] == args[1]
        if callflag
            @warn "Not tail recursive call is found."
            expr
        else
            newargs=Expr(:tuple)
            newargs.args=args[2:end]
            oldargs=Expr(:tuple)
            oldargs.args=expr.args[2:end]
            Expr(:block, Expr(:(=),newargs,oldargs), :(@goto retry) )
        end
    elseif expr.head == :block
        expr.args[end]=rewrite(expr.args[end],args, expr.head==:call || callflag)
        expr
    else
        expr.args = map(a->rewrite(a,args, expr.head==:call ||callflag),expr.args)
        expr
    end
end
end


text:
abstract type AverageMode end

module AvgMode
    import ..MLMetrics.AverageMode
    export None, Micro, Macro

    struct None <: AverageMode end
    struct Micro <: AverageMode end
    struct Macro <: AverageMode end
end

Base.convert(::Type{AverageMode}, avgmode) = throw(ArgumentError("Unknown way to specify a obsdim: $avgmode"))
Base.convert(::Type{AverageMode}, avgmode::AverageMode) = avgmode
Base.convert(::Type{AverageMode}, avgmode::String) = convert(AverageMode, Symbol(lowercase(avgmode)))
Base.convert(::Type{AverageMode}, ::Nothing) = AvgMode.None()

function Base.convert(::Type{AverageMode}, avgmode::Symbol)
    if avgmode == :micro
        AvgMode.Micro()
    elseif avgmode == :macro
        AvgMode.Macro()
    elseif avgmode == Symbol("nothing") || avgmode == :none || avgmode == :null || avgmode == :na || avgmode == :undefined
        AvgMode.None()
    else
        throw(ArgumentError("Unknown way to specify a \"avgmode\": $avgmode"))
    end
end


text:
# flexible representations of priors

# Collection of priors from -- False Discovery Rates: A new deal

ash_spiky =  MixtureModel( [Normal(0,0.25); Normal(0,0.5); Normal(0,1); Normal(0,2)] ,
                           [0.4; 0.2; 0.2; 0.2 ])

ash_nearnormal = MixtureModel([Normal(0,1); Normal(0,2)], [1/3; 2/3])

ash_flattop = MixtureModel([Normal(x,0.5) for x=-1.5:0.5:1.5])

ash_skew = MixtureModel([Normal(-2,2), Normal(-1,1.5), Normal(0,1), Normal(1,1)],
                        [1/4; 1/4; 1/3; 1/6])

ash_bignormal = Normal(0,4)

ash_bimodal = MixtureModel( [Normal(-2,1); Normal(2,1)],
                            [0.5;0.5] )

# Distributions used for simulations in this paper

hard_dist = MixtureModel([ Normal(-0.2,.2), Normal(0,0.9)],[0.7, 0.3])
easy_dist = MixtureModel([Normal(-1.5,.2), Normal(1.5, .2)])


text:
# location, measurements

# elipsis in print_range

import Base: Pair,
    pairs, getindex, length, getproperty, map, map!,
    get, keys, values, keytype, eltype, ==,
    vcat, push!, copy

include("unroll1.jl")

"""
    AbstractPairedArray

Simple data structure pairing two arrays, with the first array holding the `keys`
and the second the corresponding `values`.

For example an array of locations with an array of
corresponding measurements or a vector of times with
a vector of locations.
"""
abstract type  AbstractPairedArray
end

const APA = AbstractPairedArray

#=
struct PairedArray{T,S} <: AbstractPairedArray
    x::S
    x::T
end
=#

"""
    Trajectory{S,T} <: AbstractPairedArray

Pairs a linearly ordered (abstract) vector of `keys` (typically times) with
a vector of `values` (typically locations).
"""
struct Trajectory{S,T} <: AbstractPairedArray
    t::S
    x::T
end

keys(X::Trajectory) = X.t
values(X::Trajectory) = X.x

keytype(X::Trajectory) = eltype(X.t)
eltype(X::Trajectory) = eltype(X.x)

vcat(Xs::Trajectory...) = Trajectory(vcat(keys.(Xs)...), vcat(values.(Xs)...))

function push!(X::Trajectory, (t,x)::Pair)
    push!(keys(X), t)
    push!(values(X), x)
    X
end
copy(X::Trajectory) = map(copy, X)

function _find(r::AbstractRange, x)
    n = round(Integer, (x - first(r)) / step(r)) + 1
    if n >= 1 && n <= length(r) && r[n] == x
        return n
    else
        error("index error")
    end
end

function get(X::Trajectory{<:AbstractVector}, key)
    i = searchsorted(keys(X), key)
    isempty(i) && error("key not found")
    first(i) != last(i) && error("key not unique")
    return values(X)[first(i)]
end

trajectory(t, x) = Trajectory(t, x)
function trajectory(itr)
    local t, x
    @unroll1 for (tᵢ, xᵢ) in itr
        if $first
            t = [tᵢ]
            x = [xᵢ]
        else
            push!(t, tᵢ)
            push!(x, xᵢ)
        end
    end
    Trajectory(t, x)
end

# Iteration and indexing spared out, see issue #1
length(X::Trajectory) = length(values(X)) # seems useful enough
#=

iterate(X::Trajectory) = X.t, (X.x, nothing)
iterate(X::Trajectory, state) = state
length(X) = 2

eachindex(X::Trajectory) = eachindex(X.x)
getindex(X::Trajectory{<:AbstractVector}, i) = getindex(X.x, i)
=#

Pair(X::Trajectory) = Pair(X.t, X.x)

pairs(X::Trajectory) = (t => x for (t, x) in zip(X.t, X.x))

# sic! keys compare with `isequal`
X::Trajectory == Y::Trajectory = isequal(X.t, Y.t) && X.x == Y.x

"""
    issynchron(X, Y) = isequal(keys(X), keys(Y))
"""
issynchron(X, Y) = isequal(keys(X), keys(Y))
checksynchron(X, Y) = issynchron(X, Y) || error("checksynchron: asynchronous trajectories")

function mapeach!((f,g)::Pair, X::Trajectory)
    map!(f, X.t)
    map!(g, X.x)
    X
end
function mapeach((f,g)::Pair, X::Trajectory)
    Trajectory(map(f,X.t), map(g, X.x))
end

function map(f, X::Trajectory)
    Trajectory(f(X.t), f(X.x))
end
function map((f,g)::Pair, X::Trajectory)
    Trajectory(f(X.t), g(X.x))
end

function mapvalues(f, X::Trajectory, Y::Trajectory)
    checksynchron(X, Y)
    Trajectory(X.t, map(f, X.x, Y.x))
end

function mapvalues!(f, X::Trajectory, Y::Trajectory)
    checksynchron(X, Y)
    map!(f, X.x, Y.x)
    X
end