text:
export closure_type_counts3,
    closure_type_counts4,
    open_types3,
    open_types4,
    STRONG,
    WEAK,
    OPEN,
    simplex_key3,
    simplex_key4,
    newly_closed_types3,
    newly_closed_types4,
    initialize_type_counter3,
    initialize_type_counter4,
    triangle_counts

const STRONG = 2
const WEAK   = 1
const OPEN   = 0

simplex_key3(wij::Int64, wik::Int64, wjk::Int64) =
    sorted_tuple(min(wij, STRONG), min(wik, STRONG), min(wjk, STRONG))

function simplex_key4(wijk::Int64, wijl::Int64, wikl::Int64, wjkl::Int64,
                      num_edges::Int64)
    wijk = min(wijk, STRONG)
    wijl = min(wijl, STRONG)
    wikl = min(wikl, STRONG)
    wjkl = min(wjkl, STRONG)
    key = sort([wijk, wijl, wikl, wjkl], alg=InsertionSort)
    if     num_edges == 5; for j = 1:2; key[j] = -1; end
    elseif num_edges == 4; for j = 1:3; key[j] = -2; end
    elseif num_edges == 3; for j = 1:3; key[j] = -3; end
    end
    return tuple(key...)
end

function initialize_type_counter3()
    type_counts = counter(NTuple{3, Int64})
    for i = 0:2, j = i:2, k = j:2
        push!(type_counts, (i, j, k), 0)
    end
    return type_counts
end

function initialize_type_counter4()
    type_counts = counter(NTuple{4, Int64})
    for w1 = 0:2, w2 = w1:2, w3 = w2:2, w4 = w3:2
        push!(type_counts, (w1, w2, w3, w4), 0)  # All 6 edges
        push!(type_counts, (-1, -1, w3, w4), 0)  # Only 5 edges
        push!(type_counts, (-2, -2, -2, w4), 0)  # Only 4 edges
        push!(type_counts, (-3, -3, -3, w4), 0)  # Only 3 edges
    end
    return type_counts
end

function triangle_counts(simplices::Vector{Int64}, nverts::Vector{Int64})
    counts = counter(NTuple{3, Int64})
    curr_ind = 1
    for nvert in nverts
        simp = simplices[curr_ind:(curr_ind + nvert - 1)]
        curr_ind += nvert
        for (i, j, k) in combinations(simp, 3)
            if i == j || i == k || j == k; continue; end
            push!(counts, sorted_tuple(i, j, k))
        end
    end
    return counts
end

"""
triangle_counts
---------------

Counts the number of times each simplicial triangle appears in the data. This
function is expensive and requires storage on the order of the number of
triangles. Returns a counter that takes a sorted triple of indices as a key and
returns the count of that triangle.

triangle_counts(data::HONData)

Input parameters:
- data::HONData: The data.
"""
triangle_counts(data::HONData) =
    triangle_counts(data.simplices, data.nverts)

# Determine the type of every newly closed 4-node simplex.
function newly_closed_types4(old_simplices::Vector{Int64}, old_nverts::Vector{Int64},
                             new_simplices::Vector{Int64}, new_nverts::Vector{Int64})
    A_old, A_old_t, B_old = basic_matrices(old_simplices, old_nverts)
    simplex_order_old = simplex_degree_order(A_old_t)
    closed(i::Int64, j::Int64, k::Int64, l::Int64) =
        tetrahedron_closed(A_old, A_old_t, simplex_order_old, i, j, k, l)
    n = size(B_old, 1)
    is_new_node = vec(sum(A_old_t, dims=1)) .== 0
    degs = vec(sum(make_sparse_ones(B_old), dims=1))

    new_tetrahedra = Set{NTuple{4, Int64}}()
    type_counts = initialize_type_counter4()

    counts3 = triangle_counts(old_simplices, old_nverts)
    weight3(i::Int64, j::Int64, k::Int64) = counts3[sorted_tuple(i, j, k)]

    curr_ind = 1
    for nvert in new_nverts
        simp = new_simplices[curr_ind:(curr_ind + nvert - 1)]
        curr_ind += nvert
        for (i, j, k, l) in combinations(simp, 4)
            # SKIP if all 4 nodes are not unique
            if i == j || i == k || i == l || j == k || j == l || k == l; continue; end

            # SKIP if a node is new
            if max(i, j, k, l) > n; continue; end
            if any(is_new_node[[i, j, k, l]]); continue; end

            # SKIP if not at least one triangle
            Bo_ij, Bo_ik, Bo_il = B_old[i, j], B_old[i, k], B_old[i, l]
            Bo_jk, Bo_jl, Bo_kl = B_old[j, k], B_old[j, l], B_old[k, l]
            num_edges = sum([Bo_ij, Bo_ik, Bo_il, Bo_jk, Bo_jl, Bo_kl] .> 0)
            if num_edges < 3; continue; end
            ijk_tri = Bo_ij > 0 && Bo_ik > 0 && Bo_jk > 0
            ijl_tri = Bo_ij > 0 && Bo_il > 0 && Bo_jl > 0
            ikl_tri = Bo_ik > 0 && Bo_il > 0 && Bo_kl > 0
            jkl_tri = Bo_jk > 0 && Bo_jl > 0 && Bo_kl > 0
            if (ijk_tri + ijl_tri + ikl_tri + jkl_tri) < 1; continue; end

            # SKIP if already closed
            if num_edges == 6 && closed(i, j, k, l); continue; end

            # SKIP if we already processed the new closure
            key = sorted_tuple(i, j, k, l)
            if key in new_tetrahedra; continue; end
            push!(new_tetrahedra, key)

            wijk = weight3(i, j, k)
            wijl = weight3(i, j, l)
            wikl = weight3(i, k, l)
            wjkl = weight3(j, k, l)
            push!(type_counts, simplex_key4(wijk, wijl, wikl, wjkl, num_edges))
        end
    end

    return type_counts
end

