text:
using PyPlot

function fig_elementgeometry()
    fontsize = 10

    srcx = [-2, 2]
    srcy = [0, 0]
    obsx = -1.5
    obsy = -1.5

    angle = -35
    dx = 3.0
    dy = 2.3
    rotmat = [cosd(angle) -sind(angle); sind(angle) cosd(angle)]
    obsxglobal, obsyglobal = rotmat * [obsx ; obsy] + [dx ; dy]
    tx1, ty1 = rotmat * [srcx[1] ; srcy[1]] + [dx ; dy]
    tx2, ty2 = rotmat * [srcx[2] ; srcy[2]] + [dx ; dy]
    srcxglobal = [tx1, tx2]
    srcyglobal = [ty1, ty2]

    close("all")
    fig = figure(figsize=(5, 5))
    text(5.00, -0.30, L"$x$", fontsize=fontsize, ha="center", va="center")
    text(-0.30, 4.00, L"$y$", fontsize=fontsize, ha="center", va="center")

    plot(srcx, srcy, "-r", linewidth = 7, zorder=30, solid_capstyle="round")
    plot(srcx, srcy, ".w", linewidth = 0, zorder=30, solid_capstyle="round")
    text(-2.00, -0.5, L"$(x^{src}_1=-a, y^{src}_1=0)$", color="red", fontsize=fontsize, ha="center", va="center")
    text(2.00, -0.5, L"$(x^{src}_2=a, y^{src}_2=0)$", color="red", fontsize=fontsize, ha="center", va="center")
    plot(obsx, obsy, "ro", linewidth = 4, zorder=30)
    text(obsx, obsy-0.4, L"$(x^{obs}, y^{obs})$", color="red", fontsize=fontsize, ha="center", va="center")

    plot(srcxglobal, srcyglobal, "-b", linewidth = 7, zorder=30, solid_capstyle="round")
    plot(srcxglobal, srcyglobal, ".w", linewidth = 0, zorder=30, solid_capstyle="round")
    text(srcxglobal[1], srcyglobal[1] + 0.50, L"$(x^{src}_1, y^{src}_1)$", color="blue", fontsize=fontsize, ha="center", va="center")
    text(srcxglobal[2], srcyglobal[2] - 0.50, L"$(x^{src}_2, y^{src}_2)$", color="blue", fontsize=fontsize, ha="center", va="center")
    plot(obsxglobal, obsyglobal, "bo", linewidth = 4, zorder=30)
    text(obsxglobal, obsyglobal-0.5, L"$(\hat{x}^{obs}, \hat{y}^{obs})$", color="blue", fontsize=fontsize, ha="center", va="center")

    ax = gca()
    ax.set_xlim([-5, 5])
    ax.set_ylim([-4, 4])
    ax.set_xticks([])
    ax.set_yticks([])
    ax.spines["top"].set_visible(false) # Hide the top edge of the axis
    ax.spines["right"].set_visible(false) # Hide the right edge of the axis
    ax.spines["left"].set_position("center") # Move the right axis to the center
    ax.spines["bottom"].set_position("center") # Most the bottom axis to the center
    gca().set_aspect("equal")
    show()
end
fig_elementgeometry()


text:
using DataStructures

abstract type DiscretizeMODL <: DiscretizationAlgorithm end

struct DiscretizeMODL_Optimal <: DiscretizeMODL end
struct DiscretizeMODL_Greedy <: DiscretizeMODL end
struct DiscretizeMODL_PostGreedy <: DiscretizeMODL
    max_bin_count :: Int
    DiscretizeMODL_PostGreedy(max_bin_count::Integer = 0) = new(round(Int, max_bin_count))
end

function binedges(
    alg             :: DiscretizeMODL,
    data_continuous :: AbstractArray{N},
    data_class      :: AbstractArray{I},
    ) where {N<:Real, I<:Integer}

    @assert(length(data_continuous)==length(data_class))
    p = sortperm(data_continuous)
    data_continuous = data_continuous[p]
    data_class = data_class[p]

    if isa(alg, DiscretizeMODL_Optimal)
        return optimal_result(data_continuous,data_class)
    elseif isa(alg, DiscretizeMODL_Greedy)
        return greedy_merge(data_continuous,data_class)
    else
        @assert(isa(alg, DiscretizeMODL_PostGreedy))
        return post_greedy_result(data_continuous,data_class,alg.max_bin_count)
    end
end

function MODL_value2_oneintval(
    continuous_ij  :: AbstractArray{T},
    class_value_ij :: AbstractArray{S},
    class_uniq     :: AbstractArray{S},
    ) where {T<:AbstractFloat, S<:Integer}

    n = length(continuous_ij)
    J = length(class_uniq)

    first_part = lfactorial(n+J-1) - lfactorial(J-1) - lfactorial(n)
    second_part = lfactorial(n)
    # Note(Yi-Chun): lfactorial(n) is log(n!)

    for J_index = 1:J
            n_J = count(x->x==class_uniq[J_index],class_value_ij)
            second_part = second_part - lfactorial(n_J)
    end

    return first_part+second_part
end

function optimal_result(
    continuous      :: AbstractArray{T},
    discrete_target :: AbstractArray{S}
    ) where {T<:AbstractFloat, S<:Integer}

    n = length(continuous)

    Disc_ijk_retval = Array{Array}(undef,n,n)
    Disc_ijk_MODL_value = Array{AbstractFloat}(undef,n,n)
    for j = 1:n
            for k = 1:n
                    if k > j
                            Disc_ijk_MODL_value[j,k] = Inf
                            Disc_ijk_retval[j,k] = [0]
                    end
            end
    end

    class_uniq = unique(discrete_target)

    for k = 1:n
        for j = k:n
            if k == 1
                Disc_ijk_retval[j,k] = [j]
                Disc_ijk_MODL_value[j,k] = MODL_value2_oneintval(
                                      continuous[1:j], discrete_target[1:j],class_uniq)
            else
                MODL_value = Inf
                select_intval = 0

                for i = 1:j-1
                    second_MODL = MODL_value2_oneintval(
                                  continuous[i+1:j],discrete_target[i+1:j],class_uniq)
                    current_MODL = Disc_ijk_MODL_value[i,k-1] + second_MODL
                    if current_MODL < MODL_value
                        MODL_value = current_MODL
                        select_intval = i
                    end
                end

                Disc_ijk_retval[j,k] = append!(copy(Disc_ijk_retval[select_intval,k-1]),[j])
                Disc_ijk_MODL_value[j,k] = MODL_value
            end

        end
    end

    full_length_k_intval = copy(Disc_ijk_MODL_value[n,:])

    for l = 1:n
        full_length_k_intval[l] += lfactorial(n+l-1) - lfactorial(l-1) - lfactorial(n)
    end

    desired_intval_number = argmin(full_length_k_intval)

    bin_edges_index = Disc_ijk_retval[n,desired_intval_number]
    bin_edges = append!([1],bin_edges_index)

    bin_edge_value = Vector{AbstractFloat}(undef,length(bin_edges))

    for index = 1 : length(bin_edge_value)
        if index == 1
                bin_edge_value[index] = continuous[1]
        elseif index == length(bin_edge_value)
                bin_edge_value[index] = continuous[end]
        else
                bin_edge_value[index] = 0.5*(continuous[bin_edges[index]]+
                                             continuous[bin_edges[index]+1])
        end
    end
    return bin_edge_value
end

function merge_adj_intval(
    A_distr :: AbstractArray{S}, # Note(Yi-Chun): Distribution of class values in A
    B_distr :: AbstractArray{S}  # Note(Yi-Chun): Distribution of class values in B
    ) where {S<:Integer}

    @assert(length(A_distr) == length(B_distr))
    J = length(B_distr)

    n_A = sum(A_distr)
    n_B = sum(B_distr)

    Delta = lfactorial(n_A+n_B+J-1) + lfactorial(J-1) - lfactorial(n_A+J-1) - lfactorial(n_B+J-1)

    for j = 1:J
        n_A_J = A_distr[j]
        n_B_J = B_distr[j]
        Delta = Delta - lfactorial(n_A_J+n_B_J) + lfactorial(n_A_J) + lfactorial(n_B_J)
    end
    Delta
end

