text:
# Random forests regression

"""
Finds the threshold to split `features` with that minimizes the
mean-squared-error loss over `labels`.
Returns (best_val, best_thresh), where `best_val` is -MSE

Variable name changes from original code
  labels   -> y
  features -> x
  nl       -> n_left
  nr       -> n_right
  su       -> sum_y
  su2      -> sum_y2
  s_l      -> sum_y_left
  s2_l     -> sum_y2_left
  s_r      -> sum_y_right
  s2_r     -> sum_y2_right

"""
function _best_mse_loss{T<:Float64, U<:Real}(y::Vector{T}, x::DataArray{U, 1}, domain)
    best_val = -Inf
    best_thresh = 0.0
    n = length(y)

    sum_y_left = sum_y2_left = zero(T)            # scalar values of 0

    sum_y = sum(y)                             # scalar sum of all y_i
    sum_y2 = zero(T);                             # scalar value of 0

    # get sum of squares
    for i = 1:n
        sum_y2 += y[i]^2
    end

    n_left = 0
    i = 1

    # Since`x` is sorted, below is an O(n) algorithm for finding the optimal
    # threshold in `domain`. We iterate through the array and update sum_y_left
    # and sum_y_right (= sum(y) - sum_y_left) as we go. - @cstjean
    for thresh in domain

        # this loop checks which side of the split this x_i is on
        while i <= n && x[i] < thresh
            sum_y_left += y[i]
            sum_y2_left += y[i]^2
            n_left += 1
            i += 1
        end
        sum_y_right = sum_y - sum_y_left
        sum_y2_right = sum_y2 - sum_y2_left
        n_right = n - n_left

        # This check is necessary I think because in theory all y could
        # be the same, then either n_left or n_right would be 0. - @cstjean
        if n_right > 0 && n_left > 0

            # This isn't really squared-error loss. We're computing the sum of variance
            # in the daughter nodes times the N for that daughter node. This last point
            # is why the two paranthetical terms here look a bit diffent than mere variance.
            # We have canceled a constant term N for each term.
            loss = (sum_y2_left - sum_y_left^2/n_left) + (sum_y2_right - sum_y_right^2/n_right)

            # update best value and threshold
            if -loss > best_val
                best_val = -loss
                best_thresh = thresh
            end
        end
    end
    return best_val, best_thresh
end

"""
This function returns a tuple with the column index and the threshold
value of the predictor in the matrix X that minimizes MSE. The function
assumes complete data in both y and X.

Variable name changes from original code
  labels     -> y
  features   -> x
  nr         -> n
  nf         -> p
  i          -> j
  features_i -> x_j
  labels_i   -> y_ord
  domain_i   -> domain_j
  inds       -> col_indcs

"""

function _split_mse_df{T<:Float64}(y::Vector{T}, X::DataFrame, mtry::Int)
    n, p = size(X)
    best = (0, 0.0)
    best_val = -Inf
    if mtry > 0
        r = randperm(p)
        col_indcs = r[1:mtry]
    else
        col_indcs = collect(1:p)
    end
    for j in col_indcs
        keep_row = !isna(X[:, j])

        # x_obs = convert(Vector, X[keep_row, j])

        if !any(keep_row)
            continue
        end

        # y_obs = y[keep_row]
        # ord::Array{Int64, 1} = sortperm(x_obs)
        # x_j = convert(Array{typeof(x_obs[1])}, x_obs[ord])
        # y_ord = y_obs[ord]

        x_j, y_ord = prepare_xy_obs(X[keep_row, j], y[keep_row])

        if n > 100
            if VERSION >= v"0.4.0-dev"
                domain_j = quantile(x_j, linspace(0.01, 0.99, 99); sorted=true)
            else  # sorted=true isn't supported on StatsBase's Julia 0.3 version
                domain_j = quantile(x_j, linspace(0.01, 0.99, 99))
            end
        else
            domain_j = x_j
        end
        value, thresh = _best_mse_loss(y_ord, x_j, domain_j)
        if value > best_val
            best_val = value
            best = (j, thresh)
        end
    end
    return best
end

function prepare_xy_obs{T}(x_obs::DataArray{T, 1}, y_obs)
    ord = sortperm(x_obs)
    x_j::DataArray{T, 1} = x_obs[ord]
    y_ord = y_obs[ord]
    return (x_j, y_ord)
end

function find_na_cols(dat)
    cols_with_na = Array{Int, 1}(0)

    for j in 1:ncol(dat)
        if sum(isna(dat[:, j])) ≠ 0
            push!(cols_with_na, j)
        end
    end
    return cols_with_na
end

# function find_surrogates{T::BitArray}(left_node::T, y::Vector, X::Matrix, mtry::Int)

function build_stump{T <: Float64, U<:Real}(y::Vector{T}, X::Matrix{U})
    S = _split_mse(y, X, 0)

    if S == NO_BEST
        return Leaf(mean(y), y)
    end

    col_idx, thresh = S
    split = X[:, col_idx] .< thresh

    return Node(col_idx,
                thresh,
                Leaf(mean(y[split]), y[split]),
                Leaf(mean(y[!split]), y[!split]))
end

function build_tree_df{T <: Float64}(y::Vector{T}, X::DataFrame, maxlabels = 5, mtry = 0, maxdepth = -1, n_surrogates = 5)
    n, p = size(X)
    if mtry == 0
        mtry = p
    end
    if maxdepth < -1
        error("Unexpected value for maxdepth: $(maxdepth) (expected: maxdepth >= 0, or maxdepth = -1 for infinite depth)")
    end
    if length(y) <= maxlabels || maxdepth == 0          # stopping rules
        return Leaf(mean(y), y)
    end

    S = _split_mse_df(y, X, mtry)                       # get [complete data] optimal split point

    if S == NO_BEST
        return Leaf(mean(y), y)
    end

    col_idx, thresh = S

    if n_surrogates ≠ 0
        cols_with_na = find_na_cols(X)                  # needed for all recursive steps since rows of X change
        surrogate_vars = Vector{Tuple{Int, Real}}(n_surrogates)
    end

    # if col_idx in cols_with_na
    if col_idx in 1:p                             # find surrogates even when no missing data
        na_rows = isna(X[:, col_idx])
        split_with_na = Array{Any, 1}(n)                # vector of Bools with some NA values

        for i = 1:n
            split_with_na[i] = isna(X[i, col_idx]) ? NA : X[i, col_idx] < thresh
        end
        row_indcs = collect(1:n)[!na_rows]
        col_indcs = deleteat!(collect(1:p), col_idx)

        # Here we need a function that splits so as to optimize agreement
        # with the `split_with_na` result for each observed `row_indcs`.
        surrogate_vars = surrogate_splits(split_with_na, X, row_indcs, col_indcs, 5)
        split = apply_surrogates(split_with_na, X, surrogate_vars)
    else
        split = X[:, col_idx] .< thresh
    end

    return Node(col_idx,
                thresh,
                surrogate_vars,
                build_tree_df(y[split], X[split,:], maxlabels, mtry, max(maxdepth-1, -1)),
                build_tree_df(y[!split], X[!split,:], maxlabels, mtry, max(maxdepth-1, -1)))
end

function build_forest_df{T <: Real}(y::Vector{T}, X::DataFrame, mtry, ntrees, maxlabels = 5, maxdepth = -1; nthreads = 1, n_surrogates = 5, oob_measure = "rsq")

    n = length(y)
    tree_arr = Array{Node, 1}(ntrees)

    if nthreads ≥ 2
        yhat_mat = fill(-Inf, (n, ntrees))              # use to compute oob-score
        mse_inc_dicts = Array{Dict{Int, Float64}, 1}(ntrees)

        @threads for t in 1:ntrees
            inds = sample(1:n, n)
            tree_arr[t] = build_tree_df(y[inds], X[inds, :], maxlabels, mtry, maxdepth, n_surrogates)

            # get OOB indices to calculate OOB score
            oob_indcs = setdiff(collect(1:n), inds)
            yhat, mse_inc = apply_tree(tree_arr[t], X[oob_indcs, :], y[oob_indcs])

            yhat_mat[oob_indcs, t] = yhat
            mse_inc_dicts[t] = mse_inc              # dict with MSE increases for OOB permutation importance
        end

        oob_score = par_oob_score(yhat_mat, y, oob_measure)
        mse_inc_dict = par_mse_inc(mse_inc_dicts)
    else
        yhat_mat = zeros(n, 2)                      # col 1 is numerator, col 2 is denominator
        yhat_mat[:, 1] = fill(-Inf, n)
        mse_inc_dict = Dict{Int, Float64}()

        for t in 1:ntrees
            inds = sample(1:n, n)
            tree_arr[t] = build_tree_df(y[inds], X[inds, :], maxlabels, mtry, maxdepth, n_surrogates)

            # get OOB indices to calculate OOB score
            oob_indcs = setdiff(collect(1:n), inds)
            yhat, mse_inc = apply_tree(tree_arr[t], X[oob_indcs, :], y[oob_indcs])

            if !all(isfinite(yhat_mat[oob_indcs, 1]))
                for idx in oob_indcs
                    if !isfinite(yhat_mat[idx, 1])
                        yhat_mat[oob_indcs, 1] = 0.0
                    end
                end
            end

            yhat_mat[oob_indcs, 1] += yhat
            yhat_mat[oob_indcs, 2] += 1
            update_mse_inc!(mse_inc_dict, mse_inc)
        end

        if all(isfinite(yhat_mat[:, 1]))

            if oob_measure == "rsq"
                oob_score = rsq(y, yhat_mat[:, 1] ./ yhat_mat[:, 2])
            elseif oob_measure == "mse"
                oob_score = mean_squared_error(y, yhat_mat[:, 1] ./ yhat_mat[:, 2])
            end
        else
            warn("Not enough trees built to compute OOB score. Some rows were never out-of-bag")
            return nothing
        end
    end
    return (Ensemble(tree_arr), oob_score, mse_inc_dict)
end

function build_forest_df_novarimp{T <: Real}(y::Vector{T}, X::DataFrame, mtry, ntrees, maxlabels = 5, maxdepth = -1; nthreads = 1, n_surrogates = 5, oob_measure = "rsq")

    n = length(y)
    tree_arr = Array{Node, 1}(ntrees)

    if nthreads ≥ 2
        yhat_mat = fill(-Inf, (n, ntrees))       # use to compute oob-score

        @threads for t in 1:ntrees
            inds = sample(1:n, n)
            tree_arr[t] = build_tree_df(y[inds], X[inds, :], maxlabels, mtry, maxdepth, n_surrogates)

            # get OOB indices to calculate OOB score
            oob_indcs = setdiff(collect(1:n), inds)
            yhat = apply_tree(tree_arr[t], X[oob_indcs, :])

            yhat_mat[oob_indcs, t] = yhat
        end
        oob_score = par_oob_score(yhat_mat, y, oob_measure)
    else
        yhat_mat = zeros(n, 2)                  # col 1 is numerator, col 2 is denominator
        yhat_mat[:, 1] = fill(-Inf, n)

        for t in 1:ntrees
            inds = sample(1:n, n)
            tree_arr[t] = build_tree_df(y[inds], X[inds, :], maxlabels, mtry, maxdepth, n_surrogates)

            # get OOB indices to calculate OOB score
            oob_indcs = setdiff(collect(1:n), inds)
            yhat = apply_tree(tree_arr[t], X[oob_indcs, :])

            if !all(isfinite(yhat_mat[oob_indcs, 1]))
                for idx in oob_indcs
                    if !isfinite(yhat_mat[idx, 1])
                        yhat_mat[oob_indcs, 1] = 0.0
                    end
                end
            end

            yhat_mat[oob_indcs, 1] += yhat
            yhat_mat[oob_indcs, 2] += 1
        end

        if all(isfinite(yhat_mat[:, 1]))

            if oob_measure == "rsq"
                oob_score = rsq(y, yhat_mat[:, 1] ./ yhat_mat[:, 2])
            elseif oob_measure == "mse"
                oob_score = mean_squared_error(y, yhat_mat[:, 1] ./ yhat_mat[:, 2])
            end
        else
            warn("Not enough trees built to compute OOB score. Some rows were never out-of-bag")
            return nothing
        end
    end
    return (Ensemble(tree_arr), oob_score)
end

# This function is a helper used for getting oob-score for
# parallel version of random forest regression. The parallel case
# needs to be treated a bit differently to avoid threads overwriting
# each other. NOTE: There might be a more elegant solution.
function par_oob_score(yhat_mat, y, measure = "rsq")
    n, ntrees = size(yhat_mat)
    yhat = fill(-Inf, n)
    num_preds = zeros(Int, n)       # this will be our denomiator

    for t = 1:ntrees
        for i = 1:n
            if isfinite(yhat_mat[i, t])
                if !isfinite(yhat[i])
                    yhat[i] = 0.0
                end
                yhat[i] += yhat_mat[i, t]
                num_preds[i] += 1
            end
        end
    end
    if all(isfinite(yhat))
        yhat ./= num_preds
    else
        warn("Not enough trees built to compute OOB score. Some rows never out-of-bag")
        return nothing
    end
    if measure == "rsq"
        oob_score = rsq(y, yhat)
    elseif measure == "mse"
        oob_score = mean_squared_error(y, yhat)
    end
    return oob_score
end

function update_mse_inc!(aggregate_dict, tree_dict)
    for k in keys(tree_dict)
        aggregate_dict[k] = get(aggregate_dict, k, 0.0) + tree_dict[k]
    end
end

function par_mse_inc(mse_inc_dicts)
    out = Dict{Int, Float64}()
    ntrees = length(mse_inc_dicts)
    for i = 1:ntrees
        update_mse_inc!(out, mse_inc_dicts[i])
    end
    return out
end


text:
#=
Dimensionality Reduction methods. Each method needs
1) Type with fields necessary for algorithm
2) function "reduction" to operate on sort with type field defined above
3) any other necessary functions for alignment algorithm
=#