# Get the number of open 4-node configurations of each type containing at least
# one triangle.
function open_types4(simplices::Vector{Int64}, nverts::Vector{Int64})
    A, At, B = basic_matrices(simplices, nverts)
    simplex_order  = simplex_degree_order(At)
    triangle_order = proj_graph_degree_order(B)
    n = size(B, 1)
    type_counts = initialize_type_counter4()
    counts3 = triangle_counts(simplices, nverts)
    weight3(i::Int64, j::Int64, k::Int64) = counts3[sorted_tuple(i, j, k)]

    # map of key to integer index
    tetra_index_map = Dict{NTuple{4, Int64}, Int64}()
    for w1 = 0:2, w2 = w1:2, w3 = w2:2, w4 = w3:2
        ind = length(tetra_index_map) + 1
        tetra_index_map[(w1, w2, w3, w4)] = ind
    end
    # thread-local counters
    nthreads = Threads.nthreads()
    all_tetra_counts_arr = Vector{Vector{Int64}}(undef, nthreads)
    open_tetra_counts_arr = Vector{Vector{Int64}}(undef, nthreads)
    Threads.@threads for i = 1:nthreads
        all_tetra_counts_arr[i] = zeros(Int64, length(tetra_index_map))
        open_tetra_counts_arr[i] = zeros(Int64, length(tetra_index_map))
    end

    # 1. Get 4-node, 6-edge counts, which must be induced by looping over all
    # 2-skeletons in the data.
    print(stdout, "tetrahedral counts...\n")
    shuffled_inds = collect(1:n)
    shuffle!(shuffled_inds)
    Threads.@threads for ii = 1:n
        i = shuffled_inds[ii]
        weight_arr = zeros(Int64, 4)
        if Threads.threadid() == 1
            print(stdout, "$(ii) of $n \r")
            flush(stdout)
        end
        nbrs = sort(neighbors(B, triangle_order, i), by=v->triangle_order[v])
        nnbr = length(nbrs)
        for jj in 1:nnbr
            j = nbrs[jj]
            for kk in (jj + 1):nnbr
                k = nbrs[kk]
                if B[j, k] == 0; continue; end
                weight_arr[1] = min(weight3(i, j, k), STRONG)
                for ll in (kk + 1):nnbr
                    l = nbrs[ll]
                    if B[j, l] > 0 && B[k, l] > 0
                        weight_arr[2] = min(weight3(i, j, l), STRONG)
                        weight_arr[3] = min(weight3(i, k, l), STRONG)
                        weight_arr[4] = min(weight3(j, k, l), STRONG)
                        key = NTuple{4,Int64}(sort(weight_arr, alg=InsertionSort))
                        index = tetra_index_map[key]
                        tid = Threads.threadid()
                        all_tetra_counts_arr[tid][index] += 1
                        if minimum(weight_arr) == 0 ||
                            !tetrahedron_closed(A, At, simplex_order, i, j, k, l)
                            open_tetra_counts_arr[tid][index] += 1
                        end
                    end
                end
            end
        end
    end

    all_tetra_counts  = sum(all_tetra_counts_arr)
    open_tetra_counts = sum(open_tetra_counts_arr)
    for (key, val) in tetra_index_map
        push!(type_counts, key, open_tetra_counts[val])
    end

    index_map = Dict{NTuple{4, Int64}, Int64}()
    for (ind, (key, _)) in enumerate(type_counts)
        index_map[key] = ind
    end

    # 2. Get 5-edge, 4-node (2 triangle) counts
    #
    # 2A. Form matrices whose (i, j) entry is number of triangles with a
    # particular weight containing i and j
    W_base = triu(make_sparse_ones(B))
    W0_all = Vector{SpIntMat}(undef, nthreads)
    W1_all = Vector{SpIntMat}(undef, nthreads)
    W2_all = Vector{SpIntMat}(undef, nthreads)
    Threads.@threads for i = 1:nthreads
        W0_all[i] = copy(W_base)
        W1_all[i] = copy(W_base)
        W2_all[i] = copy(W_base)
    end

    W_all = [W0_all, W1_all, W2_all]
    Threads.@threads for i = 1:n
        for (j, k) in neighbor_pairs(B, triangle_order, i)
            if B[j, k] > 0
                weight = min(weight3(i, j, k), STRONG)
                a, b, c = sort([i, j, k], alg=InsertionSort)
                tid = Threads.threadid()
                W_all[weight + 1][tid][a, b] += 1
                W_all[weight + 1][tid][a, c] += 1
                W_all[weight + 1][tid][b, c] += 1
            end
        end
    end

    # Every edge starts with a count of 1 to get the nonzero pattern correct.
    # Here, we subtract off that count.
    Threads.@threads for i = 1:nthreads
        W0_all[i].nzval .-= 1
        W1_all[i].nzval .-= 1
        W2_all[i].nzval .-= 1
    end
    W0, W1, W2 = sum(W0_all), sum(W1_all), sum(W2_all)
    # 2B. Get non-induced counts. Each edge contributes to a count for every
    # pair of triangles to which it is adjacent.
    #
    # Same triangle weight
    for (i, W) in [(0, W0), (1, W1), (2, W2)]
        for k in 1:n, val in nz_row_vals(W, k)
            push!(type_counts, (-1, -1, i, i), binomial(val, 2))
        end
    end
    # Different triangle weights
    for (i, j, Wi, Wj) in [(0, 1, W0, W1), (0, 2, W0, W2), (1, 2, W1, W2)]
        for k in 1:n, (val, row) in zip(nz_row_vals(Wi, k), nz_row_inds(Wi, k))
            push!(type_counts, (-1, -1, i, j), val * Wj[row, k])
        end
    end

    for w1 = 0:2, w2 = w1:2, w3 = w2:2, w4 = w3:2
        count = all_tetra_counts[tetra_index_map[(w1, w2, w3, w4)]]
        for (wi, wj) in [(w1, w2), (w1, w3), (w1, w4), (w2, w3), (w2, w4), (w3, w4)]
            push!(type_counts, (-1, -1, wi, wj), -count)
        end
    end

    # 3. Triangles with an adjacent edge
    tri_adj_edge_counts_arr  = zeros(length(index_map), nthreads)
    degs = vec(sum(make_sparse_ones(B), dims=2))
    Threads.@threads for i = 1:n
        for (j, k) in neighbor_pairs(B, triangle_order, i)
            if B[j, k] > 0
                key = (-2, -2, -2, min(weight3(i, j, k), STRONG))
                count = degs[i] + degs[j] + degs[k] - 6
                tid = Threads.threadid()
                tri_adj_edge_counts_arr[index_map[key], tid] += count
            end
        end
    end
    tri_adj_edge_counts = sum(tri_adj_edge_counts_arr, dims=2)
    for (key, index) in index_map
        push!(type_counts, key, tri_adj_edge_counts[index])
    end
    # 3A. Each triangle in every 5-edge pattern contributed 2 too many
    for w1 = 0:2, w2 = w1:2
        count = type_counts[(-1, -1, w1, w2)]
        push!(type_counts, (-2, -2, -2, w1), -2 * count)
        push!(type_counts, (-2, -2, -2, w2), -2 * count)
    end
    # 3B. Each triangle in every tetrahedron contributes 3 too many
    for w1 = 0:2, w2 = w1:2, w3 = w2:2, w4 = w3:2
        count = all_tetra_counts[tetra_index_map[(w1, w2, w3, w4)]]
        for w in [w1, w2, w3, w4]
            push!(type_counts, (-2, -2, -2, w), -3 * count)
        end
    end

    # 4. Triangles with an isolated node
    tri_iso_node_counts_arr  = zeros(length(index_map), nthreads)
    num_verts = sum(sum(At, dims=1) .> 0)
    Threads.@threads for i = 1:n
       for (j, k) in neighbor_pairs(B, triangle_order, i)
           if B[j, k] > 0
               key = (-3, -3, -3, min(weight3(i, j, k), STRONG))
               tid = Threads.threadid()
               tri_iso_node_counts_arr[index_map[key], tid] += (num_verts - 3)
           end
       end
    end
    tri_iso_node_counts = sum(tri_iso_node_counts_arr, dims=2)
    for (key, index) in index_map
        push!(type_counts, key, tri_iso_node_counts[index])
    end

    # Subtract off other induced patterns
    #   --tetrahedra
    for w1 = 0:2, w2 = w1:2, w3 = w2:2, w4 = w3:2
        count = all_tetra_counts[tetra_index_map[(w1, w2, w3, w4)]]
        for w in [w1, w2, w3, w4]
            push!(type_counts, (-3, -3, -3, w), -count)
        end
    end
    #   --5-edge
    for w1 = 0:2, w2 = w1:2
        count = type_counts[(-1, -1, w1, w2)]
        push!(type_counts, (-3, -3, -3, w1), -count)
        push!(type_counts, (-3, -3, -3, w2), -count)
    end
    #   --Triangle plus adjacent edge
    for w = 0:2
        count = type_counts[(-2, -2, -2, w)]
        push!(type_counts, (-3, -3, -3, w), -count)
    end

    return type_counts
end

# Determine the previous topology of newly filled triangles.
function newly_closed_types3(old_simplices::Vector{Int64}, old_nverts::Vector{Int64},
                             new_simplices::Vector{Int64}, new_nverts::Vector{Int64})
    A_old, A_old_t, B_old = basic_matrices(old_simplices, old_nverts)
    simplex_order_old = simplex_degree_order(A_old_t)
    closed(i::Int64, j::Int64, k::Int64) =
        triangle_closed(A_old, A_old_t, simplex_order_old, i, j, k)
    n = size(B_old, 1)
    is_new_node = vec(sum(A_old_t, dims=1)) .== 0

    type_counts = initialize_type_counter3()
    new_triangles = Set{NTuple{3, Int64}}()

    # Loop over simplices for closed triangles
    curr_ind = 1
    for nvert in new_nverts
        simp = new_simplices[curr_ind:(curr_ind + nvert - 1)]
        curr_ind += nvert
        for (i, j, k) in combinations(simp, 3)
            # SKIP if any duplicates
            if i == j || i == k || j == k; continue; end
            # SKIP if 0-simplex is missing
            if max(i, j, k) > n || any(is_new_node[[i, j, k]]); continue; end
            # SKIP if the triangle was closed before
            wij, wik, wjk = B_old[i, j], B_old[i, k], B_old[j, k]
            if min(wij, wik, wjk) > 0 && closed(i, j, k); continue; end
            # SKIP if we have already looked at this triangle
            tri_key = sorted_tuple(i, j, k)
            if tri_key in new_triangles; continue; end
            push!(new_triangles, tri_key)
            # update type counts
            push!(type_counts, simplex_key3(wij, wik, wjk))
        end
    end

    return type_counts
end

# Get the number of open triangles for each type.
function open_types3(simplices::Vector{Int64}, nverts::Vector{Int64})
    A, At, B = basic_matrices(simplices, nverts)
    simplex_order  = simplex_degree_order(At)
    triangle_order = proj_graph_degree_order(B)
    n = size(B, 1)

    # Setup counts data structure
    type_counts = initialize_type_counter3()

    index_map = Dict{NTuple{3, Int64}, Int64}()
    for (ind, (key, _)) in enumerate(type_counts)
        index_map[key] = ind
    end

    # Get triangle counts
    all_tris_counts_arr  = zeros(length(index_map), n)
    open_tris_counts_arr = zeros(length(index_map), n)
    Threads.@threads for i = 1:n
        curr_triangle_order = triangle_order[i]
        nbrs = Int64[]
        weights = Int64[]
        for (nbr, weight) in zip(nz_row_inds(B, i), nz_row_vals(B, i))
            if triangle_order[nbr] > curr_triangle_order
                push!(nbrs, nbr)
                push!(weights, weight)
            end
        end
        nnbr = length(nbrs)
        for jj = 1:nnbr, kk = (jj+1):nnbr
            j, k = nbrs[[jj, kk]]
            wij, wik = weights[[jj, kk]]
            wjk = B[j, k]
            if wjk > 0
                tri_ind = index_map[simplex_key3(wij, wik, wjk)]
                all_tris_counts_arr[tri_ind, i] += 1
                if !triangle_closed(A, At, simplex_order, i, j, k)
                    open_tris_counts_arr[tri_ind, i] += 1
                end
            end
        end
    end
    all_tris_counts  = sum(all_tris_counts_arr, dims=2)
    open_tris_counts = sum(open_tris_counts_arr, dims=2)
    triad_types = [(WEAK, WEAK, WEAK), (WEAK, WEAK, STRONG),
                   (WEAK, STRONG, STRONG), (STRONG, STRONG, STRONG)]
    for key in triad_types
        push!(type_counts, key, open_tris_counts[index_map[key]])
    end

    # strong and weak degree counts
    d1, d2 = Int64[], Int64[]
    for i = 1:n
        weights = nz_row_vals(B, i)
        push!(d1, sum(weights .== 1))
        push!(d2, sum(weights .>  1))
    end

    # wedge counts
    push!(type_counts, (OPEN, WEAK,   WEAK),   sum(binomial.(d1, 2)))
    push!(type_counts, (OPEN, STRONG, STRONG), sum(binomial.(d2, 2)))
    push!(type_counts, (OPEN, WEAK,   STRONG), sum(d1 .* d2))
    for triad_key in triad_types
        count = all_tris_counts[index_map[triad_key]]
        for (w1, w2) in combinations(triad_key, 2)
            push!(type_counts, simplex_key3(OPEN, w1, w2), -count)
        end
    end

    # Edge + isolated node counts
    num_verts = sum(sum(At, dims=1) .> 0)
    push!(type_counts, (OPEN, OPEN, WEAK),   (sum(d1) / 2) * (num_verts - 2))
    push!(type_counts, (OPEN, OPEN, STRONG), (sum(d2) / 2) * (num_verts - 2))
    for triad_key in triad_types
        count = all_tris_counts[index_map[triad_key]]
        for w in triad_key
            push!(type_counts, (OPEN, OPEN, w), -count)
        end
    end
    for w1 = 1:2, w2 = w1:2
        count = type_counts[(OPEN, w1, w2)]
        for w in [w1, w2]
            push!(type_counts, (OPEN, OPEN, w), -count)
        end
    end

    # Empty subgraphs
    push!(type_counts, (OPEN, OPEN, OPEN), binomial(num_verts, 3))
    push!(type_counts, (OPEN, OPEN, OPEN), -sum(all_tris_counts))
    for key in [(OPEN, OPEN, WEAK), (OPEN, OPEN, STRONG),
                (OPEN, WEAK, WEAK), (OPEN, WEAK, STRONG), (OPEN, STRONG, STRONG)]
        push!(type_counts, (OPEN, OPEN, OPEN), -type_counts[key])
    end
    return type_counts