function greedy_merge_index(
    continuous  :: AbstractArray{T},
    class_value :: AbstractArray{S},
    ) where {T<:AbstractFloat,S<:Integer}

    @assert(length(continuous) == length(class_value))

    n = length(continuous)
    all_class_value = unique(class_value)
    J = length(all_class_value)

    if n==1
        return [1]
    end

    pq = DataStructures.PriorityQueue()
    # Note(Yi-Chun): The PriorQueue sorts differences of MODL after merging of adjacent intervals
    distr_in_intval = Dict()
    # Note(Yi-Chun): This dictionary contains the distribution in each interval
    adj_for_intval = Dict()
    # Note(Yi-Chun): This dictionary contains the two adjacent interval for interval i

    MODL = log(n) + lfactorial(2n-1) - lfactorial(n-1) - lfactorial(n) + n*log(J-1)
    # Note(Yi-Chun): Intial MODL when each attribute is treated as a single bin

    for i = 1:n-1
        A_Distrib_in_intval = zeros(Int64,J)
        A_value_index = something(findfirst(isequal(class_value[i]), all_class_value), 0)

        A_Distrib_in_intval[A_value_index] = 1
        # Note(Yi-Chun): Set up distribution in each bin as [0,0,..0,1,0,..0]
        distr_in_intval[i] = A_Distrib_in_intval
        # Note(Yi-Chun): Here we label intervals by the index of their first boundary.

        # Note(Yi-Chun): For each interval we assign 4 adjacent intervals for it.
        #                Two in the front and two in the behind.
        if i==1
            adj_for_intval[i] = [2,3]
        elseif i==2
            adj_for_intval[i] = [1,3,4]
        elseif i==n-1
            adj_for_intval[i] = [i-2,i-1,i+1]
        else
            adj_for_intval[i] = [i-2,i-1,i+1,i+2]
        end

        B_Distrib_in_intval = zeros(Int64,J)
        B_value_index = something(findfirst(isequal(class_value[i+1]), all_class_value), 0)

        B_Distrib_in_intval[B_value_index] = 1
        if i == (n-1) # Note(Yi-Chun): The adjacents of the last interval are n-2 and n-1
            distr_in_intval[n] = B_Distrib_in_intval
            adj_for_intval[i+1] = [i-1,i]
        end

        X = merge_adj_intval(A_Distrib_in_intval,B_Distrib_in_intval)
        pq[i+1] = X
    end

    n_intval = n
    first_could_remove = 2
    last_could_remove = n

    while n_intval > 4
        min_diff = DataStructures.peek(pq)[2]
        first_term_in_delta = log((n_intval-1)/(n+n_intval-1))

        if (min_diff + first_term_in_delta) < 0
            MODL = MODL + min_diff + first_term_in_delta
            removed = DataStructures.dequeue!(pq)
            if removed == first_could_remove
                # Note(Yi-Chun): Interval labels: y=1,(removed),z,w,s
                @assert(length(adj_for_intval[removed]) == 3)
                @assert(adj_for_intval[removed][1] == 1)
                z = adj_for_intval[removed][2]
                w = adj_for_intval[removed][3]
                first_could_remove = z
                # Note(Yi-Chun): Merge the distribution
                Merge_distr = distr_in_intval[1] + distr_in_intval[removed]
                # Note(Yi-Chun): Update the Delta value for the only one adjacent
                pq[z] = merge_adj_intval(Merge_distr,distr_in_intval[z])
                # Note(Yi-Chun): Update the structure
                distr_in_intval[1] = Merge_distr

                # Note(Yi-Chun): Update the adjacent structure
                adj_for_intval[1][1] = z
                adj_for_intval[1][2] = w
                s = adj_for_intval[z][4]
                adj_for_intval[z] = [1,w,s]
                adj_for_intval[w][1] = 1

                # Note(Yi-Chun): Delete the information of removed one
                delete!(adj_for_intval,removed)
                delete!(distr_in_intval,removed)

            elseif removed == last_could_remove
                @assert(length(adj_for_intval[removed]) == 2)
                x = adj_for_intval[removed][1]
                y = adj_for_intval[removed][2]
                last_could_remove = y
                Merge_distr = distr_in_intval[y] + distr_in_intval[removed]
                pq[y] = merge_adj_intval(distr_in_intval[x],Merge_distr)
                distr_in_intval[y] = Merge_distr
                pop!(adj_for_intval[x])
                pop!(adj_for_intval[y])
                delete!(adj_for_intval,removed)
                delete!(distr_in_intval,removed)

            else
                # Note(Yi-Chun): Interval labels: x,y,(removed),z,w
                # Note(Yi-Chun): Merge the distribution
                y = adj_for_intval[removed][2]
                Merge_distr = distr_in_intval[y] + distr_in_intval[removed]
                # Note(Yi-Chun): Update the Delta value for two adjacents
                x = adj_for_intval[removed][1]
                pq[y] = merge_adj_intval(distr_in_intval[x],Merge_distr)
                z= adj_for_intval[removed][3]
                pq[z] = merge_adj_intval(Merge_distr,distr_in_intval[z])
                # Note(Yi-Chun): Update the structure
                distr_in_intval[y] = Merge_distr

                # Note(Yi-Chun): Update the adjacent structure
                if length(adj_for_intval[removed]) == 3
                    adj_for_intval[y][3] = z
                    pop!(adj_for_intval[y])
                    adj_for_intval[x][end] = z
                    adj_for_intval[z][1] = x
                    adj_for_intval[z][2] = y
                else
                    w = adj_for_intval[removed][4]
                    adj_for_intval[x][end] = z
                    adj_for_intval[y][end-1] = z
                    adj_for_intval[y][end] = w
                    adj_for_intval[z][1] = x
                    adj_for_intval[z][2] = y
                    adj_for_intval[w][1] = y
                end

                # Note(Yi-Chun): Delete the information of removed one
                delete!(adj_for_intval,removed)
                delete!(distr_in_intval,removed)
            end
            n_intval = n_intval - 1
        else
            break
        end
    end
    # Note(Yi-Chun): Here we manage the case that when bin number goes under 4,
    #                since we have to use different label for intervals

    if n_intval > 4
        binedges = sort(collect(keys(distr_in_intval)))
        return binedges
    else
        binedges = sort(collect(keys(distr_in_intval)))
        while ((DataStructures.peek(pq)[2] + log((n_intval-1)/(n+n_intval-1)))<0) && (n_intval>2)
            MODL = MODL + DataStructures.peek(pq)[2] + log((n_intval-1)/(n+n_intval-1))
            n_intval = n_intval-1
            removed = DataStructures.dequeue!(pq)
            prior = binedges[something(findfirst(isequal(removed), binedges), 0)-1]

            Merge_distr = distr_in_intval[removed] + distr_in_intval[prior]
            distr_in_intval[prior] = Merge_distr
            delete!(adj_for_intval,removed)
            delete!(distr_in_intval,removed)
            binedges = sort(collect(keys(distr_in_intval)))

            for i = 1 : n_intval-1
                X = distr_in_intval[binedges[i]]
                Y = distr_in_intval[binedges[i+1]]
                pq[binedges[i+1]] = merge_adj_intval(X,Y)
            end
        end

        if (DataStructures.peek(pq)[2] + log((n_intval-1)/(n+n_intval-1)))<0
            return [1]
        else
            return binedges
        end
    end
end

function greedy_merge(
    continuous  :: AbstractArray{T},
    class_value :: AbstractArray{S},
    ) where {T<:AbstractFloat,S<:Integer}

    binedges = greedy_merge_index(continuous,class_value)
    append!(binedges,[length(continuous)])
    n = length(binedges)
    bin_edges = Vector{Float64}(undef,n)
    for i = 2: n-1
        binedges[i] = binedges[i]-1
    end

    for j = 1: n
        if j == 1
            bin_edges[1] = continuous[binedges[1]]
        elseif j==n
            bin_edges[n] = continuous[binedges[n]]
        else
            bin_edges[j] = 0.5*(continuous[binedges[j]] + continuous[binedges[j]+1])
        end
    end

    return bin_edges
end

# Note(Yi-Chun): The following 4 functions are used in post_greedy methods.

function methods_split(continuous,class_value,distr,i,j,uniq_class,I)

    N = length(continuous)
    continuous_ij = continuous[i:j-1]
    class_ij = class_value[i:j-1]
    distr_ij = distr[i]

    n = length(continuous_ij)
    J = length(uniq_class)
    distr_A = zeros(Int64,J)

    current_MODL = Inf
    current_split = 0
    current_distr_A = zeros(Int64,J)
    current_distr_B = distr_ij

    for n_A_end = 1:n-1
        distr_A[something(findfirst(isequal(class_ij[n_A_end]), uniq_class),0)] += 1
        distr_B = distr_ij - distr_A

        Del = lfactorial(n_A_end+J-1)+lfactorial(n-n_A_end+J-1)-lfactorial(n+J-1)-lfactorial(J-1)

        for j = 1:J
            Del += lfactorial(distr_A[j]+distr_B[j])-lfactorial(distr_A[j])-lfactorial(distr_B[j])
        end

        if Del < current_MODL
            current_MODL = Del
            current_split = n_A_end
            current_distr_A = deepcopy(distr_A)
            current_distr_B = deepcopy(distr_B)
        end
    end

    current_MODL += log((N+I)/(I))
    current_split += i

    return (current_MODL,current_split,current_distr_A,current_distr_B)
end

function methods_merge(
    A_distr   :: AbstractArray{S},
    B_distr   :: AbstractArray{S},
    N,
    I,
    ) where {S<:Integer}

    # Note(Yi-Chun): This function is similar to the previous function merge_adj_intval.
    #                The only different part is the difference of MODL value after merging.
    @assert(length(A_distr) == length(B_distr))
    J = length(B_distr)
    n_A = sum(A_distr)
    n_B = sum(B_distr)
    Delta = lfactorial(n_A+n_B+J-1) + lfactorial(J-1) - lfactorial(n_A+J-1) - lfactorial(n_B+J-1)
    for j = 1:J
        n_A_J = A_distr[j]
        n_B_J = B_distr[j]
        Delta = Delta - lfactorial(n_A_J+n_B_J) + lfactorial(n_A_J) + lfactorial(n_B_J)
    end
    Delta += log((I-1)/(N+I-1))
    return Delta
end

function methods_merge_split(continuous,class_value,distr,i,j,k,uniq_class,I)

    continuous_ij = continuous[i:j-1]
    class_ij = class_value[i:j-1]
    continuous_jk = continuous[j:k-1]
    class_jk = class_value[j:k-1]
    distr_ij = distr[i]
    distr_jk = distr[j]
    N = length(continuous)
    Del_1 = methods_merge(distr_ij,distr_jk,N,I)

    pseudo_dict = Dict()
    pseudo_dict[i] = distr_ij + distr_jk
    split_result = methods_split(continuous,class_value,pseudo_dict,i,k,uniq_class,I-1)

    Del = Del_1 + split_result[1]
    return (Del,split_result[2],split_result[3],split_result[4])
end

function methods_merge2_split(continuous,class_value,distr,i,j,k,l,uniq_class,I)

    N = length(continuous)
    continuous_ij = continuous[i:j-1]
    continuous_jk = continuous[j:k-1]
    continuous_kl = continuous[k:l-1]
    class_ij = class_value[i:j-1]
    class_jk = class_value[j:k-1]
    class_kl = class_value[k:l-1]
    distr_ij = distr[i]
    distr_jk = distr[j]
    distr_kl = distr[k]
    Del_1 = methods_merge(distr_ij,distr_jk,N,I)
    Del_2 = methods_merge(distr_ij + distr_jk,distr_kl,N,I-1)
    pseudo_dict = Dict()
    pseudo_dict[i] = distr_ij + distr_jk + distr_kl
    split_result = methods_split(continuous,class_value,pseudo_dict,i,l,uniq_class,I-2)

    Del = Del_1 + Del_2 + split_result[1]

    return (Del,split_result[2],split_result[3],split_result[4])
end

function PQ_methods_updata(continuous,class_value,uniq_class,binedges,distr,post_pq,effects,index,method)
    # Note(Yi-Chun): (index,method) is the index of the edge that we want to update its method in priorityqueue
    I = length(binedges) - 1
    # Note(Yi-Chun): Since in post_greedy part, we add the index N+1 as ending of bin edges for convenience
    N = length(continuous)

    if method == "M"
        A_distr = distr[binedges[index]]
        B_distr = distr[binedges[index+1]]
        post_pq[(binedges[index],method)] = methods_merge(A_distr,B_distr,N,I)
    elseif method == "S"
        method_result = methods_split(continuous,class_value,distr,binedges[index],binedges[index+1],uniq_class,I)
        post_pq[(binedges[index],method)] = method_result[1]
        effects[(binedges[index],method)] = method_result[2:end]
    elseif method == "MS"
        method_result = methods_merge_split(continuous,class_value,distr,binedges[index],
                                            binedges[index+1],binedges[index+2],uniq_class,I)
        post_pq[(binedges[index],method)] = method_result[1]
        effects[(binedges[index],method)] = method_result[2:end]
    else
        method_result = methods_merge2_split(continuous,class_value,distr,binedges[index],binedges[index+1],
                                             binedges[index+2],binedges[index+3],uniq_class,I)
        post_pq[(binedges[index],method)] = method_result[1]
        effects[(binedges[index],method)] = method_result[2:end]
    end
end

function remove_methods_on_index(index,methods_on_index,post_pq,effects)
    meth = methods_on_index[index]
    for i = 1: length(meth)
        post_pq[(index,meth[i])] = Inf
        delete!(effects,(index,meth[i]))
    end
end

insert_and_dedup!(v::Vector, x) = (splice!(v, searchsorted(v,x), [x]); v)