export ReductionNone

function reductionprepare(r::Reduction,sort::Sorting)
    nothing
end

#=
No Dimensionality Reduction
=#
mutable struct ReductionNone <: Reduction
end

#=
Maximum Difference

Gibson et al 2010
=#
#=
type ReductionMD <: Reduction
    mydims::Int64
    maximum_difference::Array{Int64,1}
    local_difference::Array{Float64,1}
    spike_old::Array{Float64,1}
    Dc::Array{Int64,1}
end

ReductionMD()=ReductionMD(10,4)

ReductionMD(dims::Int64)=ReductionMD(10,dims)

ReductionMD(N::Int64,dims::Int64)=ReductionMD(dims,zeros(Int64,N),zeros(Float64,N),zeros(Float64,N),zeros(Int64,dims))

function reductionprepare(r::ReductionMD,sort::Sorting)

    max3ind=zeros(Int64,3)
    for i=1:length(sort.fullfeature)
        sort.r.local_difference[i]=abs(sort.r.spike_old[i]-sort.fullfeature[i])

        if sort.r.local_difference[i]>max3ind[1]
            if sort.r.local_difference[i]>max3ind[2]
                if sort.r.local_difference[i]>max3ind[3]
                    max3ind[1]=max3ind[2]
                    max3ind[2]=max3ind[3]
                    max3ind[3]=i
                end
                max3ind[1]=max3ind[2]
                max3ind[2]=i
            end
            max3ind[1]=i
        end
    end

    sort.r.maximum_difference[max3ind]+=1
    sort.r.spike_old[:]=sort.fullfeature[:]

    for i in max3ind
        for j=1:length(sort.r.Dc)
            if sort.r.maximum_difference[i]>sort.r.Dc[j]
                if i != sort.dims[j]
                    sort.dims[j]=i
                end
                sort.r.Dc[j]=sort.r.maximum_difference[i]
                break
            end
        end
    end
    nothing
end
=#
#=
Lilliefor's Test
=#
#=
type ReductionLT <: Reduction
end
=#
#=
Hartigan's Dip
=#
#=
type ReductionHD <: Reduction
end
=#
#=
Uniform Sampling
=#
#=
type ReductionUS <: Reduction
end
=#

ext: jl


text:
try
    Pkg.add("Coverage")
    using Coverage
    Coveralls.submit_token(Coveralls.process_folder())
    Codecov.submit_local(Codecov.process_folder())
catch e
    println(STDERR, e)
end

max_stars_repo_name: SimiDCI/GLVisualizeCI

max_issues_repo_name: SimiDCI/GLVisualizeCI

max_forks_repo_name: SimiDCI/GLVisualizeCI


text:
# Defines data types representing a purely geometrical mesh.

##########
# Points #
##########
immutable Point2<: FixedVector{Float64, 2}
    x::Float64
    y::Float64
end
typealias Vector2 Point2

immutable Point3 <: FixedVector{Float64, 3}
    x::Float64
    y::Float64
    z::Float64
end
typealias Vector3 Point3

#########
# Nodes #
#########
abstract AbstractGeoNode

immutable GeoNode2 <: AbstractGeoNode
    n::Int
    coords::Point2
end
# Check size here?
GeoNode2(n::Int, c::Vector{Float64}) = GeoNode2(n, Point2(c[1], c[2]))
GeoNode2(n::Int, c::Vector{Int}) = GeoNode2(n, convert(Vector{Float64}, c))
get_coord(node::GeoNode2) = Point3(node.coords.x, node.coords.y, 0.0)

immutable GeoNode3 <: AbstractGeoNode
    n::Int
    coords::Point3
end

# Check size here?
GeoNode3(n::Int, c::Vector{Float64}) = GeoNode3(n, Point3(c[1], c[2], c[3]))
GeoNode3(n::Int, c::Vector{Int}) = GeoNode3(n, convert(Vector{Float64}, c))

############
# Vertices #
############
abstract AbstractVertex{L} <: FixedVector{Int, L}

immutable Vertex2 <: AbstractVertex{2}
    v1::Int
    v2::Int
end

immutable Vertex3 <: AbstractVertex{3}
    v1::Int
    v2::Int
    v3::Int
end

immutable Vertex4 <: AbstractVertex{4}
    v1::Int
    v2::Int
    v3::Int
    v4::Int
end

immutable Vertex5 <: AbstractVertex{5}
    v1::Int
    v2::Int
    v3::Int
    v4::Int
end

immutable Vertex6 <: AbstractVertex{6}
    v1::Int
    v2::Int
    v3::Int
    v4::Int
    v5::Int
    v6::Int
end

###############
# GeoElements #
###############
abstract AbstractGeoElem
getindex(geoelem::AbstractGeoElem, i0::Real) = getindex(geoelem.vertices, i0)

abstract AbstractGeoElem2 <: AbstractGeoElem
abstract AbstractGeoElem3 <: AbstractGeoElem

immutable GeoTrig <: AbstractGeoElem2
    n::Int
    vertices::Vertex3
end
GeoTrig(n::Int, v::Vector{Int}) = GeoTrig(n, Vertex3(v[1], v[2], v[3]))

immutable GeoQTrig <: AbstractGeoElem2
    n::Int
    vertices::Vertex6
end
GeoQTrig(n::Int, v::Vector{Int}) = GeoQTrig(n, Vertex6(v[1], v[2], v[3], v[4], v[5], v[6]))

immutable GeoQuad <: AbstractGeoElem2
    n::Int
    vertices::Vertex4
end
GeoQuad(n::Int, v::Vector{Int}) = GeoQuad(n, Vertex4(v[1], v[2], v[3], v[4]))

immutable GeoTetra <: AbstractGeoElem3
    n::Int
    vertices::Vertex4
end
GeoTetra(n::Int, v::Vector{Int}) = GeoTetra(n, Vertex3(v[1], v[2], v[3], v[4]))

########
# Sets #
########
immutable ElementSet
    name::ASCIIString
    element_ids::Set{Int}
end
ElementSet(s::ASCIIString, elems) = ElementSet(s, Set{Int64}(elems))
Base.append!(ns1::ElementSet, elems) = union!(ns1.element_ids, Set{Int64}(element_ids))
Base.append!(es1::ElementSet, es2::ElementSet) = union!(es1.element_ids, ne2.element_ids)

immutable NodeSet
    name::ASCIIString
    node_ids::Set{Int}
end
NodeSet(s, nodes) = NodeSet(s, Set{Int64}(nodes))

Base.append!(ns1::NodeSet, nodes) = union!(ns1.node_ids, Set{Int64}(nodes))
Base.append!(ns1::NodeSet, ns2::NodeSet) = union!(ns1.node_ids, ns2.node_ids)

immutable EdgeSet
    name::String
    # (Element_id, edge_id) tuple
    edges::Vector{Tuple{Int, Int}}
end

immutable SurfaceSet
    name::String
    # (Element_id, surface_id) tuple
    surfaces::Vector{Tuple{Int, Int}}
end

###########
# GeoMesh #
###########
abstract AbstractGeoMesh

immutable GeoMesh <: AbstractGeoMesh
    nodes::Vector{GeoNode2}
    elements::Vector{AbstractGeoElem2}
    element_sets::Dict{ASCIIString, ElementSet}
    node_sets::Dict{ASCIIString, NodeSet}
end

function GeoMesh()
    nodes = Array(GeoNode2, 0)
    elements = Array(AbstractGeoElem2, 0)
    element_sets = Dict{ASCIIString, ElementSet}()
    node_sets = Dict{ASCIIString, NodeSet}()
    GeoMesh(nodes, elements, element_sets, node_sets)
end

# Push functions
push!(mesh::GeoMesh, node::AbstractGeoNode) = push!(mesh.nodes, node)

function push!(mesh::GeoMesh, nodes::Vector{GeoNode2})
    for n in nodes
        push!(mesh, n)
    end
end

push!(mesh::GeoMesh, elem::AbstractGeoElem) = push!(mesh.elements, elem)

function push!(mesh::GeoMesh, elem_set::ElementSet)
    mesh.element_sets[elem_set.name] = elem_set
end

function push!(mesh::GeoMesh, node_set::NodeSet)
     mesh.node_sets[node_set.name] = node_set
end

# Set generation functions
function gennodeset(f::Function, name::ASCIIString, nodes::Vector{GeoNode2})
    node_ids = Set{Int}()
    for node in nodes
        if f(node)
            push!(node_ids, node.n)
        end
    end
    return NodeSet(name, node_ids)
end


text:
using JLD
using Plots
using Colors
using LaTeXStrings

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

#Extracting stored varibles
ld = load("./out/jld/KLVarySigma.jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end
for (i, sm) in enumerate(["12", "32", "52"])

  newplot = plot(xlabel=L"$\log_{10}{\sigma}$",
                 ylabel=L"$\log_{10}(\mathrm{symmetrized} \ \mathrm{KL})$",
                 legend=:bottomleft,
                 ylim=(-6,4))

  plot!(newplot,
        log10.(σList),
        log10.(KLTrueNaive[i, :] + KLNaiveTrue[i, :]),
        label=L"\mathrm{Naive}",
        markershape=:diamond,
        markersize=15,
        color=lightblue,
        linestyle=:dash,
        linewidth=5)

  plot!(newplot,
        log10.(σList),
        log10.(KLTrueApp[i, :] + KLAppTrue[i, :]),
        label=L"\mathrm{IC}, \mathrm{nonzeros}(L)",
        markershape=:dtriangle,
        markersize=15,
        color=orange,
        linestyle=:dash,
        linewidth=5)

  plot!(newplot,
        log10.(σList),
        log10.(KLTrueLarge[i, :] + KLLargeTrue[i, :]),
        label=L"\mathrm{IC}, \mathrm{nonzeros}(L L^{\top})",
        markershape=:utriangle,
        markersize=15,
        color=rust,
        linestyle=:dash,
        linewidth=5)

  plot!(newplot,
        log10.(σList),
        log10.(KLTrueExact[i, :] + KLExactTrue[i, :]),
        label=L"\mathrm{Exact}",
        markershape=:circle,
        markersize=15,
        color=silver,
        linestyle=:dot,
        linewidth=5)

  savefig(newplot, "./out/plots/KLVarySigma$(sm).pdf")
  display(newplot)
end


text:
__precompile__()

module Immerse

using GtkUtilities, Colors, Reexport, Compat, REPL, Graphics
@reexport using Gadfly
import Gtk   # because both Gadfly and Gtk define draw
import Compose, Measures
using Cairo

@eval Compose begin import Cairo end
Compose.link_cairo()

export
    Figure,
    closefig,
    closeall,
    figure,
    gcf,
    scf,
    getproperty,
    setproperty!,
    hit,
    lasso_initialize

# Stuff for Gadfly/Compose
include("compose.jl")
# using .ImmerseCompose
include("display_gadfly.jl")
# using .DisplayGadfly

# Generic (?)
include("hit_test.jl")
# using .HitTest
include("select.jl")

end # module


text:
mutable struct Orthographic
    xmag::Cfloat
    ymag::Cfloat
    zfar::Cfloat
    znear::Cfloat
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Orthographic()
        obj = new()
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Orthographic, sym::Symbol, x)
    if sym === :zfar
        x > 0 || throw(ArgumentError("the distance to the far clipping plane should be > 0"))
    elseif sym === :znear
        x ≥ 0 || throw(ArgumentError("the distance to the near clipping plane should be ≥ 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Orthographic}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Orthographic}) = (:extensions, :extras)

mutable struct Perspective
    aspectRatio::Union{Nothing,Cfloat}
    yfov::Cfloat
    zfar::Union{Nothing,Cfloat}
    znear::Cfloat
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Perspective()
        obj = new()
        obj.aspectRatio = nothing
        obj.zfar = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Perspective, sym::Symbol, x)
    if sym === :aspectRatio && x !== nothing
        x > 0 || throw(ArgumentError("the aspect ratio of the field of view should be > 0"))
    elseif sym === :yfov
        x > 0 || throw(ArgumentError("the vertical field of view in radians should be > 0"))
    elseif sym === :zfar && x !== nothing
        x > 0 || throw(ArgumentError("the distance to the far clipping plane should be > 0"))
    elseif sym === :znear
        x > 0 || throw(ArgumentError("the distance to the near clipping plane should be > 0"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Perspective}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Perspective}) = (:aspectRatio, :zfar, :extensions, :extras)

mutable struct Camera
    orthographic::Union{Nothing,Orthographic}
    perspective::Union{Nothing,Perspective}
    type::String
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Camera()
        obj = new()
        obj.orthographic = nothing
        obj.perspective = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Camera, sym::Symbol, x)
    if sym === :type
        x == "perspective" || x == "orthographic" ||
            throw(ArgumentError("""type should be "perspective" or "orthographic" """))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Camera}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Camera}) = (:orthographic, :perspective, :name, :extensions, :extras)


text:
##############################################################################
##
## Iteration: eachrow, eachcol
##
##############################################################################

# Iteration by rows
immutable DFRowIterator
    df::AbstractDataFrame
end
eachrow(df::AbstractDataFrame) = DFRowIterator(df)

Base.start(itr::DFRowIterator) = 1
Base.done(itr::DFRowIterator, i::Int) = i > size(itr.df, 1)
Base.next(itr::DFRowIterator, i::Int) = (DataFrameRow(itr.df, i), i + 1)
Base.size(itr::DFRowIterator) = (size(itr.df, 1), )
Base.length(itr::DFRowIterator) = size(itr.df, 1)
Base.getindex(itr::DFRowIterator, i::Any) = DataFrameRow(itr.df, i)
Base.map(f::Function, dfri::DFRowIterator) = [f(row) for row in dfri]

# Iteration by columns
immutable DFColumnIterator
    df::AbstractDataFrame
end
eachcol(df::AbstractDataFrame) = DFColumnIterator(df)