end

"""
closure_type_counts3
--------------------

Computes the closure probabilities of all 3-node configurations. The closure
probability is the fraction of instances of open configurations appearing in the
first 80% of the timestamped simplices that appear in a simplex in the final
20%.

closure_type_counts3(dataset::HONData, initial_cutoff::Int64=100)

Input parameters:
- dataset::String: The dataset name.
- initial_cutoff::Int64=100: Initial cutoff of the simplices. If this is set to less than 100, then the data is first preprocessed to only consider the first initial_cutoff percentage of the data.
"""
function closure_type_counts3(dataset::HONData, initial_cutoff::Int64=100)
    simps = dataset.simplices
    nverts = dataset.nverts
    times = dataset.times
    if initial_cutoff < 100
        cutoff_percentile = convert(Int64, round(percentile(times, initial_cutoff)))
        simps, nverts, times = window_data(minimum(times), cutoff_percentile,
                                           simps, nverts, times)
    end
    old_simps, old_nverts, new_simps, new_nverts =
        split_data(simps, nverts, times, 80, 100)
    closed_type_counts =
        newly_closed_types3(old_simps, old_nverts, new_simps, new_nverts)
    open_type_counts = open_types3(old_simps, old_nverts)
    type_counts = [(k..., cnt, closed_type_counts[k]) for (k, cnt) in open_type_counts]
    writedlm("$(dataset.name)-3-node-closures-$(initial_cutoff).txt", sort(type_counts))
end

"""
closure_type_counts4
--------------------

Computes the closure probabilities of all 4-node configurations that contain at
least one triangle. The closure probability is the fraction of instances of open
configurations appearing in the first 80% of the timestamped simplices that
appear in a simplex in the final 20%.

closure_type_counts4(dataset::HONData, initial_cutoff::Int64=100)

Input parameters:
- dataset::HONData: The dataset.
- initial_cutoff::Int64=100: Initial cutoff of the simplices. If this is set to less than 100, then the data is first preprocessed to only consider the first initial_cutoff percentage of the data.
"""
function closure_type_counts4(dataset::HONData, initial_cutoff::Int64=100)
    simps  = dataset.simplices
    nverts = dataset.nverts
    times  = dataset.times
    if initial_cutoff < 100
        cutoff_percentile = convert(Int64, round(percentile(times, initial_cutoff)))
        simps, nverts, times = window_data(minimum(times), cutoff_percentile,
                                           simps, nverts, times)

    end
    old_simps, old_nverts, new_simps, new_nverts =
        split_data(simps, nverts, times, 80, 100)
    closed_type_counts =
        newly_closed_types4(old_simps, old_nverts, new_simps, new_nverts)
    open_type_counts = open_types4(old_simps, old_nverts)
    type_counts = [(k..., cnt, closed_type_counts[k]) for (k, cnt) in open_type_counts]
    writedlm("$(dataset.name)-4-node-closures-$(initial_cutoff).txt", sort(type_counts))
end


text:
abstract type Parameterization end
struct Density{T} <: Parameterization
	ρ::T
end
density(parm::Density) = parm.ρ
permitivity(parm::Density,ϵmin, ϵmax) = ϵmin + density(parm) * (ϵmax - ϵmin)
apply!(f::DensityFilter,ρ::Density) = apply!(f,ρ)
apply!(p::Projector,ρ::Density)  = apply!(p,ρ)

ext: jl

max_stars_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_issues_repo_name: MKAbdElrahman/TopologyOptimizationTools

max_forks_repo_name: MKAbdElrahman/TopologyOptimizationTools


text:
using Test
using Classes

@class Foo(getter_prefix="") <: Class begin
   foo::Int

   Foo() = Foo(0)

   # Although Foo is immutable, subclasses might not be,
   # so it's still useful to define this method.
   function Foo(self::absclass(Foo))
        self.foo = 0
    end
end

@class mutable Bar <: Foo begin
    bar::Int

    # Mutable classes can use this pattern
    function Bar(self::Union{Nothing, absclass(Bar)}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Bar)(self)
        Bar(self, 0)
    end
end

@class Baz(mutable=true) <: Bar begin
   baz::Int

   function Baz(self::Union{Nothing, absclass(Baz)}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Baz)(self)
        Baz(self, 0)
    end
end

# emitted by class Foo
# @method foo(obj::Foo) = obj.foo

@method function sum(obj::Bar)
    return foo(obj) + get_bar(obj)
end

@test Set(subclasses(Foo)) == Set(Any[Bar, Baz])

x = Foo(1)
y = Bar(10, 11)
z = Baz(100, 101, 102)

@test fieldnames(Foo) == (:foo,)
@test fieldnames(Bar) == (:foo, :bar)
@test fieldnames(Baz) == (:foo, :bar, :baz)

@test foo(x) == 1
@test foo(y) == 10
@test foo(z) == 100

@test sum(y) == 21
@test get_bar(y) == 11
@test get_bar(z) == 101

@test_throws Exception get_bar(x)

# Mutable
set_foo!(z, 1000)
@test foo(z) == 1000

# Immutable
@test_throws Exception x.foo = 1000

# test that where clause is amended properly
@method zzz(obj::Foo, bar::T) where {T} = T

@test zzz(x, :x) == Symbol
@test zzz(y, 10.6) == Float64

# Test the corner case of conflicting mutability declarations.
# The eval(parse(string)) form is necessary; otherwise the error is raised
# when the macro is expanded, before @test_throws can catch it. Better way?
@test_throws(LoadError, eval(Meta.parse("@class mutable XX(mutable=false)")))

#
# Test generation of custom accessors
#
@class Buzz(setter_prefix="SET_", getter_prefix="GET_") begin
  s::Symbol
end

@test hasmethod(SET_s!, (_Buzz_, Symbol))
@test hasmethod(GET_s,  (Buzz,))

# Verify that default method doesn't exist
b = Buzz(:x)
@test_throws UndefVarError get_s(b)

# Test suppressing setters and getter prefix
@class Blub(setters=false, getter_prefix="") begin
    my_integer::Int
    my_float::Float64
end

b = Blub(10, 20.)

@test my_integer(b) == 10

@test_throws UndefVarError set_i!(b, 0)


text:
"""
    @loop begin
        body
    end

Repeatedly evaluate the body of the loop.

# Examples
```jldoctest
julia> i = 1
1

julia> @loop begin
           println(i)
           i += 1
           if i > 4
               break
           end
       end
1
2
3
4
```
"""
macro loop(body)
    quote
        while true
            $body
        end
    end |> esc
end


text:
using InvertedIndices
# needs only Base.to_indices in struct.jl to work,
# plus this to work when used in round brackets:

findindex(not::InvertedIndex, r::AbstractVector) = Base.unalias(r, not)

using IntervalSets

findindex(int::Interval, r::AbstractVector) =
    findall(x -> x in int, r)
findindex(int::Interval, r::AbstractRange{T}) where {T<:Union{Number,Char}} =
    findall(in(int), r)

# Since that is now efficient for ranges, comparisons can go there:

findindex(eq::Base.Fix2{typeof(<=)}, r::AbstractRange{T}) where {T<:Union{Number,Char}} =
    findall(in(Interval(typemin(T), eq.x)), r)

findindex(eq::Base.Fix2{typeof(>=)}, r::AbstractRange{T}) where {T<:Union{Number,Char}} =
    findall(in(Interval(eq.x, typemax(T))), r)

findindex(eq::Base.Fix2{typeof(<)}, r::AbstractRange{T}) where {T<:Union{Number,Char}} =
    findall(in(Interval{:closed, :open}(typemin(T), eq.x)), r)

findindex(eq::Base.Fix2{typeof(>)}, r::AbstractRange{T}) where {T<:Union{Number,Char}} =
    findall(in(Interval{:open, :closed}(eq.x, typemax(T))), r)

"""
    Near(val)
    Interval(lo, hi)

These selectors modify lookup using `axiskeys(A)`:
`B(time = Near(3))` matches one entry with minimum `abs(t-3)` of named dimension `:time`.
`C("cat", Interval(10,20))` matches all entries with `10 <= iter <= 20`).

`Interval` is from IntervalSets.jl, and using that you may also write `lo .. hi`,
as well as `mid ± δ`.

    ==(val)
    <(val)

Any functions can be used similarly, like C(!=("dog"), <=(33)).
They ultimately call `findall(==(val), axiskeys(A,d))`.

Functions of type `Base.Fix2`, and `Selector`s, also allow a dimension
to be chosen by type: `A(<=(3.1))` will work provided that only one of
`map(eltype, axiskeys(A))` matches `typeof(3.1)`.

# Examples
```jldoctest
julia> v = KeyedArray(Symbol.('a':'e'), 10:10:50)
1-dimensional KeyedArray(...) with keys:
↓   5-element StepRange{Int64,...}
And data, 5-element Vector{Symbol}:
 (10)  :a
 (20)  :b
 (30)  :c
 (40)  :d
 (50)  :e

julia> v[Near(33)]
:c

julia> v[==(30)]  # all matching this key
1-dimensional KeyedArray(...) with keys:
↓   1-element Vector{Int64}
And data, 1-element Vector{Symbol}:
 (30)  :c

julia> v[Interval(17, 31)]
1-dimensional KeyedArray(...) with keys:
↓   2-element StepRange{Int64,...}
And data, 2-element Vector{Symbol}:
 (20)  :b
 (30)  :c

julia> m = wrapdims(hcat(v,v), x=nothing, y=[:left, :right]);

julia> m(<(30))  # selects 1st dim by type, and makes a view
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   x ∈ 2-element StepRange{Int64,...}
→   y ∈ 2-element Vector{Symbol}
And data, 2×2 view(::Matrix{Symbol}, 1:2, :) with eltype Symbol:
       (:left)  (:right)
 (10)    :a       :a
 (20)    :b       :b
```
"""
abstract type Selector{T} end

Base.eltype(::Type{<:Selector{T}}) where {T} = T

@doc @doc(Selector)
struct Near{T} <: Selector{T}
    val::T
end

Base.show(io::IO, s::Near) = print(io, "Near(",s.val,")")
Base.show(io::IO, ::MIME"text/plain", s::Near{T}) where {T} =
    print(io, "Near(",s.val,") ::Selector{",T,"}")