function uncondi_greedy_merge_index(
    continuous  :: AbstractArray{T},
    class_value :: AbstractArray{S}
    ) where {T<:AbstractFloat,S<:Integer}
    # Note(Yi-Chun): This function processes the unconditional bin merge, Adapted from the function greedy_merge
    @assert(length(continuous) == length(class_value))

    n = length(continuous)
    all_class_value = unique(class_value)
    J = length(all_class_value)
    length_of_Greedy_merge = length(greedy_merge_index(continuous,class_value))

    pq = DataStructures.PriorityQueue()
    distr_in_intval = Dict()
    adj_for_intval = Dict()

    MODL = log(n) + lfactorial(2n-1) - lfactorial(n-1) - lfactorial(n) + n*log(J-1)

    total_distrib = zeros(Int64,J)

    for i = 1:n-1
        A_Distrib_in_intval = zeros(Int64,J)
        A_value_index = something(findfirst(isequal(class_value[i]), all_class_value), 0)
        A_Distrib_in_intval[A_value_index] = 1
        distr_in_intval[i] = A_Distrib_in_intval

        if i==1
            adj_for_intval[i] = [2,3]
        elseif i==2
            adj_for_intval[i] = [1,3,4]
        elseif i==n-1
            adj_for_intval[i] = [i-2,i-1,i+1]
        else
            adj_for_intval[i] = [i-2,i-1,i+1,i+2]
        end

        total_distrib += A_Distrib_in_intval

        B_Distrib_in_intval = zeros(Int64,J)
        B_value_index = something(findfirst(isequal(class_value[i+1]), all_class_value), 0)
        B_Distrib_in_intval[B_value_index] = 1
        if i == (n-1)
            distr_in_intval[n] = B_Distrib_in_intval
            adj_for_intval[i+1] = [i-1,i]
            total_distrib += B_Distrib_in_intval
        end

        X = merge_adj_intval(A_Distrib_in_intval,B_Distrib_in_intval)
        pq[i+1] = X
    end

    one_bin_distrib = Dict()
    one_bin_distrib[1] = total_distrib

    if n==1|J==1
        best_for_far = [1]
        best_for_far_MODL = -Inf
        return (best_for_far,one_bin_distrib)
    end

    n_intval = n
    first_could_remove = 2
    last_could_remove = n

    # Note(Yi-Chun): This best_so_far part will record the best dicretization we meet while merging unconditionally
    best_so_far = collect(1:n)
    best_so_far_MODL = Inf
    best_adj = Dict()
    best_distr = Dict()

    while n_intval > 4
        min_diff = DataStructures.peek(pq)[2]
        first_term_in_delta = log((n_intval-1)/(n+n_intval-1))

        MODL = MODL + min_diff + first_term_in_delta
        removed = DataStructures.dequeue!(pq)
        if removed == first_could_remove
            @assert(length(adj_for_intval[removed]) == 3)
            @assert(adj_for_intval[removed][1] == 1)
            z = adj_for_intval[removed][2]
            w = adj_for_intval[removed][3]
            first_could_remove = z
            # Note(Yi-Chun): Merge the distribution
            Merge_distr = distr_in_intval[1] + distr_in_intval[removed]
            # Note(Yi-Chun): Update the Delta value for the only one adjacent
            pq[z] = merge_adj_intval(Merge_distr,distr_in_intval[z])
            # Note(Yi-Chun): Update the structure
            distr_in_intval[1] = Merge_distr
            # Note(Yi-Chun): Update the pointers
            adj_for_intval[1][1] = z
            adj_for_intval[1][2] = w
            s = adj_for_intval[z][4]
            adj_for_intval[z] = [1,w,s]
            adj_for_intval[w][1] = 1
            # Note(Yi-Chun): Delete the information of removed one
            delete!(adj_for_intval,removed)
            delete!(distr_in_intval,removed)
        elseif removed == last_could_remove
            @assert(length(adj_for_intval[removed]) == 2)
            x = adj_for_intval[removed][1]
            y = adj_for_intval[removed][2]
            last_could_remove = y
            Merge_distr = distr_in_intval[y] + distr_in_intval[removed]
            pq[y] = merge_adj_intval(distr_in_intval[x],Merge_distr)
            distr_in_intval[y] = Merge_distr
            pop!(adj_for_intval[x])
            pop!(adj_for_intval[y])
            delete!(adj_for_intval,removed)
            delete!(distr_in_intval,removed)
        else
            # Note(Yi-Chun): Interval labels: x,y,(removed),z,w
            # Note(Yi-Chun): Merge the distribution
            y = adj_for_intval[removed][2]
            Merge_distr = distr_in_intval[y] + distr_in_intval[removed]
            # Note(Yi-Chun): Update the Delta value for two adjacents
            x = adj_for_intval[removed][1]
            pq[y] = merge_adj_intval(distr_in_intval[x],Merge_distr)
            z= adj_for_intval[removed][3]
            pq[z] = merge_adj_intval(Merge_distr,distr_in_intval[z])
            # Note(Yi-Chun): Update the structure
            distr_in_intval[y] = Merge_distr
            if length(adj_for_intval[removed]) == 3
                adj_for_intval[y][3] = z
                pop!(adj_for_intval[y])
                adj_for_intval[x][end] = z
                adj_for_intval[z][1] = x
                adj_for_intval[z][2] = y
            else
                w = adj_for_intval[removed][4]
                adj_for_intval[x][end] = z
                adj_for_intval[y][end-1] = z
                adj_for_intval[y][end] = w
                adj_for_intval[z][1] = x
                adj_for_intval[z][2] = y
                adj_for_intval[w][1] = y
            end
            # Note(Yi-Chun): Delete the information of removed one
            delete!(adj_for_intval,removed)
            delete!(distr_in_intval,removed)
        end
        n_intval = n_intval - 1

        if (MODL < best_so_far_MODL) && (n_intval<length_of_Greedy_merge+5)
            best_so_far_MODL = MODL
            best_so_far = deepcopy(collect(keys(distr_in_intval)))
            best_distr = deepcopy(distr_in_intval)
        end
    end

    if n_intval > 4
        binedges = sort(collect(keys(distr_in_intval)))
        return binedges
    else
        binedges = sort(collect(keys(distr_in_intval)))
        while (n_intval>2)
            MODL = MODL + DataStructures.peek(pq)[2] + log((n_intval-1)/(n+n_intval-1))
            n_intval = n_intval-1
            removed = DataStructures.dequeue!(pq)
            prior = binedges[something(findfirst(isequal(removed), binedges),0)-1]
            Merge_distr = distr_in_intval[removed] + distr_in_intval[prior]
            distr_in_intval[prior] = Merge_distr
            delete!(adj_for_intval,removed)
            delete!(distr_in_intval,removed)
            binedges = sort(collect(keys(distr_in_intval)))

            for i = 1 : n_intval-1
                X = distr_in_intval[binedges[i]]
                Y = distr_in_intval[binedges[i+1]]
                pq[binedges[i+1]] = merge_adj_intval(X,Y)
            end

            if MODL < best_so_far_MODL
                best_so_far_MODL = MODL
                best_so_far = deepcopy(collect(keys(distr_in_intval)))
                best_distr = deepcopy(distr_in_intval)
            end
        end

        if (DataStructures.peek(pq)[2] + log((n_intval-1)/(n+n_intval-1)))<0
            best_for_far = [1]
            return (best_for_far,one_bin_distrib)
        else
            return (sort(best_so_far),best_distr)
        end
    end
end

function post_greedy_index(continuous,class_value,upper_bound = 0)

    N = length(continuous)
    uniq_class = unique(class_value)
    J = length(uniq_class)

    # Note(Yi-Chun): Check whether we have restriction on bin size
    if upper_bound == 0 # Note(Yi-Chun): the case for no upper bound of bin size
        (bins,distr) = uncondi_greedy_merge_index(continuous,class_value)
        @assert(bins == sort(collect(keys(distr))))
        I = length(bins)
        methods_on_index = Dict()
        for i = 1:I
            methods_on_index[bins[i]] = ["S"]
        end
        append!(bins,[N+1]) # Note(Yi-Chun): for convenience
    else # Note(Yi-Chun): the case for having upper bound of bin size
        I = upper_bound-1
        reminder = N%I
        new_edge = 1
        bins = [1]
        for i = 1 : I-1
            new_edge += div(N,I)
            if reminder > 0
                new_edge += 1
                reminder -= 1
            end
            append!(bins,[new_edge])
        end
        append!(bins,[N+1])
        distr = Dict()
        for bin_index = 1 : I
            distr[bins[bin_index]] = Vector{Int64}(undef,J)
            for j_index = 1 : J
                n_J = count(x->x==uniq_class[j_index],class_value[bins[bin_index]:bins[bin_index+1]-1])
                distr[bins[bin_index]][j_index] = n_J
            end
        end
        methods_on_index = Dict()
        for i = 1:I
            methods_on_index[bins[i]] = ["S"]
        end
    end

        # Note(Yi-Chun): Add methods on each bin into PriorityQueue
        post_pq = DataStructures.PriorityQueue()
        # Note(Yi-Chun): keys in post_pq are represented as: (x,y), where x is the index of bin, y are methods.
        # Note(Yi-Chun): If "S"->split, "M"->merge, "MS"->merge_split, "MMS"->merge_merge_split

        effects = Dict()
        # Note(Yi-Chun): This dictionary record the effects of the methods

        # Note(Yi-Chun): Here we initialize the methods in PriorityQueue
        # Note(Yi-Chun): Methods on bin[1],... ,bin[I-2]. Each bin has 4 methods
        for i = 1 : I-2
            # Note(Yi-Chun): splitting in bin[i]
            method_result = methods_split(continuous,class_value,distr,bins[i],bins[i+1]-1,uniq_class,I)
            post_pq[(bins[i],"S")] = method_result[1]
            effects[(bins[i],"S")] = method_result[2:end]
            # Note(Yi-Chun): merge bin[i] with bin[i+1]
            method_result = methods_merge(distr[bins[i]],distr[bins[2]],N,I)
            post_pq[(bins[i],"M")] = method_result
            append!(methods_on_index[bins[i]],["M"])
            # Note(Yi-Chun): merge_split on bin[i]
            method_result = methods_merge_split(continuous,class_value,distr,
                                                bins[i],bins[i+1],bins[i+2],uniq_class,I)
            post_pq[(bins[i],"MS")] = method_result[1]
            effects[(bins[i],"MS")] = method_result[2:end]
            append!(methods_on_index[bins[i]],["MS"])
            # merge_merge_split on bin[i]
            method_result = methods_merge2_split(continuous,class_value,distr,
                                                 bins[i],bins[i+1],bins[i+2],bins[i+3],uniq_class,I)
            post_pq[(bins[i],"MMS")] = method_result[1]
            effects[(bins[i],"MMS")] = method_result[2:end]
            append!(methods_on_index[bins[i]],["MMS"])
        end

        # Note(Yi-Chun): Methods on bin[I-1]. Only Split, Merge, Merge_Split.
        if I>1
            method_result = methods_split(continuous,class_value,distr,bins[I-1],bins[I]-1,uniq_class,I)
            post_pq[(bins[I-1],"S")] = method_result[1]
            effects[(bins[I-1],"S")] = method_result[2:end]
            method_result = methods_merge(distr[bins[I-1]],distr[bins[I]],N,I)
            post_pq[(bins[I-1],"M")] = method_result
            append!(methods_on_index[bins[I-1]],["M"])
            method_result = methods_merge_split(continuous,class_value,distr,bins[I-1],bins[I],N+1,uniq_class,I)
            post_pq[(bins[I-1],"MS")] = method_result[1]
            effects[(bins[I-1],"MS")] = method_result[2:end]
            append!(methods_on_index[bins[I-1]],["MS"])
        end

        # Note(Yi-Chun): Method on bin[I]. Only Split.
        method_result = methods_split(continuous,class_value,distr,bins[I],N+1,uniq_class,I)
        post_pq[(bins[I],"S")] = method_result[1]
        effects[(bins[I],"S")] = method_result[2:end]

        # Note(Yi-Chun): Here we start the dequeue process
        current_I = I
        time_dequeue = 0
        while (DataStructures.peek(post_pq)[2] < -0.00001 )
            MODL = DataStructures.peek(post_pq)[2]
            time_dequeue += 1

            (removed_index,removed_method) = DataStructures.dequeue!(post_pq)
            removed_posit_bins = searchsorted(bins,removed_index)[1]
            # Note(Yi-Chun): Ittells the position of removed_index in binedges

            # Note(Yi-Chun): This part updates the binedges and distribution of each bin after the dequeue
            if removed_method == "M"
                merged_index = bins[removed_posit_bins+1]
                distr[removed_index] = distr[removed_index] + distr[merged_index]
                delete!(distr,merged_index)
                remove_methods_on_index(merged_index,methods_on_index,post_pq,effects)
                current_I -= 1
                splice!(bins,removed_posit_bins+1)
                methods_on_index[merged_index] = ["S"]

            elseif removed_method == "S"
                effect = effects[(removed_index,removed_method)]
                delete!(effects,(removed_index,removed_method))
                distr[removed_index] = effect[2]
                distr[effect[1]] = effect[3]
                methods_on_index[effect[1]] = ["S"]
                current_I += 1
                insert_and_dedup!(bins, effect[1])

            elseif removed_method == "MS"
                delete!(distr,bins[removed_posit_bins+1])
                effect = effects[(removed_index,removed_method)]
                delete!(effects,(removed_index,removed_method))
                distr[removed_index] = effect[2]
                distr[effect[1]] = effect[3]
                methods_on_index[effect[1]] = ["S"]
                remove_methods_on_index(bins[removed_posit_bins+1],methods_on_index,post_pq,effects)
                methods_on_index[bins[removed_posit_bins+1]] = ["S"]
                splice!(bins,removed_posit_bins+1)
                insert_and_dedup!(bins, effect[1])

            else # Note(Yi-Chun): removed_method == "MMS"
                delete!(distr,bins[removed_posit_bins+1])
                delete!(distr,bins[removed_posit_bins+2])
                remove_methods_on_index(bins[removed_posit_bins+1],methods_on_index,post_pq,effects)
                remove_methods_on_index(bins[removed_posit_bins+2],methods_on_index,post_pq,effects)
                effect = effects[(removed_index,removed_method)]
                delete!(effects,(removed_index,removed_method))
                distr[removed_index] = effect[2]
                distr[effect[1]] = effect[3]
                methods_on_index[effect[1]] = ["S"]
                current_I -= 1
                methods_on_index[bins[removed_posit_bins+1]] = ["S"]
                methods_on_index[bins[removed_posit_bins+2]] = ["S"]
                splice!(bins,removed_posit_bins+1)
                splice!(bins,removed_posit_bins+1)
                insert_and_dedup!(bins, effect[1])
            end

            # Note(Yi-Chun): This part updates the priorityqueue and the effect of each method
            for be_modified = (removed_posit_bins - 2) : (removed_posit_bins + 1)
                 # Note(Yi-Chun): Update the method "Split"
                if (be_modified > 0)
                    PQ_methods_updata(continuous,class_value,uniq_class,bins,distr,
                                        post_pq,effects,be_modified,"S")
                    methods_on_index[bins[be_modified]] = ["S"]
                end
                # Note(Yi-Chun): Update the method "Merge"
                if (be_modified > 0) && (be_modified+1 < current_I+1)
                    PQ_methods_updata(continuous,class_value,uniq_class,bins,distr,
                                        post_pq,effects,be_modified,"M")

                    append!(methods_on_index[bins[be_modified]],["M"])
                end
                # Note(Yi-Chun): Update the method "MergeSplit"
                if (be_modified > 0) && (be_modified+1 < current_I+1)
                    PQ_methods_updata(continuous,class_value,uniq_class,bins,distr,
                                        post_pq,effects,be_modified,"MS")
                    append!(methods_on_index[bins[be_modified]],["MS"])
                end
                # Note(Yi-Chun): Update the method "MergeMergeSplit"
                if (be_modified > 0) && (be_modified+2 < current_I+1)
                    PQ_methods_updata(continuous,class_value,uniq_class,bins,distr,
                                        post_pq,effects,be_modified,"MMS")
                    append!(methods_on_index[bins[be_modified]],["MMS"])
                end
            end

            if upper_bound != 0
                if current_I >= upper_bound
                    for i = 1 : current_I
                        post_pq[(bins[i],"S")] = Inf
                    end
                end
            end
        end
        return bins