Base.start(itr::DFColumnIterator) = 1
Base.done(itr::DFColumnIterator, j::Int) = j > size(itr.df, 2)
Base.next(itr::DFColumnIterator, j::Int) = ((names(itr.df)[j], itr.df[j]), j + 1)
Base.size(itr::DFColumnIterator) = (size(itr.df, 2), )
Base.length(itr::DFColumnIterator) = size(itr.df, 2)
Base.getindex(itr::DFColumnIterator, j::Any) = itr.df[:, j]
function Base.map(f::Function, dfci::DFColumnIterator)
    # note: `f` must return a consistent length
    res = DataFrame()
    for (n, v) in eachcol(dfci.df)
        res[n] = f(v)
    end
    res
end


text:
using DynamicGrids, FieldDefaults, FieldMetadata, Test
import DynamicGrids: applyrule, applyinteraction!, applyinteraction,
                     Interaction, @Image, @Graphic, @Output
import FieldMetadata: @description, description, @limits, limits,
                      @flattenable, flattenable, default

struct Double <: CellRule end
applyrule(rule::Double, data, state, index) = state * 2

struct Predation{Keys} <: CellInteraction{Keys} end
Predation(; prey=:prey, predator=:predator) = Predation{(prey, predator)}()
applyinteraction(::Predation, data, (prey, predators), index) = begin
    caught = prey * 0.02 * predators
    conversion = 0.2
    mortality = 0.1
    max(prey - caught, zero(prey)), predators + caught * 0.1 - predators * mortality
end

preyarray = rand(300, 300) .* 20
predatorarray = rand(300, 300) .* 2
init = (prey=preyarray, predator=predatorarray)

output = ArrayOutput(init, 20)
rulesets=(prey=Ruleset(Double()), predator=Ruleset());
interactions=(Predation(prey=:prey, predator=:predator),)
ruleset = MultiRuleset(rulesets=rulesets, interactions=interactions; init=init)
sim!(output, ruleset; init=init, tspan=(1, 20))

# Color processor runs
@Image @Graphic @Output mutable struct TestImageOutput{} <: ImageOutput{T} end
processor=DynamicGrids.ThreeColor(colors=(DynamicGrids.Blue(), DynamicGrids.Red()))
imageoutput = TestImageOutput(init; processor=processor, minval=(0, 0), maxval=(1000, 100), store=true)
DynamicGrids.showframe(::TestImageOutput, ::AbstractSimData, args...) = nothing
sim!(imageoutput, ruleset; init=init, tspan=(1, 20))

output[20] == imageoutput[20]


