text:
module TinyJOS

type AbstractClass

# Stadnard instance

type StdInstance
	class :: AbstractClass
	slots :: Any
end

is_instance_slot(slot) = slot_definition_allocation(slot) === :instance

std_instance_class(x :: StdInstance) = x.class
std_instance_slots(x :: StdInstance) = x.slots

const secrect_unbound_value = Symbol('unbound slot value')

allocate_std_instance(class) =
	StdInstance(
		class,
		allocate_slot_storage(
			count(is_instance_slot, class_slots(class),
			secrect_unbound_value)))

allocate_slot_storage(num, init) = fill(init, num)

the_slots_of_standard_class = Nullable()
the_class_standard_class = Nullable()

function slot_location(class, slot_name)
	if (slot_name === :effective_slots) && (class === the_class_standard_class)
		findfirst(s -> slot_definition_name(s) === :effective_slots, the_slots_of_standard_class)
	else
		pos = findfirst(s -> slot_definition_name(s) === slot_name, class_slots(class))
		if pos == 0
			error("Slot $slot_name not found in class $class")
		end
		slot = class_slots(class)[pos]
		instpos = findfirst(
			s -> slot_definition_name(s) === slot_name,
			filter(is_instance_slot, class_slots(class)))
		if instpos == 0
			error("Slot $slot_name is not an instance slot of class $class")
		end
		instpos
end

# Not sure if needed: can do with getindex and setindex! on slots type
slot_contents(slots, loc) = slots[loc]

set_slot_contents!(slots, loc, newval) = slots[loc] = newval

function std_slot_value(instance, slot_name)
	loc = slot_location(classof(instance), slot_name)
	slots = std_instance_slots(instance)
	val = slots[loc]
	if val === secrect_unbound_value
		error("Slot $slot_name is unbound in $instance")
	end
	val
end

function slot_value(instance, slot_name)
	if classof(classof(instance)) === the_class_standard_class
		std_slot_value(instance, slot_name)
	else
		slot_value_using_class(classof(instance), instance, slot_name)
	end
end

function set_std_slot_value!(instance, slot_name, newval)
	loc = slot_location((classof instance), slot_name)
	slots = std_instance_slots(instance)
	slots[loc] = newval
end

function set_slot_value!(instance, slot_name, newval)
	if classof(classof(instance)) === the_class_standard_class
		set_std_slot_value!(instance, slot_name, newval)
	else
		set_slot_value_using_class!(classof(instance), instance, slot_name, newval)
	end
end

function is_std_slot_bound(instance, slot_name)
	loc = slot_location(classof(instance), slot_name)
	slots = std_instance_slots(instance)
	stlos[loc] != secrect_unbound_value
end

function is_slot_bound(instance, slot_name)
	if classof(classof(instnace)) === the_class_standard_class
		is_std_slot_bound(instance, slot_name)
	else
		is_slot_bound_using_class(classof(instance), instnace, slot_name)
	end
end

function std_slot_make_unbound(instance, slot_name)
	loc = slot_location(classof(instnace), slot_name)
	slots = std_instance_slots(instance)
	slots[loc] = secrect_unbound_value
	nothing
end

function slot_make_unbound(instance, slot_name)
	if classof(classof(instnace)) === the_class_standard_class
		std_slot_make_unbound(instance, slot_name)
	else
		slot_make_unbound_using_class(classof(instance), instance, slot_name)
	end
end

function is_std_slot_exists(instance, slot_name)
	findfirst(
		s -> slot_definition_name(s) === slot_name,
		class_slots(classof(instance))) > 0
end

function is_slot_exists(instance, slot_name)
	if classof(classof(instance)) === the_class_standard_class
		is_std_slot_exists(instance, slot_name)
	else
		is_slot_exists_using_class(classof(instance), instance, slot_name)
	end
end

function classof(x)
	if typeof(x) == StdInstance
		std_instance_class(x)
	else
		built_in_class(x)
	end
end

function built_in_class(x)
	if typeof(x) == Bool
		findclass(:Bool)
	elseif typeof(x) == Int64
		findclass(:Int64)
	elseif typeof(x) == Float64
		findclass(:Float64)
	elseif typeof(x) == Symbol
		findclass(:Symbol)
	elseif typeof(x) == Char
		findclass(:Char)
	elseif typeof(x) == String
		findclass(:String)
	elseif typeof(x) <: Function
		findclass(:Function)
	else
		error("Not implemented for $(typeof(x))")
	end
end

function issubclass(c1, c2)
	findfirst(x -> x == c2, class_precedence_list(c1)) > 0
end

function is_subspecializer(c1, c2, cmain)
	cpl = class_precedence_list(cmain)
	c1loc = findfirst(x -> x == c1, cpl)
	if c1loc <= 0
		return false
	end
	findfirst(x -> x == c2, cpl[c1loc:end]) > 0
end

end # module