end

function post_greedy_result(continuous,class_value,upper_bound=0)

    binedges = post_greedy_index(continuous,class_value,upper_bound)
    bins = Vector{Float64}(undef,length(binedges))
    for i = 1 : length(binedges)
        if i==1
            bins[1] = continuous[1]
        elseif i == length(binedges)
            bins[i] = continuous[end]
        else
            bins[i] = 0.5*(continuous[binedges[i]-1] + continuous[binedges[i]])
        end
    end
    return bins
end

ext: jl


text:
using IntervalUnions
using Test
using Logging

# Testing tuples
gl = global_logger()
global_logger(ConsoleLogger(gl.stream, Logging.Error))

@testset "OrderedPair Tests" begin
    @info "Check operations on ordered pairs"
    include("orderedpairs.jl")
end

@testset "Simple Intervals Tests" begin
    @info "Check operations on simple intervals"
    include("intervals.jl")
end

@testset "IntervalUnionss Tests" begin
    @info "Check operations on unions of intervals"
    include("intervalunions.jl")
end


text:
using StatsBase

"""
    smd_out = isolatesmd(smd::SMD, indices::Union{Vector, UnitRange})

Isolate the requested localizations from `smd`.

# Description.
Make a copy of `smd` and isolate only those localizations corresponding to the
indices.
"""
function isolatesmd(smd::SMD, indices::Union{Vector, UnitRange})
    # Determine the length of vector fields and ensure consistency with the
    # provided indices.
    arraylength = length(smd.framenum)
    if maximum(indices) > arraylength
        error(string("Input `indices` contains values greater than the ",
            "length of `smd.framenum`"))
    end

    # Loop through the fields in `smd` and isolate the appropriate indices,
    # leaving fields with a different length than `arraylength` unchanged.
    smd_out = deepcopy(smd)
    for field in fieldnames(SMD)
        currentfield = getfield(smd_out, field)
        if length(currentfield) == arraylength
            setfield!(smd_out, field, currentfield[indices])
        end
    end

    return smd_out
end

"""
    smd = isolatesmd(smd::SMD, keepbool::Vector{Bool})

Isolate the requested localizations from `smd`.

# Description.
Make a copy of `smd` and isolate only those localizations corresponding to the
true values of `keepbool`.
"""
function isolatesmd(smd::SMD, keepbool::AbstractArray{Bool})
    # Determine the length of vector fields and ensure consistency with the
    # provided `keepbool`.
    arraylength = length(smd.framenum)
    if length(keepbool) > arraylength
        error(string("Input `keepbool` is longer than the ",
            "length of `smd.framenum`"))
    end

    # Loop through the fields in `smd` and isolate the appropriate keepbool,
    # leaving fields with a different length than `arraylength` unchanged.
    smd_out = deepcopy(smd)
    for field in fieldnames(SMD)
        currentfield = getfield(smd_out, field)
        if length(currentfield) == arraylength
            setfield!(smd_out, field, currentfield[keepbool])
        end
    end

    return smd_out
end


text:
module GridWorld

include("environment.jl")
include("distributions.jl")
include("interaction.jl")
include("solver.jl")

end # module

max_stars_repo_name: Joshua-Robison/GridWorld

max_issues_repo_name: Joshua-Robison/GridWorld

max_forks_repo_name: Joshua-Robison/GridWorld


text:
module TestConstructors
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType

@testset "Type parameter constraints" begin
    # cannot use categorical value as level type
    @test_throws TypeError CategoricalPool{CategoricalValue{Int,UInt8}, UInt8, CategoricalValue{CategoricalValue{Int,UInt8},UInt8}}(
            Dict{CategoricalValue{Int,UInt8}, UInt8}(), false)
    @test_throws TypeError CategoricalPool{CategoricalValue{Int,UInt8}, UInt8, CategoricalValue{CategoricalValue{Int,UInt8},UInt8}}(
                CategoricalValue{Int,UInt8}[], false)
    # cannot use non-categorical value as categorical value type
    @test_throws ArgumentError CategoricalPool{Int, UInt8, Int}(Int[], false)
    @test_throws ArgumentError CategoricalPool{Int, UInt8, Int}(Dict{Int, UInt8}(), false)
    # level type of the pool and categorical value must match
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{String, UInt8}}(Int[], false)
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{String, UInt8}}(Dict{Int, UInt8}(), false)
    # reference type of the pool and categorical value must match
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt16}}(Int[], false)
    @test_throws ArgumentError CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt16}}(Dict{Int, UInt8}(), false)
    # correct types combination
    @test CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt8}}(Int[], false) isa CategoricalPool
    @test CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt8}}(Dict{Int, UInt8}(), false) isa CategoricalPool
end

@testset "empty CategoricalPool{String}" begin
    pool = CategoricalPool{String}()

    @test isa(pool, CategoricalPool{String})

    @test isa(pool.levels, Vector{String})
    @test length(pool.levels) == 0

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 0
end

@testset "empty CategoricalPool{Int}" begin
    pool = CategoricalPool{Int, UInt8}()

    @test isa(pool, CategoricalPool{Int, UInt8, CategoricalValue{Int, UInt8}})

    @test isa(pool.levels, Vector{Int})
    @test length(pool.levels) == 0

    @test isa(pool.invindex, Dict{Int, UInt8})
    @test length(pool.invindex) == 0
end

@testset "CategoricalPool{String, DefaultRefType}(a b c)" begin
    pool = CategoricalPool(["a", "b", "c"])

    @test isa(pool, CategoricalPool{String, UInt32, CategoricalValue{String, UInt32}})

    @test isa(pool.levels, Vector{String})
    @test pool.levels == ["a", "b", "c"]

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)
end

@testset "CategoricalPool{String, UInt8}(a b c)" begin
    pool = CategoricalPool{String, UInt8}(["a", "b", "c"])

    @test isa(pool, CategoricalPool)

    @test isa(pool.levels, Vector{String})
    @test pool.levels == ["a", "b", "c"]

    @test isa(pool.invindex, Dict{String, UInt8})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === UInt8(1)
    @test pool.invindex["b"] === UInt8(2)
    @test pool.invindex["c"] === UInt8(3)
end

@testset "CategoricalPool(a b c) with invindex" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.levels, Vector{String})
    @test pool.levels == ["a", "b", "c"]

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with invindex" begin
    pool = CategoricalPool(
        Dict(
            "a" => 1,
            "b" => 2,
            "c" => 3,
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.levels, Vector{String})
    @test pool.levels == ["a", "b", "c"]

    @test isa(pool.invindex, Dict{String, Int})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === 1
    @test pool.invindex["b"] === 2
    @test pool.invindex["c"] === 3
end

@testset "CategoricalPool(c b a)" begin
    pool = CategoricalPool(["c", "b", "a"])

    @test isa(pool, CategoricalPool)

    @test pool.levels == ["c", "b", "a"]

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)
end

@testset "CategoricalPool(a b c) with ref codes not matching the natural order" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(3),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(1),
        )
    )

    @test isa(pool, CategoricalPool)

    @test pool.levels == ["c", "b", "a"]

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)
end

@testset "CategoricalPool{Float64, UInt8}()" begin
    pool = CategoricalPool{Float64, UInt8}([1.0, 2.0, 3.0])

    @test isa(pool, CategoricalPool{Float64, UInt8, CategoricalValue{Float64, UInt8}})
    @test CategoricalValue(1, pool) isa CategoricalValue{Float64, UInt8}
end

@testset "Invalid arguments" begin
    @test_throws ArgumentError CategoricalPool(Dict("a" => 1, "b" => 3))
    @test_throws ArgumentError CategoricalPool(["a", "a"])
end

end