findindex(sel::Near, range::AbstractArray) = argmin(map(x -> abs(x-sel.val), range))

function findindex(sel::Near, range::AbstractRange)
    iplus = searchsortedfirst(range, sel.val)
    # "index of the first value in a greater than or equal to x"
    if abs(range[iplus]-sel.val) < abs(range[iplus-1]-sel.val)
        return iplus
    else
        return iplus-1
    end
end

_index_key_doc = """
    Index[i]

This exists to let you mix in square-bracket indexing,
like `A(:b, Near(3.14), Index[4:5], "f")`.
You may also write `Index[end]`, although not yet `Index[end-2]`.

    Key(val)

This exists to perform lookup inside indexing, to allow
`A[Key(:b), Near(3.14), 4:5, Key("f")]`.

Writing `Key(isequal(:b))` is equivalent to just `isequal(:b)`, and will find all matches,
while `Key(:b)` finds only the first (and drops the dimension).

# Examples
```jldoctest
julia> v = KeyedArray(Symbol.('a':'e'), 10:10:50)
1-dimensional KeyedArray(...) with keys:
↓   5-element StepRange{Int64,...}
And data, 5-element Vector{Symbol}:
 (10)  :a
 (20)  :b
 (30)  :c
 (40)  :d
 (50)  :e

julia> v[Key(30)] == v(30) == v(Index[3]) == v[3]
true

julia> v[==(30)] == v(Index[3:3]) == v[3:3]
true
```
"""

@doc _index_key_doc
struct Index{T} <: Selector{T}
    ind::T
end

Base.show(io::IO, s::Index{T}) where {T} = print(io, "Index[",s.ind, "]")

Base.getindex(::Type{Index}, i) = Index(i)

findindex(sel::Index, range::AbstractArray) = sel.ind

struct LastIndex <: Selector{Int} end

Base.lastindex(::Type{Index}) = LastIndex()

Index(::LastIndex) = LastIndex()

findindex(sel::LastIndex, range::AbstractArray) = lastindex(range)

if VERSION >= v"1.4" # same thing for Index[begin]
    struct FirstIndex <: Selector{Int} end
    Base.firstindex(::Type{Index}) = FirstIndex()
    Index(::FirstIndex) = FirstIndex()
    findindex(sel::FirstIndex, range::AbstractArray) = firstindex(range)
end

@doc _index_key_doc
struct Key{T} <: Selector{T}
    val::T
end

Base.show(io::IO, s::Key) = print(io, "Key(",s.val,")")
Base.show(io::IO, ::MIME"text/plain", s::Key{T}) where {T} =
    print(io, "Key(",s.val,") ::Selector{",T,"}")

findindex(sel::Key, range::AbstractArray) = findindex(sel.val, range)

"""
    Base.to_indices(A, axes, inds)
    select_to_indices(A, axes, inds)

This recursively peels off the indices & axes, `select_to_indices` gets called
when the first remaining index is a Selector, Interval, or a Function.
"""
@inline function select_to_indices(A, axes, inds)
    d = ndims(A) - length(axes) + 1 # infer how many have been peeled off?
    i = findindex(first(inds), axiskeys(A, d))
    (i, Base.to_indices(A, Base.tail(axes), Base.tail(inds))...)
end

@inline Base.to_indices(A, ax, inds::Tuple{Selector, Vararg}) = select_to_indices(A, ax, inds)

# For the rest I don't own these types... which then creates ambiguities...
@inline Base.to_indices(A::Union{KeyedArray,NdaKa}, ax, inds::Tuple{IntervalSets.Domain, Vararg}) =
    select_to_indices(A, ax, inds)
@inline Base.to_indices(A::Union{KeyedArray,NdaKa}, ax, inds::Tuple{Function, Vararg}) =
    select_to_indices(A, ax, inds)

using Base: to_indices, tail, _maybetail, uncolon

@inline Base.to_indices(A::Union{KeyedArray,NdaKa}, inds, I::Tuple{Colon, Vararg{Any}}) =
    (uncolon(inds, I), to_indices(A, _maybetail(inds), tail(I))...)


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type

# using Random

"""
    LinSysVecRowOneRandCyclic <: LinSysVecRowSelect

A mutable structure with a field to store a cycling order. Randomly specifies a cycling
order over the equations of a linear system. Once this ordering is specified, the ordering
is kept fixed.

# Fields
- `order::Union{Vector{Int64},Nothing}`

Calling `LinSysVecRowOneRandCyclic()` defaults to setting `order` to `nothing`. The `sample`
function will handle the re-initialization of the fields once the system is provided.
"""
mutable struct LinSysVecRowOneRandCyclic <: LinSysVecRowSelect
    order::Union{Vector{Int64},Nothing}
end
LinSysVecRowOneRandCyclic() = LinSysVecRowOneRandCyclic(nothing)

# Common sample interface for linear systems
function sample(
    type::LinSysVecRowOneRandCyclic,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)
    if iter == 1
        type.order = randperm(length(b))
    end

    eqn_ind = type.order[mod(iter, 1:length(b))]
    return A[eqn_ind, :], b[eqn_ind]
end

#export LinSysVecRowOneRandCyclic


text:
"""
    Op

An MPI reduction operator, for use with [Reduce/Scan collective operations](@ref Reduce/Scan) to wrap binary operators. MPI.jl will perform this conversion automatically.

## Usage

    Op(op, T=Any; iscommutative=false)

Wrap the Julia reduction function `op` for arguments of type `T`. `op` is assumed to be
associative, and if `iscommutative` is true, assumed to be commutative as well.

## See also

- [`Reduce!`](@ref)/[`Reduce`](@ref)
- [`Allreduce!`](@ref)/[`Allreduce`](@ref)
- [`Scan!`](@ref)/[`Scan`](@ref)
- [`Exscan!`](@ref)/[`Exscan`](@ref)
"""
@mpi_handle Op MPI_Op fptr

const OP_NULL = _Op(MPI_OP_NULL, nothing)
const BAND    = _Op(MPI_BAND, nothing)
const BOR     = _Op(MPI_BOR, nothing)
const BXOR    = _Op(MPI_BXOR, nothing)
const LAND    = _Op(MPI_LAND, nothing)
const LOR     = _Op(MPI_LOR, nothing)
const LXOR    = _Op(MPI_LXOR, nothing)
const MAX     = _Op(MPI_MAX, nothing)
const MIN     = _Op(MPI_MIN, nothing)
const PROD    = _Op(MPI_PROD, nothing)
const REPLACE = _Op(MPI_REPLACE, nothing)
const SUM     = _Op(MPI_SUM, nothing)

if @isdefined(MPI_NO_OP)
    const NO_OP   = _Op(MPI_NO_OP, nothing)
end

Op(::typeof(min), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint}} = MIN
Op(::typeof(max), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint}} = MAX
Op(::typeof(+), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint,MPIComplex}} = SUM
Op(::typeof(*), ::Type{T}; iscommutative=true) where {T<:Union{MPIInteger,MPIFloatingPoint,MPIComplex}} = PROD
Op(::typeof(&), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BAND
Op(::typeof(|), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BOR
Op(::typeof(⊻), ::Type{T}; iscommutative=true) where {T<:MPIInteger} = BXOR

function free(op::Op)
    if op.val != OP_NULL.val && !Finalized()
        @mpichk ccall(:MPI_Op_free, Cint, (Ptr{MPI_Op},), op)
    end
    return nothing
end

struct OpWrapper{F,T}
    f::F
end

function (w::OpWrapper{F,T})(_a::Ptr{Cvoid}, _b::Ptr{Cvoid}, _len::Ptr{Cint}, t::Ptr{MPI_Datatype}) where {F,T}
    len = unsafe_load(_len)
    @assert isconcretetype(T)
    a = Ptr{T}(_a)
    b = Ptr{T}(_b)
    for i = 1:len
        unsafe_store!(b, w.f(unsafe_load(a,i), unsafe_load(b,i)), i)
    end
    return nothing
end

function Op(f, T=Any; iscommutative=false)
    if MPI_LIBRARY == MicrosoftMPI && Sys.WORD_SIZE == 32
        error("User-defined reduction operators are not supported on 32-bit Windows.\nSee https://github.com/JuliaParallel/MPI.jl/issues/246 for more details.")
    end
    w = OpWrapper{typeof(f),T}(f)
    fptr = @cfunction($w, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}, Ptr{MPI_Datatype}))

    op = Op(OP_NULL.val, fptr)
    # int MPI_Op_create(MPI_User_function* user_fn, int commute, MPI_Op* op)
    @mpichk ccall(:MPI_Op_create, Cint,
                  (Ptr{Cvoid}, Cint, Ptr{MPI_Op}),
                  fptr, iscommutative, op)

    finalizer(free, op)
    return op
end


text:
"""
    pyimport(m)
    pyimport(m => k)
    pyimport(m => (k1, k2, ...))
    pyimport(m1, m2, ...)

Import a module `m`, or an attribute `k`, or a tuple of attributes.

If several arguments are given, return the results of importing each one in a tuple.
"""
pyimport(m) = pynew(errcheck(@autopy m C.PyImport_Import(getptr(m_))))
pyimport((m,k)::Pair) = (m_=pyimport(m); k_=pygetattr(m_,k); pydel!(m_); k_)
pyimport((m,ks)::Pair{<:Any,<:Tuple}) = (m_=pyimport(m); ks_=map(k->pygetattr(m_,k), ks); pydel!(m_); ks_)
pyimport(m1, m2, ms...) = map(pyimport, (m1, m2, ms...))
export pyimport


text:
"""
Store a diagonally-scaled Toeplitz∘Hankel matrix:

    DL(T∘H)DR

where the Hankel matrix `H` is non-negative definite. This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹.
"""
struct ToeplitzHankelPlan{S}
    T::TriangularToeplitz{S}
    C::Vector{Vector{S}}
    DL::Vector{S}
    DR::Vector{S}
    ToeplitzHankelPlan{S}(T,C,DL,DR) where {S} = new{S}(T,C,DL,DR)
end

function ToeplitzHankelPlan(T::TriangularToeplitz,C::Vector,DL::AbstractVector,DR::AbstractVector)
    S=promote_type(eltype(T),eltype(C[1]),eltype(DL),eltype(DR))
    ToeplitzHankelPlan{S}(T,C,collect(S,DL),collect(S,DR))