text:
srand(1)
using ExperimentalSVM
iris = readcsv(joinpath(dirname(@__FILE__), "iris.csv"))
irisX = float(iris[:, 1:4]')
p, n = size(irisX)
train = bitrand(n)
irisY = [species == "setosa" ? 1.0 : -1.0 for species in  iris[:, 5]]
iris_model_pegasos_optim = svm(irisX[:,train], irisY[train])
@assert (predict(iris_model_pegasos_optim, irisX[:, ~train])) == irisY[~train]
iris_model_cddual_optim = cddual(irisX[:,train], irisY[train])
@assert (predict(iris_model_cddual_optim, irisX[:, ~train])) == irisY[~train]
iris_model_cddual_shrinking_optim = cddual_shrinking(irisX[:,train], irisY[train])
@assert (predict(iris_model_cddual_shrinking_optim, irisX[:, ~train])) == irisY[~train]
rcvY, rcvX = read_svm_data(joinpath(dirname(@__FILE__), "rcv1_train.binary"))
rcv_model_pegasos_optim = svm(rcvX, rcvY)
rcv_model_cddual_optim = cddual(rcvX, rcvY)
rcv_model_cddual_shrinking_optim = cddual_shrinking(rcvX, rcvY)

ext: jl


text:
###############
# differences #
###############

# NOTE: values encountered during an update that are not wrapped in Diffed are
# assumed to have not changed (i.e. constants)

# NOTE: the requirement is that if a method takes any Diffed arguments, it must
# return a Diffed value.

# a Julia method may take non-diffed arguments (which indicate no change, in
# the update context) and return a non-diffed arguments. this means that
# constant expressions (and expressions that depend only on contsant
# expressions) can use code that is not overloaded for Diff types.

# Q: what about regular functions that operate on Diffed values. for example,
# suppose there was a method isdiffed(diff::Diffed) that returned true or
# false. can this function not be used in the model?

# A: a different Julia type (e.g. DiffedValue) would need to be used to
# represent Diffed values that do not represent the semantics of Diff. then
# different Julia methods would be implemented that take in Diffed and
# DiffedValues.

# a value with type Diffed should never appear outside of the update context.

"""
    abstract type Diff end

Abstract type for information about a change to a value.
"""
abstract type Diff end

"""
    UnknownChange

No information is provided about the change to the value.
"""
struct UnknownChange <: Diff end

"""
    NoChange

The value did not change.
"""
struct NoChange <: Diff end

struct SetDiff{V} <: Diff

    # elements that were added
    added::Set{V}

    # elements that were deleted
    deleted::Set{V}
end

struct DictDiff{K,V} <: Diff

    # keys that that were added and their values
    added::AbstractDict{K,V}

    # keys that were deleted
    deleted::AbstractSet{K}

    # map from key to diff value for that key
    updated::AbstractDict{K,Diff}
end

struct VectorDiff <: Diff
    new_length::Int
    prev_length::Int
    updated::Dict{Int,Diff}
end

struct IntDiff <: Diff
    difference::Int # new - old
end

export Diff
export get_diff, strip_diff
export UnknownChange, NoChange
export SetDiff, DictDiff, VectorDiff
export IntDiff

###############################
## differencing of Julia code #
###############################

"""
   Diffed{V,DV <: Diff}

Container for a value and information about a change to its value.
"""
struct Diffed{V,DV <: Diff}
   value::V
   diff::DV
end

# obtain the diff part of a Diffed value
get_diff(diffed::Diffed) = diffed.diff

# a value that is not wrapped in Diffed is a constant
get_diff(value) = NoChange()

strip_diff(diffed::Diffed) = diffed.value

strip_diff(value) = value

export Diffed, diff, strip_diff

# sets

function Base.in(element::Diffed{V,DV}, set::Diffed{T,DT}) where {V, T <: AbstractSet{V}, DV, DT}
    result = strip_diff(element) in strip_diff(set)
    Diffed(result, UnknownChange())
end

function Base.in(element::Diffed{V,NoChange}, set::Diffed{T,NoChange}) where {V, T <: AbstractSet{V}}
    result = strip_diff(element) in strip_diff(set)
    Diffed(result, NoChange())
end

function Base.in(element::Diffed{V,NoChange}, set::Diffed{T, SetDiff{V}}) where {V, T <: AbstractSet{V}}
    el = strip_diff(element)
    result = el in strip_diff(set)
    changed = (el in get_diff(set).added) || (el in get_diff(set).deleted)
    if changed
        Diffed(result, UnknownChange())
    else
        Diffed(result, NoChange())
    end
end

# TODO handle case where the set is itself a constant but the element is a Diffed

function Base.length(set::Diffed{T,UnknownChange}) where {V, T <: AbstractSet{V}}
    result = length(strip_diff(set))
    Diffed(result, UnknownChange())
end

function Base.length(set::Diffed{T,NoChange}) where {V, T <: AbstractSet{V}}
    result = length(strip_diff(set))
    Diffed(result, NoChange())
end

function Base.length(set::Diffed{T,SetDiff{V}}) where {V, T <: AbstractSet{V}}
    result = length(strip_diff(set))
    n_added = length(get_diff(set).added)
    n_deleted = length(get_diff(set).deleted)
    if n_added != n_deleted
        Diffed(result, IntDiff(n_added - n_deleted))
    else
        Diffed(result, NoChange())
    end
end

# dictionaries

function Base.haskey(dict::Diffed{T,DT}, key::Diffed{K,DK}) where {K, V, T <: AbstractDict{K,V}, DT, DK}
    result = haskey(strip_diff(dict), strip_diff(key))
    Diffed(result, UnknownChange())
end

function Base.haskey(dict::Diffed{T,NoChange}, key::Union{Diffed{K,NoChange},K}) where {K, V, T <: AbstractDict{K,V}}
    result = haskey(strip_diff(dict), strip_diff(key))
    Diffed(result, NoChange())
end

function Base.haskey(dict::Diffed{T,DictDiff{K,V}}, key::Union{Diffed{K,NoChange},K}) where {K, V, T <: AbstractDict{K,V}}
    result = haskey(strip_diff(dict), key)
    changed = (key in get_diff(dict).deleted) || haskey(get_diff(dict).added, key)
    if changed
        Diffed(result, UnknownChange())
    else
        Diffed(has, NoChange())
    end
end

function Base.getindex(dict::Diffed{T,DT}, key::Diffed{K,DK}) where {K, V, T  <: AbstractDict{K,V}, DT, DK}
    result = getindex(strip_diff(dict), key)
    Diffed(result, UnknownChange())
end

function Base.getindex(dict::Diffed{T,NoChange}, key::Union{Diffed{K,NoChange},K}) where {K, V, T <: AbstractDict{K,V}}
    result = getindex(strip_diff(dict), key)
    Diffed(result, NoChange())
end

function Base.getindex(dict::Diffed{T,DictDiff{K,V}}, key::Union{Diffed{K,NoChange},K}) where {K, V, T <: AbstractDict{K,V}}
    result = getindex(strip_diff(dict), key)
    changed = (key in get_diff(dict).deleted) || haskey(get_diff(dict).added, key)
    if changed
        Diffed(result, UnknownChange())
    else
        Diffed(val, NoChange())
    end
end

# TODO handle case where the dictionary is itself a constant, but the key is a Diffed

# vectors and tuples

function Base.length(vec::Union{Diffed{T,UnknownChange}, Diffed{T,UnknownChange}}) where {T <: Union{AbstractVector,Tuple}}
    result = length(strip_diff(vec))
    Diffed(result, UnknownChange())
end

function Base.length(vec::Diffed{T,NoChange}) where {T <: Union{AbstractVector,Tuple}}
    result = length(strip_diff(vec))
    Diffed(result, NoChange())
end

function Base.length(vec::Diffed{T,VectorDiff}) where {T <: Union{AbstractVector,Tuple}}
    len = length(strip_diff(vec))
    len_diff = get_diff(vec).new_length - get_diff(vec).prev_length
    if len_diff == 0
        Diffed(len, NoChange())
    else
        Diffed(len, IntDiff(len_diff))
    end
end

# TODO: we know that indeices before teh deleted/inserted idnex have not been changed

function Base.getindex(vec::Union{AbstractVector,Tuple}, idx::Diffed{U,DU}) where {U <: Integer, DU}
    result = vec[strip_diff(idx)]
    Diffed(result, UnknownChange())
end

function Base.getindex(vec::Union{AbstractVector,Tuple}, idx::Diffed{U,NoChange}) where {U <: Integer}
    result = vec[strip_diff(idx)]
    Diffed(result, NoChange())
end

function Base.getindex(vec::Diffed{T,DT}, idx::Union{Diffed{U,DU},Integer}) where {T <: Union{AbstractVector,Tuple}, U <: Integer, DT, DU}
    result = strip_diff(vec)[strip_diff(idx)]
    Diffed(result, UnknownChange())
end

function Base.getindex(vec::Diffed{T,NoChange}, idx::Union{Diffed{U,NoChange},Integer}) where {T <: Union{AbstractVector,Tuple}, U <: Integer}
    result = strip_diff(vec)[strip_diff(idx)]
    Diffed(result, NoChange())
end

function Base.getindex(vec::Diffed{T,VectorDiff}, idx::Union{Diffed{U,NoChange},Integer}) where {T <: Union{AbstractVector,Tuple}, U <: Integer}
    v = strip_diff(vec)
    i = strip_diff(idx)
    d = get_diff(vec)
    result = v[i]
    if i > d.prev_length
        Diffed(result, UnknownChange())
    elseif haskey(d.updated, i)
        Diffed(result, d.updated[i])
    else
        Diffed(result, NoChange())
    end
end

# TODO handle case where the vector or tuple is itself a constant, but the index is a Diffed

# built-in mathematical operators and functions

macro diffed_unary_function(fn)
    quote
        function $(fn)(value::Diffed{T,NoChange}) where {T}
            Diffed($(fn)(strip_diff(value)), NoChange())
        end

        function $(fn)(value::Diffed{T,UnknownChange}) where {T}
            Diffed($(fn)(strip_diff(value)), UnknownChange())
        end
    end
end

macro diffed_binary_function(fn)
    quote
        function $(fn)(a::Diffed{T,NoChange}, b::Diffed{U,NoChange}) where {T,U}
            result = $(fn)(strip_diff(a), strip_diff(b))
            Diffed(result, NoChange())
        end

        function $(fn)(a, b::Diffed{U,NoChange}) where {U}
            result = $(fn)(a, strip_diff(b))
            Diffed(result, NoChange())
        end

        function $(fn)(a::Diffed{T,NoChange}, b) where {T}
            result = $(fn)(strip_diff(a), b)
            Diffed(result, NoChange())
        end

        function $(fn)(a, b::Diffed{U,DU}) where {U,DU}
            result = $(fn)(a, strip_diff(b))
            Diffed(result, UnknownChange())
        end

        function $(fn)(a::Diffed{T,DT}, b) where {T,DT}
            result = $(fn)(strip_diff(a), b)
            Diffed(result, UnknownChange())
        end

        function $(fn)(a::Diffed{T,DT}, b::Diffed{U,DU}) where {T,U,DT,DU}
            result = $(fn)(strip_diff(a), strip_diff(b))
            Diffed(result, UnknownChange())
        end

        function $(fn)(a::Diffed{T,NoChange}, b::Diffed{U,DU}) where {T,U,DU}
            result = $(fn)(strip_diff(a), strip_diff(b))
            Diffed(result, UnknownChange())
        end

        function $(fn)(a::Diffed{T,DT}, b::Diffed{U,NoChange}) where {T,U,DT}
            result = $(fn)(strip_diff(a), strip_diff(b))
            Diffed(result, UnknownChange())
        end
    end
end

@diffed_binary_function Base.:+
@diffed_binary_function Base.:*
@diffed_binary_function Base.:-
@diffed_binary_function Base.:/
@diffed_binary_function Base.atan

@diffed_unary_function Base.exp
@diffed_unary_function Base.log
@diffed_unary_function Base.sqrt
@diffed_unary_function Base.sin
@diffed_unary_function Base.cos
@diffed_unary_function Base.tan

# fill

function Base.fill(value::Diffed{V,NoChange}, n::Integer) where {V}
    result = fill(strip_diff(value), strip_diff(n))
    Diffed(result, NoChange())
end

function Base.fill(value::V, n::Diffed{U,NoChange}) where {V,U <: Integer}
    result = fill(value, strip_diff(n))
    Diffed(result, NoChange())
end

function Base.fill(value::Diffed{V,NoChange}, n::Diffed{U,NoChange}) where {V,U <: Integer}
    result = fill(strip_diff(value), strip_diff(n))
    Diffed(result, NoChange())
end

function Base.fill(value::Diffed{V,DV}, n::Diffed{U,NoChange}) where {V,U <: Integer,DV}
    result = fill(strip_diff(value), strip_diff(n))
    Diffed(result, UnknownChange())
end

function Base.fill(value::Diffed{V,DV}, n::Diffed{U,DU}) where {V,U <: Integer,DU,DV}
    result = fill(strip_diff(value), strip_diff(n))
    Diffed(result, UnknownChange())
end

# TODO filter, map, reduce, foldl, etc.

# NOTE: just handle the case where the function argument is a constant (a Function not a Diffed{Function,Diff})

# broadcasting (?)

function Base.broadcast(f, a::Diffed{T,NoChange}, b::Diffed{U,NoChange}) where {T,U}
    result = broadcast(f, strip_diff(a), strip_diff(b))
    Diffed(result, NoChange())
end

function Base.broadcast(f, a::Diffed{T,NoChange}, b) where {T}
    result = broadcast(f, strip_diff(a), b)
    Diffed(result, NoChange())
end

function Base.broadcast(f, a, b::Diffed{U,NoChange}) where {U}
    result = broadcast(f, a, strip_diff(b))
    Diffed(result, NoChange())
end

function Base.broadcast(f, a::Diffed{T,DT}, b::Diffed{U,NoChange}) where {T,U,DT}
    result = broadcast(f, strip_diff(a), strip_diff(b))
    Diffed(result, UnknownChange())
end

function Base.broadcast(f, a::Diffed{T,NoChange}, b::Diffed{U,DU}) where {T,U,DU}
    result = broadcast(f, strip_diff(a), strip_diff(b))
    Diffed(result, UnknownChange())
end

function Base.broadcast(f, a::Diffed{T,DT}, b::Diffed{U,DU}) where {T,U,DT,DU}
    result = broadcast(f, strip_diff(a), strip_diff(b))
    Diffed(result, UnknownChange())
end

# control flow

function ifelse(c::Union{Bool,Diffed{Bool,NoChange}}, x::Diffed{T,NoChange}, y::Diffed{U,NoChange}) where {T,U}
    result = strip_diff(c) ? strip_diff(x) : strip_diff(y)
    Diffed(result, NoChange())
end

function ifelse(c::Diffed{Bool,DC}, x, y) where {DC}
    result = strip_diff(c) ? strip_diff(x) : strip_diff(y)
    Diffed(result, UnknownChange())
end

function ifelse(c::Bool, x, y)
    c ? x : y
end

export ifelse


text:
module LSystems

export lindenmayer

function can_transition{T<:String, U<:Char}(transitions::Dict{U, T}, c::U)
    return haskey(transitions, c)
end

function transition{T<:String, U<:Char}(transitions::Dict{U, T}, c::U)
    return transitions[c]
end

function lindenmayer{T<:String, U<:Char, V<:Int}(start::T, transitions::Dict{U, T}, level::V)
    return lindenmayer(start, c -> can_transition(transitions, c), c -> transition(transitions, c), level)
end

function lindenmayer{T<:String, V<:Int}(start::T, transition_checker::Function, transition_function::Function, level::V)
    if level <= 0
        for c in start
            produce(c)
        end
    else
        for c in start
            if transition_checker(c)
                for v in @task lindenmayer(transition_function(c), transition_checker, transition_function, level - 1)
                    produce(v)
                end
            else
                produce(c)
            end
        end
    end
end

end # module


text:
include("roblib.jl")

function find_border!(r::Robot, mov::HorizonSide, bord::HorizonSide)
    while isborder(r, bord) == false
        if isborder(r, mov) == false
            move!(r, mov)
        else
            move!(r, bord)
            mov = inverse(mov)
        end
    end
    return mov
end

function fin!(r::Robot)
    num_y1 = moves!(r, Sud)
    num_x = moves!(r, West)
    num_y2 = moves!(r, Sud)
    side = find_border!(r, Nord, Ost)
    mark_perimetr!(r, side)
    moves!(r, West)
    moves!(r, Sud)
end

fin!(r)


text:
# ---
# title: 1664. Ways to Make a Fair Array
# id: problem1664
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/ways-to-make-a-fair-array/description/>
# hidden: true
# ---
#
# You are given an integer array `nums`. You can choose **exactly one** index (
# **0-indexed** ) and remove the element. Notice that the index of the elements
# may change after the removal.
#
# For example, if `nums = [6,1,7,4,1]`:
#
#   * Choosing to remove index `1` results in `nums = [6,7,4,1]`.
#   * Choosing to remove index `2` results in `nums = [6,1,4,1]`.
#   * Choosing to remove index `4` results in `nums = [6,1,7,4]`.
#
# An array is **fair** if the sum of the odd-indexed values equals the sum of
# the even-indexed values.
#
# Return the _**number** of indices that you could choose such that after the
# removal, _`nums` ___is **fair**. _
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,1,6,4]
#     Output: 1
#     Explanation:
#     Remove index 0: [1,6,4] -> Even sum: 1 + 4 = 5. Odd sum: 6. Not fair.
#     Remove index 1: [2,6,4] -> Even sum: 2 + 4 = 6. Odd sum: 6. Fair.
#     Remove index 2: [2,1,4] -> Even sum: 2 + 4 = 6. Odd sum: 1. Not fair.
#     Remove index 3: [2,1,6] -> Even sum: 2 + 6 = 8. Odd sum: 1. Not fair.
#     There is 1 index that you can remove to make nums fair.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,1,1]
#     Output: 3
#     Explanation:  You can remove any index and the remaining array is fair.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: 0
#     Explanation:  You cannot make a fair array after removing any index.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using SparseArrays
using FastGaussQuadrature