text:
@doc raw"""
    Entropy <: AbstractImageBinarizationAlgorithm
    Entropy()

    binarize([T,] img, f::Entropy)
    binarize!([out,] img, f::Entropy)

An algorithm for finding the binarization threshold value using
the entropy of the image histogram.

# Output

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

# Details

This algorithm uses the entropy of a one-dimensional histogram to produce a threshold
value.

Let ``f_1, f_2, \ldots, f_I`` be the frequencies in the various bins of the
histogram and ``I`` the number of bins. With ``N = \sum_{i=1}^{I}f_i``, let
``p_i = \frac{f_i}{N}`` (``i = 1, \ldots, I``) denote the probability
distribution of gray levels. From this distribution one derives two additional
distributions. The first defined for discrete values ``1`` to ``s`` and the
other, from ``s+1`` to ``I``. These distributions are

```math
A: \frac{p_1}{P_s}, \frac{p_2}{P_s}, \ldots, \frac{p_s}{P_s}
\quad \text{and} \quad
B: \frac{p_{s+1}}{1-P_s}, \ldots, \frac{p_n}{1-P_s}
\quad \text{where} \quad
P_s = \sum_{i=1}^{s}p_i.
```

The entropies associated with each distribution are as follows:

```math
H(A) = \ln(P_s) + \frac{H_s}{P_s}
```
```math
H(B) = \ln(1-P_s) + \frac{H_n-H_s}{1-P_s}
```
```math
\quad \text{where} \quad
H_s = -\sum_{i=1}^{s}p_i\ln{p_i}
\quad \text{and} \quad
H_n = -\sum_{i=1}^{I}p_i\ln{p_i}.
```

Combining these two entropy functions we have

```math
\psi(s) = \ln(P_s(1-P_s)) + \frac{H_s}{P_s} + \frac{H_n-H_s}{1-P_s}.
```
Finding the discrete value ``s`` which maximises the function ``\psi(s)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

See Section 4 of [1] for more details on the derivation of the entropy.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

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

# Example

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

```julia
using TestImages, ImageBinarization

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

# References
1. J. N. Kapur, P. K. Sahoo, and A. K. C. Wong, “A new method for gray-level picture thresholding using the entropy of the histogram,” *Computer Vision, Graphics, and Image Processing*, vol. 29, no. 1, p. 140, Jan. 1985.[doi:10.1016/s0734-189x(85)90156-2](https://doi.org/10.1016/s0734-189x%2885%2990156-2)
"""
struct Entropy <: AbstractImageBinarizationAlgorithm end

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

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


text:
__precompile__()

# HamNNNNclassifier utility functions

"""
Reads a dataset file and massages it to return a datatable suitable for training or testing,
and a vector containing values of the class attribute.
minInteger is a parameter that determines the cutoff for considering an integer attribute
as being a categorical (ie discrete) attribute. If the number of uniques for an attribute
is equal to or greater than minInteger; treat as a continuous attribute; else discrete.
"""
function generatetraintesttable(dataFile)
	codes, names, table = readdatafile(dataFile)
	i, classAttributeValues = extractclassattribute(table, codes)
	# @show codes
	for i in 1:length(names)
		# take out null values to get type purity for type dispatching
		d = dropnull(table[i])
		# for missing codes, fill in "C" or "D"
		typeCode = generateattributecode(d, codes[i])
		# convert "C" to "D" for integer attributes in range 0 to 9
		if typeCode == "C"
			typeCode = changeattributetype(d)
		end
		# add the new codes to the existing codes, in a newCodes array
		codes[i] = typeCode
	end
	delete!(table, unusedcolumns(codes))
	# @show table unusedcodes(codes) classAttributeValues
	return table, unusedcodes(codes), classAttributeValues
end

"""
Takes a vector of orange format codes and returns a vector of indices for delete codes.
"""
unusedcolumns(codes, deletes=["c", "i", "m", "w"]) = [i for i in eachindex(codes) if codes[i] in deletes]
unusedcodes(codes, deletes=["c", "i", "m", "w"]) = [codes[i] for i in eachindex(codes) if ∉(codes[i], deletes)]

function extractclassattribute(table, codes)
	# returns an index for the class attribute within the data table, and
	# a vector with the values for the class attribute
	for i in eachindex(codes)
		if codes[i] == "c"
			return i, table[i]
		end
	end
end

"""
Round to the number of decimal places desired.
"""
roundpercentage(a::Float64, decimals::Integer) = a < 10.0 ? signif(a,decimals+1) : signif(a,decimals+2)

"""
The following two function methods are used to differentiate numeric (ie continuous)
attributes from non-numeric, ie strings (discrete).
"""
numericattribute{T<:Real}(a::AbstractArray{T}) = true
numericattribute(a::AbstractArray) = false

"""
The first method returns true for integer attributes; the second returns false for
non-integer attributes
"""
integerattribute{T<:Integer}(a::AbstractArray{T}) = true
integerattribute(a::AbstractArray) = false

"""
the first method, for numeric attributes, returns the minimum, maximum, mean, and median.
the second method, for string (ie discrete) attributes, returns the number of unique values
and an array with the unique values.
"""
describecontinuousattribute{T<:Real}(a::AbstractArray{T}) = return extrema(a), mean(a), median(a)

function describediscreteattribute(a::AbstractArray)
	u = unique(a)
	return length(u), u
end

"""
Provide codes for continuous or discrete attributes ("C" or "D" respectively) if no code
was in the data file, based on whether the attribute values are numbers or not.
"""
function generateattributecode(vector, code)
	if code != ""
		return code
	else
		# test for whether attribute values are from type Real
		return numericattribute(vector) ? "C" : "D"
	end
end

"""
Change a continuous attribute to discrete, if the attribute is an integer and is a single digit.
"""
function changeattributetype(vector)
	min, max = extrema(vector)
	if integerattribute(vector) && min >= 0 && max < 10
		return "D"
	else
		return "C"
	end
end

"""
Return an array of the attributes sorted into rank order for relative ability to separate the classes.
"""
function generateattributeranking(dataFile)
	cs(x::NullableArrays.NullableArray{Int64,1}) = convert(Array, x, 0)
	cs(x::NullableArrays.NullableArray{Float64,1}) = convert(Array, x, 0.)
	cs(x::NullableArrays.NullableArray{String,1}) = convert(Array, x, "")
	table, codes, classVals = generatetraintesttable(dataFile)
	classValues = cs(classVals)
	classes = unique(classValues)
	labels = names(table)
	A = falses(length(codes),length(classValues),length(classes),length(classValues))
	B = zeros(Int32,length(codes),length(classValues),length(classes))
	for i in indices(A,1) # attributes
		uniques = unique(dropnull(table[i]))
		# @show uniques
		# @show table[i]
		values = cs(table[i])
		# @show i labels[i] values uniques
		for j in eachindex(uniques) # unique values for each attribute
			# @show j uniques[j]
			for k in indices(A,3) # classes
				# @show k classes[k]
				for m in indices(A,4) # cases
					# @show m classValues[m] values[m]
					# set an element for a particular class k and unique value j of an attribute i
					# to true if the class for its case m is the same as class k
					# @show classes[k] classValues[m] values[m] uniques[j]
					A[i,j,k,m] = classes[k] == classValues[m] && values[m] == uniques[j]
					# B[i,j,k] = sum(classes[k] == classValues[m] && values[m] == uniques[j])

				end
				# @show A[i,j,k,:] sum(A[i,j,k,:])
				# sum along the m dimension; put the result in B
				B[i,j,k] = sum(A[i,j,k,:])
			end
		end
	end
	# @show size(A) A
	# @show size(B,1)
	# for each attribute, get the absolute values of the differences between each class
	C = Array{Int64,2}(length(codes), 2)
	for i in indices(B,1) # attributes
		rank = 0
		for j in indices(B,2) # cases
			for k in indices(B,3) # classes
				for m in k:size(B,3) # other classes
					# @show i j k m
					rank += abs(B[i,j,k] - B[i,j,m])
				end
			end
		end
		# @show i labels[i] rank
		C[i,1] = rank
		C[i,2] = i
	end
	D = sortrows(C,rev=true)
	# @show labels table D
	return labels, codes, table, D
end


text:
printstyled("  \"splat\" structure creation\n", color=:light_green)
U = randSeisData()
C = randSeisChannel()
W = randSeisEvent()
TD = convert(EventTraceData, randSeisData())
EC = convert(EventChannel, randSeisChannel())

S = SeisData( U, C, EC, TD, W)
T = EventTraceData( U, C, EC, TD, W)

for f in datafields
  @test getfield(S, f) == getfield(T, f)
end


text:
@testset "Partitions" begin
    trials = 100

    for m = 1:100, K = 1:4, _ = 1:trials
        Φ1 = MapPartition(K, rand(1:K, m))
        Φ2 = convert(DomainPartition, Φ1)
        @test convert(MapPartition, Φ2) == Φ1

        Φ1 = DomainPartition(K, randperm(m), [1, sort(rand(1:(m + 1), K - 1))..., m + 1])
        for k = 1:K
            sort!(@view Φ1.prm[Φ1.spl[k] : Φ1.spl[k + 1] - 1])
        end
        Φ2 = convert(MapPartition, Φ1)
        @test convert(DomainPartition, Φ2) == Φ1
    end
end


text:
struct Normal <: Distribution{Float64} end
struct BroadcastedNormal <: Distribution{Array{Float64}} end

function broadcast_shapes_or_crash(args...)
    # This expression produces a DimensionMismatch exception if the shapes are
    # incompatible.  There seems to currently be no ready-made exception-free
    # way to check this, and I think it makes more sense to wait until that
    # capability exists in Julia proper than to re-implement it here.
    s = Base.Broadcast.broadcast_shape(map(size, args)...)
end

function assert_has_shape(x, expected_shape; msg="Shape assertion failed")
    if size(x) != expected_shape
        throw(DimensionMismatch(string(msg,
                                       " Expected shape: $expected_shape",
                                       " Actual shape: $(size(x))")))
    end
    nothing
end

"""
    normal(mu::Real, std::Real)

Samples a `Float64` value from a normal distribution.
"""
const normal = Normal()

"""
    broadcasted_normal(mu::AbstractArray{<:Real, N1},
                       std::AbstractArray{<:Real, N2}) where {N1, N2}

Samples an `Array{Float64, max(N1, N2)}` of shape
`Broadcast.broadcast_shapes(size(mu), size(std))` where each element is
independently normally distributed.  This is equivalent to (a reshape of) a
multivariate normal with diagonal covariance matrix, but its implementation is
more efficient than that of the more general [`mvnormal`](@ref) for this case.

The shapes of `mu` and `std` must be broadcast-compatible.

If all args are 0-dimensional arrays, then sampling via
`broadcasted_normal(...)` returns a `Float64` rather than properly returning an
`Array{Float64, 0}`.  This is consistent with Julia's own inconsistency on the
matter:

```jldoctest
julia> typeof(ones())
Array{Float64,0}

julia> typeof(ones() .* ones())
Float64
```
"""
const broadcasted_normal = BroadcastedNormal()

function logpdf(::Normal, x::Real, mu::Real, std::Real)
    z = (x - mu) / std
    - (abs2(z) + log(2π))/2 - log(std)
end

function logpdf(::BroadcastedNormal,
                x::Union{AbstractArray{<:Real}, Real},
                mu::Union{AbstractArray{<:Real}, Real},
                std::Union{AbstractArray{<:Real}, Real})
    assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                     msg="Shape of `x` does not agree with the sample space")
    z = (x .- mu) ./ std
    sum(- (abs2.(z) .+ log(2π)) / 2 .- log.(std))
end

function logpdf_grad(::Normal, x::Real, mu::Real, std::Real)
    z = (x - mu) / std
    deriv_x = - z / std
    deriv_mu = -deriv_x
    deriv_std = -1. / std + abs2(z) / std
    (deriv_x, deriv_mu, deriv_std)
end

function logpdf_grad(::BroadcastedNormal,
                     x::Union{AbstractArray{<:Real}, Real},
                     mu::Union{AbstractArray{<:Real}, Real},
                     std::Union{AbstractArray{<:Real}, Real})
    assert_has_shape(x, broadcast_shapes_or_crash(mu, std);
                     msg="Shape of `x` does not agree with the sample space")
    z = (x .- mu) ./ std
    deriv_x = - z ./ std
    deriv_mu = -deriv_x
    deriv_std = -1. ./ std .+ abs2.(z) ./ std
    (_unbroadcast_like(x, deriv_x),
     _unbroadcast_like(mu, deriv_mu),
     _unbroadcast_like(std, deriv_std))