end
ToeplitzHankelPlan(T::TriangularToeplitz,C::Matrix) =
    ToeplitzHankelPlan(T,C,ones(size(T,1)),ones(size(T,2)))
ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,DL::AbstractVector,DR::AbstractVector) =
    ToeplitzHankelPlan(T,partialchol(H),DL,DR)
ToeplitzHankelPlan(T::TriangularToeplitz,H::Hankel,D::AbstractVector,DL::AbstractVector,DR::AbstractVector) =
    ToeplitzHankelPlan(T,partialchol(H,D),DL,DR)

*(P::ToeplitzHankelPlan,v::AbstractVector)=P.DL.*toeplitzcholmult(P.T,P.C,P.DR.*v)

function partialchol(H::Hankel)
    # Assumes positive definite
    σ = eltype(H)[]
    n = size(H,1)
    C = Vector{eltype(H)}[]
    v = [H[:,1]; vec(H[end,2:end])]
    d = diag(H)
    @assert length(v) ≥ 2n-1
    reltol = maximum(abs,d)*eps(eltype(H))*log(n)
    for k=1:n
        mx,idx = findmax(d)
        if mx ≤ reltol break end
        push!(σ, inv(mx))
        push!(C, v[idx:n+idx-1])
        for j=1:k-1
            nCjidxσj = -C[j][idx]*σ[j]
            Compat.LinearAlgebra.axpy!(nCjidxσj, C[j], C[k])
        end
        @simd for p=1:n
            @inbounds d[p] -= C[k][p]^2/mx
        end
    end
    for k=1:length(σ) rmul!(C[k],sqrt(σ[k])) end
    C
end

function partialchol(H::Hankel, D::AbstractVector)
    # Assumes positive definite
    T = promote_type(eltype(H),eltype(D))
    σ = T[]
    n = size(H,1)
    C = Vector{T}[]
    v = [H[:,1];vec(H[end,2:end])]
    d = diag(H).*D.^2
    @assert length(v) ≥ 2n-1
    reltol = maximum(abs,d)*eps(T)*log(n)
    for k=1:n
        mx,idx=findmax(d)
        if mx ≤ reltol break end
        push!(σ,inv(mx))
        push!(C,v[idx:n+idx-1].*D.*D[idx])
        for j=1:k-1
            nCjidxσj = -C[j][idx]*σ[j]
            Compat.LinearAlgebra.axpy!(nCjidxσj, C[j], C[k])
        end
        @simd for p=1:n
            @inbounds d[p]-=C[k][p]^2/mx
        end
    end
    for k=1:length(σ) rmul!(C[k],sqrt(σ[k])) end
    C
end

function toeplitzcholmult(T,C,v)
    n,K = length(v),length(C)
    ret,temp1,temp2 = zero(v),zero(v),zero(v)
    un,ze = one(eltype(v)),zero(eltype(v))
    broadcast!(*, temp1, C[K], v)
    mul!(temp2, T, temp1, un, ze)
    broadcast!(*, ret, C[K], temp2)
    for k=K-1:-1:1
        broadcast!(*, temp1, C[k], v)
        mul!(temp2, T, temp1, un, ze)
        broadcast!(*, temp1, C[k], temp2)
        broadcast!(+, ret, ret, temp1)
    end
    ret
end

# Diagonally-scaled Toeplitz∘Hankel polynomial transforms

function leg2chebTH(::Type{S},n) where S
    λ = Λ.(0:half(S):n-1)
    t = zeros(S,n)
    t[1:2:end] = λ[1:2:n]
    T = TriangularToeplitz(2t/π,:U)
    H = Hankel(λ[1:n], λ[n:end])
    DL = ones(S,n)
    DL[1] /= 2
    T,H,DL
end

function cheb2legTH(::Type{S},n) where S
    t = zeros(S,n-1)
    t[1:2:end] = Λ.(0:one(S):div(n-2,2), -half(S), one(S))
    T = TriangularToeplitz(t,:U)
    h = Λ.(1:half(S):n-1, zero(S), 3half(S))
    H = Hankel(h[1:n-1],h[n-1:end])
    D = 1:one(S):n-1
    DL = (3half(S):n-half(S))./D
    DR = -(one(S):n-one(S))./4D
    T,H,D,DL,DR
end

function leg2chebuTH(::Type{S},n) where S
    λ = Λ.(0:half(S):n-1)
    t = zeros(S,n)
    t[1:2:end] = λ[1:2:n]./(((1:2:n).-2))
    T = TriangularToeplitz(-2t/π,:U)
    H = Hankel(λ[1:n]./((1:n).+1),λ[n:end]./((n:2n-1).+1))
    T,H
end

function ultra2ultraTH(::Type{S},n,λ₁,λ₂) where S
    @assert abs(λ₁-λ₂) < 1
    DL = (zero(S):n-one(S)) .+ λ₂
    jk = 0:half(S):n-1
    t = zeros(S,n)
    t[1:2:n] = Λ.(jk,λ₁-λ₂,one(S))[1:2:n]
    T = TriangularToeplitz(lmul!(inv(gamma(λ₁-λ₂)),t),:U)
    h = Λ.(jk,λ₁,λ₂+one(S))
    lmul!(gamma(λ₂)/gamma(λ₁),h)
    H = Hankel(h[1:n],h[n:end])
    DR = ones(S,n)
    T,H,DL,DR
end

function jac2jacTH(::Type{S},n,α,β,γ,δ) where S
    if β == δ
        @assert abs(α-γ) < 1
        @assert α+β > -1
        jk = zero(S):n-one(S)
        DL = (2jk .+ γ .+ β .+ one(S)).*Λ.(jk,γ+β+one(S),β+one(S))
        T = TriangularToeplitz(Λ.(jk,α-γ,one(S)),:U)
        H = Hankel(Λ.(jk,α+β+one(S),γ+β+two(S)),Λ.(jk.+n.-one(S),α+β+one(S),γ+β+two(S)))
        DR = Λ.(jk,β+one(S),α+β+one(S))./gamma(α-γ)
        T,H,DL,DR
    elseif α == γ
        T,H,DL,DR = jac2jacTH(S,n,β,α,δ,γ)
        ve = T.ve
        @inbounds for k = 2:2:length(ve)
            ve[k] *= -1
        end
        TriangularToeplitz(ve, :U),H,DL,DR
    else
        throw(ArgumentError("Cannot create Toeplitz dot Hankel, use a sequence of plans."))
    end
end

struct ChebyshevToLegendrePlanTH{TH}
    toeplitzhankel::TH
end

ChebyshevToLegendrePlanTH(::Type{S},n) where {S} = ChebyshevToLegendrePlanTH(th_cheb2legplan(S,n))

function *(P::ChebyshevToLegendrePlanTH,v::AbstractVector)
    w = zero(v)
    S,n = eltype(v),length(v)
    w[1:2:end] = -one(S)./(one(S):two(S):n)./(-one(S):two(S):n-two(S))
    [dot(w,v);P.toeplitzhankel*view(v,2:n)]
end

th_leg2chebplan(::Type{S},n) where {S} = ToeplitzHankelPlan(leg2chebTH(S,n)...,ones(S,n))
th_cheb2legplan(::Type{S},n) where {S} = ChebyshevToLegendrePlanTH(ToeplitzHankelPlan(cheb2legTH(S,n)...))
th_leg2chebuplan(::Type{S},n) where {S} = ToeplitzHankelPlan(leg2chebuTH(S,n)...,1:n,ones(S,n))
th_ultra2ultraplan(::Type{S},n,λ₁,λ₂) where {S} = ToeplitzHankelPlan(ultra2ultraTH(S,n,λ₁,λ₂)...)
th_jac2jacplan(::Type{S},n,α,β,γ,δ) where {S} = ToeplitzHankelPlan(jac2jacTH(S,n,α,β,γ,δ)...)

th_leg2cheb(v) = th_leg2chebplan(eltype(v),length(v))*v
th_cheb2leg(v) = th_cheb2legplan(eltype(v),length(v))*v
th_leg2chebu(v) = th_leg2chebuplan(eltype(v),length(v))*v
th_ultra2ultra(v,λ₁,λ₂) = th_ultra2ultraplan(eltype(v),length(v),λ₁,λ₂)*v
th_jac2jac(v,α,β,γ,δ) = th_jac2jacplan(eltype(v),length(v),α,β,γ,δ)*v