text:
# function calculate how close a matrix is to identity (up to permutations and rescaling)
using Combinatorics, LinearAlgebra
"""
    err_pd(Ie)

Calculate the error (up to permuataion and rescaling) between Ie and identity matrix:

```math
err = \\mathrm{min}_{P, D} ||PDI_{e} - I||_{F},
```

where ``P`` runs over all permuataion matrices and `` D `` runs over all diagonal non-singular matrices.
``||~||_F`` denotes the Frobeneous norm.

# Arguments
-   `Ie`: the input matrix, supposed to be close to identity matrix.

# Outputs
-   `err`: the error between `Ie` and identity matrix defined as above.
"""
function err_pd(Ie)
    # get dimension
    s = size(Ie, 1)

    # generate all permutations in Sym(s)
    perm = permutations(1:s) |> collect

    # identity matrix
    idm = Matrix{Float64}(I, s, s)

    err = []
    for i in 1:factorial(s)
        Iep = Ie[perm[i],:]
        errtemp = []
        for j in 1: s
            eie = Iep[:,j]
            ei = idm[:,j]
            push!(errtemp, norm((I - eie*eie'/norm(eie)^2)*ei))
        end
        push!(err, norm(errtemp))
    end

    return minimum(err)
end


text:
module MPIStructuredMesh
    using MPI

    include("mesh.jl")
    export DistributedMesh, update_halo!, create_buffers, finalize_mesh, get_rank, get_nelts, nelts_with_halo, local_indices
end


text:
"""
Abstract supertype for quantum error correcting codes.

Subtypes should include the following fields:

    stabilizers::Vector{Vector{Int}}
    logicals::Vector{Vector{Int}}
    pure_errors::Vector{Vector{Int}}

where `stabilizers` are a linearly independent collection of stabilizer generators,
`logicals` are representatives of logical operators ordered as ``X``-type, ``Z``-type, ...
for each logical qubit, and `pure_errors` (or destabilizers) are operators the anticommute
with the corresponding stabilizer generator and commute with all other stabilizer
generators. Each of `stabilizers`, `logicals` and `pure_errors` are collections of
multi-qubit Pauli operators, see [Pauli functions](@ref). The relationships between the
operators can be checked using [`verify`](@ref).

!!! note

    Transformation and contraction functions do not modify `stabilizers`, `logicals` or
    `pure_errors` fields directly but rather return new codes. Therefore these fields can be
    treated as immutable.

"""
abstract type QuantumCode end

"""
    SimpleCode <: QuantumCode

    SimpleCode(name, stabilizers, logicals, pure_errors)

Simple named implementation of a quantum code.

Fields in addition to [`QuantumCode`](@ref):

    name::String

This code contains no information on physical qubit layout but it can be used as a seed code
in a [`TensorNetworkCode`](@ref).
"""
struct SimpleCode <: QuantumCode
    name::String
    stabilizers::Vector{Vector{Int}}
    logicals::Vector{Vector{Int}}
    pure_errors::Vector{Vector{Int}}
end

"""
    SimpleCode()

Construct an empty simple code with empty string name.
"""
SimpleCode() = SimpleCode("", Vector{Int}[], Vector{Int}[], Vector{Int}[])

"""
    SimpleCode(name, stabilizers, logicals)

Construct a simple code with pure errors evaluated using [`find_pure_errors`](@ref).

An `ErrorException` is thrown if the pure errors cannot be evaluated; for example, if the
stabilizers are not linearly independent.

# Examples
```jldoctest
julia> code = SimpleCode(
           "5-qubit code",
           [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]],
           [[1, 1, 1, 1, 1], [3, 3, 3, 3, 3]]
       );

julia> code.pure_errors
4-element Vector{Vector{Int64}}:
 [0, 1, 0, 0, 0]
 [1, 3, 0, 0, 0]
 [3, 1, 0, 0, 0]
 [1, 0, 0, 0, 0]

julia> verify(code)
true
```
"""
function SimpleCode(name, stabilizers, logicals)
    # convert stabilizers for find_pure_errors
    typed_stabilizers = convert(Vector{Vector{Int}}, stabilizers)
    return SimpleCode(name, stabilizers, logicals, find_pure_errors(typed_stabilizers))
end

"""
    CodeGraph

Type aggregating all geometrical data used by [`TensorNetworkCode`](@ref), including
coordinates and `ITensor` indices.

See also [`TensorNetworkCode`](@ref).

!!! note

    `CodeGraph` is not typically accessed directly but rather created, read and updated
    using [Transformation](@ref), [Contraction](@ref) and [Code graph functions](@ref).
"""
struct CodeGraph
    coords::Dict{Int,Vector{T}} where T <: Real
    node_types::Dict{Int,String}
    edge_types::Dict{Set{Int},String}
    node_indices::Dict{Int,Vector{Index{Int}}}
    edge_indices::Dict{Set{Int},Vector{Index{Int}}}
end

"""
    TensorNetworkCode <: QuantumCode

    TensorNetworkCode(stabilizers, logicals, pure_errors, code_graph, seed_codes)

Tensor-network implementation of a quantum code, with associated code graph and seed codes.

Fields in addition to [`QuantumCode`](@ref):

    code_graph::CodeGraph
    seed_codes::Dict{String,SimpleCode}

The `code_graph` contains information on physical qubit layout. The `seed_codes` are simple
codes that are used to build the tensor-network code.

See also [`CodeGraph`](@ref).

!!! note

    The `code_graph` and `seed_codes` fields are not typically accessed directly but rather
    read and updated using [Transformation](@ref), [Contraction](@ref) and
    [Code graph functions](@ref).
"""
struct TensorNetworkCode <: QuantumCode
    stabilizers::Vector{Vector{Int}}
    logicals::Vector{Vector{Int}}
    pure_errors::Vector{Vector{Int}}
    code_graph::CodeGraph
    seed_codes::Dict{String,SimpleCode}
end

# non-public constructor
function TensorNetworkCode(
    code::SimpleCode,
    code_graph::CodeGraph,
    seed_codes::Dict{String,SimpleCode}
)
    return TensorNetworkCode(
        code.stabilizers, code.logicals, code.pure_errors, code_graph, seed_codes
    )
end

"""
    TensorNetworkCode(code::SimpleCode)

Construct a tensor-network code from a seed simple code and initialize the code graph with
with appropriate defaults.

The constructed code shares the `stabilizers`, `logicals` and `pure_errors` with the given
code.

See also [`SimpleCode`](@ref).

# Examples
```jldoctest
julia> simple_code = five_qubit_code();

julia> tn_code = TensorNetworkCode(simple_code);

julia> verify(tn_code)
true
```
"""
function TensorNetworkCode(code::SimpleCode)
    n = num_qubits(code)
    # initialize code graph dictionaries
    coords = Dict{Int,Vector{Float64}}()
    node_types = Dict{Int,String}()
    edge_types = Dict{Set{Int},String}()
    node_indices = Dict{Int,Vector{Index{Int}}}()
    edge_indices = Dict{Set{Int},Vector{Index{Int}}}()
    # assign coords to qubits
    coords[-1] = [0.0, 0.0]
    vec = [0.0, -1.0]
    θ = 2 * pi / n
    R = [cos(θ) sin(θ); -sin(θ) cos(θ)]
    for node in 1:n
        coords[node] = vec
        vec = R * vec
    end
    # label nodes
    node_types[-1] = code.name
    for label in 1:n
        node_types[label] = "physical"
    end
    # label edges
    for node in 1:n
        edge = Set([-1,node])
        edge_types[edge] = "physical"
    end
    # assign indices
    indices = [Index(4, "physical") for _ in 1:n]
    node_indices[-1] = indices
    for label in 1:n
        edge = Set([-1,label])
        edge_indices[edge] = [indices[label]]
    end
    # create code_graph
    code_graph = CodeGraph(coords, node_types, edge_types, node_indices, edge_indices)
    # return new code
    return TensorNetworkCode(
        code.stabilizers,
        code.logicals,
        code.pure_errors,
        code_graph,
        Dict([(code.name, code)])
    )
end

"""
     SimpleCode(code::TensorNetworkCode)

Construct a simple code from a tensor-network code with an empty string name.

The constructed code shares the `stabilizers`, `logicals` and `pure_errors` with the given
code.

See also [`TensorNetworkCode`](@ref).
"""
function SimpleCode(code::TensorNetworkCode)
    #TODO: add name="" as parameter?
    return SimpleCode("", code.stabilizers, code.logicals, code.pure_errors)
end


text:
# MLStyle patches
@active Argument(x) begin
    if x isa Argument
        Some(x.n)
    else
        nothing
    end
end

@active SSAValue(x) begin
    if x isa SSAValue
        Some(x.id)
    else
        nothing
    end
end

@active SlotNumber(x) begin
    if x isa SlotNumber
        Some(x.id)
    else
        nothing
    end
end

@active NewvarNode(x) begin
    if x isa NewvarNode
        Some(x.slot)
    else
        nothing
    end
end

@active NewSSAValue(x) begin
    if x isa NewSSAValue
        Some(x.id)
    else
        nothing
    end
end

@active ReturnNode(x) begin
    if x isa ReturnNode && isdefined(x, :val)
        Some(x.val)
    else
        nothing
    end
end

@active Const(x) begin
    if x isa Const
        Some(x.val)
    else
        nothing
    end
end

@active GotoNode(x) begin
    if x isa GotoNode
        Some(x.label)
    else
        nothing
    end
end

@active GotoIfNot(x) begin
    if x isa GotoIfNot
        x.cond, x.dest
    else
        nothing
    end
end

@active Signature(x) begin
    if x isa Signature && isdefined(x, :atype)
        x.f, x.ft, x.atypes, x.atype
    elseif x isa Signature
        x.f, x.ft, x.atypes
    else
        nothing
    end
end

# move this to Base?
Base.iterate(ic::Core.Compiler.IncrementalCompact) = Core.Compiler.iterate(ic)
Base.iterate(ic::Core.Compiler.IncrementalCompact, st) = Core.Compiler.iterate(ic, st)
Base.getindex(ic::Core.Compiler.IncrementalCompact, idx) = Core.Compiler.getindex(ic, idx)
Base.setindex!(ic::Core.Compiler.IncrementalCompact, v, idx) = Core.Compiler.setindex!(ic, v, idx)

Base.getindex(ic::Core.Compiler.Instruction, idx) = Core.Compiler.getindex(ic, idx)
Base.setindex!(ic::Core.Compiler.Instruction, v, idx) = Core.Compiler.setindex!(ic, v, idx)

Base.getindex(ir::Core.Compiler.IRCode, idx) = Core.Compiler.getindex(ir, idx)
Base.setindex!(ir::Core.Compiler.IRCode, v, idx) = Core.Compiler.setindex!(ir, v, idx)

Base.getindex(ref::UseRef) = Core.Compiler.getindex(ref)
Base.iterate(uses::UseRefIterator) = Core.Compiler.iterate(uses)
Base.iterate(uses::UseRefIterator, st) = Core.Compiler.iterate(uses, st)

@static if VERSION < v"1.8-"
Base.iterate(p::Core.Compiler.Pair) = Core.Compiler.iterate(p)
Base.iterate(p::Core.Compiler.Pair, st) = Core.Compiler.iterate(p, st)
end

Base.getindex(m::Core.Compiler.MethodLookupResult, idx::Int) = Core.Compiler.getindex(m, idx)

Base.length(x::Core.Compiler.BitSet) = Core.Compiler.length(x)
Base.iterate(x::Core.Compiler.BitSet) = Core.Compiler.iterate(x)
Base.iterate(x::Core.Compiler.BitSet, st) = Core.Compiler.iterate(x, st)


text:
@testset "Solver_Wrappers" begin

    include("nls_solver.jl")

end


text:
using DataStructures

module LM
  using CxxWrap
  @wrapmodule(joinpath("LanguageModel","lib","liblanguageModel.so"))
  function __init__()
    @initcxx
  end
end

function modelPrior(next_letter::String,lmModel,lmState)
  alphabet = "abcdefghijklmnopqrstuvwxyz "
  next_letter = lowercase(next_letter)

  logprobs = LM.model(lmModel,lmState,next_letter,alphabet)
  prior = OrderedDict([(alphabet[i], 10^(logprobs[i])) for i in 1:length(alphabet)])
  total = sum(values(prior))
  return OrderedDict([(k==' ' ? :SPACE : Symbol(uppercase(k)), v/total) for (k,v) in prior])
end

keyboardStrings = OrderedDict(:A => "A",
                              :B => "B",
                              :C => "C",
                              :D => "D",
                              :E => "E",
                              :F => "F",
                              :G => "G",
                              :H => "H",
                              :I => "I",
                              :J => "J",
                              :K => "K",
                              :L => "L",
                              :M => "M",
                              :N => "N",
                              :O => "O",
                              :P => "P",
                              :Q => "Q",
                              :R => "R",
                              :S => "S",
                              :T => "T",
                              :U => "U",
                              :V => "V",
                              :W => "W",
                              :X => "X",
                              :Y => "Y",
                              :Z => "Z",
                              :SPACE => "SPACE",
                              :UNDO => "UNDO")

# history of color selections (Int) along with the assignment at the time of selection
ColorSelections = Vector{Tuple{Int,Dict{Symbol,Int}}}

mutable struct Belief
  b::OrderedDict{Symbol,Float64} # belief distribution over all keys
  # these two counts represent beta distribution over user's error rate
  right_color_count::Int # count of how many times the user chose the right color for their letter
  wrong_color_count::Int # count of how many times the user chose the wrong color for their letter
  selections::ColorSelections # history of what colors users chose
end

function Belief(b::OrderedDict{Symbol,Float64},right_color_count::Int,wrong_color_count::Int)
  return Belief(b,right_color_count,wrong_color_count,ColorSelections())
end

BeliefHistory = Vector{Tuple{Symbol,Belief}}

# likelihood - probability of user choosing a color given they want to choose the key l
function colorProbability(color::Int,l::Symbol,belief::Belief,assignment::Dict{Symbol,Int})
  p_right_color = belief.right_color_count / (belief.right_color_count + belief.wrong_color_count)
  if assignment[l] == color
    return p_right_color
  else
    return 1 - p_right_color
  end
end

# probability user will choose a color given our current belief and the current assignment
function colorProbability(color::Int,belief::Belief,assignment::Dict{Symbol,Int})
  return sum(belief.b[l]*colorProbability(color,l,belief,assignment) for l in keys(belief.b))
end

# entropy of user's color choice
function colorEntropy(belief::Belief, assignment::Dict{Symbol,Int})
  p(c) = colorProbability(c,belief,assignment)
  return -sum(p(c)log(2,p(c)) for c in 1:2)
end

function randomAssignments(m)
  assignments = [Dict([l => c for (l,c) in zip(keys(keyboardStrings),digits(n,base=2,pad=28).+1)]) for n in rand(1:((2^27)-1), m)]
end

function updateBelief(belief::Belief, color::Int, assignment::Dict{Symbol,Int})
  push!(belief.selections,(color,copy(assignment)))

  for l in keys(belief.b)
    belief.b[l] = colorProbability(color,l,belief,assignment)*belief.b[l]
  end
  total = sum(values(belief.b))
  map!(x->x/total,values(belief.b))
end

function huffmanAssignment(belief::Belief,assignment::Dict{Symbol,Int})
  vals = unique(values(assignment))
  if length(vals) == 2
    # convert assignment to colors (1 or 2)
    for l in keys(assignment)
      assignment[l] = assignment[l] == vals[1] ? 1 : 2
    end
    return assignment
  else
    probs = OrderedDict{Int,Float64}(v => 0 for v in vals)
    for l in keys(assignment)
      probs[assignment[l]] += belief.b[l]
    end
    sorted_vals = map(first, sort(collect(probs),by=last))
    # merge two least probable values in assignment
    for l in keys(assignment)
      if assignment[l] == sorted_vals[2]
        assignment[l] = sorted_vals[1]
      end
    end

    return huffmanAssignment(belief,assignment)
  end
end

function huffmanAssignment(belief::Belief)
  assignment = Dict{Symbol,Int}()
  c = 1
  for s in keys(belief.b)
    assignment[s] = c
    c += 1
  end
  return huffmanAssignment(belief,assignment)
end

function greedyPartition(belief::Belief)
  sorted_belief = sort(collect(belief.b),rev=true,by=x->x[2])
  assignment = Dict{Symbol,Int}()
  for (sym,_) in sorted_belief
    redSum = sum(vcat([belief.b[k] for (k,v) in assignment if v==1],0))
    blueSum = sum(vcat([belief.b[k] for (k,v) in assignment if v==2],0))
    assignment[sym] = redSum < blueSum ? 1 : 2
  end
  return assignment
end

function changeAssignment(belief::Belief, assignment::Dict{Symbol,Int})
  best = greedyPartition(belief)

  #best = huffmanAssignment(belief)

  for k in keys(assignment)
    assignment[k] = best[k]
  end
end

function getUniformPrior()
  nChoices = length(keyboardStrings)
  prior = OrderedDict{Symbol,Float64}([letter => 1.0/nChoices for letter in keys(keyboardStrings)])
  return prior
end

function getPrior(lmModel,lmState)
  #return getUniformPrior()
  prior = modelPrior("",lmModel,lmState)
  prior[:UNDO] = 0
  return prior
end

function getPrior(commString::String, belief::Belief,selected_letter::Symbol,lmModel,lmState)
  #return getUniformPrior()
  if isempty(commString) #can't do UNDO
    return getPrior(lmState)
  else # use past belief to inform
    nextLetter = selected_letter == :SPACE ? " " : keyboardStrings[selected_letter]
    prior = modelPrior(nextLetter,lmModel,lmState)

    prior[:UNDO] = 1 - belief.b[selected_letter]
    #normalize other probablities
    for key in keys(prior)
      if key != :UNDO
        prior[key] = prior[key] / (1 - prior[:UNDO])
      end
    end
  end
  return prior
end

function learnLikelihood(belief::Belief,selected_letter::Symbol)
    selections = belief.selections
    right_count = sum(color == assignment[selected_letter] for (color,assignment) in selections)
    wrong_count = sum(color != assignment[selected_letter] for (color,assignment) in selections)
    belief.right_color_count += right_count
    belief.wrong_color_count += wrong_count
end

function undo(belief::Belief,history::BeliefHistory)
  (prev_letter,prev_belief) = pop!(history)
  belief.right_color_count = prev_belief.right_color_count
  belief.wrong_color_count = prev_belief.wrong_color_count
  belief.selections = prev_belief.selections

  new_belief = deepcopy(prev_belief.b)
  new_belief[prev_letter] = 1 - belief.b[:UNDO] # the probability that I still wanted that previous letter is the probability my current undo was wrong
  #normalize other probablities
  for key in keys(new_belief)
    if key != prev_letter
      new_belief[key] = new_belief[key] / (1 - new_belief[prev_letter])
    end
  end
  belief.b = new_belief
end

# choose letter IF we are confident, and update belief
function chooseLetter(belief::Belief, commString::String, certaintyThreshold::Float64, history::BeliefHistory,lmModel, lmState)
  selected_letter = findfirst(prob->prob>=certaintyThreshold,belief.b)
  if !isnothing(selected_letter)
    if selected_letter == :UNDO
      commString = commString[1:end-1]
      undo(belief,history)
    else
      nextLetter = selected_letter == :SPACE ? " " : keyboardStrings[selected_letter]
      commString = commString * nextLetter
      prior = getPrior(commString,belief,selected_letter,lmModel,lmState)
      push!(history,(selected_letter,deepcopy(belief)))
      learnLikelihood(belief,selected_letter)
      belief.b = prior
      belief.selections = ColorSelections()
    end
  end
  return commString
end

max_stars_repo_name: mrdaly/ColorCode

max_issues_repo_name: mrdaly/ColorCode

max_forks_repo_name: mrdaly/ColorCode


text:
"""
    function get_thermal_nonzeroslope(args..., kwargs...)

Identify generators which cost function has a non zero slope, locate their buses and return
the corresponding PTDF matrix. Skips Slack generators, generators which Active Power is
binding and Generators with constant price. Note: For a quadratic cost function of the form:
αPg²+ βPg + γ, a non zero slope generator will have the parameter α ≠ 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

"""
function get_thermal_nonzeroslope(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_gen_tol::Float64
    )

    # Define elements
    gen_loc = 0
    gen_bindPg = 0
    gen_nonzero = 0
    gens_thermal_nonzeroslope = Array{Any}(undef, 0)
    gens_thermal_nonzeroslope_id = Array{Int64}(undef, 0)
    gens_thermal_nonzeroslope_name = Array{String}(undef, 0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc + 1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1, gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                gen_bindPg = gen_bindPg +1
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                if α ≠ 0 # Only gens with α(Pg^2)
                    gen_nonzero = gen_nonzero +1
                    resize!(gens_thermal_nonzeroslope, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_id, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_name, gen_nonzero)
                    gens_thermal_nonzeroslope[gen_nonzero] =
                        get_component(ThermalStandard, sys, gen_thermal.name)
                    gens_thermal_nonzeroslope_id[gen_nonzero] = gen_loc
                    gens_thermal_nonzeroslope_name[gen_nonzero] = gen_thermal.name
                end
            end
        end
    end

    return (
        gens_thermal_nonzeroslope,
        gens_thermal_nonzeroslope_id,
        gens_thermal_nonzeroslope_name
    )
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
using aoc2019.computer

input = joinpath(@__DIR__, "input")
p = computer.load_program(input)

struct NIC
    stdin::Channel{Int64}
    stdout::Channel{Int64}
end

function init(p::Array{Int,1}, stdout::Channel{Int64})
    nics = Dict{Int, Channel{Int64}}()
    for i in 0:49
        stdin = Channel{Int64}(Inf)
        nics[i] = stdin
        @async computer.run(copy(p), stdin, stdout)
        put!(stdin, i)
        ## hack to get nics started
        put!(stdin, -1)
    end
    return nics
end

function run(listen::Channel{Int64}, nics::Dict{Int, Channel{Int64}})
    nat = nothing
    first = nothing
    last = nothing
    idle = false
    while true
        yield()
        while isready(listen)
            address = take!(listen)
            x = take!(listen)
            y = take!(listen)
            if address == 255
                nat = (x, y)
                if first === nothing first = y end
            else
                dst = nics[address]
                put!(dst, x)
                put!(dst, y)
            end
            idle = false
        end
        if idle && nat !== nothing
            origin = nics[0]
            x, y = nat
            put!(origin, x)
            put!(origin, y)
            if y === last break end
            last = y
        end
        idle = true
    end
    return (first, last)
end

stdout = Channel{Int64}(Inf)
nics = init(p, stdout)
p1, p2 = run(stdout, nics)
@show p1, p2
@assert p1 == 26779
@assert p2 == 19216

print("-----------------------------------------------------------------------\n")
print("category six -- part one\n    first y : $p1\n")
print("category six -- part two\n    last y  : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
export simulatemodel

function modify_constraint(m::JuMP.AbstractModel, consname::Symbol, data::Array{Float64,2})

    !(size(m[consmane]) == size(data)) ? @error("The data and the constraint are size inconsistent") : true

    for (n, c) in enumerate(IndexCartesian(), data)

        JuMP.setRHS(m[consname], data[n])

    end

    return m

end

function run_simulations(simulation::PowerSimulationsModel{S}, solver, ps_dict::Dict; args...) where {S<:AbstractOperationsModel}
    # TODO: refactor system to be mutable and remove ps_dict from arguments

    # CheckPowerModel(m::PowerSimulationsModel{T}) where T<:AbstractPowerSimulationType
    # AssignSolver(m::PowerSimulationsModel{T}) where T<:AbstractPowerSimulationType
    # WarmUpModel(m::PowerSimulationsModel{T}) where T<:AbstractPowerSimulationType

    # Precalculate PTDF and initialize generator statuses
    sys = simulation.model.system

    if simulation.model.transmission <: StandardPTDF
        PTDF,  A = PowerSystems.buildptdf(sys.branches, sys.buses)
    else
        PTDF = nothing
    end

    name_index = [gen.name for gen in sys.generators.thermal];

    initialpowerdict = :initialpowerdict in keys(args) ? args[:initialpowerdict] : Dict(zip(name_index,[sys.generators.thermal[i].tech.activepower for i in 1:length(name_index)]))
    initialstatusdict = :initialstatusdict in keys(args) ? args[:initialstatusdict] : Dict(zip(name_index,ones(length(name_index))));
    initialondurationdict = :initialondurationdict in keys(args) ? args[:initialondurationdict] : Dict(zip(name_index,ones(length(name_index))*100));
    initialoffdurationdict = :initialoffdurationdict in keys(args) ? args[:initialoffdurationdict] : Dict(zip(name_index,zeros(length(name_index))));

    simulation_results = Dict();

    # run a PCM
    for (step, step_ts) in sort(simulation.timeseries)
        # assign TS to ps_dict
        ps_dict = PowerSystems.assign_ts_data(ps_dict,step_ts);

        # build sys
        sys = PowerSystems.PowerSystem(ps_dict; args...);

        # make model
        tmp_model = PS.PowerOperationModel(simulation.model.psmodel,
            simulation.model.generation,
            simulation.model.demand,
            simulation.model.storage,
            simulation.model.branches,
            simulation.model.transmission,
            simulation.model.services,
            sys,
            Model(), false, PTDF)

        # build model
        println( "Building model for $step ...")
        buildmodel!(sys, tmp_model; PTDF = PTDF,
            initialpower = initialpowerdict,
            initialstatus = initialstatusdict,
            initialonduration = initialondurationdict,
            initialoffduration = initialoffdurationdict)

        # solve model
        println("Solving model for $step ...")
        optimize!(tmp_model.model,with_optimizer(solver))

        status = :termination_status in fieldnames(typeof(tmp_model.model.moi_backend.model.optimizer)) ? tmp_model.model.moi_backend.model.optimizer.termination_status : nothing

        if status == JuMP.MOI.Success
            println( "Problem solved successfully...")
        elseif status == nothing
            @info "$solver does not provide a termination status"
        else
            @warn "Problem solve unsuccessful, solver returned with status: $status"
        end

        # TODO: Subset model results when a lookahead is provided

        # extract results
        res = PS.get_model_result(tmp_model)

        simulation_results[step] = copy(res)

        # update initial...
        initialpowerdict = :p_th in keys(res) ? PS.get_previous_value(res[:p_th]) : nothing
        initialstatusdict = :on_th in keys(res) ? PS.get_previous_value(res[:on_th]) : nothing
        initialondurationdict = :start_th in keys(res) ? PS.commitment_duration(res,initialondurationdict,:start_th) : nothing
        initialoffdurationdict = :stop_th in keys(res) ? PS.commitment_duration(res,initialoffdurationdict,:stop_th) : nothing

    end

    return simulation_results

end;


text:
#
# basics.jl --
#
# Basic types and method for 2-D points and bounding-boxes.
#
#-------------------------------------------------------------------------------
#
# This file is part of the TwoDimensional Julia package licensed under the MIT
# license (https://github.com/emmt/TwoDimensional.jl).
#
# Copyright (C) 2019-2022, Éric Thiébaut.
#

using Base: @propagate_inbounds

"""
    AbstractPoint{T}

is the abstract type of objects with at least 2 properties: `x` and `y`, their
respective abscissa and ordinate, both of type `T`.

See also: [`Point`](@ref), [`WeightedPoint`](@ref).

"""
abstract type AbstractPoint{T<:Real} end

"""
    Point(x,y)

yields an instance of a 2D point of coordinates `(x,y)`.

A point may be multiplied or divided by a scalar to scale its coordinates.  The
addition (resp. subtraction) of two points adds (resp. subtracts) their
coordinates.

Coordinates can be specified by keywords:

    Point(x=xval, y=yval)

There are no default values for keywords `x` and `y` so both must be specified.

The coordinates of a `Point`, say `pnt`, can be retrieved as follows:

    pnt.x  or  pnt[1]  ->  x
    pnt.y  or  pnt[2]  ->  y

or:

    x, y = pnt

See also: [`WeightedPoint`](@ref), [`AbstractPoint`](@ref).

""" Point

struct Point{T} <: AbstractPoint{T}
    x::T
    y::T
end

Base.getindex(obj::Point, i::Integer) =
    (i == 1 ? obj.x :
     i == 2 ? obj.y :
     error("out of range index for `Point` object"))

"""
    WeightedPoint{T}(w,x,y)

yields a weighted point which has 3 fields: `w` its weight, `x` its abscissa
and `y` its ordinate, all of type `T`.  By convention `w ≥ 0` but this is not
checked for efficiency reasons.

See also: [`Point`](@ref), [`AbstractPoint`](@ref).

""" WeightedPoint

struct WeightedPoint{T<:AbstractFloat}  <: AbstractPoint{T}
    w::T # weight
    x::T # abscissa
    y::T # ordinate
end

"""
    BoundingBox(xmin,xmax,ymin,ymax)

yields an instance of a 2D rectangular bounding-box whose sides are aligned
with the coordinate axes and containing points of coordinates `(x,y)` such that
`xmin ≤ x ≤ xmax` and `ymin ≤ y ≤ ymax`.  The box is *empty* if `xmin > xmax`
or `ymin > ymax`.

A bounding-box can be constructed from the first and last points (i.e. at the
lower-left and upper right opposite corners) of the box:

    BoundingBox(P0::Point, P1::Point)
    BoundingBox(I0::CartesianIndex{2}, I1::CartesianIndex{2})

Coordinates can be specified by keywords:

    BoundingBox(xmin=x0, ymin=y0, xmax=x1, ymax=y1)

There are no default values for keywords `xmin`, `xmax`, `ymin` and `ymax` so
all must be specified.

The coordinates of a `BoundingBox`, say `box`, can be retrieved as follows:

    box.xmin  or  box[1]  ->  xmin
    box.xmax  or  box[2]  ->  xmax
    box.ymin  or  box[3]  ->  ymin
    box.ymax  or  box[4]  ->  ymax

or:

    xmin, xmax, ymin, ymax = box

See also [`Point`](@ref), [`interior`](@ref), [`exterior`](@ref).

""" BoundingBox

struct BoundingBox{T<:Real}
    xmin::T
    xmax::T
    ymin::T
    ymax::T
end

Base.getindex(obj::BoundingBox, i::Integer) =
    (i == 1 ? obj.xmin :
     i == 2 ? obj.xmax :
     i == 3 ? obj.ymin :
     i == 4 ? obj.ymax :
     error("out of range index for `BoundingBox` object"))

# Allowed types to construct (or convert to) a point.
const PointTypes = Union{AbstractPoint,NTuple{2,Real},CartesianIndex{2}}

# Constructors of points and conversion to/from a Cartesian index.
Point(P::Point) = P
Point{T}(P::Point{T}) where {T} = P
Point(P::AbstractPoint) = Point(P.x, P.y)
Point{T}(P::AbstractPoint) where {T} = Point{T}(P.x, P.y)
Point(x::Tx, y::Ty) where {Tx<:Real,Ty<:Real} = Point{promote_type(Tx,Ty)}(x, y)
Point(; x::Real, y::Real) = Point(x, y)
Point{T}(; x::Real, y::Real) where {T} = Point{T}(x, y)
Point(I::CartesianIndex{2}) = Point(I[1], I[2])
Point{T}(I::CartesianIndex{2}) where {T} = Point{T}(I[1], I[2])
Point(P::NTuple{2,Real}) = Point(P[1], P[2])
Point{T}(P::NTuple{2,Real}) where {T} = Point{T}(P[1], P[2])

# Conversion to points (rely on constructors).
Base.convert(::Type{T}, arg::T) where {T<:Point} = arg
Base.convert(::Type{T}, arg::PointTypes) where {T<:Point} = T(arg)

# Other basic methods.
Base.CartesianIndex(P::Point{<:Integer}) = CartesianIndex(P.x, P.y)
Base.Tuple(P::Point) = (P.x, P.y)
Base.eltype(::AbstractPoint{T}) where {T} = T
Broadcast.broadcasted(::Type{T}, obj::Point) where {T<:Real} =
    Point{T}(obj)
Base.promote_type(::Type{Point{T}}, ::Type{Point{U}}) where {T,U} =
    Point{promote_type(T,U)}
Base.promote_type(::Type{Point{T}}, ::Type{Point{T}}) where {T} =
    Point{T}
Base.zero(::Type{Point{T}}) where {T} = Point(zero(T),zero(T))
Base.one(::Type{Point{T}}) where {T} = Point(one(T),one(T))
Base.oneunit(::Type{Point{T}}) where {T} = Point(oneunit(T),oneunit(T))

# Constructors of weighted points.
WeightedPoint(P::WeightedPoint) = P
WeightedPoint{T}(P::WeightedPoint{T}) where {T<:AbstractFloat} = P
WeightedPoint{T}(P::WeightedPoint) where {T} = WeightedPoint{T}(P.w, P.x, P.y)
WeightedPoint(w::Tw, x::Tx, y::Ty) where {Tw<:Real,Tx<:Real,Ty<:Real} =
    WeightedPoint{float(promote_type(Tw,Tx,Ty))}(w, x, y)
WeightedPoint(; w::Real, x::Real, y::Real) = WeightedPoint(w, x, y)
WeightedPoint(P::Point{T}) where {T} = WeightedPoint(one(T), P.x, P.y)
WeightedPoint(P::NTuple{3,Real}) = WeightedPoint(P[1], P[2], P[3])
Base.Tuple(P::WeightedPoint) = (P.w, P.x, P.y)
Broadcast.broadcasted(::Type{T}, obj::WeightedPoint) where {T<:AbstractFloat} =
    WeightedPoint{T}(obj)
Base.promote_type(::Type{WeightedPoint{T}}, ::Type{WeightedPoint{U}}) where {T,U} =
    WeightedPoint{promote_type(T,U)}
Base.promote_type(::Type{WeightedPoint{T}}, ::Type{WeightedPoint{T}}) where {T} =
    WeightedPoint{T}

"""
    round([T,] obj::Union{Point,BoundingBox})

yields the object that is the nearest to `obj` by rounding its coordinates to
the nearest integer.  Argument `T` can be the type of the result (a point or a
bounding-box) or the type of the coordinates of the result.

For points, see also: [`floor(::Point)`](@ref), [`ceil(::Point)`](@ref).

For bounding-boxes, see also: [`interior`](@ref), [`exterior`](@ref).

"""
Base.round(obj::Point{T}) where {T} = round(T, obj)
Base.round(::Type{Point{T}}, obj::Point) where {T} = round(T, obj)
Base.round(::Type{T}, obj::Point{T}) where {T<:Integer} = obj
Base.round(::Type{T}, obj::Point{T}) where {T<:Real} =
    Point(round(obj.x, RoundNearest),
          round(obj.y, RoundNearest))
Base.round(::Type{T}, obj::Point{<:Integer}) where {T<:Integer} =
    Point{T}(obj)
Base.round(::Type{T}, obj::Point{<:Real}) where {T<:Integer} =
    Point(round(T, obj.x, RoundNearest),
          round(T, obj.y, RoundNearest))
Base.round(::Type{T}, obj::Point{<:Integer}) where {T<:Real} =
    Point{T}(obj)
Base.round(::Type{T}, obj::Point{U}) where {T<:Real,U<:Real} =
    Point{T}(round(U, obj))

# Extend round for bounding-boxes.
Base.round(obj::BoundingBox{T}) where {T} = round(T, obj)
Base.round(::Type{BoundingBox{T}}, obj::BoundingBox) where {T} = round(T, obj)
Base.round(::Type{T}, obj::BoundingBox{T}) where {T<:Integer} = obj
Base.round(::Type{T}, obj::BoundingBox{T}) where {T<:Real} =
    BoundingBox(round(obj.xmin, RoundNearest),
                round(obj.xmax, RoundNearest),
                round(obj.ymin, RoundNearest),
                round(obj.ymax, RoundNearest))
Base.round(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Integer} =
    BoundingBox{T}(obj)
Base.round(::Type{T}, obj::BoundingBox{<:Real}) where {T<:Integer} =
    BoundingBox(round(T, obj.xmin, RoundNearest),
                round(T, obj.xmax, RoundNearest),
                round(T, obj.ymin, RoundNearest),
                round(T, obj.ymax, RoundNearest))
Base.round(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Real} =
    BoundingBox{T}(obj)
Base.round(::Type{T}, obj::BoundingBox{U}) where {T<:Real,U<:Real} =
    BoundingBox{T}(round(U, obj))

"""
    floor([T,] P::Point)

yields the point with the largest integer coordinates smaller or equal those of
the point `P`.  Argument `T` can be the type of the result or the type of the
coordinates of the result.

See also: [`round`](@ref), [`ceil`](@ref).

"""
Base.floor(obj::Point{T}) where {T} = floor(T, obj)
Base.floor(::Type{Point{T}}, obj::Point) where {T} = floor(T, obj)
Base.floor(::Type{T}, obj::Point{T}) where {T<:Integer} = obj
Base.floor(::Type{T}, obj::Point{T}) where {T<:Real} =
    Point(floor(obj.x), floor(obj.y))
Base.floor(::Type{T}, obj::Point{<:Integer}) where {T<:Integer} =
    Point{T}(obj)
Base.floor(::Type{T}, obj::Point{<:Real}) where {T<:Integer} =
    Point(floor(T, obj.x), floor(T, obj.y))
Base.floor(::Type{T}, obj::Point{<:Integer}) where {T<:Real} =
    Point{T}(obj)
Base.floor(::Type{T}, obj::Point{U}) where {T<:Real,U<:Real} =
    Point{T}(floor(U, obj))

"""
    ceil([T,] P::Point)

yields the point with the smallest integer coordinates larger or equal those of
the point `P`.  Argument `T` can be the type of the result or the type of the
coordinates of the result.

See also: [`round`](@ref), [`floor`](@ref).

"""
Base.ceil(obj::Point{T}) where {T} = ceil(T, obj)
Base.ceil(::Type{Point{T}}, obj::Point) where {T} = ceil(T, obj)
Base.ceil(::Type{T}, obj::Point{T}) where {T<:Integer} = obj
Base.ceil(::Type{T}, obj::Point{T}) where {T<:Real} =
    Point(ceil(obj.x), ceil(obj.y))
Base.ceil(::Type{T}, obj::Point{<:Integer}) where {T<:Integer} =
    Point{T}(obj)
Base.ceil(::Type{T}, obj::Point{<:Real}) where {T<:Integer} =
    Point(ceil(T, obj.x), ceil(T, obj.y))
Base.ceil(::Type{T}, obj::Point{<:Integer}) where {T<:Real} =
    Point{T}(obj)
Base.ceil(::Type{T}, obj::Point{U}) where {T<:Real,U<:Real} =
    Point{T}(ceil(U, obj))

function Base.clamp(P::Point{T1}, B::BoundingBox{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    return Point{T}(clamp(T(P.x), T(B.xmin), T(B.xmax)),
                    clamp(T(P.y), T(B.ymin), T(B.ymax)))
end

# Methods hypot() and atan() yield the polar coordinates of a point.
Base.hypot(P::Point) = hypot(P.x, P.y)
Base.atan(P::Point) = atan(P.y, P.x)

"""
    distance(A, B)

yields the Euclidean distance between the 2 points `A` and `B`.

"""
distance(A::Point, B::Point) =
    distance(promote(A, B)...)
distance(A::Point{T}, B::Point{T}) where {T<:Unsigned} =
    hypot(ifelse(A.x > B.x, A.x - B.x, B.x - A.x),
          ifelse(A.y > B.y, A.y - B.y, B.y - A.y))
distance(A::Point{T}, B::Point{T}) where {T<:Real} =
    hypot(B.x - A.x, B.y - A.y)

# Allowed types to construct (or convert to) a bounding box.
const BoundingBoxTypes = Union{BoundingBox,NTuple{2,AbstractPoint},
                               NTuple{4,Real},NTuple{2,CartesianIndex{2}},
                               NTuple{2,AbstractUnitRange{<:Integer}},
                               CartesianIndices{2}}

# Constructors of bounding-boxes and conversion.
function BoundingBox(xmin::Txmin, xmax::Txmax,
                     ymin::Tymin, ymax::Tymax) where {Txmin,Txmax,Tymin,Tymax}
    T = promote_type(Txmin, Txmax, Tymin, Tymax)
    return BoundingBox{T}(xmin, xmax, ymin, ymax)
end

BoundingBox(; xmin::Real, xmax::Real, ymin::Real, ymax::Real) =
    BoundingBox(xmin, xmax, ymin, ymax)
BoundingBox{T}(; xmin::Real, xmax::Real, ymin::Real, ymax::Real) where {T} =
    BoundingBox{T}(xmin, xmax, ymin, ymax)

BoundingBox(B::BoundingBox) = B
BoundingBox{T}(B::BoundingBox{T}) where {T<:Real} = B
BoundingBox{T}(B::BoundingBox) where {T<:Real} =
    BoundingBox{T}(B.xmin, B.xmax, B.ymin, B.ymax)

BoundingBox(arg::NTuple{2,CartesianIndex{2}}) = BoundingBox(arg...)
BoundingBox{T}(arg::NTuple{2,CartesianIndex{2}}) where {T} =
    BoundingBox{T}(arg...)

BoundingBox(arg::NTuple{2,AbstractPoint}) = BoundingBox(arg...)
BoundingBox{T}(arg::NTuple{2,AbstractPoint}) where {T} = BoundingBox{T}(arg...)

BoundingBox(I0::CartesianIndex{2}, I1::CartesianIndex{2}) =
    BoundingBox{Int}(I0, I1)
BoundingBox{T}(I0::CartesianIndex{2}, I1::CartesianIndex{2}) where {T} =
    BoundingBox{T}(I0[1], I1[1], I0[2], I1[2])

BoundingBox(P0::AbstractPoint, P1::AbstractPoint) =
    BoundingBox(P0.x, P1.x, P0.y, P1.y)
BoundingBox{T}(P0::AbstractPoint, P1::AbstractPoint) where {T} =
    BoundingBox{T}(P0.x, P1.x, P0.y, P1.y)

BoundingBox(P0::NTuple{2,Real}, P1::NTuple{2,Real}) =
    BoundingBox(P0[1], P1[1], P0[2], P1[2])
BoundingBox{T}(P0::NTuple{2,Real}, P1::NTuple{2,Real}) where {T} =
    BoundingBox{T}(P0[1], P1[1], P0[2], P1[2])

BoundingBox(P::NTuple{4,Real}) = BoundingBox(P[1], P[2], P[3], P[4])
BoundingBox{T}(P::NTuple{4,Real}) where {T} =
    BoundingBox{T}(P[1], P[2], P[3], P[4])

BoundingBox(inds::NTuple{2,AbstractUnitRange{<:Integer}}) =
    BoundingBox(inds[1], inds[2])
BoundingBox{T}(inds::NTuple{2,AbstractUnitRange{<:Integer}}) where {T} =
    BoundingBox{T}(inds[1], inds[2])

BoundingBox(X::AbstractUnitRange{<:Integer}, Y::AbstractUnitRange{<:Integer}) =
    BoundingBox{Int}(X, Y)
BoundingBox{T}(X::AbstractUnitRange{<:Integer}, Y::AbstractUnitRange{<:Integer}) where {T} =
    BoundingBox{T}(first(X), last(X), first(Y), last(Y))

BoundingBox(R::CartesianIndices{2}) = BoundingBox(first(R), last(R))
BoundingBox{T}(R::CartesianIndices{2}) where {T} = BoundingBox{T}(first(R), last(R))

@deprecate BoundingBox(A::AbstractMatrix) BoundingBox(axes(A))
BoundingBox(A::AbstractMatrix{Bool}) = BoundingBox(identity, A)

# Conversion to bounding-boxes (rely on constructors).
Base.convert(::Type{T}, arg::T) where {T<:BoundingBox} = arg
Base.convert(::Type{T}, arg::BoundingBoxTypes) where {T<:BoundingBox} = T(arg)

# Other basic methods.
Base.eltype(::BoundingBox{T}) where {T} = T
Base.Tuple(B::BoundingBox) = (B.xmin, B.xmax, B.ymin, B.ymax)
Broadcast.broadcasted(::Type{T}, obj::BoundingBox) where {T<:Real} =
    BoundingBox{T}(obj)
Base.promote_type(::Type{BoundingBox{T}}, ::Type{BoundingBox{U}}) where {T,U} =
    BoundingBox{promote_type(T,U)}
Base.promote_type(::Type{BoundingBox{T}}, ::Type{BoundingBox{T}}) where {T} =
    BoundingBox{T}

# See
# https://stackoverflow.com/questions/9852159/calculate-bounding-box-of-arbitrary-pixel-based-drawing
# for the basic ideas under the following algorithm.
function BoundingBox(f::Function, A::AbstractMatrix)
    I, J = axes(A)
    i0, i1 = getaxisbounds(I)
    j0, j1 = getaxisbounds(J)
    imin = jmin = typemax(Int)
    imax = jmax = typemin(Int)
    # Assuming column-major order, first start by scanning rows to narrow the
    # subsequent searches along columns.
    #
    # 1. Find bottom bound `jmin` by scanning rows from bottom to top.
    flag = false
    @inbounds for j in j0:j1, i in i0:i1
        if f(A[i,j])
            # This definitively set the value of `jmin` and gives limits for
            # the other bounds.
            imin = imax = i
            jmin = jmax = j
            flag = true
            break
        end
    end
    if flag
        # 2. Find top bound `jmax` by scanning rows from top to bottom.  No
        #    needs to go beyond `jmax+1`.
        @inbounds for j in j1:-1:jmax+1, i in i0:i1
            if f(A[i,j])
                jmax = j
                imin = min(imin, i)
                imax = max(imax, i)
                break
            end
        end
        # 3. Find leftmost bound `imin` by scanning columns from left to right.
        #    No needs to go beyond `imin-1`.
        @inbounds for i in i0:imin-1, j in jmin:jmax
            if f(A[i,j])
                imin = i
                imax = max(imax, i)
                break
            end
        end
        # 4. Find rightmost bound `imax` by scanning columns from right to
        #    left.  No needs to go beyond `imax+1`.
        @inbounds for i in i1:-1:imax+1, j in jmin:jmax
            if f(A[i,j])
                imax = i
                break
            end
        end
    end
    return BoundingBox(imin, imax, jmin, jmax)
end

"""
    getaxisbounds(I) = (i0,i1)

yields the bounds `i0` and `i1` of index range `I` as a 2-tuple of `Int`'s and
such that `i0:i1` represents the same indices as `I` (although not in the same
order if `step(I) < 0`).  If `step(I)` is not equal to ±1, an `ArgumentError`
exception is thrown.

"""
@inline function getaxisbounds(I::AbstractRange{<:Integer})
    i0, i1, s = Int(first(I)), Int(last(I)), step(I)
    return (s == +1 ? (i0,i1) :
            s == -1 ? (i1,i0) :
            throw(ArgumentError("expecting a range with a step equal to ±1")))
end

# Empty bounding and unlimited boxes.
BoundingBox{T}(::Nothing) where {T<:Real} = typemin(BoundingBox{T})
Base.typemin(::Type{BoundingBox{T}}) where {T<:Real} =
    BoundingBox(typemax(T), typemin(T), typemax(T), typemin(T))
Base.typemax(::Type{BoundingBox{T}}) where {T<:Real} =
    BoundingBox(typemin(T), typemax(T), typemin(T), typemax(T))

# Conversion of bounding-boxes to/from Cartesian indices.
Base.CartesianIndices(B::BoundingBox{<:Integer}) =
    CartesianIndices((Int(B.xmin):Int(B.xmax), Int(B.ymin):Int(B.ymax)))

# Lower left and upper right corners of a bounding-box.
Base.first(B::BoundingBox) = Point(B.xmin, B.ymin)
Base.last(B::BoundingBox) = Point(B.xmax, B.ymax)

Base.isempty(B::BoundingBox) = ((B.xmin > B.xmax)|(B.ymin > B.ymax))

Base.size(B::BoundingBox{Int}) =
    (isempty(B) ? (0, 0) : (B.xmax - B.xmin + 1,
                            B.ymax - B.ymin + 1))
Base.size(B::BoundingBox{Int}, k::Integer) =
    (isempty(B) ? (k ≥ 1 ? 0 : throw_bad_dimension_index()) :
     k == 1 ? B.xmax - B.xmin + 1 :
     k == 2 ? B.ymax - B.ymin + 1 :
     k > 2 ? 1 : throw_bad_dimension_index())

Base.size(B::BoundingBox{<:Integer}) =
    (isempty(B) ? (0, 0) : (Int(B.xmax) - Int(B.xmin) + 1,
                            Int(B.ymax) - Int(B.ymin) + 1))

Base.size(B::BoundingBox{<:Integer}, k::Integer) =
    (isempty(B) ? (k ≥ 1 ? 0 : throw_bad_dimension_index()) :
     k == 1 ? Int(B.xmax) - Int(B.xmin) + 1 :
     k == 2 ? Int(B.ymax) - Int(B.ymin) + 1 :
     k > 2 ? 1 : throw_bad_dimension_index())

Base.axes(B::BoundingBox{Int}) = (B.xmin:B.xmax, B.ymin:B.ymax)
Base.axes(B::BoundingBox{Int}, k::Integer) =
    (k == 1 ? (B.xmin:B.xmax) :
     k == 2 ? (B.ymin:B.ymax) :
     k > 2 ? (1:1) : throw_bad_dimension_index())

Base.axes(B::BoundingBox{<:Integer}) = (Int(B.xmin):Int(B.xmax),
                                        Int(B.ymin):Int(B.ymax))
Base.axes(B::BoundingBox{<:Integer}, k::Integer) =
    (k == 1 ? (Int(B.xmin):Int(B.xmax)) :
     k == 2 ? (Int(B.ymin):Int(B.ymax)) :
     k > 2 ? (1:1) : throw_bad_dimension_index())

@noinline throw_bad_dimension_index() =
    error("invalid dimension index")

function Base.iterate(itr::Union{AbstractPoint,BoundingBox})
    vals = Tuple(itr)
    return (vals[1], (vals, 2))
end

function Base.iterate(itr::Union{AbstractPoint,BoundingBox},
                      state::Tuple{Tuple,Int})
    vals, i = state
    return (i ≤ length(vals) ? (vals[i], (vals, i+1)) : nothing)
end

Base.in(pnt::AbstractPoint, box::BoundingBox) =
    ((box.xmin ≤ pnt.x ≤ box.xmax)&
     (box.ymin ≤ pnt.y ≤ box.ymax))

# To deprecate `A ∈ B` in favor of `A ⊆ B` for bounding boxes `A` and `B`, the
# `Base.in` method must be imported.
import Base: in
@deprecate in(A::BoundingBox, B::BoundingBox) (A ⊆ B) false

# Extend ⊆ operator.
Base.issubset(A::BoundingBox, B::BoundingBox) =
    (isempty(A)|((A.xmin ≥ B.xmin)&(A.xmax ≤ B.xmax)&
                 (A.ymin ≥ B.ymin)&(A.ymax ≤ B.ymax)))
# Union of bounding-boxes:
Base.:(∪)(A::BoundingBox, B::BoundingBox) =
    BoundingBox(min(A.xmin, B.xmin), max(A.xmax, B.xmax),
                min(A.ymin, B.ymin), max(A.ymax, B.ymax))

# Intersection of bounding-boxes:
Base.:(∩)(A::BoundingBox, B::BoundingBox) =
    BoundingBox(max(A.xmin, B.xmin), min(A.xmax, B.xmax),
                max(A.ymin, B.ymin), min(A.ymax, B.ymax))

# Use bounding-boxes to extract a sub-array or a view.
@propagate_inbounds function Base.getindex(A::AbstractMatrix,
                                           B::BoundingBox{<:Integer})
    A[B.xmin:B.xmax, B.ymin:B.ymax]
end

Base.view(A::AbstractMatrix, B::BoundingBox{<:Integer}) =
    view(A, B.xmin:B.xmax, B.ymin:B.ymax)

"""
    interior([T,] B)

yields the largest bounding-box with integer valued bounds and which is
contained by the bounding-box `B`.  Optional argument `T` is to specify the
type of the result or of the coordinates of the result which is the same as `B`
by default.

See also: [`exterior`](@ref), [`round`](@ref).

"""
interior(obj::BoundingBox{T}) where {T} = interior(T, obj)
interior(::Type{BoundingBox{T}}, obj::BoundingBox) where {T} = interior(T, obj)
interior(::Type{T}, obj::BoundingBox{T}) where {T<:Integer} = obj
interior(::Type{T}, obj::BoundingBox{T}) where {T<:Real} =
    BoundingBox(ceil(obj.xmin), floor(obj.xmax),
                ceil(obj.ymin), floor(obj.ymax))
interior(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Integer} =
    BoundingBox{T}(obj)
interior(::Type{T}, obj::BoundingBox{<:Real}) where {T<:Integer} =
    BoundingBox(ceil(T, obj.xmin), floor(T, obj.xmax),
                ceil(T, obj.ymin), floor(T, obj.ymax))
interior(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Real} =
    BoundingBox{T}(obj)
interior(::Type{T}, obj::BoundingBox{U}) where {T<:Real,U<:Real} =
    BoundingBox{T}(interior(U, obj))

"""
    exterior([T,] B)

yields the smallest bounding-box with integer valued bounds and which contains
the bounding-box `B`.  Optional argument `T` is to specify the type of the
result or of the coordinates of the result which is the same as `B` by default.

See also: [`interior`](@ref), [`round`](@ref).

"""
exterior(obj::BoundingBox{T}) where {T} = exterior(T, obj)
exterior(::Type{BoundingBox{T}}, obj::BoundingBox) where {T} = exterior(T, obj)
exterior(::Type{T}, obj::BoundingBox{T}) where {T<:Integer} = obj
exterior(::Type{T}, obj::BoundingBox{T}) where {T<:Real} =
    BoundingBox(floor(obj.xmin), ceil(obj.xmax),
                floor(obj.ymin), ceil(obj.ymax))
exterior(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Integer} =
    BoundingBox{T}(obj)
exterior(::Type{T}, obj::BoundingBox{<:Real}) where {T<:Integer} =
    BoundingBox(floor(T, obj.xmin), ceil(T, obj.xmax),
                floor(T, obj.ymin), ceil(T, obj.ymax))
exterior(::Type{T}, obj::BoundingBox{<:Integer}) where {T<:Real} =
    BoundingBox{T}(obj)
exterior(::Type{T}, obj::BoundingBox{U}) where {T<:Real,U<:Real} =
    BoundingBox{T}(exterior(U, obj))

"""
    center(B::BoundingBox) -> c::Point

yields the central point of the bounding-box `B`.

"""
center(B::BoundingBox{<:Integer}) =
    Point(0.5*(B.xmin + B.xmax), 0.5*(B.ymin + B.ymax))
center(B::BoundingBox{T}) where {T<:AbstractFloat} =
    Point(half(T)*(B.xmin + B.xmax), half(T)*(B.ymin + B.ymax))

half(::Type{T}) where {T<:AbstractFloat} = one(T)/convert(T, 2)

"""
    area(B)

yields the area of the bounding-box `B`.

"""
area(B::BoundingBox{T}) where {T<:Real} =
    max(B.xmax - B.xmin, zero(T))*max(B.ymax - B.ymin, zero(T))

# Scaling of a point coordinates.
Base.:(*)(P::Point, α::Real) = α*P
Base.:(*)(α::Real, P::Point) = Point(α*P.x, α*P.y)
Base.:(/)(P::Point, α::Real) = Point(P.x/α, P.y/α)
Base.:(\)(α::Real, P::Point) = P/α

# Unary plus does nothing.
Base.:(+)(P::Point) = P
Base.:(+)(B::BoundingBox) = B

# Unary minus applied to a point negate coordinates.
Base.:(-)(P::Point{<:Union{AbstractFloat,Signed,Irrational}}) =
    Point(-P.x, -P.y)

# Unary minus applied to a bounding-box.
Base.:(-)(B::BoundingBox{<:Union{AbstractFloat,Signed,Irrational}}) =
    BoundingBox(-B.xmax, -B.xmin, -B.ymax, -B.ymin)

# Addition and subtraction of point coordinates.
Base.:(+)(A::Point, B::Point) = Point(A.x + B.x, A.y + B.y)
Base.:(-)(A::Point, B::Point) = Point(A.x - B.x, A.y - B.y)

# Scaling of bounding-box bounds (e.g. to change units).
Base.:(*)(B::BoundingBox, α::Real) = α*B
Base.:(*)(α::Real, B::BoundingBox) = BoundingBox(α*B.xmin, α*B.xmax,
                                                 α*B.ymin, α*B.ymax)
Base.:(/)(B::BoundingBox{T}, α::Real) where {T} = (one(T)/α)*B
Base.:(\)(α::Real, B::BoundingBox) = B/α

# Add or remove a margin δ to a bounding-box B.
Base.:(+)(B::BoundingBox, δ::Real) =
    BoundingBox(B.xmin - δ, B.xmax + δ,
                B.ymin - δ, B.ymax + δ)
Base.:(-)(B::BoundingBox, δ::Real) =
    BoundingBox(B.xmin + δ, B.xmax - δ,
                B.ymin + δ, B.ymax - δ)

# Translate a bounding-box.
Base.:(+)(B::BoundingBox, P::Point) =
    BoundingBox(B.xmin + P.x, B.xmax + P.x,
                B.ymin + P.y, B.ymax + P.y)
Base.:(-)(B::BoundingBox, P::Point) =
    BoundingBox(B.xmin - P.x, B.xmax - P.x,
                B.ymin - P.y, B.ymax - P.y)


text:
ontol = parseontology("c:\\home\\resource\\geneontology\\go-basic.obo")
annot = load_annot_long("c:\\home\\resource\\Xt\\9.1\\go_b2g_lf.map", ontol)
gographs = ontologygraphs(ontol)
exannot = extendannotation(annot, gographs)

### enrichment tests
function loadclusters(file="c:\\home\\projects\\pot\\gc\\pot_xen_tpm_gc_corrected_gp_smoothed.tsv")
    CSV.read(file)
end
data = loadclusters()
#
background = data.Gene
genesel = data.ClusterBothBIC .== 6
enrichtest(data.ClusterBothBIC .== 6, background, exannot, ontol) |> display

rexannot = reverseannotation(exannot)
rannot = reverseannotation(annot)

Z = enrichtest(data.ClusterBothBIC .== 6, background, exannot, ontol)

for row in eachrow(Z)
    id = goid_to_num(row.GT)
    for g in row.Genes
        @assert id ∈ exannot["P"][g]
    end
end


text:
using Test
using FinancialToolbox, DualNumbers

#Test Parameters
testToll = 1e-14;

#print_colored("Starting Standard Test\n",:green)
spot = 10;
K = 10;
r = 0.02;
T = 2.0;
sigma = 0.2;
d = 0.01;
spotDual = dual(spot, 1.0);

toll = 1e-6

#EuropeanCall Option
PriceCall = blsprice(spot, K, r, T, sigma, d);
PriceCallBlack = blkprice(spot, K, r, T, sigma);

df(f, x, dx) = (f(x + dx) - f(x)) / dx
dx = 1e-8;
assert_(value, toll) = @test abs(value) < toll

sigma1 = blsimpv(spot, K, r, T, PriceCall, d);
ResDual = blsimpv(spotDual, K, r, T, PriceCall, d);
DerDF_ = df(spot -> blsimpv(spot, K, r, T, PriceCall, d), spot, dx);

assert_(sigma1 - ResDual, toll)
assert_(DerDF_ - ResDual.epsilon, toll)

sigma1 = blkimpv(spot, K, r, T, PriceCallBlack);
ResDual = blkimpv(spotDual, K, r, T, PriceCallBlack);
DerDF_ = df(spot -> blkimpv(spot, K, r, T, PriceCallBlack), spot, dx);

assert_(sigma1 - ResDual, toll)
assert_(DerDF_ - ResDual.epsilon, toll)

print_colored("Starting DualNumbers Test\n", :green)

#Test Parameters
testToll = 1e-14;
spot = 10;
K = 10;
r = 0.02;
T = 2.0;
sigma = 0.2;
d = 0.01;

#EuropeanCall Option
PriceCall = blsprice(spot, K, r, T, sigma, d);
DeltaCall = blsdelta(spot, K, r, T, sigma, d);
ThetaCall = blstheta(spot, K, r, T, sigma, d);
RhoCall = blsrho(spot, K, r, T, sigma, d);
LambdaCall = blslambda(spot, K, r, T, sigma, d);
VannaCall = blsvanna(spot, K, r, T, sigma, d);

#EuropeanPut Option
PricePut = blsprice(spot, K, r, T, sigma, d, false);
DeltaPut = blsdelta(spot, K, r, T, sigma, d, false);
ThetaPut = blstheta(spot, K, r, T, sigma, d, false);
RhoPut = blsrho(spot, K, r, T, sigma, d, false);
LambdaPut = blslambda(spot, K, r, T, sigma, d, false);
VannaPut = blsvanna(spot, K, r, T, sigma, d, false);

#Equals for both Options
Gamma = blsgamma(spot, K, r, T, sigma, d);
Vega = blsvega(spot, K, r, T, sigma, d);

########DUAL NUMBERS
DerToll = 1e-13;
#Function definition
#Call
FcallDual(spot) = blsprice(spot, K, r, T, sigma, d);
GcallDual(r) = blsprice(spot, K, r, T, sigma, d);
HcallDual(T) = blsprice(spot, K, r, T, sigma, d);
LcallDual(sigma) = blsprice(spot, K, r, T, sigma, d);
PcallDual(spot) = blsprice(spot, K, r, T, sigma, d) * spot.value / blsprice(spot.value, K, r, T, sigma, d);
VcallDual(sigma) = blsdelta(spot, K, r, T, sigma, d);
#Put
FputDual(spot) = blsprice(spot, K, r, T, sigma, d, false);
GputDual(r) = blsprice(spot, K, r, T, sigma, d, false);
HputDual(T) = blsprice(spot, K, r, T, sigma, d, false);
PputDual(spot) = blsprice(spot, K, r, T, sigma, d, false) * spot.value / blsprice(spot.value, K, r, T, sigma, d, false);
VPutDual(sigma) = blsdelta(spot, K, r, T, sigma, d, false);
#Input
SpotDual = Dual(spot, 1.0);
rDual = Dual(r, 1.0);
TDual = Dual(T, 1.0);
SigmaDual = Dual(sigma, 1.0);

#Automatic Differentiation Test
#TEST
print_colored("--- European Call Sensitivities: DualNumbers\n", :yellow)
print_colored("-----Testing Delta\n", :blue);
assert_(FcallDual(SpotDual).epsilon - DeltaCall, DerToll)
print_colored("-----Testing Rho\n", :blue);
assert_(GcallDual(rDual).epsilon - RhoCall, DerToll)
print_colored("-----Testing Theta\n", :blue);
assert_(-HcallDual(TDual).epsilon - ThetaCall, DerToll)
print_colored("-----Testing Lambda\n", :blue);
assert_(PcallDual(SpotDual).epsilon - LambdaCall, DerToll)
print_colored("-----Testing Vanna\n", :blue);
assert_(VcallDual(SigmaDual).epsilon - VannaCall, DerToll)
#TEST
print_colored("--- European Put Sensitivities: DualNumbers\n", :yellow)
print_colored("-----Testing Delta\n", :blue);
assert_(FputDual(SpotDual).epsilon - DeltaPut, DerToll)
print_colored("-----Testing Rho\n", :blue);
assert_(GputDual(rDual).epsilon - RhoPut, DerToll)
print_colored("-----Testing Theta\n", :blue);
assert_(-HputDual(TDual).epsilon - ThetaPut, DerToll)
print_colored("-----Testing Lambda\n", :blue);
assert_(PputDual(SpotDual).epsilon - LambdaPut, DerToll)
print_colored("-----Testing Vanna\n", :blue);
assert_(VPutDual(SigmaDual).epsilon - VannaPut, DerToll)

print_colored("-----Testing Vega\n", :blue);
assert_(LcallDual(SigmaDual).epsilon - Vega, DerToll)
print_colored("Dual Numbers Test Passed\n", :green)
println("")

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Dual\n", :magenta)
print_colored("----Testing Negative  Spot Price \n", :cyan)
@test_throws(ErrorException, blsprice(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blkprice(-SpotDual, K, r, T, sigma))
@test_throws(ErrorException, blsdelta(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blsgamma(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blstheta(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blsrho(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blspsi(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blslambda(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blsvanna(-SpotDual, K, r, T, sigma, d))
@test_throws(ErrorException, blsvega(-SpotDual, K, r, T, sigma, d))

print_colored("----Testing Negative  Strike Price \n", :cyan)
KK = Dual(K, 0);
@test_throws(ErrorException, blsprice(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blkprice(spot, -KK, r, T, sigma))
@test_throws(ErrorException, blsdelta(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blsgamma(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blstheta(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blsrho(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blspsi(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blslambda(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blsvanna(spot, -KK, r, T, sigma, d))
@test_throws(ErrorException, blsvega(spot, -KK, r, T, sigma, d))

print_colored("----Testing Negative  Time to Maturity \n", :cyan)
@test_throws(ErrorException, blsprice(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blkprice(spot, K, r, -TDual, sigma))
@test_throws(ErrorException, blsdelta(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blsgamma(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blstheta(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blsrho(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blspsi(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blslambda(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blsvanna(spot, K, r, -TDual, sigma, d))
@test_throws(ErrorException, blsvega(spot, K, r, -TDual, sigma, d))

print_colored("----Testing Negative  Volatility \n", :cyan)
@test_throws(ErrorException, blsprice(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blkprice(spot, K, r, T, -SigmaDual))
@test_throws(ErrorException, blsdelta(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blsgamma(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blstheta(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blsrho(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blspsi(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blslambda(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blsvanna(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blsvega(spot, K, r, T, -SigmaDual, d))
@test_throws(ErrorException, blsimpv(spotDual, K, r, T, -0.2, d))

print_colored("Dual Input Validation Test Passed\n", :magenta)


text:
@testset "Maximal Independent Set" begin

    g3 = StarGraph(5)
    for g in testgraphs(g3)
        z = @inferred(independent_set(g, MaximalIndependentSet()))
        @test (length(z)== 1 || length(z)== 4)
    end

    g4 = CompleteGraph(5)
    for g in testgraphs(g4)
        z = @inferred(independent_set(g, MaximalIndependentSet()))
        @test length(z)== 1 #Exactly one vertex
    end

    g5 = PathGraph(5)
    for g in testgraphs(g5)
        z = @inferred(independent_set(g, MaximalIndependentSet()))
        sort!(z)
        @test (z == [2, 4] || z == [2, 5] || z == [1, 3, 5] || z == [1, 4])
    end
end


text:
using LazySets.Approximations: overapproximate

for N in [Float64, Float32, Rational{Int}]
    # Cartesian Product of a centered 1D BallInf and a centered 2D BallInf
    # Here a 3D BallInf
    b1 = BallInf(N[0], N(1))
    b2 = BallInf(N[0, 0], N(1))
    # Test Construction
    cp = CartesianProduct(b1, b2)
    @test cp.X == b1
    @test cp.Y == b2
    # Test Dimension
    @test dim(cp) == 3
    # Test Support Vector
    d = N[1, 1, 1]
    @test σ(d, cp) == N[1, 1, 1]
    d = N[1, 1, -1]
    @test σ(d, cp) == N[1, 1, -1]
    d = N[1, -1, 1]
    @test σ(d, cp) == N[1, -1, 1]
    d = N[1, -1, -1]
    @test σ(d, cp) == N[1, -1, -1]
    d = N[-1, 1, 1]
    @test σ(d, cp) == N[-1, 1, 1]
    d = N[-1, 1, -1]
    @test σ(d, cp) == N[-1, 1, -1]
    d = N[-1, -1, 1]
    @test σ(d, cp) == N[-1, -1, 1]
    d = N[-1, -1, -1]
    @test σ(d, cp) == N[-1, -1, -1]

    # boundedness
    @test isbounded(cp)
    @test !isbounded(Singleton(N[1]) * HalfSpace(N[1], N(1)))

    # isempty
    @test !isempty(cp)

    # Cartesian Product of a not-centered 1D BallInf and a not-centered 2D BallInf
    # Here a Hyperrectangle where c = [1, -3, 4] and r = [3, 2, 2]
    b1 = BallInf(N[1], N(3))
    b2 = BallInf(N[-3, 4], N(2))
    # Test Construction
    c = CartesianProduct(b1, b2)
    # Test Dimension
    @test dim(c) == 3
    # Test Support Vector
    d = N[1, 1, 1]
    @test σ(d, c) == N[4, -1, 6]
    d = N[1, 1, -1]
    @test σ(d, c) == N[4, -1, 2]
    d = N[1, -1, 1]
    @test σ(d, c) == N[4, -5, 6]
    d = N[1, -1, -1]
    @test σ(d, c) == N[4, -5, 2]
    d = N[-1, 1, 1]
    @test σ(d, c) == N[-2, -1, 6]
    d = N[-1, 1, -1]
    @test σ(d, c) == N[-2, -1, 2]
    d = N[-1, -1, 1]
    @test σ(d, c) == N[-2, -5, 6]
    d = N[-1, -1, -1]
    @test σ(d, c) == N[-2, -5, 2]

    # Test Cartesian Product with EmptySet
    s = Singleton(N[1])
    E = EmptySet{N}()
    cs1 = E * s
    cs2 = s * E
    @test cs1 isa EmptySet
    @test cs2 isa EmptySet

    # Test Cartesian Product of an array
    # 0-elements
    as = LazySet{N}[]
    cs = CartesianProduct(as)
    @test cs isa EmptySet
    # 1-element
    as = [Singleton(N[1])]
    cs = CartesianProduct(as)
    @test cs.element == N[1]
    # 3-elements
    as = [Singleton(N[1]), Singleton(N[2]), Singleton(N[3])]
    cs = CartesianProduct(as)
    @test cs.X.element == N[1]
    @test cs.Y.X.element == N[2]
    @test cs.Y.Y.element == N[3]

    # Test containment with respect to CartesianProduct
    p1 = HPolygon{N}()
    addconstraint!(p1, LinearConstraint(N[2, 2], N(12)))
    addconstraint!(p1, LinearConstraint(N[-3, 3], N(6)))
    addconstraint!(p1, LinearConstraint(N[-1, -1], N(0)))
    addconstraint!(p1, LinearConstraint(N[2, -4], N(0)))
    p2 = HPolygon{N}()
    addconstraint!(p2, LinearConstraint(N[1, 0], N(1)))
    addconstraint!(p2, LinearConstraint(N[-1, 0], N(1)))
    addconstraint!(p2, LinearConstraint(N[0, 1], N(1)))
    addconstraint!(p2, LinearConstraint(N[0, -1], N(1)))
    cp = CartesianProduct(p1, p2)

    @test ∈(N[0, 0, 0, 0], cp)
    @test ∈(N[4, 2, 1, 0], cp)
    @test ∈(N[2, 4, -1, 0], cp)
    @test ∈(N[-1, 1, 0.5, 0.7], cp)
    @test ∈(N[2, 3, -0.8, 0.9], cp)
    @test ∈(N[1, 1, -1, 0], cp)
    @test ∈(N[3, 2, 0, 1], cp)
    @test ∈(N[5 / 4, 7 / 4, 1, 1], cp)
    @test !∈(N[4, 1, 0, 0], cp)
    @test !∈(N[5, 2, 0, 0], cp)
    @test !∈(N[3, 4, 0, 0], cp)
    @test !∈(N[-1, 5, 0, 0], cp)
    @test !∈(N[4, 2, 3, 1], cp)
    @test !∈(N[2, 3, 3, 1], cp)
    @test !∈(N[0, 0, 3, 1], cp)
    @test !∈(N[1, 1, 3, 1], cp)

    # vertices_list
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    vlist = vertices_list(CartesianProduct(i1, i2))
    @test ispermutation(vlist, [N[0, 2], N[0, 3], N[1, 2], N[1, 3]])

    #constraints_list
    hlist = constraints_list(CartesianProduct(i1, i2))
    @test ispermutation(hlist,
        [LinearConstraint(sparsevec(N[1], N[1], 2), N(1)),
        LinearConstraint(sparsevec(N[1], N[-1], 2), N(0)),
        LinearConstraint(sparsevec(N[2], N[1], 2),  N(3)),
        LinearConstraint(sparsevec(N[2], N[-1], 2), N(-2))])
    @test all(H -> dim(H) == 2, hlist)

    # ==================================
    # Conversions of Cartesian Products
    # ==================================

    # convert the cartesian product of two hyperrectangles to one hyperrectangle
    h1 = Hyperrectangle(N[1/2],  N[1/2])
    h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
    H = convert(Hyperrectangle, h1 × h2)
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    # make one of the sets a BallInf, to test that the conversion does
    # not depend on the Hyperrectangle type
    b1 = BallInf(N[1/2],  N(1/2))
    H = convert(Hyperrectangle, b1 × h2)
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    # convert a hyperrectangle to a cartesian product of intervals
    H = Hyperrectangle(N[0, 0], N[1, 1])
    Hcp = convert(CartesianProduct{N, Interval{N}, Interval{N}}, H)
    @test Hcp isa CartesianProduct &&
          Hcp.X == Interval(N(-1), N(1)) &&
          Hcp.Y == Interval(N(-1), N(1))

    # checking that the order is correct
    H = Hyperrectangle(N[0, 1], N[1, 3])
    Hcp = convert(CartesianProduct{N, Interval{N}, Interval{N}}, H)
    @test Hcp isa CartesianProduct &&
          Hcp.X == Interval(N(-1), N(1)) &&
          Hcp.Y == Interval(N(-2), N(4))

    # =====================
    # CartesianProductArray
    # =====================

    # relation to base type (internal helper functions)
    @test LazySets.array_constructor(CartesianProduct) == CartesianProductArray
    @test LazySets.is_array_constructor(CartesianProductArray)

    # standard constructor
    v = Vector{LazySet{N}}()
    push!(v, Singleton(N[1, 2]))
    push!(v, Singleton(N[3, 4]))
    cpa = CartesianProductArray(v)

    # constructor with size hint and type
    CartesianProductArray(10, N)

    # array getter
    @test array(cpa) ≡ v

    # boundedness
    @test isbounded(cpa)
    @test !isbounded(CartesianProductArray([Singleton(N[1]), HalfSpace(N[1], N(1))]))

    # membership
    @test ∈(N[1, 2, 3, 4], cpa)
    @test !∈(N[3, 4, 1, 2], cpa)

    # in-place modification
    b = BallInf(N[0, 0], N(2))
    cpa = CartesianProductArray(LazySet{N}[])
    @test CartesianProduct!(b, b) isa CartesianProduct &&
          length(array(cpa)) == 0
    res = CartesianProduct!(b, cpa)
    @test res isa CartesianProductArray && length(array(cpa)) == 1
    res = CartesianProduct!(cpa, b)
    @test res isa CartesianProductArray && length(array(cpa)) == 2
    res = CartesianProduct!(cpa, cpa)
    @test res isa CartesianProductArray && length(array(cpa)) == 4

    # vertices_list
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    i3 = Interval(N[4, 5])
    vlist = vertices_list(CartesianProductArray([i1, i2, i3]))
    @test ispermutation(vlist, [N[0, 2, 4], N[0, 3, 4], N[1, 2, 4],
        N[1, 3, 4], N[0, 2, 5], N[0, 3, 5], N[1, 2, 5], N[1, 3, 5]])

    #constraints_list
    hlist = constraints_list(CartesianProductArray([i1, i2, i3]))
    @test ispermutation(hlist,
        [LinearConstraint(sparsevec(N[1], N[1], 3), N(1)),
        LinearConstraint(sparsevec(N[1], N[-1], 3), N(0)),
        LinearConstraint(sparsevec(N[2], N[1], 3), N(3)),
        LinearConstraint(sparsevec(N[2], N[-1], 3), N(-2)),
        LinearConstraint(sparsevec(N[3], N[1], 3), N(5)),
        LinearConstraint(sparsevec(N[3], N[-1], 3), N(-4)),
        ])
    @test all(H -> dim(H) == 3, hlist)

    # ========================================
    # Conversions of Cartesian Product Arrays
    # ========================================

    # conversion of the cartesian product array of intervals to a hyperrectangle
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    i3 = Interval(N[4, 5])
    H = convert(Hyperrectangle, CartesianProductArray([i1, i2, i3]))
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    M = hcat(N(1))
    h1 = Hyperrectangle(N[1/2],  N[1/2])
    h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
    H = convert(Hyperrectangle, CartesianProductArray([h1, h2]))
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    # same block structure
    aX = [Ball1(N[2], N(1)), Ball1(N[0, 0], N(2))]
    aY = [Interval(N(1), N(3)), BallInf(N[0, 0], N(3))]
    @test LazySets.same_block_structure(aX, aY)

    # inclusion
    cpa1 = CartesianProductArray(aX)
    cpa2 = CartesianProductArray(aY)
    res, w = ⊆(cpa1, cpa2, true)
    @test cpa1 ⊆ cpa2 && res && w == N[]
    res, w = ⊆(cpa2, cpa1, true)
    @test !⊆(cpa2, cpa1) && !res && w ∈ cpa2 && w ∉ cpa1

    # convert a hyperrectangle to the cartesian product array of intervals
    # convert a hyperrectangle to a cartesian product of intervals
    H = Hyperrectangle(N[0, 0, 0], N[1, 1, 1])
    Hcp = convert(CartesianProductArray{N, Interval{N}}, H)
    @test Hcp isa CartesianProductArray &&
          all([Hcpi == Interval(N(-1), N(1)) for Hcpi in array(Hcp)])

    # check that the order is correct
    H = Hyperrectangle(N[0, 1, 2], N[1, 2, 3])
    Hcp = convert(CartesianProductArray{N, Interval{N}}, H)
    @test Hcp isa CartesianProductArray &&
          all([array(Hcp)[i] == Interval(N(-1), N(2*i-1)) for i in 1:3])

    # ================
    # common functions
    # ================

    # absorbing element
    e = EmptySet{N}()
    b = BallInf(N[0, 0], N(2))
    @test absorbing(CartesianProduct) == absorbing(CartesianProductArray) ==
          EmptySet
    @test b × e == e × b == cpa × e == e × cpa == e × e == e
end

for N in [Float64, Float32]
    # is_intersection_empty
    i1 = Interval(N[0, 1])
    i2 = Interval(N[2, 3])
    h1 = Hyperrectangle(low=N[3, 4], high=N[5, 7])
    h2 = Hyperrectangle(low=N[5, 5], high=N[6, 8])
    cpa1 = CartesianProductArray([i1, i2, h1])
    cpa2 = CartesianProductArray([i1, i2, h2])
    G = HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(1))])
    G_empty = HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(-1))])
	cpa1_box = overapproximate(cpa1)
	cpa2_box = overapproximate(cpa2)

    @test is_intersection_empty(cpa1, cpa2) == is_intersection_empty(cpa1_box, cpa2_box) == false
    @test is_intersection_empty(cpa1, G_empty) == is_intersection_empty(cpa1_box, G_empty) == true
    @test is_intersection_empty(cpa1, G) == is_intersection_empty(Approximations.overapproximate(cpa1), G) == false
end


text:
module GDALUtils

    import GDAL

    include("types.jl")
    include("metadata.jl")
    include("driver.jl")
    include("gcp.jl")
    include("spatialref.jl")
    include("dataset.jl") # depends: types.jl, driver.jl
    include("raster/rasterband.jl") # depends: dataset.jl
    include("raster/rasterio.jl")
    include("ogr/geometry.jl")
    include("ogr/featurelayer.jl")
    include("ogr/featuredefn.jl")
    include("ogr/fielddefn.jl")
    include("base/display.jl")
    include("base/iteration.jl")

end # module