end

_unbroadcast_like(::Real, full_arr) = sum(full_arr)
_unbroadcast_like(::AbstractArray{<:Real, 0}, full_arr::Real) = fill(full_arr)
function _unbroadcast_like(a::AbstractArray{<:Real, N},
                           full_arr::AbstractArray{T}
                          )::AbstractArray{T, N} where {N,T}
    if size(a) == size(full_arr)
        return full_arr
    end
    return _unbroadcast_to_shape(size(a), full_arr)
end

"""
"Unbroadcasts" `full_arr` to have shape `target_shape` by:

* Summing over all dims that would be increased by a broadcast from shape
  `target_shape` to shape `size(full_arr)`
* Then dropping trailing dims (which will all be 1's) as needed so that the
  result has shape `target_shape`.

Requires that `size(full_arr)` is "strictly bigger" than `target_shape`, in the
sense that

    Broadcast.broadcast_shapes(target_shape, size(full_arr)) == size(full_arr)
"""
function _unbroadcast_to_shape(target_shape::NTuple{target_ndims, Int},
                               full_arr::AbstractArray{T, full_ndims}
                         ) where {T, target_ndims, full_ndims}
    @assert full_ndims >= target_ndims
    should_sum_dim(i) = (i > target_ndims) || (target_shape[i] == 1 &&
                                               size(full_arr, i) > 1)
    dropdims(sum(full_arr; dims=filter(should_sum_dim, 1:full_ndims));
             dims=Dims(target_ndims + 1 : full_ndims))
end

random(::Normal, mu::Real, std::Real) = mu + std * randn()
is_discrete(::Normal) = false

function random(::BroadcastedNormal,
                mu::Union{AbstractArray{<:Real}, Real},
                std::Union{AbstractArray{<:Real}, Real})
    broadcast_shape = broadcast_shapes_or_crash(mu, std)
    mu .+ std .* randn(broadcast_shape)
end

(::Normal)(mu, std) = random(Normal(), mu, std)
(::BroadcastedNormal)(mu, std) = random(BroadcastedNormal(), mu, std)

has_output_grad(::Normal) = true
has_argument_grads(::Normal) = (true, true)

has_output_grad(::BroadcastedNormal) = true
has_argument_grads(::BroadcastedNormal) = (true, true)

export normal
export broadcasted_normal


text:
#Calculating the denominator in Racah Formula

function fung(tt,j1,j2,j3,J1,J2,J3)

  return r = factorial(tt-j1-j2-j3)*factorial(tt-j1-J2-J3)*factorial(tt-J1-j2-J3)*factorial(tt-J1-J2-j3)*factorial(j1+j2+J1+J2-tt)*
   factorial(j2+j3+J2+J3-tt)*factorial(j3+j1+J3+J1-tt)
 end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
import CSV: Rows
import Statistics: mean
using Plots

"""
    write_circ(file_name,ps)

Takes PorousStructure `ps` and writes `circ` file
(text file) with name `file_name`
"""
function write_circ(file_name::String,ps::PorousStructure)
    filepath = string(file_name)
    open(filepath, "w") do io
        println(io, "# nbods")
        println(io, ps.param.n_objects)
        println(io, "# data below: radius, xc, yc for each body.")
        for ti=1:ps.param.n_objects
            println(io, ps.olist[ti].radius)
            println(io, ps.olist[ti].center.x)
            println(io, ps.olist[ti].center.y)
        end
    end
end

"""
    csv_to_polygon(file_name)

Reads in csv file containing (x,y) coordinates of points. Returns
arrays of `Point` type.
"""
function csv_to_polygon(file_name::String)
    n_lines = 0
    for row in Rows(file_name;datarow=1)
        n_lines += 1
    end

    xArr = zeros(n_lines)
    yArr = zeros(n_lines)

    ind = 1
    for row in Rows(file_name;datarow=1)
        xArr[ind] = parse(Float64,row[1])
        yArr[ind] = parse(Float64,row[2])
        ind += 1
    end

    rescale_polygon!(xArr,yArr)
    center_polygon!(xArr,yArr)

    pList = [Point(xArr[ti],yArr[ti]) for ti=1:n_lines]

    return pList
end

"""
    center_polygon!(xArr,yArr)

Takes in list of points as two arrays, and centers them. Especially
useful for taking in points obtains from, for example, Inkscape which
often is translated to random points in the plane
"""
function center_polygon!(xArr::Vector{T},yArr::Vector{T}) where T<:Real
    xCenter = 0.5*(maximum(xArr) + minimum(xArr))
    yCenter = 0.5*(maximum(yArr) + minimum(yArr))
    for ti=1:length(xArr)
        xArr[ti] -= xCenter
        yArr[ti] -= yCenter
    end
    nothing
end

"""
    rescale_polygon!(xArr,yArr)

Takes in list of points as two arrays, and rescales them so that the
minimum axis (either horizontal or vertical) has length 1. Notably,
this function keeps the same aspect ratio.
"""
function rescale_polygon!(xArr::Vector{T},yArr::Vector{T}) where T<:Real
    xlength = maximum(xArr) - minimum(xArr)
    ylength = maximum(yArr) - minimum(yArr)
    len = minimum([xlength,ylength])
    for ti=1:length(xArr)
        xArr[ti] /= len
        yArr[ti] /= len
    end
    nothing
end

function translate!(ps::PorousStructure,xval::T,yval::T) where T<:Real
    for ti=1:length(ps.olist)
        ps.olist[ti].center.x += xval
        ps.olist[ti].center.y += yval
    end
    copyCentersToArrays!(ps)
    nothing
end

function translate!(outline::OutlineCircle,xval::T,yval::T) where T<:Real
    outline.center.x += xval;
    outline.center.y += yval;
    nothing
end
function translate!(outline::OutlineRectangle,xval::T,yval::T) where T<:Real
    outline.p1.x += xval;
    outline.p1.y += yval;
    outline.p2.x += xval;
    outline.p2.y += yval;
    nothing
end
function translate!(outline::OutlinePolygon,xval::T,yval::T) where T<:Real
    for ti=1:length(outline.pList)
        outline.pList[ti].x += xval;
        outline.pList[ti].y += yval;
    end
    nothing
end

"""
    save_image(output_name,ps;fill=:blue)

Saves image of filled in objects to file
with name output_name (recommended file type: svg)
with option to choose fill color
"""
function save_image(output_name::String,ps::PorousStructure;fill=:blue)
    # initialize figure
    p = plot(color=:black,aspect_ratio=1)

    # plot circles
    for ti=1:ps.param.n_objects
        plot_circle!(p,ps.radiiArr[ti],ps.xArr[ti],ps.yArr[ti])
    end

    savefig(p,output_name)
end
"""
    save_image(output_name,ps,outline)

Saves image of outline and filled in objects to file
with name output_name (recommended file type: svg)
with option to choose fill color
"""
function save_image(output_name::String,ps::PorousStructure,outline::O;fill=:blue) where O<:AbstractOutlineObject
    # initialize figure
    p = plot(color=:black,aspect_ratio=1)

    # plot circles
    for ti=1:ps.param.n_objects
        plot_circle!(p,ps.radiiArr[ti],ps.xArr[ti],ps.yArr[ti])
    end

    # plot outline
    plot_outline!(p,outline)

    savefig(p,output_name)
end
"""
    save_image(output_name,circ_file,outline)

Saves image of outline and filled in objects from
circ file to file with name output_name (recommended file
type: svg) with option to choose fill color
"""
function save_image(output_name::String,circ_file::String,outline::O;fill=:blue) where O<:AbstractOutlineObject
    # load in circ object
    radiusArr, xArr, yArr = read_in_circ(circ_file)

    # initialize figure
    p = plot(color=:black,aspect_ratio=1)

    # plot circles
    n_objects = length(radiusArr)
    for ti=1:n_objects
        plot_circle!(p,radiusArr[ti],xArr[ti],yArr[ti];fill=fill)
    end
    plot_outline!(p,outline)

    savefig(p,output_name)
end

function save_image(output_name::String,psArr::Vector{PorousStructure{T}},cArr;alpha=0.2) where T<:Real
    # initialize figure
    p = plot(color=:black,aspect_ratio=1)

    # plot circles
    for tk=1:length(psArr)
        for ti=1:psArr[tk].param.n_objects
            plot_circle!(
                p,psArr[tk].radiiArr[ti],psArr[tk].xArr[ti],psArr[tk].yArr[ti];fill = cArr[tk],alpha=alpha)
        end
    end

    savefig(p,output_name)
end

"""
    read_in_circ(file_input)

Takes circ file with name `file_input` and returns arrays
of radii, x, and y coordinates.
"""
function read_in_circ(file_input::String)
    # read in circ file
    circfile = string(file_input)
    f         = open(circfile,"r")
    lines     = readlines(f)
    num_lines = length(lines)
    close(f)

    # obtain pre-processing parameters (# bodies, radii, x/y centers)
    nbodies = parse(Int,lines[2])
    radiusArr = zeros(nbodies)
    xArr = zeros(nbodies)
    yArr = zeros(nbodies)

    tind = 1
    for ti = 4:3:num_lines-1
        radiusArr[tind] = parse(Float64,lines[ti])
        xArr[tind] = parse(Float64,lines[ti+1])
        yArr[tind] = parse(Float64,lines[ti+2])
        tind += 1
    end
    return radiusArr, xArr, yArr
end

"""
    plot_circle!(p,radius,x,y;color=:blue)

Adds a circle to plot object `p`. Circle has center (x,y). Has option
to select fill-in color
"""
function plot_circle!(p::Plots.Plot{Plots.GRBackend},radius::T,x::T,y::T;fill=:blue,alpha=0.2) where T<:Real
    nθ = 40
    θarr = 0:2*pi/(nθ-1):2*pi

    plot!(p,
        x .+ radius*cos.(θarr),y .+ radius*sin.(θarr),
        seriestype=[:shape,],lw=0.0,c=fill,
        linecolor=:black,legend=false,
        fillalpha=alpha)

    nothing
end

"""
    plot_outline!(p,outline)

Adds a outline to plot object `p`. Contains different methods
for different outline types (circle, rectangle, polygon).
"""
function plot_outline!(p::Plots.Plot{Plots.GRBackend},outline::OutlineCircle{T}) where T<:Real
    nθ = 100
    θarr = 0:2*pi/(nθ-1):2*pi

    plot!(p,
        outline.center.x .+ outline.radius*cos.(θarr),
        outline.center.y .+ outline.radius*sin.(θarr),
        label="",color=:black,linestyle=:dash)

    nothing
end
function plot_outline!(p::Plots.Plot{Plots.GRBackend},outline::OutlineRectangle{T}) where T<:Real
    nn = 3

    # plot top
    x = range(outline.p1.x,outline.p2.x,length=nn)
    y = outline.p2.y*ones(nn)
    plot!(p,x,y,label="",color=:black,linestyle=:dash)
    # plot bottom
    y = outline.p1.y*ones(nn)
    plot!(p,x,y,label="",color=:black,linestyle=:dash)
    # plot left
    x = outline.p1.x*ones(nn)
    y = range(outline.p1.y,outline.p2.y,length=nn)
    plot!(p,x,y,label="",color=:black,linestyle=:dash)
    # plot right
    x = outline.p2.x*ones(nn)
    plot!(p,x,y,label="",color=:black,linestyle=:dash)

    nothing