text:
@testset "Scripts" begin

    # julia bin
    julia = joinpath(Base.Sys.BINDIR, "julia")
    # ../../
    mitos_folder = splitdir(splitdir(dirname(@__FILE__))[1])[1]

    @testset "Distances.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "Distances.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "small.pdb")
        out_intra = read(`$julia $path_script $path_file -o STDOUT`, String)
        out_inter = read(`$julia $path_script $path_file --inter -o STDOUT`, String)

        intra = Float64[parse(Float64,split(line,',')[end]) for line in
            split(out_intra,'\n') if !startswith(line,'#') && occursin(r"\d\.\d+$",line)]
        @test sum(intra) == 2.0*3.7836265671971385
        @test length(collect((m.match for m = eachmatch(r"\n1,A,", out_intra)))) == 1
        @test length(collect((m.match for m = eachmatch(r"\n1,B,", out_intra)))) == 1

        inter = Float64[parse(Float64,split(line,',')[end]) for line in
            split(out_inter,'\n') if !startswith(line,'#') && occursin(r"\d\.\d+$",line)]
        @test sum(inter) == 4.0*3.7836265671971385
        @test length(collect((m.match for m = eachmatch(r"\n1,A,", out_inter)))) == 5
        @test length(collect((m.match for m = eachmatch(r"\n1,B,", out_inter)))) == 1
    end

    @testset "AlignedColumns.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "AlignedColumns.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "PF09645_full.stockholm")
        output = read(`$julia $path_script $path_file -o STDOUT`, String)

        @test read(path_file, Stockholm, generatemapping=true, useidcoordinates=true,
            deletefullgaps=true) == parse(output, Stockholm)
    end

    @testset "BLMI.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "BLMI.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "simple.fasta")
        output = read(`$julia $path_script $path_file -o STDOUT --samples 0 --format FASTA --apc`, String)

        @test occursin(r"1,2,0.0,0.17",output)
    end

    @testset "Buslje09.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "Buslje09.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "simple.fasta")
        output = read(`$julia $path_script $path_file -o STDOUT --samples 0 --format FASTA --apc`, String)

        @test occursin(r"1,2,0.0,0.13",output)
    end

    @testset "Conservation.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "Conservation.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "simple.fasta")
        output = read(`$julia $path_script $path_file -o STDOUT -c --format FASTA`, String)

        @test occursin(r"1,0.6931471805599453,2.0901394960274127",output)
        @test occursin(r"2,0.6931471805599453,2.0901394960274127",output)
    end

    @testset "DownloadPDB.jl" begin
        # Helper functions
        _get_files() = filter!(f -> startswith(f, "3NIR"),readdir())
        function _delete_files()
            for file in _get_files()
                rm(file)
            end
        end
        # Test
        _delete_files()
        path_script = joinpath(mitos_folder, "scripts", "DownloadPDB.jl")
        output = read(`$julia $path_script --code 3NIR`, String)
        files = _get_files()
        @test length(files) == 1
        _delete_files()
    end

    @testset "MSADescription.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "MSADescription.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "PF09645_full.stockholm")
        output = read(`$julia $path_script $path_file -o STDOUT`, String)

        @test occursin(r"clusters,number,,4", output)
        @test occursin(r"gapfraction,quantile,0.75,0.25", output)
        @test occursin(r"coverage,mean,,0.80", output)
    end

    @testset "PairwiseGapPercentage.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "PairwiseGapPercentage.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "PF09645_full.stockholm")
        output = read(`$julia $path_script $path_file -o STDOUT`, String)

        @test occursin(r"115,115,75.0,75.0", output)
        @test occursin(r"40,70,0.0,0.0", output)
        @test occursin(r"115,115,75.0,75.0", output)
    end

    @testset "PercentIdentity.jl" begin
        path_script = joinpath(mitos_folder, "scripts", "PercentIdentity.jl")
        path_file   = joinpath(mitos_folder, "test", "data", "PF09645_full.stockholm")
        output = read(`$julia $path_script $path_file -o STDOUT`, String)

        @test occursin(
            r"110,4,29.5,15.21,14.13,19.75[0-9]*,26.3[0-9]+,33.78[0-9]*,56.38", output)
    end
end


text:
# Active patterns
module Active
using MLStyle
using MLStyle.Qualification
using MLStyle.AbstractPatterns

export @active, active_def
@nospecialize
function active_def(P, body, mod::Module, line::LineNumberNode)

    @switch P begin
        @case Expr(:call,
            Expr(:curly, t, type_args...) || t && let type_args = [] end,
            arg
        ) && if t isa Symbol end

        @case _
            error("malformed active pattern definition: $P")
    end

    definition = if isdefined(mod, t)
        line
    else
        :(struct $t end)
    end
    parametric = isempty(type_args) ? [] : type_args
    prepr = "$P"
    token = gensym(prepr)
    v_ty = Val{(view, token)}
    v_val = Val((view, token))

    quote
        $definition
        (::$v_ty)($(parametric...), ) = $arg -> $body
        $line
        function $MLStyle.pattern_uncall(t::($t isa Function ? typeof($t) : Type{$t}), self::Function, type_params, type_args, args)
            $line
            isempty(type_params) || error("A ($t) pattern requires no type params.")
            parametric = isempty(type_args) ? [] : type_args
            n_args = length(args)

            function trans(expr)
                Expr(:call, Expr(:call, $v_val, parametric...), expr)
            end

            function guard2(expr)
                if n_args === 0
                    :($expr isa Bool && $expr)
                elseif n_args === 1
                    expr_s = "$t(x)"
                    msg = "invalid use of active patterns: " *
                          "1-ary view pattern($expr_s) should accept Union{Some{T}, Nothing} " *
                          "instead of Union{T, Nothing}! " *
                           "A simple solution is:\n" *
                           "  (@active $expr_s ex) =>\n  (@active $expr_s let r=ex; r === nothing? r : Some(r)) end"

                    :($expr !== nothing && ($expr isa $Some || begin
                        $error($msg)
                    end))
                else
                    :($expr isa $Tuple && length($expr) === $n_args)
                end
            end

            extract = if n_args <= 1
                function (expr::Any, i::Int, ::Any, ::Any)
                    expr
                end
            else
                function (expr::Any, i::Int, ::Any, ::Any)
                    :($expr[$i])
                end
            end

            type_infer(_...) = Any

            comp = $PComp(
                $prepr, type_infer;
                view=$SimpleCachablePre(trans),
                guard2=$NoncachablePre(guard2)
            )
            ps = if n_args === 0
                []
            elseif n_args === 1
                [self(Expr(:call, Some, args[1]))]
            else
                [self(e) for e in args]
            end
            $decons(comp, extract, ps)
        end
    end
end

"""
Simple active pattern implementation.
You can give a qualifier in the first argument of `@active` to customize its visibility in other modules.
```julia
    @active F(x) begin
        if x > 0
            nothing
        else
            :ok
        end
    end

    @match -1 begin
        F(:ok) => false
        _ => true
    end # true

    @active public IsEven(x) begin
        x % 2 === 0
    end

    @match 4 begin
        IsEven() => :ok
        _ => :err
    end # :ok
```
"""
macro active(qualifier, case, active_body)
    deprecate_qualifier_macro(qualifier, __source__)
    active_def(case, active_body, __module__, __source__) |> esc
end

macro active(case, active_body)
    active_def(case, active_body, __module__, __source__) |> esc
end
@specialize

end