include("SplineModel.jl")
include("MeshUtils.jl")

struct Boundary2D
    type::String #can be "Dirichlet", "Neumann" or "Robin"
    side::String #can be "Down" (v=0), "Right" (u=1), "Up" (v=1) or "Left" (u=0)
    α::Number #parameter value for Robin boundary condition (αu+u')
    op_val::Function #boundary operator value (as a function of (x,y))
    #To Do: Partial boundary conditions, Lagrange multipliers (?)
end

struct BoundaryIndex2D
    down::Array{Int64, 1}
    up::Array{Int64, 1}
    left::Array{Int64, 1}
    right::Array{Int64, 1}
end

struct Material2D
    Emod::Float64
    nu::Float64
    Cmat::Array{Float64,2}
end

"""
Constructor for Boundary2D which sets α=0
"""
function Boundary2D(type::String, side::String, op_val::Function)
    return Boundary2D(type, side, 0., op_val)
end

struct Problem2D
    rhs_f::Function
    a0::Function
    a1::Function
    boundary_cond::Array{Boundary2D}
    domain::NURBS
end

struct GaussQuad
    nodes::Array{Float64,1}
    weights::Array{Float64,1}
end

"""
Generates the Gauss-Legendre points of order n
"""
function genGaussLegendre(n)
    nodes, weights = gausslegendre(n)
    gauss_rule = GaussQuad(nodes, weights)
    return gauss_rule
end

"""
Generates the 2D Bernstein polynomials at points (pts_u, pts_v)
"""
function bernsteinBasis2D(pts_u::Array{Float64,1}, pts_v::Array{Float64,1},
                            degP::Array{Int64,1})
    B_u, dB_u = bernsteinBasis(pts_u, degP[1])
    B_v, dB_v = bernsteinBasis(pts_v, degP[2])

    numPtsU = length(pts_u)
    numPtsV = length(pts_v)
    basisCounter = 0
    Buv = zeros(numPtsU, numPtsV, (degP[1]+1)*(degP[2]+1))
    dBdu = zeros(numPtsU, numPtsV, (degP[1]+1)*(degP[2]+1))
    dBdv = zeros(numPtsU, numPtsV, (degP[1]+1)*(degP[2]+1))
    for j=1:degP[2]+1
        for i=1:degP[1]+1
            basisCounter += 1
            Buv[:,:,basisCounter] = B_u[:,i]*B_v[:,j]'
            dBdu[:,:,basisCounter] = dB_u[:,i]*B_v[:,j]'
            dBdv[:,:,basisCounter] = B_u[:,i]*dB_v[:,j]'
        end
    end
    return Buv, dBdu, dBdv
end