end
function plot_outline!(p::Plots.Plot{Plots.GRBackend},outline::OutlinePolygon{T}) where T<:Real
    x = zeros(2)
    y = zeros(2)

    # plot first nn lines
    nn = length(outline.pList)
    for ti=1:(nn-1)
        x[1] = outline.pList[ti].x; x[2] = outline.pList[ti+1].x
        y[1] = outline.pList[ti].y; y[2] = outline.pList[ti+1].y
        plot!(p,x,y,label="",color=:black,linestyle=:dash)
    end
    # plot last one
    x[1] = outline.pList[nn].x; x[2] = outline.pList[1].x
    y[1] = outline.pList[nn].y; y[2] = outline.pList[1].y
    plot!(p,x,y,label="",color=:black,linestyle=:dash)

    nothing
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
using AbstractLogic
using Test
using Suppressor

@test (@capture_out abstractlogic("clearall")) == "Clearing Everything!\n"

# Basic logicparser features
@suppress abstractlogic("a = 1 [clear]"); @test replerror

@suppress abstractlogic("range abc"); @test replerror
@suppress abstractlogic("clear;; a,b,c in 1:2;; range a b"); @test !replerror
@suppress abstractlogic("a,b,c in 1:2;; range a [clear]"); @test !replerror

@suppress abstractlogic("check a==b"); @test !replerror
@suppress abstractlogic("a in 1; check z=1[clear]"); @test replerror

@suppress abstractlogic("clear; a,b,c in 0:1; clear"); @test nfeasible(replset) != 1
@suppress abstractlogic("clearall"); @test sessionhistory.current == 1
@suppress abstractlogic("a,b,c ∈ 0:1; a=b [clear]"); @test activehistory.current == 2

# Navigation
@suppress abstractlogic("back [clear]"); @test replerror
@suppress abstractlogic("next [clear]"); @test replerror
@suppress abstractlogic("a,b,c ∈ 0:1;; back"); @test !replerror
@suppress abstractlogic("next"); @test !replerror
@test @suppress abstractlogic("", returnactive = true) |> nfeasible == 8

# Show
@suppress abstractlogic("show"); @test !replerror
@suppress abstractlogic("clear;; show"); @test replerror

# Import/Export
@suppress abstractlogic("clear;; a,b in 1:2;; export testset"); @test testset |> nfeasible == 4
@test @suppress abstractlogic("clear;; import testset", returnactive = true) |> nfeasible == 4

# REPL Commands
@suppress abstractlogic("compare noname"); @test replerror
@suppress abstractlogic("keys"); @test !replerror
@suppress abstractlogic("d"); @test !replerror
@suppress abstractlogic("discover"); @test !replerror

x = @capture_out(abstractlogic("clear;; a,b in 1:2;; compare testset"))[(end-30):end]
  @test x == "lues. replset[:] == testset[:]\n"
x = @capture_out(abstractlogic("a = 1;; compare testset"))[(end-40):end]
  @test x == " nfeasible(replset) < nfeasible(testset)\n"
x = @capture_out(abstractlogic("d, e in 1:5;; compare testset"))[(end-40):end]
  @test x == " nfeasible(replset) > nfeasible(testset)\n"
x = @capture_out(abstractlogic("clear;; a ∈ 1:2;; export t2;; clear;; b ∈ 1:3; b>1 ;; compare t2 [silent]"))[(end-40):end]
  @test x == "ues. nfeasible(replset) == nfeasible(t2)\n"

@suppress abstractlogic("clear; a in 1; b = 1"); @test replerror

@suppress abstractlogic("clear;; a,b,c in 1:2;; search {{i}}>{{!i}}"); @test !replerror
@suppress abstractlogic("clear;; a,b,c in 1:2;; search z>z"); @test replerror

@suppress abstractlogic("h"); @test !replerror
@suppress abstractlogic("history"); @test !replerror

@suppress abstractlogic("clear ;; a in 1;; check aa = 1"); @test replerror

printcleaner(x) = replace(x, r"( |\n|–|\"|\t|Feasible|Perceived|Outcomes)"=>"")

@test (@capture_out abstractlogic("H")) == (@capture_out abstractlogic("History"))
@test (@capture_out abstractlogic("clearall;; H")) |> printcleaner ==
  "ClearingEverything!Command#feasible#SessionStarted0<<present>>..."

# Silence and noisy
@test @suppress abstractlogic("clear;; a,b in 1:2;; preserve;; clear;; restore", returnactive = true) |> nfeasible == 4
@suppress abstractlogic("silence"); @test !replerror
@suppress abstractlogic("noisy"); @test !replerror
@suppress abstractlogic("dash;; dashboard"); @test !replerror

@suppress abstractlogic("clear;; a,b,c in 1")
output = (@capture_out abstractlogic("history")) |> printcleaner
@test  output == "Command#feasible#SessionStarted0a,b,c∈11<<present>>..."
@test @capture_out(abstractlogic("keys")) == "a, b, c\n"

output = (@capture_out abstractlogic("b")) |> printcleaner
@test output == "#SessionStarted-:0"

@suppress abstractlogic("n")
output = (@capture_out abstractlogic("a=b")) |> printcleaner
@test output == "a=b:1:1✓✓:111"

@suppress abstractlogic("silence")
@test (@capture_out abstractlogic("a ∈ 1,2 [clear]")) == ""

@suppress abstractlogic("noisy")
@test (@capture_out abstractlogic("a ∈ 1,2 [clear]")) != ""

@test (@capture_out abstractlogic("clear;; a,b ∈ 1,2; a=b [silent]")) == ""
@capture_out(AbstractLogic.parse_to_expr("clear")) == "Clearing Activeset\n"

# Running all examples: examples 4 and 8 have been simplified to improve run speed
@test length(@capture_out abstractlogic("t(:)")) > 1000
@suppress abstractlogic("t(undefinedtest)"); @test replerror

# @suppress abstractlogic("t(1)")  ; @test replset |> nfeasible == 1
# @suppress abstractlogic("t(2)")  ; @test replset |> nfeasible == 6
# @suppress abstractlogic("t(3)")  ; @test replset |> nfeasible == 3
# @suppress abstractlogic("t(4)")  ; @test replset |> nfeasible == 6
# # @suppress abstractlogic("t(5)")  ; @test replset |> nfeasible == 12


text:
# Tests for Cost Hölder DEA Models
@testset "CostHölderDEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Cost Hölder L1
    costholderl1 = deacostholder(X, Y, W, l = 1)

    @test typeof(costholderl1) == CostHolderDEAModel

    @test nobs(costholderl1) == 8
    @test ninputs(costholderl1) == 2
    @test noutputs(costholderl1) == 1
    @test ismonetary(costholderl1) == false

    @test efficiency(costholderl1) == efficiency(costholderl1, :Economic)
    @test efficiency(costholderl1, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 6.0; 3.0; 3.0; 5.6]
    @test efficiency(costholderl1, :Technical)  ≈ [0.0; 0.0; 0.0; 2.0; 4.0; 0.0; 1.0; 0.6]
    @test efficiency(costholderl1, :Allocative) ≈ [0.0; 1.0; 1.0; 1.0; 2.0; 3.0; 2.0; 5.0]
    @test normfactor(costholderl1) ≈ [1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]

    # With monetary = true
    costholderl1mom = deacostholder(X, Y, W, l = 1, monetary = true)

    @test ismonetary(costholderl1mom) == true

    @test efficiency(costholderl1mom) ≈ efficiency(costholderl1) .* normfactor(costholderl1)
    @test efficiency(costholderl1mom, :Technical) ≈ efficiency(costholderl1, :Technical) .* normfactor(costholderl1)
    @test efficiency(costholderl1mom, :Allocative) ≈ efficiency(costholderl1, :Allocative) .* normfactor(costholderl1)

    # With CRS (Same result in this example)
    costholderl1crs = deacostholder(X, Y, W, l = 1, rts = :CRS)

    @test efficiency(costholderl1crs, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 6.0; 3.0; 3.0; 5.6]
    @test efficiency(costholderl1crs, :Technical)  ≈ [0.0; 0.0; 0.0; 2.0; 4.0; 0.0; 1.0; 0.6]
    @test efficiency(costholderl1crs, :Allocative) ≈ [0.0; 1.0; 1.0; 1.0; 2.0; 3.0; 2.0; 5.0]

    # Cost Hölder L1 Weighted (weakly)
    costholderl1weight= deacostholder(X, Y, W, l = 1, weight = true)

    @test efficiency(costholderl1weight, :Economic)   ≈ [0.0; 0.25; 0.25; 0.75; 1.2; 0.5; 0.6; 0.7]
    @test efficiency(costholderl1weight, :Technical)  ≈ [0.0; 0.0; 0.0; 0.625; 0.8; 0.0; 0.5; 0.375]
    @test efficiency(costholderl1weight, :Allocative) ≈ [0.0; 0.25; 0.25; 0.125; 0.4; 0.5; 0.1; 0.325]

    # Cost Hölder L2
    logs, value = Test.collect_test_logs() do
        costholderl2 = deacostholder(X, Y, W, l = 2, optimizer = DEAOptimizer(:NLP))

        @test normfactor(costholderl2) ≈ sqrt(2) .* ones(8)
        show(IOBuffer(), costholderl2)
    end

    # Cost Hölder L2 Weighted (weakly)
    logs, value = Test.collect_test_logs() do
        costholderl2weight = deacostholder(X, Y, W, l = 2, weight = true, optimizer = DEAOptimizer(:NLP))

        @test normfactor(costholderl2weight) ≈ sqrt.(sum((W .* X).^2, dims = 2))
        show(IOBuffer(), costholderl2weight)
    end

    # Cost Hölder LInf
    costholderlInf = deacostholder(X, Y, W, l = Inf)

    @test efficiency(costholderlInf, :Economic)   ≈ [0.0; 0.5; 0.5; 1.5; 3.0; 1.5; 1.5; 2.8]
    @test efficiency(costholderlInf, :Technical)  ≈ [0.0; 0.0; 0.0; 4/3; 3.0; 0.0; 1.0; 0.6]
    @test efficiency(costholderlInf, :Allocative) ≈ [0.0; 0.5; 0.5; 1/6; 0.0; 1.5; 0.5; 2.2]
    @test normfactor(costholderlInf) ≈ [2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0]

    # Cost Hölder LInf Weighted (weakly)
    costholderlInfweight = deacostholder(X, Y, W, l = Inf, weight = true)

    @test efficiency(costholderlInfweight, :Economic)   ≈ [0.0; 0.2; 0.2; 3/7; 0.6; 3/7; 3/7; 7/12]
    @test efficiency(costholderlInfweight, :Technical)  ≈ [0.0; 0.0; 0.0; 0.4; 0.6; 0.0; 1/3; 3/8]
    @test efficiency(costholderlInfweight, :Allocative) ≈ [0.0; 0.2; 0.2; 1/35; 0.0; 3/7; 2/21; 5/24]
    @test normfactor(costholderlInfweight) ≈ [4.0; 5.0; 5.0; 7.0; 10.0; 7.0; 7.0; 9.6]

    # Print
    show(IOBuffer(), costholderl1)
    show(IOBuffer(), costholderl1weight)
    show(IOBuffer(), costholderlInf)
    show(IOBuffer(), costholderlInfweight)

    # Test errors
    @test_throws DimensionMismatch deacostholder([1; 2 ; 3], [4 ; 5], [1; 1; 1], l = 1) #  Different number of observations
    @test_throws DimensionMismatch deacostholder([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], l = 1) # Different number of observation in prices
    @test_throws DimensionMismatch deacostholder([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], l = 1) # Different number of input prices and inputs
    @test_throws ArgumentError deacostholder([1; 2; 3], [4; 5; 6], [1; 2; 3], l = 5) # Invalid l
    @test_throws ArgumentError deacostholder([1; 2; 3], [4; 5; 6], [1; 2; 3], l = 1, rts = :Error) # Invalid returns to scale