text:
function get_parameters_drivers(obj::EstimationModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::EstimationModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
using ElectrochemicalKinetics
using Test

@testset "ElectrochemicalKinetics.jl" begin
    # Write your tests here.
end


text:
if haskey(ENV, "CI")
    ENV["PLOTS_TEST"] = "true"
    ENV["GKSwstype"] = "100" # gr segfault workaround
end

using FluxOptTools, Optim, Zygote, Flux, Plots, Test, Statistics, Random
##
@testset "FluxOptTools" begin
@info "Testing FluxOptTools"
@testset "copy" begin
@info "Testing copy"

m = Chain(Dense(1,5,tanh), Dense(5,5,tanh) , Dense(5,1))
x = collect(LinRange(-pi,pi,100)')
y = sin.(x)
sp = sortperm(x[:])

loss() = mean(abs2, m(x) .- y)
Zygote.refresh()
pars = Flux.params(m)
pars0 = deepcopy(pars)
npars = veclength(pars)
@test npars == 46

copy!(pars, zeros(pars))
@test all(all(iszero, p) for p in pars)

p = zeros(pars)
copy!(pars, 1:npars)
copy!(p, pars)
@test p == 1:npars

grads = Zygote.gradient(loss, pars)
grads0 = deepcopy(grads)

copy!(grads, zeros(grads))
@test all(all(iszero,grads[k]) for k in keys(grads.grads))

p = zeros(grads)
copy!(grads, 1:npars)
copy!(p, grads)
@test p == 1:npars
end

## Test optimization ============================================

end

# NOTE: tests below fail if they are in a testset, probably Zygote's fault

m = Chain(Dense(1,5,tanh), Dense(5,5,tanh) , Dense(5,1))
x = collect(LinRange(-pi,pi,100)')
y = sin.(x)
sp = sortperm(x[:])

loss() = mean(abs2, m(x) .- y)
@show loss()
Zygote.refresh()
pars = Flux.params(m)

opt = ADAM(0.01)
@show loss()
for i = 1:500
    grads = Zygote.gradient(loss, pars)
    Flux.Optimise.update!(opt, pars, grads)
end
@show loss()
@test loss() < 1e-1
plot(x[sp], [y[sp] m(x)[sp]])  |> display
plot(loss, pars, l=0.5, npoints=50, seriestype=:contour) |> display

lossfun, gradfun, fg!, p0 = optfuns(loss, pars)
res = Optim.optimize(Optim.only_fg!(fg!), p0, BFGS())
@test loss() < 1e-3
plot(loss, pars, l=0.1, npoints=50) |> display

plot(x[sp], [y[sp] m(x)[sp]]) |> display

## Benchmark Optim vs ADAM
losses_adam = map(1:10) do i
    @show i
    Random.seed!(i)
    m = Chain(Dense(1,5,tanh), Dense(5,5,tanh) , Dense(5,1))
    x = collect(LinRange(-pi,pi,100)')
    y = sin.(x)
    loss() = mean(abs2, m(x) .- y)
    Zygote.refresh()
    pars = Flux.params(m)
    opt = Flux.ADAM(0.2)
    trace = [loss()]
    for i = 1:500
        l,back = Zygote.pullback(loss, pars)
        push!(trace, l)
        grads = back(l)
        Flux.Optimise.update!(opt, pars, grads)
    end
    trace
end

res_lbfgs = map(1:10) do i
    @show i
    Random.seed!(i)
    m = Chain(Dense(1,5,tanh), Dense(5,5,tanh) , Dense(5,1))
    x = LinRange(-pi,pi,100)'
    y = sin.(x)
    loss() = mean(abs2, m(x) .- y)
    Zygote.refresh()
    pars = Flux.params(m)
    lossfun, gradfun, fg!, p0 = optfuns(loss, pars)
    res = Optim.optimize(Optim.only_fg!(fg!), p0, BFGS(), Optim.Options(iterations=500, store_trace=true))
    res
end

losses_SLBFGS = map(1:10) do i
    @show i
    Random.seed!(i)
    m = Chain(Dense(1,5,tanh), Dense(5,5,tanh) , Dense(5,1))
    x = LinRange(-pi,pi,100)'
    y = sin.(x)
    loss() = mean(abs2, m(x) .- y)
    Zygote.refresh()
    pars = Flux.params(m)
    lossfun, gradfun, fg!, p0 = optfuns(loss, pars)
    opt = SLBFGS(lossfun,p0; m=3, ᾱ=1., ρ=false, λ=.0001, κ=0.1)
    function train(opt, p0, iters=20)
        p = copy(p0)
        g = zeros(veclength(pars))
        trace = [loss()]
        for i = 1:iters
            g = gradfun(g,p)
            p = apply(opt, g, p)
            push!(trace, opt.fold)
        end
        trace
    end
    trace = train(opt,p0, 500)
end

##
valuetrace(r) = getfield.(r.trace, :value)
valuetraces = valuetrace.(res_lbfgs)
plot(valuetraces, yscale=:log10, xscale=:identity, lab="", c=:red)
plot!(losses_adam, lab="", c=:blue, xlabel="Epochs", ylabel="Loss")
plot!(losses_SLBFGS, lab="", c=:green)


text:
# Tree.jl: Based on Sedgewick's "Left-leaning Red-Black Tree".

# Red-black trees are binary search trees obeying two key invariants:
# 1) All paths from root to leaf have same number of black nodes.
# 2) Red nodes always have black children.

# (c) 2016 E. Adrián Garro S. Costa Rica Institute of Technology.

import Base.isempty

RED = true
BLACK = false

"""
Tree helper node data type.
"""
type Node{Key, Value}
	# key
	key::Key
	# associated data
	value::Value
	# color of parent link
	color::Bool
	# link to left subtree
	left::Nullable{Node{Key, Value}}
	# link to right subtree
	right::Nullable{Node{Key, Value}}
	# constructor
	function Node(key::Key, value::Value, color::Bool)
		new(
			key, value, color,
			Nullable{Node{Key, Value}}(),
			Nullable{Node{Key, Value}}()
		)
	end
end

"""
Tree structure.
"""
type Tree{Key, Value}
	# root of the Tree
	root::Nullable{Node{Key, Value}}
	function Tree()
	  new(Nullable{Node{Key, Value}}())
	end
end

"""
Check if color of node x is red.
"""
function isred{Key, Value}(x::Nullable{Node{Key, Value}})
	if isnull(x)
		false
	else
		get(x).color == RED
	end
end

"""
Is tree t empty?
"""
function isempty{Key, Value}(t::Tree{Key, Value})
	isnull(t.root)
end

"""
Returns the value associated with the given key.
return error if the search returns a null node.
"""
function search{Key, Value}(t::Tree{Key, Value}, key::Key)
	result = search(t.root, key)
	if isnull(result)
		error("Key not found.")
	end
	get(result)
end

"""
Value associated with the given key in subtree rooted at x;
null value if no such key.
"""
function search{Key, Value}(x::Nullable{Node{Key, Value}}, key::Key)
	while !isnull(x) && key != get(x).key
		if key < get(x).key
			x = get(x).left
		else
			x = get(x).right
		end
	end
	if isnull(x)
		return Nullable{Value}()
	else
		return Nullable(get(x).value)
	end
end

"""
Does this tree contain the given key?
"""
function contains{Key, Value}(t::Tree{Key, Value}, key::Key)
	!isnull(search(t.root, key))
end

"""
Inserts the specified key-value pair into the tree,
overwriting the old value with the new value if the
tree already contains the specified key.
"""
function insert!{Key, Value}(t::Tree{Key, Value}, key::Key, value::Value)
	t.root = insert!(t.root, key, value)
	get(t.root).color = BLACK
end

"""
Insert the key-value pair in the subtree rooted at h.
"""
function insert!{Key, Value}(h::Nullable{Node{Key, Value}}, key::Key, value::Value)
	if isnull(h)
		Nullable(Node{Key, Value}(key, value, RED))
	else
		if key < get(h).key
			get(h).left = insert!(get(h).left, key, value)
		elseif key > get(h).key
			get(h).right = insert!(get(h).right, key, value)
		else
			get(h).value = value
		end
		# fix-up any right-leaning links
		if isred(get(h).right) && !isred(get(h).left)
			h = rotate_left(h)
		end
		if isred(get(h).left) && isred(get(get(h).left).left)
			h = rotate_right(h)
		end
		if isred(get(h).left) && isred(get(h).right)
			flip_colors!(h)
		end
		h
	end
end

"""
Removes the smallest key and associated value from the tree.
error if the tree is empty.
"""
function delete_min!{Key, Value}(t::Tree{Key, Value})
	if isempty(t)
		error("BST underflow.")
	end
	# if both children of root are black, set root to red
	if !isred(get(t.root).left) && !isred(get(t.root).right)
		get(t.root).color = RED
	end
	t.root = delete_min!(t.root)
	if !isempty(t)
		get(t.root).color = BLACK
	end
end

"""
delete! the key-value pair with the minimum key rooted at h.
"""
function delete_min!{Key, Value}(h::Nullable{Node{Key, Value}})
	if !isnull(get(h).left)
		if !isred(get(h).left) && !isred(get(get(h).left).left)
			h = move_red_left(h)
		end
		get(h).left = delete_min!(get(h).left)
		balance(h)
	end
end

"""
Removes the specified key and its associated value
from this tree (if the key is in this tree).
error if key is not in tree.
"""
function delete!{Key, Value}(t::Tree{Key, Value}, key::Key)
	if !contains(t, key)
		error("Key not found.")
	end
	# if both children of root are black, set root to red
	if !isred(get(t.root).left) && !isred(get(t.root).right)
		get(t.root).color = RED
	end
	t.root = delete!(t.root, key)
	if !isempty(t)
		get(t.root).color = BLACK
	end
end

"""
delete! the key-value pair with the given key rooted at h.
"""
function delete!{Key, Value}(h::Nullable{Node{Key, Value}}, key::Key)
	if key < get(h).key
		if !isred(get(h).left) && !isred(get(get(h).left).left)
			h = move_red_left(h)
		end
		get(h).left = delete!(get(h).left, key)
	else
		if isred(get(h).left)
			h = rotate_right(h)
		end
		if key == get(h).key && isnull(get(h).right)
			return Nullable{Node{Key, Value}}()
		end
		if !isred(get(h).right) && !isred(get(get(h).right).left)
			h = move_red_right(h)
		end
		if key == get(h).key
			x = min(get(h).right)
			get(h).key = get(x).key
			get(h).value = get(x).value
			get(h).right = delete_min!(get(h).right)
		else
			get(h).right = delete!(get(h).right, key)
		end
	end
	balance(h)
end

"""
Delete all nodes from tree.
"""
function clear!{Key, Value}(t::Tree{Key, Value})
	clear!(t.root)
	t.root = Nullable{Node{Key, Value}}()
end

function clear!{Key, Value}(x::Nullable{Node{Key, Value}})
	if !isnull(x)
		# recurse: visit all nodes in the two subtrees.
		clear!(get(x).left)
		clear!(get(x).right)
		# after both subtrees have been visited,
		# set links of this node to None
		get(x).left = Nullable{Node{Key, Value}}()
		get(x).right = Nullable{Node{Key, Value}}()
	end
end

"""
Make a left-leaning link lean to the right.
"""
function rotate_right{Key, Value}(h::Nullable{Node{Key, Value}})
	x = get(h).left
	get(h).left = get(x).right
	get(x).right = h
	get(x).color = get(get(x).right).color
	get(get(x).right).color = RED
	x
end

"""
Make a right-leaning link lean to the left.
"""
function rotate_left{Key, Value}(h::Nullable{Node{Key, Value}})
	x = get(h).right
	get(h).right = get(x).left
	get(x).left = h
	get(x).color = get(get(x).left).color
	get(get(x).left).color = RED
	x
end

"""
Flip the colors of a node and its two children,
h must have opposite color of its two children.
"""
function flip_colors!{Key, Value}(h::Nullable{Node{Key, Value}})
	get(h).color = !get(h).color
	get(get(h).left).color = !get(get(h).left).color
	get(get(h).right).color = !get(get(h).right).color
end

"""
Assuming that h is red and both h.left and h.left.left
are black, make h.left or one of its children red.
"""
function move_red_left{Key, Value}(h::Nullable{Node{Key, Value}})
	flip_colors!(h)
	if isred(get(get(h).right).left)
		get(h).right = rotate_right(get(h).right)
		h = rotate_left(h)
		flip_colors!(h)
	end
	h
end

"""
Assuming that h is red and both h.right and h.right.left
are black, make h.right or one of its children red.
"""
function move_red_right{Key, Value}(h::Nullable{Node{Key, Value}})
	flip_colors!(h)
	if isred(get(get(h).left).left)
		h = rotate_right(h)
		flip_colors!(h)
	end
	h
end

"""
Restore red-black tree invariant.
"""
function balance{Key, Value}(h::Nullable{Node{Key, Value}})
	if isred(get(h).right)
		h = rotate_left(h)
	end
	if isred(get(h).left) && isred(get(get(h).left).left)
		h = rotate_right(h)
	end
	if isred(get(h).left) && isred(get(h).right)
		flip_colors!(h)
	end
	h
end

"""
Returns the smallest key in the tree,
error if the tree is empty.
"""
function min{Key, Value}(t::Tree{Key, Value})
	if isempty(t)
		error("Called min with empty tree.")
	end
	get(min(t.root)).key
end

"""
Node with smallest key in subtree rooted at x;
"""
function min{Key, Value}(x::Nullable{Node{Key, Value}})
	while !isnull(get(x).left)
		x = get(x).left
	end
	x
end

"""
Returns the biggest key in the tree,
error if the tree is empty.
"""
function max{Key, Value}(t::Tree{Key, Value})
	if isempty(t)
		error("Called max with empty tree.")
	end
	get(max(t.root)).key
end

"""
Node with biggest key in subtree rooted at x;
"""
function max{Key, Value}(x::Nullable{Node{Key, Value}})
	while !isnull(get(x).right)
		x = get(x).right
	end
	x
end

function inorderwalk{Key, Value}(t::Tree{Key, Value})
	inorderwalk(t.root)
end

function inorderwalk{Key, Value}(x::Nullable{Node{Key, Value}})
	s = []
	while !isempty(s) || !isnull(x)
		if !isnull(x)
			push!(s, x)
			x = get(x).left
		else
			x = pop!(s)
			produce(get(x).value)
			x = get(x).right
		end
	end
end

function preorderwalk{Key, Value}(t::Tree{Key, Value})
	preorderwalk(t.root)
end

function preorderwalk{Key, Value}(x::Nullable{Node{Key, Value}})
	if !isnull(x)
		depth = 0
		s = []
		push!(s, (x, depth))
		while !isempty(s)
			node_and_depth = pop!(s)
			x = node_and_depth[1]
			depth = node_and_depth[2]
			produce(node_and_depth)
			if !isnull(get(x).right)
				push!(s, (get(x).right, depth + 1))
			end
			if !isnull(get(x).left)
				push!(s, (get(x).left, depth + 1))
			end
		end
	end
end

function prettyprint{Key, Value}(t::Tree{Key, Value})
	prettyprint(t.root, 0)
end

function prettyprint{Key, Value}(h::Nullable{Node{Key, Value}}, depth::Int64)
	for (node, depth) in @task preorderwalk(h)
		print("  " ^ depth)
		print("+-")
		if get(node).color == RED
			println(
				string("\033[91m", string(
						string(get(node).value), "\033[0m"
					)
				)
			)
		else
			println(
				string("\033[90m", string(
						string(get(node).value), "\033[0m"
					)
				)
			)
		end
	end
end

"""
Returns true if a binary tree is a binary search tree.
"""
function isbst{Key, Value}(t::Tree{Key, Value})
	if !isempty(t)
		minkey = min(t)
		maxkey = max(t)
		isbst(t.root, minkey, maxkey)
	else
		true
	end
end

"""
Returns true if a binary tree is a binary search tree.
"""
function isbst{Key, Value}(x::Nullable{Node{Key, Value}}, minkey::Key, maxkey::Key)
	if isnull(x)
		return true
	end
	if get(x).key < minkey
		return false
	end
	if get(x).key > maxkey
		return false
	end
	isbst(get(x).left, minkey, get(x).key) && isbst(get(x).right, get(x).key, maxkey)
end

"""
Do all paths from x to leaf have same number of black nodes?
"""
function isbalanced{Key, Value}(t::Tree{Key, Value})
	# number of black links on path from root to min.
	blacklinks = 0
	x = t.root
	while !isnull(x)
		if !isred(x)
			blacklinks += 1
		end
		x = get(x).left
	end
	isbalanced(t.root, blacklinks)
end

"""
Do all paths from x to leaf have same number of black nodes?
"""
function isbalanced{Key, Value}(x::Nullable{Node{Key, Value}}, blacklinks::Int64)
	if isnull(x)
		return blacklinks == 0
	end
	if !isred(x)
		blacklinks -= 1
	end
	isbalanced(get(x).left, blacklinks) && isbalanced(get(x).right, blacklinks)
end

max_stars_repo_name: adriangarro/llrb

max_issues_repo_name: adriangarro/llrb

max_forks_repo_name: adriangarro/llrb


text:
# Setup everything using the common setup.
include("CloudGraphSetup.jl")

@testset "Encoding and decoding of packed data types" begin
    # And check that if we encode and decode this type, it's exactly the same.
    # Make a packed data test structure.
    global fullType = DataTest(rand(10,10), "This is a test string", rand(Int32,10,10));
    typePackedRegName = string(PackedDataTest);
    typeOriginalRegName = string(DataTest);
    # Now lets encode and decode to see.
    testPackedType = convert(PackedDataTest, fullType);
    # testPackedType = cloudGraph.packedOriginalDataTypes[typeOriginalRegName].encodingFunction(PackedDataTest, fullType);

    testFullType = convert(DataTest, testPackedType)
    # testFullType = cloudGraph.packedPackedDataTypes[typePackedRegName].decodingFunction(DataTest, testPackedType);

    @test json(testFullType) == json(fullType)
end

# bsonDoc = Mongoc.BSON()
# bsonDoc["test"] = 1

@testset "BigData testing" begin

cloudVertex = CloudGraphs.CloudVertex()

@testset "Creating a CloudVertex from an ExVertex with bigdata elements..." begin

localGraph = graph(ExVertex[], ExEdge{ExVertex}[]);
#Make an ExVertex that may be encoded
v = make_vertex(localGraph, "TestVertex");
vertex = Graphs.add_vertex!(localGraph, v);
vertex.attributes["data"] = fullType;
vertex.attributes["age"] = 64;
vertex.attributes["latestEstimate"] = [0.0,0.0,0.0];
bigData = CloudGraphs.BigData();
testElementLegacy = CloudGraphs.BigDataElement("TestElement1", "Performance test dataset legacy.", rand(UInt8,100), -1); #Data element
testElementDict = CloudGraphs.BigDataElement("TestElement2", "Performance test dataset new dict type.", Dict{String, Any}("testString"=>"Test String", "randUint8"=>rand(UInt8,100)), -1); #Data element
append!(bigData.dataElements, [testElementLegacy, testElementDict]);
vertex.attributes["bigData"] = bigData;
# Now encoding the structure to CloudGraphs vertex
cloudVertex = CloudGraphs.exVertex2CloudVertex(vertex);

end

cloudNodeId = -1

@testset "Adding a vertex with bigdata elements..." begin

CloudGraphs.add_vertex!(cloudGraph, cloudVertex);
cloudNodeId = cloudVertex.neo4jNode.id
@test cloudVertex.neo4jNode.id != -1

end

@testset "Checking the big data is persisted for NeoNode $cloudNodeId..." begin

cloudVertexRet = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
@test length(cloudVertexRet.bigData.dataElements) == 2
@test cloudVertexRet.bigData.dataElements[1].data == cloudVertex.bigData.dataElements[1].data
@test json(cloudVertexRet.bigData.dataElements[2].data) == json(cloudVertex.bigData.dataElements[2].data)
@test cloudVertexRet.bigData.isRetrieved

end

@testset "Testing bigdata update method..." begin

cloudVertexOrig = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
cloudVertexUpdate = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
cloudVertexUpdate.bigData.dataElements[1].description = "Updated!"
cloudVertexUpdate.bigData.dataElements[1].data = zeros(UInt8,100)
update_vertex!(cloudGraph, cloudVertexUpdate, true)
cloudVertexRet = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
@test cloudVertexRet.bigData.dataElements[1].data != cloudVertexOrig.bigData.dataElements[1].data

end

# Saving an image as binary in a separate collection
@testset "Saving a raw string..." begin

cloudVertex = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required

sourceParams = Dict{String, Any}("database" => "TestCGRepo", "collection" => "RawImages")
testElementString = CloudGraphs.BigDataElement("TestString", "String test! Lots of lorem ipsums", "String test! Lots of lorem ipsums", -1, mimeType="application/text");
push!(cloudVertex.bigData.dataElements, testElementString);
update_vertex!(cloudGraph, cloudVertex, true)
# Reading it back to see if all is good.
bDEData = read_MongoData(cloudGraph, testElementString)
@test testElementString.data == bDEData

end

# Saving an image as binary in a separate collection
@testset "Saving an image to a custom database/collection..." begin

cloudVertex = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
fid = open(dirname(Base.source_path()) * "/IMG_1407.JPG","r")
imgBytes = read(fid)
close(fid)

sourceParams = Dict{String, Any}("database" => "TestCGRepo", "collection" => "RawImages")
testElementImg = CloudGraphs.BigDataElement("TestImage", "Image Test", imgBytes, -1, sourceParams=sourceParams); #Data element
push!(cloudVertex.bigData.dataElements, testElementImg);
update_vertex!(cloudGraph, cloudVertex, true)
# Reading it back to see if all is good
bDEData = read_MongoData(cloudGraph, testElementImg)
@test testElementImg.data == bDEData

end

@testset "Deleting nodes and all bigdata..." begin

cloudVertexGet = CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required
delete_vertex!(cloudGraph, cloudVertexGet)
@test cloudVertexGet.neo4jNode == nothing
@test cloudVertexGet.neo4jNodeId == -1

# TODO restore
# @fact_throws ErrorException CloudGraphs.get_vertex(cloudGraph, cloudNodeId, true) # fullType not required

# Going deeper for underlying read failure on bigdata
cloudVertex.bigData.dataElements[1].sourceId = "DoesntExist"
# TODO restore
# @fact_throws ErrorException CloudGraphs.read_BigData!(cloudGraph, cloudVertex)
end
end


text:
function epsilon_dominates_clear(u, v, epsilon)
  veps = v .+ epsilon
  all(u .<= veps) && any(u .< veps)
end

function epsilon_dominates_fast(u, v, epsilon)
  pareto_dominates_fast(u, v .+ epsilon)
end

function pareto_dominates_clear(u, v)
  all(u .<= v) && any(u .< v)
end

function pareto_dominates_fast(u, v)
  one_smaller = false
  for i in 1:length(u)
    if u[i] > v[i]
      return false
    elseif u[i] < v[i]
      one_smaller = true
    end
  end
  return one_smaller
end

# Returns two booleans where the first indicates e-box dominance and the 2nd indicates
# epsilon progress
function epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)
  uindex = floor(u ./ epsilon)
  vindex = floor(v ./ epsilon)
  inner_epsilon_box_dominates_and_epsilon_progress(
    u, uindex, (epsilon .* uindex),
    v, vindex, (epsilon .* vindex))
end

epsilon_box_dominates(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[1]
epsilon_box_progress(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[2]

# Inner boolean expression when the intermediate results have already been calculated.
# This is to speed up processing when there is a whole archive of solutions to compare
# a new solution to.
function inner_epsilon_box_dominates_and_epsilon_progress(
  u, uindex, uindextimesepsilon,
  v, vindex, vindextimesepsilon)

  if pareto_dominates_fast(uindex, vindex)
    return (true, true)
  elseif all(uindex .== vindex)
    smaller_distance = norm(u .- uindextimesepsilon) < norm(v .- vindextimesepsilon)
    return (smaller_distance, false)
  else
    return (false, false)
  end

end

# An EpsilonBoxArchive keeps a epsilon Pareto set of the best, non-dominated solutions
# found so far for a multi-objective search/optimization problem.
mutable struct EpsilonBoxArchive
  epsilon
  numsolutions::Int
  solutions::Matrix{Float64}       # solutions in this archive, one solution per column

  # For speedup we cache the boxes and leftpoints of all solutions. They have the same index as in
  # the solutions above.
  boxes::Matrix{Float64}      # the epsilon-boxes corresponding to each solution (same index as solutions)
  leftpoints::Matrix{Float64} # Position of "lower left" point of box (same index as solutions)

  function EpsilonBoxArchive(dimension::Integer, epsilon = 0.01, size = 100)
    new(
      epsilon,
      0,
      zeros(Float64, dimension, size),
      zeros(Float64, dimension, size),
      zeros(Float64, dimension, size)
      )
  end
end

add_if_dominates!(a::EpsilonBoxArchive)