"""
Aseembles the stiffness matrix K_ij = ∫ a0(x)∇ϕ_i(x,y)∇ϕ_j(x) dΩ
"""
function assemble_stiff2D(mesh::Mesh, a0::Function, gauss_rule::Array{GaussQuad,1})
    Buv, dBdu, dBdv = bernsteinBasis2D(gauss_rule[1].nodes, gauss_rule[2].nodes, mesh.degP)
    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    #allocate memory for the triplet arrays
    indexCounter = 0
    for iElem=1:mesh.numElem
        numNodes = length(mesh.elemNode[iElem])
        indexCounter += numNodes^2
    end
    II = zeros(Int64, indexCounter)
    JJ = zeros(Int64, indexCounter)
    S = zeros(indexCounter)

    indexCounter = 0
    domainArea = 0
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localStiff = zeros(numNodes, numNodes)

        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)

                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2
                #compute the derivatives w.r.t to the physical space
                dR = [dRdu'; dRdv']

                dxdxi = dR * cpts'
                dR = dxdxi\dR
                Jac_par_phys = det(dxdxi)
                RR /= w_sum
                phys_pt = cpts*RR

                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                localStiff += localArea * a0(phys_pt[1], phys_pt[2]) * (dR'*dR)
                domainArea += localArea
            end
        end
        II[indexCounter+1:indexCounter+numNodes^2] = repeat(curNodes, numNodes)
        JJ[indexCounter+1:indexCounter+numNodes^2] = reshape(repeat(curNodes', numNodes, 1), numNodes^2)
        S[indexCounter+1:indexCounter+numNodes^2] = reshape(localStiff, numNodes^2)
        indexCounter += numNodes^2

    end
    @show domainArea
    stiff = sparse(II, JJ, S, mesh.numBasis, mesh.numBasis)
    return stiff
end

"""
Aseembles the stiffness matrix for 2D Elasticity Kₑ = ∫ B^T*C*B dΩ
"""
function assemble_stiff_elast2D(mesh::Mesh, Cmat::Array{Float64,2}, gauss_rule::Array{GaussQuad,1})
    Buv, dBdu, dBdv = bernsteinBasis2D(gauss_rule[1].nodes, gauss_rule[2].nodes, mesh.degP)
    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    #allocate memory for the triplet arrays
    indexCounter = 0
    for iElem=1:mesh.numElem
        numNodes = length(mesh.elemNode[iElem])
        indexCounter += 4*numNodes^2
    end
    II = zeros(Int64, indexCounter)
    JJ = zeros(Int64, indexCounter)
    S = zeros(indexCounter)

    indexCounter = 0
    domainArea = 0
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        curNodesXY = reshape(hcat(2*curNodes.-1, 2*curNodes)', 2*numNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localStiff = zeros(2*numNodes, 2*numNodes)

        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)
                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2

                #compute the derivatives w.r.t to the physical space
                dR = [dRdu'; dRdv']
                dxdxi = dR * cpts'
                dR = dxdxi\dR
                Jac_par_phys = det(dxdxi)
                RR /= w_sum
                phys_pt = cpts*RR

                B = zeros(2*numNodes,3);
                B[1:2:2*numNodes-1,1] = dR[1,:]
                B[2:2:2*numNodes,2] = dR[2,:]
                B[1:2:2*numNodes-1,3] = dR[2,:]
                B[2:2:2*numNodes,3] = dR[1,:]

                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                localStiff += localArea *  (B*Cmat*B')
                domainArea += localArea
            end
        end
        II[indexCounter+1:indexCounter+4*numNodes^2] = repeat(curNodesXY, 2*numNodes)
        JJ[indexCounter+1:indexCounter+4*numNodes^2] = reshape(repeat(curNodesXY',
                                                    2*numNodes, 1), 4*numNodes^2)
        S[indexCounter+1:indexCounter+4*numNodes^2] = reshape(localStiff, 4*numNodes^2)
        indexCounter += 4*numNodes^2

    end
    @show domainArea
    stiff = sparse(II, JJ, S, 2*mesh.numBasis, 2*mesh.numBasis)
    return stiff
end

"""
Aseembles the mass matrix M_ij = ∫ a1(x)ϕ_i(x,y)ϕ_j(x) dΩ
"""
function assemble_mass2D(mesh::Mesh, a1::Function, gauss_rule::Array{GaussQuad,1})
    Buv, dBdu, dBdv = bernsteinBasis2D(gauss_rule[1].nodes, gauss_rule[2].nodes, mesh.degP)
    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    #allocate memory for the triplet arrays
    indexCounter = 0
    for iElem=1:mesh.numElem
        numNodes = length(mesh.elemNode[iElem])
        indexCounter += numNodes^2
    end
    II = zeros(Int64, indexCounter)
    JJ = zeros(Int64, indexCounter)
    S = zeros(indexCounter)

    indexCounter = 0
    domainArea = 0
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localMass = zeros(numNodes, numNodes)

        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)

                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2
                #compute the derivatives w.r.t to the physical space
                dR = [dRdu'; dRdv']

                dxdxi = dR * cpts'
                dR = dxdxi\dR
                Jac_par_phys = det(dxdxi)
                RR /= w_sum
                phys_pt = cpts*RR

                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                localMass += localArea * a1(phys_pt[1], phys_pt[2]) * (RR*RR')
                domainArea += localArea
            end
        end
        II[indexCounter+1:indexCounter+numNodes^2] = repeat(curNodes, numNodes)
        JJ[indexCounter+1:indexCounter+numNodes^2] = reshape(repeat(curNodes', numNodes, 1), numNodes^2)
        S[indexCounter+1:indexCounter+numNodes^2] = reshape(localMass, numNodes^2)
        indexCounter += numNodes^2

    end
    @show domainArea
    stiff = sparse(II, JJ, S, mesh.numBasis, mesh.numBasis)
    return stiff
end

"""
Assembles the RHS vector corresponding to the body force
RHS[i]=∫_Ω ϕ_i(x,y)*f(x,y) dΩ
"""
function assemble_rhs2D(mesh::Mesh, f::Function,  gauss_rule::Array{GaussQuad,1})
    Buv, dBdu, dBdv = bernsteinBasis2D(gauss_rule[1].nodes, gauss_rule[2].nodes, mesh.degP)
    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    domainArea = 0
    rhs = zeros(Complex{Float64}, mesh.numBasis)
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(numNodes)
        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)

                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2
                #compute the derivatives w.r.t to the physical space
                dR = [dRdu'; dRdv']

                dxdxi = dR * cpts'
                dR = dxdxi\dR
                Jac_par_phys = det(dxdxi)

                RR /= w_sum
                phys_pt = cpts*RR
                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                localRhs += localArea * f(phys_pt[1], phys_pt[2]) * RR
                domainArea += localArea
            end
        end
        rhs[curNodes] += localRhs
    end
    return rhs
end

function uniquetol(arr, tol=1e-10)
    sort!(arr)
    delete_index = Array{Int64,1}(undef, 0)
    for i=1:length(arr)-1
        if abs(arr[i+1]-arr[i])<tol
            push!(delete_index, i+1)
        end
    end
    deleteat!(arr, delete_index)
end

"""
Helper function to compute the setdifference within a tolerance
"""
function setdifftol(arrA, arrB, tol=1e-10)
    delete_index = Array{Int64,1}(undef, 0)
    for i=1:length(arrA)
        for j=1:length(arrB)
            if abs(arrA[i]-arrB[j])<tol
                push!(delete_index, i)
                break
            end
        end
    end
    deleteat!(arrA, delete_index)
end

"""
Returns the boundary (down, right, up, left) indices for a matrix of dimension
(p+1)*(q+1)
"""
function getBoundaryIndices(degP)
    indexDown = 1:degP[1]+1
    indexRight = degP[1]+1:degP[1]+1:prod(degP.+1)
    indexUp = 1+(degP[1]+1)*degP[2]:prod(degP.+1)
    indexLeft = 1:degP[1]+1:1+(degP[1]+1)*degP[2]
    index_all = BoundaryIndex2D(indexDown, indexUp, indexLeft, indexRight)
    return index_all
end

"""
Classifies the boundary nodes and boundary elements according to the side in the
parameter space (i.e. bcdof_down and elem_down for v=0, bcdof_up and elem_up
for v=1, bcdof_left and elem_left for u=0, bcdof_right and elem_right for u=1)
"""
function classifyBoundary2D(mesh::Mesh)
    bcdof_down = Array{Int64,1}(undef, 0)
    elem_down = Array{Int64, 1}(undef, 0)
    bcdof_up = Array{Int64,1}(undef, 0)
    elem_up = Array{Int64,1}(undef, 0)
    bcdof_left = Array{Int64,1}(undef, 0)
    elem_left = Array{Int64,1}(undef, 0)
    bcdof_right = Array{Int64,1}(undef, 0)
    elem_right = Array{Int64,1}(undef, 0)
    index_all = getBoundaryIndices(mesh.degP)
    tolEq = 1e-10

    for iElem = 1:mesh.numElem
        if abs(mesh.elemVertex[iElem,1])<tolEq
            #u_min = 0
            bcdof_left = vcat(bcdof_left, mesh.elemNode[iElem][index_all.left])
            elem_left = vcat(elem_left, iElem)
        end
        if abs(mesh.elemVertex[iElem,2])<tolEq
            #v_min = 0
            bcdof_down = vcat(bcdof_down, mesh.elemNode[iElem][index_all.down])
            elem_down = vcat(elem_down, iElem)
        end
        if abs(mesh.elemVertex[iElem,3]-1)<tolEq
            #u_max = 1
            bcdof_right = vcat(bcdof_right, mesh.elemNode[iElem][index_all.right])
            elem_right = vcat(elem_right, iElem)
        end
        if abs(mesh.elemVertex[iElem,4]-1)<tolEq
            #v_max = 1
            bcdof_up = vcat(bcdof_up, mesh.elemNode[iElem][index_all.up])
            elem_up = vcat(elem_up, iElem)
        end
    end
    bcdof_down = unique(bcdof_down)
    bcdof_up = unique(bcdof_up)
    bcdof_left = unique(bcdof_left)
    bcdof_right = unique(bcdof_right)
    bcdof_all = BoundaryIndex2D(bcdof_down, bcdof_up, bcdof_left, bcdof_right)
    elem_all = BoundaryIndex2D(elem_down, elem_up, elem_left, elem_right)
    return bcdof_all, elem_all
end

"""
Apply the Neumann boundary conditions to the RHS for a 2D NURBS geometry, scalar
solution field
"""
function applyNeumannScalar2D(mesh::Mesh, rhs, elem_list, direction::String,
                                quad_rule_1d::GaussQuad, op_val::Function)
    #Select the Gauss/edge points for evaluating the  1D the basis functions
    #select the scatter vector corresponding to the edge DOFs
    index_all = getBoundaryIndices(mesh.degP)
    numNodesEdge = length(quad_rule_1d.nodes)
    if direction == "right"
        pts_u = [1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.right
    elseif direction == "left"
        pts_u = [-1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.left
    elseif direction == "down"
        pts_u = quad_rule_1d.nodes
        pts_v = [-1.]
        sctr = index_all.down
    elseif direction == "up"
        pts_u = quad_rule_1d.nodes
        pts_v = [1.]
        sctr = index_all.up
    else
        error("Wrong direction given")
    end

    #Form the 2D tensor product of the basis functions
    Buv, dBdu, dBdv = bernsteinBasis2D(pts_u, pts_v, mesh.degP)

    #Evaluate the Neumann integral on each element
    domainLength = 0.
    for iElem in elem_list
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        if direction == "right" || direction == "left"
            Jac_ref_par = (vMax-vMin)/2
        else
            Jac_ref_par = (uMax-uMin)/2
        end

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(numNodesEdge)
        for iGauss = 1:length(quad_rule_1d.nodes)
            #compute the (B-)spline basis functions and derivatives with Bezier extraction
            if direction == "right" || direction == "left"
                N_mat = mesh.C[iElem] * Buv[1, iGauss, :]
                dN_du = mesh.C[iElem] * dBdu[1, iGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[1, iGauss, :] * 2/(vMax-vMin)
            else
                N_mat = mesh.C[iElem] * Buv[iGauss, 1, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, 1, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, 1, :] * 2/(vMax-vMin)
            end

            #compute the rational basis
            RR = N_mat.* wgts
            dRdu = dN_du.* wgts
            dRdv = dN_dv.* wgts
            w_sum = sum(RR)
            dw_xi = sum(dRdu)
            dw_eta = sum(dRdv)
            dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
            dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2

            #compute the derivatives w.r.t to the physical space
            dR = [dRdu'; dRdv']
            dxdxi = dR * cpts'

            #Jacobian of face mapping
            if direction=="right" || direction=="left"
                eJac = dxdxi[2,1]^2 + dxdxi[2,2]^2
            else
                eJac = dxdxi[1,1]^2 + dxdxi[1,2]^2
            end
            Jac_par_phys = sqrt(eJac)
            dR = dxdxi\dR
            RR /= w_sum
            phys_pt = cpts*RR
            gFunc = op_val(phys_pt[1], phys_pt[2])
            localLength = Jac_par_phys * Jac_ref_par * quad_rule_1d.weights[iGauss]
            localRhs += RR[sctr].*gFunc.*localLength
            domainLength += localLength
        end
        rhs[curNodes[sctr]] += localRhs
    end
    @show domainLength
    return rhs
end

"""
Apply the Neumann boundary conditions to the RHS for a 2D NURBS geometry, elasticity
solution field
"""
function applyNeumannElast2D(mesh::Mesh, rhs, elem_list, direction::String,
                                quad_rule_1d::GaussQuad, op_val::Function)
    #Select the Gauss/edge points for evaluating the  1D the basis functions
    #select the scatter vector corresponding to the edge DOFs
    index_all = getBoundaryIndices(mesh.degP)
    numNodesEdge = length(quad_rule_1d.nodes)
    if direction == "right"
        pts_u = [1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.right
    elseif direction == "left"
        pts_u = [-1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.left
    elseif direction == "down"
        pts_u = quad_rule_1d.nodes
        pts_v = [-1.]
        sctr = index_all.down
    elseif direction == "up"
        pts_u = quad_rule_1d.nodes
        pts_v = [1.]
        sctr = index_all.up
    else
        error("Wrong direction given")
    end

    #Form the 2D tensor product of the basis functions
    Buv, dBdu, dBdv = bernsteinBasis2D(pts_u, pts_v, mesh.degP)

    #Evaluate the Neumann integral on each element
    domainLength = 0.
    for iElem in elem_list
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        if direction == "right" || direction == "left"
            Jac_ref_par = (vMax-vMin)/2
        else
            Jac_ref_par = (uMax-uMin)/2
        end

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        curNodesXY = reshape(hcat(2*curNodes[sctr].-1, 2*curNodes[sctr])', 2*numNodesEdge)

        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(2*numNodesEdge)
        for iGauss = 1:length(quad_rule_1d.nodes)
            #compute the (B-)spline basis functions and derivatives with Bezier extraction
            if direction == "right" || direction == "left"
                N_mat = mesh.C[iElem] * Buv[1, iGauss, :]
                dN_du = mesh.C[iElem] * dBdu[1, iGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[1, iGauss, :] * 2/(vMax-vMin)
            else
                N_mat = mesh.C[iElem] * Buv[iGauss, 1, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, 1, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, 1, :] * 2/(vMax-vMin)
            end

            #compute the rational basis
            RR = N_mat.* wgts
            dRdu = dN_du.* wgts
            dRdv = dN_dv.* wgts
            w_sum = sum(RR)
            dw_xi = sum(dRdu)
            dw_eta = sum(dRdv)
            dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
            dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2

            #compute the derivatives w.r.t to the physical space
            dR = [dRdu'; dRdv']
            dxdxi = dR * cpts'

            #Jacobian of face mapping
            if direction=="right" || direction=="left"
                eJac = dxdxi[2,1]^2 + dxdxi[2,2]^2
            else
                eJac = dxdxi[1,1]^2 + dxdxi[1,2]^2
            end
            Jac_par_phys = sqrt(eJac)

            dR = dxdxi\dR
            RR /= w_sum
            phys_pt = cpts*RR

            gFunc = op_val(phys_pt[1], phys_pt[2])
            localLength = Jac_par_phys * Jac_ref_par * quad_rule_1d.weights[iGauss]
            localRhs[1:2:end-1] += RR[sctr].*gFunc[1].*localLength
            localRhs[2:2:end] += RR[sctr].*gFunc[2].*localLength
            domainLength += localLength

        end
        rhs[curNodesXY] += localRhs
    end
    @show domainLength
    return rhs
end

"""
Apply the Robin boundary conditions to the RHS for a 2D NURBS geometry, scalar
solution field
"""
function applyRobinScalar2D(mesh::Mesh, lhs, rhs, elem_list, direction::String,
                                quad_rule_1d::GaussQuad, α::Number, op_val::Function)
    #Select the Gauss/edge points for evaluating the  1D the basis functions
    #select the scatter vector corresponding to the edge DOFs
    index_all = getBoundaryIndices(mesh.degP)
    numNodesEdge = length(quad_rule_1d.nodes)
    if direction == "right"
        pts_u = [1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.right
    elseif direction == "left"
        pts_u = [-1.]
        pts_v = quad_rule_1d.nodes
        sctr = index_all.left
    elseif direction == "down"
        pts_u = quad_rule_1d.nodes
        pts_v = [-1.]
        sctr = index_all.down
    elseif direction == "up"
        pts_u = quad_rule_1d.nodes
        pts_v = [1.]
        sctr = index_all.up
    else
        error("Wrong direction given")
    end

    #allocate memory for the triplet arrays
    indexCounter = length(elem_list)*length(sctr)^2
    II = zeros(Int64, indexCounter)
    JJ = zeros(Int64, indexCounter)
    S = zeros(Complex{Float64}, indexCounter)
    indexCounter = 0

    #Form the 2D tensor product of the basis functions
    Buv, dBdu, dBdv = bernsteinBasis2D(pts_u, pts_v, mesh.degP)

    #Evaluate the Neumann integral on each element
    domainLength = 0.
    for iElem in elem_list
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        if direction == "right" || direction == "left"
            Jac_ref_par = (vMax-vMin)/2
        else
            Jac_ref_par = (uMax-uMin)/2
        end

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(sctr)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localBndMass = zeros(numNodesEdge,numNodesEdge)
        for iGauss = 1:length(quad_rule_1d.nodes)
            #compute the (B-)spline basis functions and derivatives with Bezier extraction
            if direction == "right" || direction == "left"
                N_mat = mesh.C[iElem] * Buv[1, iGauss, :]
                dN_du = mesh.C[iElem] * dBdu[1, iGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[1, iGauss, :] * 2/(vMax-vMin)
            else
                N_mat = mesh.C[iElem] * Buv[iGauss, 1, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, 1, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, 1, :] * 2/(vMax-vMin)
            end

            #compute the rational basis
            RR = N_mat.* wgts
            dRdu = dN_du.* wgts
            dRdv = dN_dv.* wgts
            w_sum = sum(RR)
            dw_xi = sum(dRdu)
            dw_eta = sum(dRdv)
            dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
            dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2

            #compute the derivatives w.r.t to the physical space
            dR = [dRdu'; dRdv']
            dxdxi = dR * cpts'

            #Jacobian of face mapping
            if direction=="right" || direction=="left"
                eJac = dxdxi[1,2]^2 + dxdxi[2,2]^2
            else
                eJac = dxdxi[1,1]^2 + dxdxi[2,1]^2
            end
            Jac_par_phys = sqrt(eJac)
            dR = dxdxi\dR
            RR /= w_sum
            phys_pt = cpts*RR

            localBndMass += (RR[sctr]*RR[sctr]').*α.*Jac_ref_par.*Jac_par_phys.*quad_rule_1d.weights[iGauss]
            domainLength += Jac_par_phys * Jac_ref_par * quad_rule_1d.weights[iGauss]
        end
        II[indexCounter+1:indexCounter+numNodes^2] = repeat(curNodes[sctr], numNodes)
        JJ[indexCounter+1:indexCounter+numNodes^2] = reshape(repeat(curNodes[sctr]', numNodes, 1), numNodes^2)
        S[indexCounter+1:indexCounter+numNodes^2] = reshape(localBndMass, numNodes^2)
        indexCounter += numNodes^2
    end
    @show domainLength
    bnd_mass = sparse(II, JJ, S, mesh.numBasis, mesh.numBasis)
    lhs += bnd_mass
    return lhs, rhs
end

"""
Apply the boundary conditions to a linear system for a 2D
NURBS geometry and elasticity problems
"""
function applyBCnurbsElast(mesh::Mesh, bound_cond::Array{Boundary2D, 1}, lhs, rhs,
                            bcdof_all, elem_all, quad_rule)
    #collect the dofs and values corresponding to the boundary
    bcdof = Array{Int64,1}(undef, 0)
    bcval = Array{Float64,1}(undef, 0)
    #assume homogeneous boundary conditions
    evalPt = [0.,0.]
    for i=1:length(bound_cond)
        if bound_cond[i].type=="Dirichlet"
            if bound_cond[i].side=="Down"
                #check x-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[1]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.down.-1)
                end
                #check y-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[2]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.down)
                end
            end
            if bound_cond[i].side=="Up"
                #check x-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[1]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.up.-1)
                end
                #check y-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[2]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.up)
                end
            end
            if bound_cond[i].side=="Left"
                #check x-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[1]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.left.-1)
                end
                #check y-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[2]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.left)
                end
            end
            if bound_cond[i].side=="Right"
                #check x-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[1]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.right.-1)
                end
                #check y-direction
                if bound_cond[i].op_val(evalPt[1], evalPt[2])[2]!=undef
                    bcdof = vcat(bcdof, 2*bcdof_all.right)
                end
            end
        elseif bound_cond[i].type=="Neumann"
            if bound_cond[i].side=="Down"
                rhs = applyNeumannElast2D(mesh, rhs, elem_all.down, "down",
                                        quad_rule[1], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Up"
                rhs = applyNeumannElast2D(mesh, rhs, elem_all.up, "up",
                                        quad_rule[1], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Left"
                rhs = applyNeumannElast2D(mesh, rhs, elem_all.left, "left",
                                        quad_rule[2], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Right"
                rhs = applyNeumannElast2D(mesh, rhs, elem_all.right, "right",
                                        quad_rule[2], bound_cond[i].op_val)
            end
        else
            error("Boundary condition type not implemented")
        end
    end
    unique!(bcdof)
    bcval = zero(bcdof)
    rhs = rhs - lhs[:,bcdof]*bcval
    rhs[bcdof] = bcval
    lhs[bcdof, :] .= 0.
    lhs[:, bcdof] .= 0.
    lhs[bcdof, bcdof] = sparse(I, length(bcdof), length(bcdof))
    return lhs, rhs
end

function applyBCnurbs(mesh::Mesh, bound_cond::Array{Boundary2D, 1}, lhs, rhs,
                        bcdof_all, elem_all, quad_rule)
    #collect the dofs and values corresponding to the boundary
    bcdof = Array{Int64,1}(undef, 0)
    bcval = Array{Float64,1}(undef, 0)
    for i=1:length(bound_cond)
        if bound_cond[i].type=="Dirichlet"
            if bound_cond[i].side=="Down"
                bcdof = vcat(bcdof, bcdof_all.down)
            end
            if bound_cond[i].side=="Up"
                bcdof = vcat(bcdof, bcdof_all.up)
            end
            if bound_cond[i].side=="Left"
                bcdof = vcat(bcdof, bcdof_all.left)
            end
            if bound_cond[i].side=="Right"
                bcdof = vcat(bcdof, bcdof_all.right)
            end
        elseif bound_cond[i].type=="Neumann"
            if bound_cond[i].side=="Down"
                rhs = applyNeumannScalar2D(mesh, rhs, elem_all.down, "down",
                                        quad_rule[1], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Up"
                rhs = applyNeumannScalar2D(mesh, rhs, elem_all.up, "up",
                                        quad_rule[1], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Left"
                rhs = applyNeumannScalar2D(mesh, rhs, elem_all.left, "left",
                                        quad_rule[2], bound_cond[i].op_val)
            end
            if bound_cond[i].side=="Right"
                rhs = applyNeumannScalar2D(mesh, rhs, elem_all.right, "right",
                                        quad_rule[2], bound_cond[i].op_val)
            end
        elseif bound_cond[i].type=="Robin"
            α = bound_cond[i].α
            op_val = bound_cond[i].op_val
            if bound_cond[i].side=="Down"
                lhs, rhs = applyRobinScalar2D(mesh, lhs, rhs, elem_all.down, "down",
                                        quad_rule[1], α, op_val)
            end
            if bound_cond[i].side=="Up"
                lhs, rhs = applyRobinScalar2D(mesh, lhs, rhs, elem_all.up, "up",
                                        quad_rule[1], α, op_val)
            end
            if bound_cond[i].side=="Left"
                lhs, rhs = applyRobinScalar2D(mesh, lhs, rhs, elem_all.left, "left",
                                        quad_rule[2], α, op_val)
            end
            if bound_cond[i].side=="Right"
                lhs, rhs = applyRobinScalar2D(mesh, lhs, rhs, elem_all.right, "right",
                                        quad_rule[2], α, op_val)
            end
        else
            error("Boundary condition type not implemented")
        end
    end
    unique!(bcdof)
    bcval = zero(bcdof)
    rhs = rhs - lhs[:,bcdof]*bcval
    rhs[bcdof] = bcval
    lhs[bcdof, :] .= 0.
    lhs[:, bcdof] .= 0.
    lhs[bcdof, bcdof] = sparse(I, length(bcdof), length(bcdof))
    return lhs, rhs
end

function compErrorNorm(mesh::Mesh, sol0, exactSol::Function, derivExactSol::Function,
                        a0::Function, gauss_rule::Array{GaussQuad,1})
    B_u, dB_u = bernsteinBasis(gauss_rule[1].nodes, mesh.degP[1])
    B_v, dB_v = bernsteinBasis(gauss_rule[2].nodes, mesh.degP[2])
    domainArea = 0

    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    basisCounter = 0
    Buv = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))
    dBdu = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))
    dBdv = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))

    for j=1:degP[2]+1
        for i=1:degP[1]+1
            basisCounter += 1
            Buv[:,:,basisCounter] = B_u[:,i]*B_v[:,j]'
            dBdu[:,:,basisCounter] = dB_u[:,i]*B_v[:,j]'
            dBdv[:,:,basisCounter] = B_u[:,i]*dB_v[:,j]'
        end
    end
    l2NormErr = 0.
    l2NormSol = 0.
    h1NormErr = 0.
    h1NormSol = 0.
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(numNodes)
        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)

                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2
                #compute the derivatives w.r.t. the physical space
                dR = [dRdu'; dRdv']

                dxdxi = dR * cpts'
                dR = dxdxi\dR
                Jac_par_phys = det(dxdxi)

                RR /= w_sum
                phys_pt = cpts*RR
                solVal = RR'*sol0[curNodes]
                dXsolVal = dR*sol0[curNodes]
                exSolVal = real(exactSol(phys_pt[1], phys_pt[2]))
                dXexSolVal = real(derivExactSol(phys_pt[1], phys_pt[2]))
                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                domainArea += localArea
                l2NormErr += localArea * (exSolVal-solVal)^2
                l2NormSol += localArea*exSolVal^2
                funValue = a0(phys_pt[1], phys_pt[2])
                h1NormErr += localArea * funValue *sum((dXexSolVal.-dXsolVal).^2)
                h1NormSol += localArea * funValue *sum((dXexSolVal).^2)
            end
        end
    end
    relL2Err = sqrt(l2NormErr/l2NormSol)
    relH1Err = sqrt(h1NormErr/h1NormSol)
    return real(relL2Err), real(relH1Err)
end

function compErrorNormElast(mesh::Mesh, Cmat, sol0, exactSolDisp::Function, exactSolStress::Function,
                         gauss_rule::Array{GaussQuad,1})
    invC = inv(Cmat)
    B_u, dB_u = bernsteinBasis(gauss_rule[1].nodes, mesh.degP[1])
    B_v, dB_v = bernsteinBasis(gauss_rule[2].nodes, mesh.degP[2])
    domainArea = 0

    numGaussU = length(gauss_rule[1].nodes)
    numGaussV = length(gauss_rule[2].nodes)
    basisCounter = 0
    Buv = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))
    dBdu = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))
    dBdv = zeros(numGaussU, numGaussV, (degP[1]+1)*(degP[2]+1))

    for j=1:degP[2]+1
        for i=1:degP[1]+1
            basisCounter += 1
            Buv[:,:,basisCounter] = B_u[:,i]*B_v[:,j]'
            dBdu[:,:,basisCounter] = dB_u[:,i]*B_v[:,j]'
            dBdv[:,:,basisCounter] = B_u[:,i]*dB_v[:,j]'
        end
    end
    l2NormErr = 0.
    l2NormSol = 0.
    h1NormErr = 0.
    h1NormSol = 0.
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 3]
        vMin = mesh.elemVertex[iElem, 2]
        vMax = mesh.elemVertex[iElem, 4]
        Jac_ref_par = (uMax-uMin)*(vMax-vMin)/4

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        curNodesXY = reshape(hcat(2*curNodes.-1, 2*curNodes)', 2*numNodes)
        cpts = mesh.controlPoints[1:2, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(numNodes)
        for jGauss = 1:numGaussV
            for iGauss = 1:numGaussU
                #compute the (B-)spline basis functions and derivatives with Bezier extraction
                N_mat = mesh.C[iElem] * Buv[iGauss, jGauss, :]
                dN_du = mesh.C[iElem] * dBdu[iGauss, jGauss, :] * 2/(uMax-uMin)
                dN_dv = mesh.C[iElem] * dBdv[iGauss, jGauss, :] * 2/(vMax-vMin)

                #compute the rational basis
                RR = N_mat.* wgts
                dRdu = dN_du.* wgts
                dRdv = dN_dv.* wgts
                w_sum = sum(RR)
                dw_xi = sum(dRdu)
                dw_eta = sum(dRdv)

                dRdu = dRdu/w_sum - RR*dw_xi/w_sum^2
                dRdv = dRdv/w_sum - RR*dw_eta/w_sum^2
                #compute the derivatives w.r.t to the physical space
                dR = [dRdu'; dRdv']

                dxdxi = dR * cpts'
                if abs(det(dxdxi))<1e-12
                    @warn "Singularity in mapping at $phys_pt"
                    dR = pinv(dxdxi)*dR
                else
                    dR = dxdxi\dR
                end
                Jac_par_phys = det(dxdxi)
                RR /= w_sum
                phys_pt = cpts*RR

                B = zeros(2*numNodes,3);
                B[1:2:2*numNodes-1,1] = dR[1,:]
                B[2:2:2*numNodes,2] = dR[2,:]
                B[1:2:2*numNodes-1,3] = dR[2,:]
                B[2:2:2*numNodes,3] = dR[1,:]
                solValX = RR'*sol0[2*curNodes.-1]
                solValY = RR'*sol0[2*curNodes]
                stressVect = Cmat*B'*sol0[curNodesXY]

                exSolVal = exactSolDisp(phys_pt[1], phys_pt[2])
                exStressVal = exactSolStress(phys_pt[1], phys_pt[2])
                localArea = Jac_par_phys * Jac_ref_par *
                    gauss_rule[1].weights[iGauss] * gauss_rule[2].weights[jGauss]
                domainArea += localArea
                l2NormErr += localArea * ((exSolVal[1]-solValX)^2 + (exSolVal[2]-solValY)^2)
                l2NormSol += localArea * (exSolVal[1]^2 + exSolVal[2]^2)

                h1NormErr += localArea * (exStressVal-stressVect)'*invC*(exStressVal-stressVect)
                h1NormSol += localArea * exStressVal'*invC*exStressVal
            end
        end
    end
    relL2Err = sqrt(l2NormErr/l2NormSol)
    relH1Err = sqrt(h1NormErr/h1NormSol)
    return real(relL2Err), real(relH1Err)
end


text:
@testset "bellman_ford" begin

    @testset "small_network" begin
        file = "data/networks/distance/weighted_elist_small.txt"
        graph = read_graph(file, "distance")
        distances, pred = bellman_ford(graph, 1)
        @test distances[5] == 36.0
        @test distances[10] == 19.0
    end

    @testset "large_network" begin
        file = "data/networks/distance/weighted_elist_large.txt"
        graph = read_graph(file, "distance")
        distances, pred = bellman_ford(graph, 1)
        @test distances[130] == 49.0
        @test distances[180] == 25.0
        @test distances[200] == 46.0
    end

    @testset "vlarge_network" begin
        file = "data/networks/distance/weighted_elist_vlarge.txt"
        graph = read_graph(file, "distance")
        distances, pred = bellman_ford(graph, 1)
        @test distances[250] == 19.0
        @test distances[1000] == 20.0
        @test distances[1500] == 16.0
    end

end


text:
using AlgebraOfGraphics, Makie, Test

@testset "utils" begin
    v1 = [1, 2, 7, 11]
    v2 = [1, 3, 4, 5.1]
    @test AlgebraOfGraphics.mergesorted(v1, v2) == [1, 2, 3, 4, 5.1, 7, 11]
    @test AlgebraOfGraphics.mergesorted(v2, v1) == [1, 2, 3, 4, 5.1, 7, 11]
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([2, 1], [1, 3])
    @test_throws ArgumentError AlgebraOfGraphics.mergesorted([1, 2], [10, 3])

    e1 = (-3, 11)
    e2 = (-5, 10)
    @test AlgebraOfGraphics.extend_extrema(e1, e2) == (-5, 11)

    @test AlgebraOfGraphics.midpoints(1:10) == 1.5:9.5
end

@testset "layers" begin
    df = (x = rand(1000), y = rand(1000), c = rand(["a", "b", "c"], 1000))
    d = mapping(:x, :y, color=:c)
    s = visual(color=:red) + mapping(markersize=:c)
    layers = data(df) * d * s
    @test layers[1].transformations[1] isa AlgebraOfGraphics.Visual
    @test layers[1].transformations[1].attributes[:color] == :red
    @test layers[1].positional == (:x, :y)
    @test layers[1].named == (color=:c,)
    @test layers[1].data == df
end

@testset "process_mappings" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), c=rand(["a", "b", "c"], 1000))
    d = mapping(:x => exp, [:y, :z], color=:c, marker = dims(1) => renamer(["a", "b"]))
    layer = data(df) * d
    entry = AlgebraOfGraphics.process_mappings(layer)
    @test entry.positional[1] == fill(map(exp, df.x))
    @test entry.positional[2] == [df.y, df.z]
    @test entry.primary[:color] == fill(df.c)
    @test entry.primary[:marker] == [fill("a"), fill("b")]
    @test entry.named == (;)
    @test entry.labels[1] == fill("x")
    @test entry.labels[2] ==  ["y", "z"]
    @test entry.labels[:color] == fill("c")
    @test entry.labels[:marker] == ""
end

@testset "grouping" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), w=rand(1000), c=rand(["a", "b", "c"], 1000))
    df.c[1:3] .= ["a", "b", "c"] # ensure all three values exist
    d = mapping(:x => exp, [:y, :z], color=:c, marker=dims(1) => t -> ["1", "2"][t], markersize=:w)
    layer = data(df) * d
    e = AlgebraOfGraphics.to_entry(layer)
    entries = map(CartesianIndices(AlgebraOfGraphics.shape(e))) do c
        primary, positional, named = map((e.primary, e.positional, e.named)) do tup
            return map(v -> v[c], tup)
        end
        labels = copy(e.labels)
        map!(l -> AlgebraOfGraphics.getnewindex(l, c), values(labels))
        return Entry(e; primary, positional, named, labels)
    end
    @test length(entries) == 6
    for i in 1: 6
        @test entries[i].plottype === Any
        @test isempty(entries[i].attributes)
        @test entries[i].labels[1] == "x"
        @test entries[i].labels[2] == (i ≤ 3 ? "y" : "z")
        @test entries[i].labels[:color] == "c"
        @test entries[i].labels[:marker] == ""
        @test entries[i].labels[:markersize] == "w"
    end

    @test entries[1].primary[:color] == "a"
    @test entries[1].primary[:marker] == "1"
    @test entries[1].positional[1] == exp.(df.x[df.c .== "a"])
    @test entries[1].positional[2] == df.y[df.c .== "a"]
    @test entries[1].named[:markersize] == df.w[df.c .== "a"]

    @test entries[2].primary[:color] == "b"
    @test entries[2].primary[:marker] == "1"
    @test entries[2].positional[1] == exp.(df.x[df.c .== "b"])
    @test entries[2].positional[2] == df.y[df.c .== "b"]
    @test entries[2].named[:markersize] == df.w[df.c .== "b"]

    @test entries[3].primary[:color] == "c"
    @test entries[3].primary[:marker] == "1"
    @test entries[3].positional[1] == exp.(df.x[df.c .== "c"])
    @test entries[3].positional[2] == df.y[df.c .== "c"]
    @test entries[3].named[:markersize] == df.w[df.c .== "c"]

    @test entries[4].primary[:color] == "a"
    @test entries[4].primary[:marker] == "2"
    @test entries[4].positional[1] == exp.(df.x[df.c .== "a"])
    @test entries[4].positional[2] == df.z[df.c .== "a"]
    @test entries[4].named[:markersize] == df.w[df.c .== "a"]

    @test entries[5].primary[:color] == "b"
    @test entries[5].primary[:marker] == "2"
    @test entries[5].positional[1] == exp.(df.x[df.c .== "b"])
    @test entries[5].positional[2] == df.z[df.c .== "b"]
    @test entries[5].named[:markersize] == df.w[df.c .== "b"]

    @test entries[6].primary[:color] == "c"
    @test entries[6].primary[:marker] == "2"
    @test entries[6].positional[1] == exp.(df.x[df.c .== "c"])
    @test entries[6].positional[2] == df.z[df.c .== "c"]
    @test entries[6].named[:markersize] == df.w[df.c .== "c"]
end

@testset "helpers" begin
    r = renamer("a" => "A", "b" => "B", "c" => "C")
    @test r("a") == AlgebraOfGraphics.Sorted(1, "A")
    @test r("b") == AlgebraOfGraphics.Sorted(2, "B")
    @test r("c") == AlgebraOfGraphics.Sorted(3, "C")
    @test_throws KeyError r("d")
    @test string(r("a")) == "A"
    @test string(r("b")) == "B"
    @test string(r("c")) == "C"
    @test r("a") < r("b") < r("c")
    @test r("a") == r("a")
    @test r("a") != r("b")

    s = sorter("b", "c", "a")
    @test s("a") == AlgebraOfGraphics.Sorted(3, "a")
    @test s("b") == AlgebraOfGraphics.Sorted(1, "b")
    @test s("c") == AlgebraOfGraphics.Sorted(2, "c")
    @test_throws KeyError s("d")
    @test string(s("a")) == "a"
    @test string(s("b")) == "b"
    @test string(s("c")) == "c"
    @test s("b") < s("c") < s("a")
    @test s("a") == s("a")
    @test s("a") != s("b")

    a = AlgebraOfGraphics.Sorted(1, [1, 2])
    b = AlgebraOfGraphics.Sorted(1, [1, 2])
    c = AlgebraOfGraphics.Sorted(1, [1, 3])
    @test a == b
    @test hash(a) == hash(b)
    @test a != c
    @test hash(a) != hash(c)

    @test string(nonnumeric(1)) == "1"
    @test isless(nonnumeric(1), nonnumeric(2))
end

@testset "legend_merging" begin
    mwe_data = (; t = 1:20, y1 = ones(20), y2 = rand(20))
    mwe_names = ["f_n", "f_d"]
    plt = data(mwe_data) *
        mapping("t", ["y1", "y2"] .=> "y"; color = dims(1) => (i -> mwe_names[i]), marker=dims(1)) *
        (visual(Lines) + visual(Scatter))
    @test_throws ArgumentError draw(plt)
end


text:
function _pathway_to(id::I,
                     network::TransmissionNetwork;
                     depth::Real=Inf) where {I <: Integer}
  path = I[]
  if network.external[id] | any(network.internal[:, id])
    push!(path, id)
    next_id = findfirst(network.internal[:, id])
    while (length(path) <= depth) & (typeof(next_id) != Nothing)
      push!(path, next_id)
      next_id = findfirst(network.internal[:, path[end]])
    end
  end
  @debug "_pathway_to: transmission pathway to $id: $path"
  return path
end


text:
#module LoadTestDataIntoJulia
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..BasicStructs
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates
using ..MainOverlap, ..RandIndex , ..ProbabilisticMetrics , ..VolumeMetric ,..InformationTheorhetic
using ..CUDAAtomicUtils, ..TpfpfnKernel, ..InterClassCorrKernel,..MeansMahalinobis
using ..MainAbstractions

using Conda
using PyCall
using Pkg

Conda.pip_interop(true)
Conda.pip("install", "SimpleITK")
Conda.pip("install", "pandas")
Conda.pip("install", "pymia")

sitk = pyimport("SimpleITK")
pym = pyimport("pymia")
pymMetr = pyimport("pymia.evaluation.metric")
pymEval = pyimport("pymia.evaluation.evaluator")
pymWrite = pyimport("pymia.evaluation.writer")
np= pyimport("numpy")

data_dir = "C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\exampleForTestsData"
result_file = "C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results.csv"
result_summary_file = "C:\\GitHub\\GitHub\\NuclearMedEval\\test\\data\\pymiaOutput\\results_summary.csv"

#given directory it gives all mhd file names concateneted with path - to get full file path and second in subarray will be file name
function getListOfExampleDatasFromFolder(folderPath::String) ::Vector{Vector{AbstractString}}
    return readdir(folderPath) |>
    (arr)-> filter((str)-> occursin("Subject",str), arr) |>
    (arr)-> map(str-> [split(str,".")[1], joinpath(folderPath,str), "$(joinpath(joinpath(folderPath,str),str))_GT.mha" ],arr)
end

function getDataAndEvaluationFromPymia(examplemhaDat)
    examplemha= examplemhaDat[3]
    ground_truth = sitk.ReadImage(examplemha)
    prediction = ground_truth

   py"""
    import SimpleITK as sitk
    def getPred(prediction, numb):
        return sitk.BinaryErode(prediction, [1,1,1], sitk.sitkBall, 0,numb)
    """
    prediction= py"getPred"(prediction,1)
    prediction= py"getPred"(prediction,2)

    labels = Dict([(1,"WHITEMATTER"),(2,"GREYMATTER") ])

    # metrics = [pymMetr.DiceCoefficient()
    #             , pymMetr.DiceCoefficient()
    #             , pymMetr.JaccardCoefficient()
    #             , pymMetr.GlobalConsistencyError()
    #             , pymMetr.AdjustedRandIndex()
    #             , pymMetr.CohenKappaCoefficient()
    #             , pymMetr.MutualInformation()
    #             , pymMetr.VariationOfInformation()
    #             , pymMetr.VolumeSimilarity()
    #             ,pymMetr.TruePositive()
    #             ,pymMetr.TrueNegative()
    #             ,pymMetr.FalsePositive()
    #             ,pymMetr.FalsePositive()
    # pymMetr.CohenKappaCoefficient()
    #             ,pymMetr.AdjustedRandIndex()
    #             ,pymMetr.InterclassCorrelation()

    #             ]
    metrics = [pymMetr.MahalanobisDistance()                ]

    evaluator = pymEval.SegmentationEvaluator(metrics, labels)

    evaluator.evaluate(prediction, ground_truth, examplemhaDat[1])

    pymWrite.ConsoleWriter().write(evaluator.results)

    goldS = np.array(sitk.GetArrayViewFromImage(ground_truth))
    segmAlgo = np.array(sitk.GetArrayViewFromImage(prediction))

    return (goldS,segmAlgo )
end#getDataAndEvaluationFromPymia

exampleFiles = getListOfExampleDatasFromFolder(data_dir)
examplemhaDat = exampleFiles[2]

goldS,segmAlgo =getDataAndEvaluationFromPymia(examplemhaDat);

arrGold = CuArray(goldS)
arrAlgo = CuArray(segmAlgo)
sizz= size(goldS)
##### load tests ...
conf = ConfigurtationStruct(false,trues(11)...)
numberToLooFor = UInt8(1)

preparedDict=MainAbstractions.prepareMetrics(conf)
res = calcMetricGlobal(preparedDict,conf,arrGold,arrAlgo,numberToLooFor)
res
#numbers below taken from pymia

@test isapprox(res.dice,0.654; atol = 0.1) #4) dice
@test isapprox(res.jaccard,0.486; atol = 0.1) #5) jaccard
@test isapprox(res.gce ,0.000; atol = 0.1) #6) gce
@test isapprox(res.randInd,0.618699; atol = 0.1) #7) randInd
@test isapprox(res.kc,0.640; atol = 0.1) #8) cohen kappa
@test isapprox(res.vol,0.654; atol = 0.1) #9) volume metric
@test isapprox(res.mi,0.130; atol = 0.1) #10) mutual information
@test isapprox(res.vi,0.256; atol = 0.1) #11) variation of information
@test isapprox(res.ic,0.6381813122385622; atol = 0.1)# interclas correlation
@test isapprox(res.md,0.08;atol = 0.1 ) #Mahalinobis

# ################# icc
# argsMain, threads,blocks, totalNumbOfVoxels=InterClassCorrKernel.prepareInterClassCorrKernel(arrGold ,arrAlgo,numberToLooFor)
# globalICC= InterClassCorrKernel.calculateInterclassCorr(arrGold,arrAlgo,threads,blocks,argsMain)

# # ################ Mahalinobis
# using ..MeansMahalinobis
# args,threads ,blocks= MeansMahalinobis.prepareMahalinobisKernel()
# mahalanobisResGlob=  MeansMahalinobis.calculateMalahlinobisDistance(arrGold,arrAlgo,args,threads ,blocks,1)
# goldS3d= CuArray(goldS);
# segmS3d= CuArray(segmAlgo);
# #we will fill it after we work with launch configuration
# loopXdim = UInt32(1);loopYdim = UInt32(1) ;loopZdim = UInt32(1) ;
# sizz = size(goldS3d);maxX = UInt32(sizz[1]);maxY = UInt32(sizz[2]);maxZ = UInt32(sizz[3])
# #gold
# totalXGold= CuArray([0.0]);
# totalYGold= CuArray([0.0]);
# totalZGold= CuArray([0.0]);
# totalCountGold= CuArray([0]);
# #segm
# totalXSegm= CuArray([0.0]);
# totalYSegm= CuArray([0.0]);
# totalZSegm= CuArray([0.0]);

# varianceXGlobalGold,covarianceXYGlobalGold,covarianceXZGlobalGold,varianceYGlobalGold,covarianceYZGlobalGold,varianceZGlobalGold= CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]);
# varianceXGlobalSegm,covarianceXYGlobalSegm,covarianceXZGlobalSegm,varianceYGlobalSegm,covarianceYZGlobalSegm,varianceZGlobalSegm= CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]);

# totalCountSegm= CuArray([0]);
# totalCountGold
# #countPerZGold= CUDA.zeros(Float32,sizz[3]+1);
# #countPerZSegm= CUDA.zeros(Float32,sizz[3]+1);
# countPerZGold= CUDA.zeros(Float32,500);
# countPerZSegm= CUDA.zeros(Float32,500);

# #covariancesSliceWise= CUDA.zeros(Float32,12,sizz[3]+1);
# covariancesSliceWiseGold= CUDA.zeros(Float32,6,500);
# covariancesSliceWiseSegm= CUDA.zeros(Float32,6,500);
# covarianceGlobal= CUDA.zeros(Float32,12,1);

# mahalanobisResGlobal= CUDA.zeros(1);
# mahalanobisResSliceWise= CUDA.zeros(500);
# #mahalanobisResSliceWise= CUDA.zeros(sizz[3]);

# args = (goldS3d,segmS3d,numberToLooFor
# ,loopYdim,loopXdim,loopZdim
# ,(maxX, maxY,maxZ)
# ,totalXGold,totalYGold,totalZGold,totalCountGold
# ,totalXSegm,totalYSegm,totalZSegm,totalCountSegm,countPerZGold
# , countPerZSegm,covariancesSliceWiseGold, covariancesSliceWiseSegm,
# varianceXGlobalGold,covarianceXYGlobalGold,covarianceXZGlobalGold,varianceYGlobalGold,covarianceYZGlobalGold,varianceZGlobalGold
#     ,varianceXGlobalSegm,covarianceXYGlobalSegm,covarianceXZGlobalSegm,varianceYGlobalSegm,covarianceYZGlobalSegm,varianceZGlobalSegm
#     ,mahalanobisResGlobal, mahalanobisResSliceWise)

#     # calculate the amount of dynamic shared memory for a 2D block size
#     get_shmem(threads) = (sizeof(UInt32)*3*4)

#     function get_threads(threads)
#         threads_x = 32
#         threads_y = cld(threads,threads_x )
#         return (threads_x, threads_y)
#     end

#     kernel = @cuda launch=false MeansMahalinobis.meansMahalinobisKernel(args...)

#     config = launch_configuration(kernel.fun, shmem=threads->get_shmem(get_threads(threads)))

#    # convert to 2D block size and figure out appropriate grid size
#     threads = get_threads(config.threads)
#     blocks = UInt32(config.blocks)
#     loopXdim = UInt32(cld(maxX, threads[1]))
#     loopYdim = UInt32(cld(maxY, threads[2]))
#     loopZdim = UInt32(cld(maxZ,blocks ))

# #covariancesSliceWise= CUDA.zeros(Float32,12,sizz[3]+1);
# covariancesSliceWiseGold= CUDA.zeros(Float32,6,500);
# covariancesSliceWiseSegm= CUDA.zeros(Float32,6,500);

# #gold
# totalXGold= CuArray([0.0]);
# totalYGold= CuArray([0.0]);
# totalZGold= CuArray([0.0]);
# totalCountGold= CuArray([0]);
# #segm
# totalXSegm= CuArray([0.0]);
# totalYSegm= CuArray([0.0]);
# totalZSegm= CuArray([0.0]);
# totalCountSegm= CuArray([0]);

# varianceXGlobalGold,covarianceXYGlobalGold,covarianceXZGlobalGold,varianceYGlobalGold,covarianceYZGlobalGold,varianceZGlobalGold= CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]);
# varianceXGlobalSegm,covarianceXYGlobalSegm,covarianceXZGlobalSegm,varianceYGlobalSegm,covarianceYZGlobalSegm,varianceZGlobalSegm= CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]),CuArray([0.0]);

# args = (goldS3d,segmS3d,numberToLooFor
# ,loopYdim,loopXdim,loopZdim
# ,(maxX, maxY,maxZ)
# ,totalXGold,totalYGold,totalZGold,totalCountGold
# ,totalXSegm,totalYSegm,totalZSegm,totalCountSegm,countPerZGold
# , countPerZSegm,covariancesSliceWiseGold, covariancesSliceWiseSegm,
# varianceXGlobalGold,covarianceXYGlobalGold,covarianceXZGlobalGold,varianceYGlobalGold,covarianceYZGlobalGold,varianceZGlobalGold
#     ,varianceXGlobalSegm,covarianceXYGlobalSegm,covarianceXZGlobalSegm,varianceYGlobalSegm,covarianceYZGlobalSegm,varianceZGlobalSegm
#     ,mahalanobisResGlobal, mahalanobisResSliceWise)

#     @cuda cooperative=true threads=threads blocks=blocks MeansMahalinobis.meansMahalinobisKernel(args...)

#     @test isapprox(mahalanobisResGlobal[1],0.08;atol = 0.1 )

# ####################### some benchmarks

# BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
# BenchmarkTools.DEFAULT_PARAMETERS.seconds =60
# BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

# function toBench(FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn)
#     CUDA.@sync TpfpfnKernel.getTpfpfnData!(FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn
#                             ,sliceMetricsTupl
#                             ,metricsTuplGlobal
#                             ,sizz[1]*sizz[2]
#                             ,sizz[3]
#                             ,UInt8(1)
#                             ,conf
#                             ,totalNumberOfVoxels)
#                         end

# bb2 = @benchmark toBench(FlattGoldGPU,FlattSegGPU,tp,tn,fp,fn)  setup=(goldBoolGPU,segmBoolGPU,tp,tn,fp,fn, tpArr,tnArr,fpArr, fnArr, blockNum , nx,ny,nz ,tpTotalTrue,tnTotalTrue,fpTotalTrue, fnTotalTrue ,tpPerSliceTrue,  tnPerSliceTrue,fpPerSliceTrue,fnPerSliceTrue ,flattG, flattSeg ,FlattGoldGPU,FlattSegGPU,intermediateResTp,intermediateResFp,intermediateResFn = getSmallTestBools())

# CUDA.@profile begin
#     TpfpfnKernel.getTpfpfnData!(arrGold,arrAlgo,tp,tn,fp,fn
#     , intermediateResTp,intermediateResFp
#     ,intermediateResFn,sizz[1]*sizz[2],sizz[3]
#     ,UInt8(1)
#     ,IndexesArray)
# end

# ## stored indexes  and now we are inspecting it
# diff = Int64(round(length(goldS)- maximum(IndexesArray)))

# diff/pixelNumberPerSlice
# cpuZZ = zeros(Int32,10000000)
# copyto!(cpuZZ, zz)
# zz = filter(it->it>0 ,cpuZZ)

# diffrenceOfIndexes = filter(pair->pair[1]!=pair[2], collect(enumerate(zz))  )
# length(zz) - length(goldS)

# sort(zz)

# pixelNumberPerSlice-512*77

# # Subject_2  GREYMATTER   0.298  10.863     74392.000   0.298
# # Subject_2  WHITEMATTER  0.654  6.000      206422.000  0.654

# #tpTotalTrue = filter(pair->pair[2]== vec(goldS)[pair[1]] ==1 ,collect(enumerate(vec(segmAlgo))))|>length

# arrOnesA = CUDA.ones(sizz);
# arrOnesB = CUDA.ones(sizz);

# goldBoolGPU,segmBoolGPU,tp,tn,fp,fn, tpArr,tnArr,fpArr, fnArr, blockNum , nx,ny,nz ,tpTotalTrue,tnTotalTrue,fpTotalTrue, fnTotalTrue ,tpPerSliceTrue,  tnPerSliceTrue,fpPerSliceTrue,fnPerSliceTrue ,flattG, flattSeg ,FlattGoldGPU,FlattSegGPU,intermediateResTp,intermediateResFp,intermediateResFn = getSmallTestBools();
# IndexesArray= CUDA.zeros(Int32,10000000)

# TpfpfnKernel.getTpfpfnData!(arrOnesA,arrOnesB,tp,tn,fp,fn, intermediateResTp,intermediateResFp,intermediateResFn,sizz[1]*sizz[2],sizz[3],Float32(1),IndexesArray)
# # sum(IndexesArray)
# tp[1]
# tp[1] -length(arrOnesA)
# # length(arrOnesA) -sum(IndexesArray)
# #length(arrOnesA) - tp[1]

# tp = 5
# fp = 5
# fn = 5
# tn = 5

# using ..RandIndex

# MainOverlap.dice(tp,fp, fn)
# MainOverlap.jaccard(tp,fp, fn)
# MainOverlap.gce(tn,tp,fp, fn)
# RandIndex.calculateAdjustedRandIndex(tn,tp,fp, fn)
# ProbabilisticMetrics.calculateCohenCappa(tn,tp,fp, fn )
# VolumeMetric.getVolumMetric(tp,fp, fn )
# InformationTheorhetic.mutualInformationMetr(tn,tp,fp, fn)
# InformationTheorhetic.variationOfInformation(tn,tp,fp, fn)

# rand(1,3,7)[19]

# end #module


text:
school = GroupedErrors.exampletable("school.csv") |> loadtable
tables = []

#Test scatter plot
processed_table = @> school begin
    @splitby _.Sx
    @across _.School
    @x _.MAch
    @y _.SSS
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t1"))

#Test different estimator
processed_table = @> school begin
    @splitby _.Sx
    @across _.School
    @x _.MAch median
    @y _.SSS median
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t2"))

#Test xy comparison
processed_table = @> school begin
    @splitby _.Minrty
    @across _.School
    @xy _.SSS
    @compare _.Sx
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t3"))

#Test cumulative
processed_table = @> school begin
    @splitby _.Sx
    @across _.School
    @x _.MAch
    @y :cumulative
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t4"))

#Test density
processed_table = @> school begin
    @splitby _.Minrty
    @across _.School
    @x _.CSES
    @y :density bandwidth = 0.2
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t5"))

#Test locreg
processed_table = @> school begin
    @splitby _.Minrty
    @across :all
    @x _.Sx :discrete
    @y _.MAch
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t6"))

#Test binning
processed_table = @> school begin
    @splitby _.Minrty
    @x _.MAch :binned 40
    @y :density
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t7"))

#test bootstrap
processed_table = @> school begin
    @splitby _.Minrty
    @bootstrap 500
    @x _.CSES
    @y :density bandwidth = 0.2
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t8"))

#test continuous locreg
processed_table = @> school begin
    @splitby _.Minrty
    @x _.MAch :continuous
    @y _.SSS
    ProcessedTable
end

push!(tables, processed_table.table)
#Dagger.save(table(processed_table.table, chunks = 1), joinpath(@__DIR__, "tables", "t9"))


text:
function run_with_timeout(
             command; log=stdout, timeout = 10*60, name = "",
             wait_time = 1, verbose = true, kill_timeout = 60, ondone=identity
         )
    out_io = IOBuffer()
    err_io = IOBuffer()
    out_file, err_file = "", ""
    if VERSION < v"1.1"
        out_file, out_io = mktemp()
        err_file, err_io = mktemp()
    end
    pipe = pipeline(command, stdout = out_io, stderr = err_io)
    process = run(pipe, wait = false)
    if VERSION < v"1.1"
        out_io = open(out_file)
        err_io = open(err_file)
    end
    timeout_start = time()
    task = @async begin
        logfallback = false
        io = try
            log isa String ? open(log, "w") : log
        catch err
            @error "Error opening logfile, falling back to stdout" error=err
            logfallback = true
            stdout
        end
        try
            tstart = time()
            verbose && @info("starting $name")
            while process_running(process)
                elapsed = (time() - timeout_start)
                if elapsed > timeout
                    verbose && @info("Terminating $name")
                    kill(process)
                    # Handle scenarios where SIGTERM is blocked/ignored/handled by the process
                    start_time = time()
                    while process_running(process)
                        if time() - start_time > kill_timeout
                            verbose && @info("Killing $name")
                            kill(process, signum = SIGKILL)
                        end
                        sleep(5)
                    end
                    break
                end
                errstr, outstr = readstr_buffer.((out_io, err_io))
                is_silent = length(errstr) == 0 && length(outstr) == 0
                isempty(outstr) || println(io, outstr)
                isempty(errstr) || println(io, errstr)

                # if something printed reset timeout
                if !is_silent
                    timeout_start = time()
                end

                sleep(wait_time)
            end
            ondone()
            verbose && @info("$name completed in $(round(time() - tstart, digits=1)) seconds")
        catch err
            @error "Error while running $(name) with timeout." error=err
        finally
            errstr, outstr = readstr_buffer.((out_io, err_io))
            isempty(outstr) || println(io, outstr)
            isempty(errstr) || println(io, errstr)

            flush(io)
            if log isa String && !logfallback
                close(io)
            end
        end
    end
    return process, task
end

function run_queued(f, arr; sleeptime = 1, processes = 3, ondone=identity)
    process_queue = []

    for (i, el) in enumerate(arr)
        while length(process_queue) >= processes
            filter!(process_running, process_queue)
            sleep(sleeptime)
        end
        args = f(el)
        process, task = run_with_timeout(first(args); args[2]..., name=el, ondone = () -> ondone(i, el), verbose=false)
        push!(process_queue, process)
    end
    for proc in process_queue
        wait(proc)
    end
end

function readstr_buffer(x::IOStream)
    return read(x, String)
end

function readstr_buffer(x::Base.GenericIOBuffer{Array{UInt8,1}})
    return String(take!(x))
end