end


text:
abstract type AbstractSoil end

@units @description mutable struct BaseSoil <: AbstractSoil
    pct_sand::Float64        | "percent sand"
    pct_clay::Float64        | "percent clay"
    s_sat::Float64           | "relative saturation"
    N_cell::Array{Float64,1} | "number of cells in a microsite"
end

function clapp_hornberger_par(pct_sand::Float64, pct_clay::Float64)
    chb = 2.91 + 0.195*pct_clay
    sat = 0.489 - 0.00126*pct_sand
    psisat = -10.0^(1.88-0.0131*pct_sand)
    ksat = 0.0070556*10.0^(-0.884+0.0153*pct_sand)
    return chb, sat, psisat, ksat
end

function moldrup_tau(por, epsi, theta, kappa)
    taug = epsi*(epsi/por)^(3.0/kappa)
    tauw = theta*((theta+1e-20)/por)^(kappa/3.0-1.0)
    return taug, tauw
end

function cosby_psi(s_sat, psisat, sat, chb)
    psi=max(psisat*(s_sat+1.e-10)^(-chb)*1.e-3*1.e4, -1e8*1e-3*1e4)
    dpsidvsm = -chb*psi/(s_sat*sat+1e-20)
    return psi, dpsidvsm
end

function cosby_hk(s_sat, ksat, chb)
    hk = ksat*(s_sat)^(2.0*chb+3.0)
end

function cosby_Dwpsi(s_sat, psisat, ksat, sat, chb)
    psi, dpsidvsm = cosby_psi(s_sat, psisat, sat, chb)
    hk = cosby_hk(s_sat, ksat, chb)
    Dwpsi = hk*1e-3*dpsidvsm*1e-4
end

function soil_affinity_properties(pct_sand::Float64, pct_clay::Float64, s_sat::Float64)
    chb, sat, psisat, ksat = clapp_hornberger_par(pct_sand, pct_clay)
    theta = s_sat*sat # phi_w
    epsi = sat-theta
    psi, dpsidvsm = cosby_psi(s_sat, psisat, sat, chb)
    film = max(exp(-13.65-0.857*log(-psi*1e-6)), 1e-7)
    τ_g, τ_w = moldrup_tau(sat, epsi, theta, chb)
    return theta, τ_g, τ_w, film
end

function aqueous_diffusivity(molecular_weight::Array{Float64,1})
    # at 298 K
    D_sucrose = 0.5216
    M_sucrose = 342.0
    D_S = @. D_sucrose*sqrt(M_sucrose/molecular_weight)
    return D_S.*1e-9
end

function effective_diffusivity(molecular_weight::Array{Float64,1}, pct_sand::Float64, pct_clay::Float64, s_sat::Float64)
    # at 298 K
    D_sucrose = 0.5216
    M_sucrose = 342.0
    D_S = @. D_sucrose*sqrt(M_sucrose/molecular_weight)
    ϕ_w, τ_g, τ_w, δ  = soil_affinity_properties(pct_sand, pct_clay, s_sat)
    D_eff = D_S*τ_w*ϕ_w
    return D_eff.*1e-9
end


text:
@info "Pandas.jl is available - including additional file formats using Python Pandas"

import .Pandas
import DataFrames

# HDF5 formats

function read_table(::TableIOInterface.HDF5Format, filename, key; kwargs...):: Pandas.DataFrame
    df_pandas = Pandas.read_hdf(filename, key; kwargs)
    return df_pandas
end

function read_table(t::TableIOInterface.HDF5Format, filename; kwargs...):: Pandas.DataFrame
    hdf = Pandas.HDFStore(filename)
    try
        table_list = _list_tables(t, hdf)
        length(table_list) > 1 && @warn "File contains more than one table, the alphabetically first one is taken"
        key = first(table_list)
        return Pandas.read_hdf(filename, key; kwargs)
    finally
        close(hdf)
    end
end

function write_table!(::TableIOInterface.HDF5Format, filename, tablename, table:: Pandas.DataFrame; kwargs...)
    # for HDF5 files, the table name may contain a path, therefore _checktablename is not executed.
    # Pandas.DataFrame has no Tables.jl compatible interface, therefore _checktable is not executed.
    Pandas.to_hdf(table, filename, tablename; kwargs...)
    nothing
end

write_table!(::TableIOInterface.HDF5Format, filename, tablename, table; kwargs...) = write_table!(TableIOInterface.HDF5Format(), filename, tablename, Pandas.DataFrame(table); kwargs...)

function list_tables(t::TableIOInterface.HDF5Format, filename:: AbstractString)
    hdf = Pandas.HDFStore(filename)
    try
        _list_tables(t, hdf)
    finally
        close(hdf)
    end
end

_list_tables(::TableIOInterface.HDF5Format, hdf) = keys(hdf) |> sort


text:
module BipedControlUtil

export Side,
    left,
    right,
    flipsign_if_right

include("side.jl")

end # module


text:
function sanitycheck(app_module)
    try getfield(app_module, :ui)
    catch
        rethrow(ErrorException("You have not defined the `ui` function in your app module"))
    end
    try getfield(app_module, :Server)
    catch
        rethrow(ErrorException("You have not defined the `Server` module in your app module"))
    end
    try getfield(app_module, :title)
    catch
        rethrow(ErrorException("You have not set a `title` in your app module"))
    end
end


text:
function slice(array :: AbstractArray, periodic :: Bool, iterators...)
    indices = zip(iterators...)
    array = CircularArray(array)

    if periodic
        # Since periodic diagonals work only on cubic arrays you can
        # call size(a, i) with any i to get side of a cube.
        stop_iter = take(indices, size(array, 1))
    else
        stop_iter = takewhile(pair -> checkbounds(Bool, array, pair...), indices)
    end

    return [@inbounds array[indices...] for indices in stop_iter]
end

function diagonals(array     :: AbstractArray{T,2},
                   periodic  :: Bool,
                   direction :: Tuple{Int, Int}) where T
    h, w = size(array)
    δx, δy = direction
    sx, sy = map((dir, dim) -> (dir == 1) ? 1 : dim, direction, (h, w))
    diagonal(x, y) = slice(array, periodic, countfrom(x, δx), countfrom(y, δy))

    part1 = (diagonal(x,sy) for x in 1:h)
    # Start from 2:w to not get the longest slice twice
    part2 = (diagonal(sx,y) for y in 2:w)

    return periodic ? part1 : flatten((part1, part2))
end

# Code for 3D diagonals is really hard for understanding.
# TODO: Try to "fuse" these 4 methods into more simple one (if
# possible at all).
function diagonals(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:xyz}) where T
    # Direction (1, 1, 1)
    h, w, d = size(array)
    diagonal(x, y, z) = slice(array, periodic,
                              countfrom(x, 1),
                              countfrom(y, 1),
                              countfrom(z, 1))
    part1 = (diagonal(1,y,z) for y in 1:w for z in 1:d)
    part2 = (diagonal(x,1,z) for x in 2:h for z in 1:d)
    part3 = (diagonal(x,y,1) for x in 2:h for y in 2:w)
    return periodic ? part1 : flatten((part1, part2, part3))
end

function diagonals(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:yxz}) where T
    # Direction (-1, 1, 1)
    h, w, d = size(array)
    diagonal(x, y, z) = slice(array, periodic,
                              countfrom(x, -1),
                              countfrom(y,  1),
                              countfrom(z,  1))
    part1 = (diagonal(h,y,z) for y in 1:w   for z in 1:d)
    part2 = (diagonal(x,1,z) for x in 1:h-1 for z in 1:d)
    part3 = (diagonal(x,y,1) for x in 1:h-1 for y in 2:w)
    return periodic ? part1 : flatten((part1, part2, part3))
end

function diagonals(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:xzy}) where T
    # Direction (1, -1, 1)
    h, w, d = size(array)
    diagonal(x, y, z) = slice(array, periodic,
                              countfrom(x,  1),
                              countfrom(y, -1),
                              countfrom(z,  1))
    part1 = (diagonal(1,y,z) for y in 1:w for z in 1:d)
    part2 = (diagonal(x,w,z) for x in 2:h for z in 1:d)
    part3 = (diagonal(x,y,1) for x in 2:h for y in 1:w-1)
    return periodic ? part1 : flatten((part1, part2, part3))
end

function diagonals(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:zyx}) where T
    # Direction (1, 1, -1)
    h, w, d = size(array)
    diagonal(x, y, z) = slice(array, periodic,
                              countfrom(x,  1),
                              countfrom(y,  1),
                              countfrom(z, -1))
    part1 = (diagonal(1,y,z) for y in 1:w for z in 1:d)
    part2 = (diagonal(x,1,z) for x in 2:h for z in 1:d)
    part3 = (diagonal(x,y,d) for x in 2:h for y in 2:w)
    return periodic ? part1 : flatten((part1, part2, part3))
end

# Slicers for "true" diagonals (when there are no zeros in the
# direction vector).
TrueDiagonal = Union{Val{:xyz}, Val{:yxz}, Val{:xzy}, Val{:zyx}}
slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, dir :: TrueDiagonal) where T =
    diagonals(array, periodic, dir)

# Slicers for other directions (3D)
slice_generators(array :: AbstractArray{T,3}, :: Bool, :: Val{:x}) where T =
    (array[:,j,k] for j in 1:size(array, 2) for k in 1:size(array, 3))

slice_generators(array :: AbstractArray{T,3}, :: Bool, :: Val{:y}) where T =
    (array[i,:,k] for i in 1:size(array, 1) for k in 1:size(array, 3))

slice_generators(array :: AbstractArray{T,3}, :: Bool, :: Val{:z}) where T =
    (array[i,j,:] for i in 1:size(array, 1) for j in 1:size(array, 2))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:xy}) where T =
    flatten(diagonals(array[:,:,k], periodic, (1, 1)) for k in 1:size(array, 3))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:xz}) where T =
    flatten(diagonals(array[:,j,:], periodic, (1, 1)) for j in 1:size(array, 2))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:yz}) where T =
    flatten(diagonals(array[i,:,:], periodic, (1, 1)) for i in 1:size(array, 1))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:yx}) where T =
    flatten(diagonals(array[:,:,k], periodic, (-1, 1)) for k in 1:size(array, 3))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:zx}) where T =
    flatten(diagonals(array[:,j,:], periodic, (-1, 1)) for j in 1:size(array, 2))

slice_generators(array :: AbstractArray{T,3}, periodic :: Bool, :: Val{:zy}) where T =
    flatten(diagonals(array[i,:,:], periodic, (-1, 1)) for i in 1:size(array, 1))

# Slicers for different directions (2D)
slice_generators(array :: AbstractArray{T,2}, :: Bool, :: Val{:x}) where T =
    (array[:,j] for j in 1:size(array, 2))

slice_generators(array :: AbstractArray{T,2}, :: Bool, :: Val{:y}) where T =
    (array[i,:] for i in 1:size(array, 1))

slice_generators(array :: AbstractArray{T,2}, periodic :: Bool, :: Val{:xy}) where T =
    diagonals(array, periodic, (1, 1))

slice_generators(array :: AbstractArray{T,2}, periodic :: Bool, :: Val{:yx}) where T =
    diagonals(array, periodic, (-1, 1))

# Trivial slicer for 1D case
slice_generators(array :: AbstractArray{T,1}, :: Bool, :: Val{:x}) where T =
    # Ugly
    (array for x in 0:0)