text:
struct StackedArea
    time_range::Array
    data_matrix::Matrix
    labels::Array

end

struct BarPlot
    time_range::Array
    bar_data::Matrix
    labels::Array

end

struct StackedGeneration
    time_range::Array
    data_matrix::Matrix
    labels::Array

end

struct BarGeneration
    time_range::Array
    bar_data::Matrix
    labels::Array

end

"""

get_stacked_plot_data(res::OperationModelResults, variable::String)

This function takes in results of struct OperationModelResult. It takes the
dataframe from whichever variable name was given and converts it to type StackedArea.
StackedArea is the type of struct that signals the plot() function to use the
StackedArea plot recipe method.

#Example

to make a single stack plot for the P_ThermalStandard variable:

P_ThermalStandard = get_stacked_plot_data(res, "P_ThermalStandard")
plot(P_ThermalStandard)

"""

function get_stacked_plot_data(res::OperationModelResults, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:,sort]
    end

    data_matrix = convert(Matrix, variable)
    labels = collect(names(variable))
    legend = string.(labels)

    return StackedArea(time_range, data_matrix, legend)

end

function get_bar_plot_data(res::OperationModelResults, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:,sort]
    end

    data = convert(Matrix, variable)
    bar_data = sum(data, dims = 1)
    labels = collect(names(variable))
    legend = string.(labels)

    return BarPlot(time_range, bar_data, legend)

end

function get_stacked_generation_data(res::OperationModelResults; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    key_name = collect(keys(res.variables))
    alphabetical = sort!(key_name)

    if !isnothing(sort)
        labels = sort
    else
        labels = alphabetical
    end

    variable = res.variables[Symbol(labels[1])]
    data_matrix = sum(convert(Matrix, variable), dims = 2)
    legend = string.(labels)

    for i in 1:length(labels)
        if i !== 1
            variable = res.variables[Symbol(labels[i])]
            data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
        end
    end

    return StackedGeneration(time_range, data_matrix, legend)

end

function get_bar_gen_data(res::OperationModelResults)

   time_range = res.time_stamp[!,:Range]
   key_name = collect(keys(res.variables))

   variable = res.variables[Symbol(key_name[1])]
   data_matrix = sum(convert(Matrix, variable), dims = 2)
   legend = string.(key_name)

    for i in 1:length(key_name)
       if i !== 1
           variable = res.variables[Symbol(key_name[i])]
           data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
       end
   end
   bar_data = sum(data_matrix, dims = 1)
   return BarGeneration(time_range, bar_data, legend)

end

"""
    sort_data(results::OperationModelResults)

This function takes in struct OperationModelResults,
sorts the generators in each variable, and outputs the sorted
results. The generic function sorts the generators alphabetically.

kwargs: 'Variables' to choose which variables to be sorted.

each variable has a kwarg, so that a specific order can
be generated, such that when plotted, the first generator is on the bottom.
if a list of generator names has fewer generators than the variable, only the
generators on the list will be outputted.

#Examples

example 1:
sorted_results = sort_data(res)

>sorted_results.variables[P_RenewableDispatch] will be in the order
    [:WindBusA :WindBusB :WindBusC] (alphabetical)

example 2:
my_order = [:WindBusC :WindBusB :WindBusA]
sorted_results = sort_data(res; P_RenewableDispatch = my_order)

>sorted_results.variables[P_RenewableDispatch] will be in the order
    [:WindBusC :WindBusB :WindBusA] (my_order)

example 3:
my_order = [:WindBusC :WindBusA]
sorted_results = sort_data(res; P_RenewableDispatch = my_order)

>sorted_results.variables[P_RenewableDispatch] will be in the order
    [:WindBusC :WindBusA] (my_order)

example 4:
my_variable_order = [:P_ThermalStandard :ON_ThermalStandard]
sorted_results = sort_data(res; Variables = my_variable_order)

>sorted_results.variables
    Dict{Symbol,DataFrames.DataFrame} with 2 entries:
    :P_ThermalStandard => 24×5 DataFrames.DataFrame…
    :ON_ThermalStandard => 24×5 DataFrames.DataFrame…

* note that only the generators included in 'my_order' will be in the
results, and consequently, only these will be plotted. This can be a nice
feature for variables with more than 5 generators.

"""
function sort_data(res::OperationModelResults; kwargs...)

    Variables = Dict()
    Variables[:P_ThermalStandard]  = get(kwargs, :P_ThermalStandard, nothing)
    Variables[:P_RenewableDispatch]  = get(kwargs, :P_RenewableDispatch, nothing)
    Variables[:START_ThermalStandard] = get(kwargs, :START_ThermalStandard, nothing)
    Variables[:STOP_ThermalStandard] = get(kwargs, :STOP_ThermalStandard, nothing)
    Variables[:ON_ThermalStandard] = get(kwargs, :ON_ThermalStandard, nothing)
    Variable_dict = get(kwargs, :Variables, nothing)

    key_name = collect(keys(res.variables))
    alphabetical = sort!(key_name)

    if !isnothing(Variable_dict)
        labels = Variable_dict
    else
        labels = alphabetical
    end

    variable_dict = Dict()

    for i in 1:length(labels)

          variable_dict[labels[i]] = res.variables[labels[i]]

    end

    for (k,v) in Variables, k in keys(variable_dict)

        variable = variable_dict[k]
        alphabetical = sort!(names(variable))
        order = Variables[k]

        if isnothing(order)
            variable = variable[:, alphabetical]
        else

            variable = variable[:, order]

        end
        variable_dict[k] = variable

    end

    res = OperationModelResults(variable_dict, res.total_cost, res.optimizer_log, res.time_stamp)

    return res
end


text:
@testset "warlpiri-division" begin

  #comprehensive negative distribution testset
  for a in WP, b in WP, c in WP, d in WP
      #println("testing $a → $b, $c → $d")
      @test -((a → b) / (c → d)) == ((-(a → b)) / (c → d))
      @test -((a → b) / (c → d)) == ((a → b) / (-(c → d)))
      @test ((a → b) / (c → d)) == ((-(a → b)) / (-(c → d)))
  end

  const warlpiri_div_matrix = [
   #W0_0 W0_1 W0_2 W0_3 W0_4 W0_5 W0_6 W0_7 W0_8   W1_1 W1_2 W1_3 W1_4 W1_5 W1_6 W1_7 W1_8   W2_2 W2_3 W2_4 W2_5 W2_6 W2_7 W2_8   W3_3 W3_4 W3_5 W3_6 W3_7 W3_8   W4_4 W4_5 W4_6 W4_7 W4_8   W5_5 W5_6 W5_7 W5_8   W6_6 W6_7 W6_8   W7_7 W7_8   W8_8
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_0 W0_0 W0_0 W0_0 W0_0 W0_0 W0_0 W0_0   W0_0 W0_0 W0_0 W0_0 W0_0 W0_0 W0_0   W0_0 W0_0 W0_0 W0_0 W0_0 W0_0   W0_0 W0_0 W0_0 W0_0 W0_0   W0_0 W0_0 W0_0 W0_0   W0_0 W0_0 W0_0   W0_0 W0_0   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_3 W0_3 W0_3 W0_3 W0_3 W0_3 W0_3   W0_3 W0_3 W0_3 W0_3 W0_3 W0_3   W0_1 W0_1 W0_1 W0_1 W0_1   W0_1 W0_1 W0_1 W0_1   W0_1 W0_1 W0_1   W0_1 W0_1   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_4 W0_4 W0_4 W0_4 W0_4 W0_4 W0_4   W0_3 W0_3 W0_3 W0_3 W0_3 W0_3   W0_2 W0_2 W0_2 W0_2 W0_2   W0_1 W0_1 W0_1 W0_1   W0_1 W0_1 W0_1   W0_1 W0_1   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_5 W0_5 W0_5 W0_5 W0_5 W0_5 W0_5   W0_5 W0_5 W0_5 W0_5 W0_5 W0_5   W0_3 W0_3 W0_3 W0_3 W0_3   W0_3 W0_3 W0_3 W0_3   W0_1 W0_1 W0_1   W0_1 W0_1   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_6 W0_6 W0_6 W0_6 W0_6 W0_6 W0_6   W0_5 W0_5 W0_5 W0_5 W0_5 W0_5   W0_4 W0_4 W0_4 W0_4 W0_4   W0_3 W0_3 W0_3 W0_3   W0_2 W0_2 W0_2   W0_1 W0_1   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_5 W0_5 W0_5 W0_5 W0_5   W0_5 W0_5 W0_5 W0_5   W0_3 W0_3 W0_3   W0_3 W0_3   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_6 W0_6 W0_6 W0_6 W0_6   W0_5 W0_5 W0_5 W0_5   W0_4 W0_4 W0_4   W0_3 W0_3   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7 W0_7   W0_7 W0_7 W0_7   W0_7 W0_7   W0_0 # * W0_0
    WALL WALL WALL WALL WALL WALL WALL WALL WALL   W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 WALL   W0_8 W0_8 W0_8 W0_8 W0_8 W0_8 WALL   W0_8 W0_8 W0_8 W0_8 W0_8 WALL   W0_8 W0_8 W0_8 W0_8 WALL   W0_8 W0_8 W0_8 WALL   W0_8 W0_8 WALL   W0_8 WALL   WALL # * W0_0

    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_3 W1_3 W1_3 W1_3 W1_3 W1_3 W0_3   W1_3 W1_3 W1_3 W1_3 W1_3 W0_3   W1_1 W1_1 W1_1 W1_1 W0_1   W1_1 W1_1 W1_1 W0_1   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # 8 W1_1
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_4 W1_4 W1_4 W1_4 W1_4 W1_4 W0_4   W1_3 W1_3 W1_3 W1_3 W1_3 W0_3   W1_2 W1_2 W1_2 W1_2 W0_2   W1_1 W1_1 W1_1 W0_1   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # 8 W1_2
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_5 W1_5 W1_5 W1_5 W1_5 W1_5 W0_5   W1_5 W1_5 W1_5 W1_5 W1_5 W0_5   W1_3 W1_3 W1_3 W1_3 W0_3   W1_3 W1_3 W1_3 W0_3   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # 8 W1_3
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_6 W1_6 W1_6 W1_6 W1_6 W1_6 W0_6   W1_5 W1_5 W1_5 W1_5 W1_5 W0_5   W1_4 W1_4 W1_4 W1_4 W0_4   W1_3 W1_3 W1_3 W0_3   W1_2 W1_2 W0_2   W1_1 W0_1   W0_0 # 8 W1_4
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_5 W1_5 W1_5 W1_5 W0_5   W1_5 W1_5 W1_5 W0_5   W1_3 W1_3 W0_3   W1_3 W0_3   W0_0 # 8 W1_5
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_6 W1_6 W1_6 W1_6 W0_6   W1_5 W1_5 W1_5 W0_5   W1_4 W1_4 W0_4   W1_3 W0_3   W0_0 # 8 W1_6
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W0_8   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W0_7   W1_7 W0_7   W0_0 # 8 W1_7
    W8_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 WALL   W1_8 W1_8 WALL   W1_8 WALL   WALL # 8 W1_8

    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_4 W3_4 W2_4 W1_4 W1_4 W1_4 W0_4   W3_3 W2_3 W1_3 W1_3 W1_3 W0_3   W2_2 W1_2 W1_2 W1_2 W0_2   W1_1 W1_1 W1_1 W0_1   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # * W2_2
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_5 W3_5 W2_5 W1_5 W1_5 W1_5 W0_5   W3_5 W2_5 W1_5 W1_5 W1_5 W0_5   W2_3 W1_3 W1_3 W1_3 W0_3   W1_3 W1_3 W1_3 W0_3   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # * W2_3
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_6 W3_6 W2_6 W1_6 W1_6 W1_6 W0_6   W3_5 W2_5 W1_5 W1_5 W1_5 W0_5   W2_4 W1_4 W1_4 W1_4 W0_4   W1_3 W1_3 W1_3 W0_3   W1_2 W1_2 W0_2   W1_1 W0_1   W0_0 # * W2_4
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W2_5 W1_5 W1_5 W1_5 W0_5   W1_5 W1_5 W1_5 W0_5   W1_3 W1_3 W0_3   W1_3 W0_3   W0_0 # * W2_5
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W2_6 W1_6 W1_6 W1_6 W0_6   W1_5 W1_5 W1_5 W0_5   W1_4 W1_4 W0_4   W1_3 W0_3   W0_0 # * W2_6
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 W0_8   W5_7 W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W4_7 W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W3_7 W2_7 W1_7 W1_7 W1_7 W0_7   W2_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W2_7
    W8_8 W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 WALL   W5_8 W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 WALL   W4_8 W3_8 W2_8 W1_8 W1_8 W1_8 WALL   W3_8 W2_8 W1_8 W1_8 W1_8 WALL   W2_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 WALL   W1_8 W1_8 WALL   W1_8 WALL   WALL # * W2_8

    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 W0_8   W5_7 W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W5_5 W3_5 W3_5 W1_5 W1_5 W1_5 W0_5   W3_5 W3_5 W1_5 W1_5 W1_5 W0_5   W3_3 W1_3 W1_3 W1_3 W0_3   W1_3 W1_3 W1_3 W0_3   W1_1 W1_1 W0_1   W1_1 W0_1   W0_0 # * W3_3
    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 W0_8   W5_7 W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W5_6 W3_6 W3_6 W1_6 W1_6 W1_6 W0_6   W3_5 W3_5 W1_5 W1_5 W1_5 W0_5   W3_4 W1_4 W1_4 W1_4 W0_4   W1_3 W1_3 W1_3 W0_3   W1_2 W1_2 W0_2   W1_1 W0_1   W0_0 # * W3_4
    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 W0_8   W5_7 W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_5 W1_5 W1_5 W1_5 W0_5   W1_5 W1_5 W1_5 W0_5   W1_3 W1_3 W0_3   W1_3 W0_3   W0_0 # * W3_5
    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 W0_8   W5_7 W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_6 W1_6 W1_6 W1_6 W0_6   W1_5 W1_5 W1_5 W0_5   W1_4 W1_4 W0_4   W1_3 W0_3   W0_0 # * W3_6
    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 W0_8   W5_7 W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W5_7 W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_7 W3_7 W1_7 W1_7 W1_7 W0_7   W3_7 W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W1_7 W0_7   W1_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W3_7
    W8_8 W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 WALL   W5_8 W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 WALL   W5_8 W3_8 W3_8 W1_8 W1_8 W1_8 WALL   W3_8 W3_8 W1_8 W1_8 W1_8 WALL   W3_8 W1_8 W1_8 W1_8 WALL   W1_8 W1_8 W1_8 WALL   W1_8 W1_8 WALL   W1_8 WALL   WALL # * W3_8

    W8_8 W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 W0_8   W7_7 W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W6_6 W5_6 W4_6 W3_6 W2_6 W1_6 W0_6   W5_5 W4_5 W3_5 W2_5 W1_5 W0_5   W4_4 W3_4 W2_4 W1_4 W0_4   W3_3 W2_3 W1_3 W0_3   W2_2 W1_2 W0_2   W1_1 W0_1   W0_0 # * W4_4
    W8_8 W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 W0_8   W7_7 W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W4_5 W3_5 W2_5 W1_5 W0_5   W3_5 W2_5 W1_5 W0_5   W2_3 W1_3 W0_3   W1_3 W0_3   W0_0 # * W4_5
    W8_8 W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 W0_8   W7_7 W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W4_6 W3_6 W2_6 W1_6 W0_6   W3_5 W2_5 W1_5 W0_5   W2_4 W1_4 W0_4   W1_3 W0_3   W0_0 # * W4_6
    W8_8 W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 W0_8   W7_7 W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W6_7 W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W5_7 W4_7 W3_7 W2_7 W1_7 W0_7   W4_7 W3_7 W2_7 W1_7 W0_7   W3_7 W2_7 W1_7 W0_7   W2_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W4_7
    W8_8 W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 WALL   W7_8 W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 WALL   W6_8 W5_8 W4_8 W3_8 W2_8 W1_8 WALL   W5_8 W4_8 W3_8 W2_8 W1_8 WALL   W4_8 W3_8 W2_8 W1_8 WALL   W3_8 W2_8 W1_8 WALL   W2_8 W1_8 WALL   W1_8 WALL   WALL # * W4_8

    W8_8 W7_8 W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 W0_8   W7_7 W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_5 W3_5 W3_5 W1_5 W0_5   W3_5 W3_5 W1_5 W0_5   W3_3 W1_3 W0_3   W1_3 W0_3   W0_0 # * W5_5
    W8_8 W7_8 W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 W0_8   W7_7 W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_6 W3_6 W3_6 W1_6 W0_6   W3_5 W3_5 W1_5 W0_5   W3_4 W1_4 W0_4   W1_3 W0_3   W0_0 # * W5_6
    W8_8 W7_8 W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 W0_8   W7_7 W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W7_7 W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_7 W5_7 W3_7 W3_7 W1_7 W0_7   W5_7 W3_7 W3_7 W1_7 W0_7   W3_7 W3_7 W1_7 W0_7   W3_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W5_7
    W8_8 W7_8 W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 WALL   W7_8 W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 WALL   W7_8 W5_8 W5_8 W3_8 W3_8 W1_8 WALL   W5_8 W5_8 W3_8 W3_8 W1_8 WALL   W5_8 W3_8 W3_8 W1_8 WALL   W3_8 W3_8 W1_8 WALL   W3_8 W1_8 WALL   W1_8 WALL   WALL # * W5_8

    W8_8 W7_8 W7_8 W7_8 W6_8 W5_8 W4_8 W1_8 W0_8   W7_7 W7_7 W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W7_7 W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W6_6 W5_6 W4_6 W1_6 W0_6   W5_5 W4_5 W1_5 W0_5   W4_4 W1_4 W0_4   W1_3 W0_3   W0_0 # * W6_6
    W8_8 W7_8 W7_8 W7_8 W6_8 W5_8 W4_8 W1_8 W0_8   W7_7 W7_7 W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W7_7 W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W7_7 W6_7 W5_7 W4_7 W1_7 W0_7   W6_7 W5_7 W4_7 W1_7 W0_7   W5_7 W4_7 W1_7 W0_7   W4_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W6_7
    W8_8 W7_8 W7_8 W7_8 W6_8 W5_8 W4_8 W1_8 WALL   W7_8 W7_8 W7_8 W6_8 W5_8 W4_8 W1_8 WALL   W7_8 W7_8 W6_8 W5_8 W4_8 W1_8 WALL   W7_8 W6_8 W5_8 W4_8 W1_8 WALL   W6_8 W5_8 W4_8 W1_8 WALL   W5_8 W4_8 W1_8 WALL   W4_8 W1_8 WALL   W1_8 WALL   WALL # * W6_8

    W8_8 W7_8 W7_8 W7_8 W7_8 W5_8 W5_8 W1_8 W0_8   W7_7 W7_7 W7_7 W7_7 W5_7 W5_7 W1_7 W0_7   W7_7 W7_7 W7_7 W5_7 W5_7 W1_7 W0_7   W7_7 W7_7 W5_7 W5_7 W1_7 W0_7   W7_7 W5_7 W5_7 W1_7 W0_7   W5_7 W5_7 W1_7 W0_7   W5_7 W1_7 W0_7   W1_7 W0_7   W0_0 # * W7_7
    W8_8 W7_8 W7_8 W7_8 W7_8 W5_8 W5_8 W1_8 WALL   W7_8 W7_8 W7_8 W7_8 W5_8 W5_8 W1_8 WALL   W7_8 W7_8 W7_8 W5_8 W5_8 W1_8 WALL   W7_8 W7_8 W5_8 W5_8 W1_8 WALL   W7_8 W5_8 W5_8 W1_8 WALL   W5_8 W5_8 W1_8 WALL   W5_8 W1_8 WALL   W1_8 WALL   WALL # * W7_8

    W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 W8_8 WALL   W8_8 W8_8 W8_8 WALL   W8_8 W8_8 WALL   W8_8 WALL   WALL # * W8_8
    ]

    warlpiri_div_test_matrix = [a / b for a in pos_warlpiri, b in pos_warlpiri]
    warlpiri_label_matrix = [(a, b) for a in pos_warlpiri, b in pos_warlpiri]

    for (expected, evaluated, l) in zip(warlpiri_div_matrix, warlpiri_div_test_matrix, warlpiri_label_matrix)
        #println(l)
        @test (expected, l) == (evaluated, l)
    end

    #next, test values when dividing by something which goes "round the
    #zero"... We'll use the interval (-1, 1], which is representable as 0b1101
    #-> 0b1100

    WRTZ = WP1101 → WP0100
    warlpiri_div_rtz_values = [
    WP1000 → WP0111, #WP0000 → WP0000   0 / (-1, 1] = all reals
    WP1000 → WP0111, #WP0000 → WP0001
    WP1000 → WP0111, #WP0000 → WP0010
    WP1000 → WP0111, #WP0000 → WP0011
    WP1000 → WP0111, #WP0000 → WP0100
    WP1000 → WP0111, #WP0000 → WP0101
    WP1000 → WP0111, #WP0000 → WP0110
    WP1000 → WP0111, #WP0000 → WP0111
    WP1000 → WP0111, #WP0000 → WP1000

    WP0001 → WP1111, #WP0001 → WP0001  0->0.5... / (-1,1] = all but zero.
    WP0001 → WP1111, #WP0001 → WP0010
    WP0001 → WP1111, #WP0001 → WP0011
    WP0001 → WP1111, #WP0001 → WP0100
    WP0001 → WP1111, #WP0001 → WP0101
    WP0001 → WP1111, #WP0001 → WP0110
    WP0001 → WP1111, #WP0001 → WP0111
    WP0001 → WP1111, #WP0001 → WP1000

    WP0010 → WP1101, #WP0010 → WP0010  0.5... / (-1, 1] = [0.5, -0.5)
    WP0010 → WP1101, #WP0010 → WP0011
    WP0010 → WP1101, #WP0010 → WP0100
    WP0010 → WP1101, #WP0010 → WP0101
    WP0010 → WP1101, #WP0010 → WP0110
    WP0010 → WP1101, #WP0010 → WP0111
    WP0010 → WP1101, #WP0010 → WP1000

    WP0011 → WP1101, #WP0011 → WP0011  (0.5... / (-1, 1] = (0.5, -0.5)
    WP0011 → WP1101, #WP0011 → WP0100
    WP0011 → WP1101, #WP0011 → WP0101
    WP0011 → WP1101, #WP0011 → WP0110
    WP0011 → WP1101, #WP0011 → WP0111
    WP0011 → WP1101, #WP0011 → WP1000

    WP0100 → WP1011, #WP0100 → WP0100  1... / (-1, 1] = [1, -1)
    WP0100 → WP1011, #WP0100 → WP0101
    WP0100 → WP1011, #WP0100 → WP0110
    WP0100 → WP1011, #WP0100 → WP0111
    WP0100 → WP1011, #WP0100 → WP1000

    WP0101 → WP1011, #WP0101 → WP0101  (1... / (-1, 1] = (1, -1)
    WP0101 → WP1011, #WP0101 → WP0110
    WP0101 → WP1011, #WP0101 → WP0111
    WP0101 → WP1011, #WP0101 → WP1000

    WP0110 → WP1001, #WP0110 → WP0110  [2... / (-1, 1] = [2, -2)
    WP0110 → WP1001, #WP0110 → WP0111
    WP0110 → WP1001, #WP0110 → WP1000

    WP0111 → WP1001, #WP0111 → WP0111  (2... / (-1, 1] = (2, -2)
    WP0111 → WP1001, #WP0111 → WP1000

    WP1000 → WP1000, #WP1000 → WP1000
    ]

    for (a,b) in zip(warlpiri_div_rtz_values,pos_warlpiri)
        @test (a, b) == (b / WRTZ, b)
    end

    #next, test values when dividing from something which goes "round the
    #zero"... We'll use the interval (-1, 1], which is representable as 0b1101
    #-> 0b1100

    warlpiri_rtz_div_values = [
    WP1000 → WP0111, #WP0000 → WP0000   (-1, 1] / 0 = all reals
    WP1000 → WP0111, #WP0000 → WP0001
    WP1000 → WP0111, #WP0000 → WP0010
    WP1000 → WP0111, #WP0000 → WP0011
    WP1000 → WP0111, #WP0000 → WP0100
    WP1000 → WP0111, #WP0000 → WP0101
    WP1000 → WP0111, #WP0000 → WP0110
    WP1000 → WP0111, #WP0000 → WP0111
    WP1000 → WP0111, #WP0000 → WP1000

    WP1001 → WP0111, #WP0001 → WP0001  (-1,1] / (0... = all but inf.
    WP1001 → WP0111, #WP0001 → WP0010
    WP1001 → WP0111, #WP0001 → WP0011
    WP1001 → WP0111, #WP0001 → WP0100
    WP1001 → WP0111, #WP0001 → WP0101
    WP1001 → WP0111, #WP0001 → WP0110
    WP1001 → WP0111, #WP0001 → WP0111
    WP1001 → WP0111, #WP0001 → WP1000

    WP1011 → WP0110, #WP0010 → WP0010  (-1, 1] / [0.5... = (-2, 2]
    WP1011 → WP0110, #WP0010 → WP0011
    WP1011 → WP0110, #WP0010 → WP0100
    WP1011 → WP0110, #WP0010 → WP0101
    WP1011 → WP0110, #WP0010 → WP0110
    WP1011 → WP0110, #WP0010 → WP0111
    WP1011 → WP0110, #WP0010 → WP1000

    WP1011 → WP0101, #WP0011 → WP0011  (-1, 1] / (0.5... = (-2, 2)
    WP1011 → WP0101, #WP0011 → WP0100
    WP1011 → WP0101, #WP0011 → WP0101
    WP1011 → WP0101, #WP0011 → WP0110
    WP1011 → WP0101, #WP0011 → WP0111
    WP1011 → WP0101, #WP0011 → WP1000

    WP1101 → WP0100, #WP0100 → WP0100  (-1, 1] / [1... = (-1, 1]
    WP1101 → WP0100, #WP0100 → WP0101
    WP1101 → WP0100, #WP0100 → WP0110
    WP1101 → WP0100, #WP0100 → WP0111
    WP1101 → WP0100, #WP0100 → WP1000

    WP1101 → WP0011, #WP0101 → WP0101  (-1, 1] / (1... = (-1, 1)
    WP1101 → WP0011, #WP0101 → WP0110
    WP1101 → WP0011, #WP0101 → WP0111
    WP1101 → WP0011, #WP0101 → WP1000

    WP1111 → WP0010, #WP0110 → WP0110  (-1, 1] / [2... = (0.5, -0.5]
    WP1111 → WP0010, #WP0110 → WP0111
    WP1111 → WP0010, #WP0110 → WP1000

    WP1111 → WP0001, #WP0111 → WP0111  (-1, 1] / (2.. = (0.5, -0.5)
    WP1111 → WP0001, #WP0111 → WP1000

    WP0000 → WP0000, #WP1000 → WP1000  (-1, 1] / Inf = 0
    ]

    for (a,b) in zip(warlpiri_rtz_div_values,pos_warlpiri)
        @test (a, b) == (WRTZ / b, b)
    end

    #test two rounding-zero
    @test (WP1110 → WP0001) / (WP1010 → WP0101) == (WP1001 → WP1000)

    #test dividing by a plain rounding-infinity [1, -1)
    WRTI = WP0100 → WP1011

    warlpiri_div_rti_values = [
    WP0000 → WP0000, #WP0000 → WP0000  0       / [1, -1) = 0
    WP1111 → WP0001, #WP0000 → WP0001  [0, 0.5) / [1, -1) = (-0.5, 0.5)
    WP1111 → WP0010, #WP0000 → WP0010  [0, 0.5] / [1, -1) = (-0.5, 0.5]
    WP1101 → WP0011, #WP0000 → WP0011  [0,   1) / [1, -1) = (-1,     1)
    WP1101 → WP0100, #WP0000 → WP0100  [0,   1] / [1, -1) = (-1,     1]
    WP1011 → WP0101, #WP0000 → WP0101  [0,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0000 → WP0110  [0,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0000 → WP0111  [0, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0000 → WP1000  [0, Inf] / [1, -1) = allreals

    WP1111 → WP0001, #WP0001 → WP0001  (0, 0.5) / [1, -1) = (-0.5, 0.5)
    WP1111 → WP0010, #WP0001 → WP0010  (0, 0.5] / [1, -1) = (-0.5, 0.5]
    WP1101 → WP0011, #WP0001 → WP0011  (0,   1) / [1, -1) = (-1,     1)
    WP1101 → WP0100, #WP0001 → WP0100  (0,   1] / [1, -1) = (-1,     1]
    WP1011 → WP0101, #WP0001 → WP0101  (0,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0001 → WP0110  (0,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0001 → WP0111  (0, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0001 → WP1000  (0, Inf] / [1, -1) = allreals

    WP1111 → WP0010, #WP0010 → WP0010  0.5        / [1, -1) = (-0.5, 0.5]
    WP1101 → WP0011, #WP0010 → WP0011  [0.5,   1) / [1, -1) = (-1,     1)
    WP1101 → WP0100, #WP0010 → WP0100  [0.5,   1] / [1, -1) = (-1,     1]
    WP1011 → WP0101, #WP0010 → WP0101  [0.5,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0010 → WP0110  [0.5,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0010 → WP0111  [0.5, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0010 → WP1000  [0.5, Inf] / [1, -1) = allreals

    WP1101 → WP0011, #WP0011 → WP0011  (0.5,   1) / [1, -1) = (-1,     1)
    WP1101 → WP0100, #WP0011 → WP0100  (0.5,   1] / [1, -1) = (-1,     1]
    WP1011 → WP0101, #WP0011 → WP0101  (0.5,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0011 → WP0110  (0.5,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0011 → WP0111  (0.5, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0011 → WP1000  (0.5, Inf] / [1, -1) = allreals

    WP1101 → WP0100, #WP0100 → WP0100  1        / [1, -1) = (-1,     1]
    WP1011 → WP0101, #WP0100 → WP0101  [1,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0100 → WP0110  [1,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0100 → WP0111  [1, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0100 → WP1000  [1, Inf] / [1, -1) = allreals

    WP1011 → WP0101, #WP0101 → WP0101  (1,   2) / [1, -1) = (-2,     2)
    WP1011 → WP0110, #WP0101 → WP0110  (1,   2] / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0101 → WP0111  (1, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0101 → WP1000  (1, Inf] / [1, -1) = allreals

    WP1011 → WP0110, #WP0110 → WP0110  2        / [1, -1) = (-2,     2]
    WP1001 → WP0111, #WP0110 → WP0111  [2, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0110 → WP1000  [2, Inf] / [1, -1) = allreals

    WP1001 → WP0111, #WP0111 → WP0111  (2, Inf) / [1, -1) = (-Inf, Inf)
    WP1001 → WP1000, #WP0111 → WP1000  (2, Inf] / [1, -1) = allreals

    WP1001 → WP1000, #WP1000 → WP1000  Inf / [1, -1) = allreals
    ]

    for (a,b) in zip(warlpiri_div_rti_values,pos_warlpiri)
        @test (a, b) == (b / WRTI, b)
    end

    #test dividing a rounding-inifinity
    warlpiri_rti_div_values = [
    WP1000 → WP1000, #WP0000 → WP0000  [1, -1) / 0         = Inf
    WP0111 → WP1001, #WP0000 → WP0001  [1, -1) / [0, 0.5)  = (2,     -2)
    WP0110 → WP1001, #WP0000 → WP0010  [1, -1) / [0, 0.5]  = [2,     -2)
    WP0101 → WP1011, #WP0000 → WP0011  [1, -1) / [0,   1)  = (1,     -1)
    WP0100 → WP1011, #WP0000 → WP0100  [1, -1) / [0,   1]  = [1,     -1)
    WP0011 → WP1101, #WP0000 → WP0101  [1, -1) / [0,   2)  = (0.5, -0.5)
    WP0010 → WP1101, #WP0000 → WP0110  [1, -1) / [0,   2]  = [0.5, -0.5)
    WP0001 → WP1111, #WP0000 → WP0111  [1, -1) / [0, Inf)  = allreals*
    WP1001 → WP1000, #WP0000 → WP1000  [1, -1) / [0, Inf]  = allreals

    WP0111 → WP1001, #WP0001 → WP0001  [1, -1) / (0, 0.5)  = (2,      -2)
    WP0110 → WP1001, #WP0001 → WP0010  [1, -1) / (0, 0.5]  = [2,      -2)
    WP0101 → WP1011, #WP0001 → WP0011  [1, -1) / (0,   1)  = (1,      -1)
    WP0100 → WP1011, #WP0001 → WP0100  [1, -1) / (0,   1]  = [1,      -1)
    WP0011 → WP1101, #WP0001 → WP0101  [1, -1) / (0,   2)  = (0.5, -`0.5)
    WP0010 → WP1101, #WP0001 → WP0110  [1, -1) / (0,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0001 → WP0111  [1, -1) / (0, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0001 → WP1000  [1, -1) / (0, Inf]  = allreal`s

    WP0110 → WP1001, #WP0010 → WP0010  [1, -1) / [0.5, 0.5]  = [2,      -2)
    WP0101 → WP1011, #WP0010 → WP0011  [1, -1) / [0.5,   1)  = (1,      -1)
    WP0100 → WP1011, #WP0010 → WP0100  [1, -1) / [0.5,   1]  = [1,      -1)
    WP0011 → WP1101, #WP0010 → WP0101  [1, -1) / [0.5,   2)  = (0.5, -`0.5)
    WP0010 → WP1101, #WP0010 → WP0110  [1, -1) / [0.5,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0010 → WP0111  [1, -1) / [0.5, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0010 → WP1000  [1, -1) / [0.5, Inf]  = allreal`s

    WP0101 → WP1011, #WP0011 → WP0011  [1, -1) / (0.5,   1)  = (1,      -1)
    WP0100 → WP1011, #WP0011 → WP0100  [1, -1) / (0.5,   1]  = [1,      -1)
    WP0011 → WP1101, #WP0011 → WP0101  [1, -1) / (0.5,   2)  = (0.5, -`0.5)
    WP0010 → WP1101, #WP0011 → WP0110  [1, -1) / (0.5,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0011 → WP0111  [1, -1) / (0.5, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0011 → WP1000  [1, -1) / (0.5, Inf]  = allreal`s

    WP0100 → WP1011, #WP0100 → WP0100  [1, -1) / [1,   1]  = [1,      -1)
    WP0011 → WP1101, #WP0100 → WP0101  [1, -1) / [1,   2)  = (0.5, -`0.5)
    WP0010 → WP1101, #WP0100 → WP0110  [1, -1) / [1,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0100 → WP0111  [1, -1) / [1, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0100 → WP1000  [1, -1) / [1, Inf]  = allreal`s

    WP0011 → WP1101, #WP0101 → WP0101  [1, -1) / (1,   2)  = (0.5, -`0.5)
    WP0010 → WP1101, #WP0101 → WP0110  [1, -1) / (1,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0101 → WP0111  [1, -1) / (1, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0101 → WP1000  [1, -1) / (1, Inf]  = allreal`s

    WP0010 → WP1101, #WP0101 → WP0110  [1, -1) / [2,   2]  = (0.5, -`0.5)
    WP0001 → WP1111, #WP0101 → WP0111  [1, -1) / [2, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0101 → WP1000  [1, -1) / [2, Inf]  = allreal`s

    WP0001 → WP1111, #WP0111 → WP0111  [1, -1) / (2, Inf)  = allreal`s*
    WP1001 → WP1000, #WP0111 → WP1000  [1, -1) / (2, Inf]  = allreal`s

    WP1001 → WP1000, #WP1000 → WP1000  [1, -1) / Inf = allreals
    ]

    for (a,b) in zip(warlpiri_rti_div_values,pos_warlpiri)
        @test (a, b) == (WRTI / b, b)
    end

    #test dividing by values that round both zero and infinity...  In this  case (2, 0.5]
    #values copied straight from multiplication with [2, 0.5)
    WRTB = WP0111 → WP0010

    warlpiri_div_rtb_values = [
    WP1001 → WP1000, #WP0000 → WP0000  #anything that contains zero will trigger all reals.
    WP1001 → WP1000, #WP0000 → WP0001
    WP1001 → WP1000, #WP0000 → WP0010
    WP1001 → WP1000, #WP0000 → WP0011
    WP1001 → WP1000, #WP0000 → WP0100
    WP1001 → WP1000, #WP0000 → WP0101
    WP1001 → WP1000, #WP0000 → WP0110
    WP1001 → WP1000, #WP0000 → WP0111
    WP1001 → WP1000, #WP0000 → WP1000

    WP1001 → WP1000, #WP0001 → WP0001
    WP1001 → WP1000, #WP0001 → WP0010
    WP1001 → WP1000, #WP0001 → WP0011
    WP1001 → WP1000, #WP0001 → WP0100
    WP1001 → WP1000, #WP0001 → WP0101
    WP1001 → WP1000, #WP0001 → WP0110
    WP1001 → WP1000, #WP0001 → WP0111
    WP1001 → WP1000, #WP0001 → WP1000

    WP0100 → WP0001, #WP0010 → WP0010
    WP0100 → WP0001, #WP0010 → WP0011
    WP0100 → WP0001, #WP0010 → WP0100
    WP1001 → WP1000, #WP0010 → WP0101
    WP1001 → WP1000, #WP0010 → WP0110
    WP1001 → WP1000, #WP0010 → WP0111
    WP1001 → WP1000, #WP0010 → WP1000

    WP0101 → WP0001, #WP0011 → WP0011
    WP0101 → WP0001, #WP0011 → WP0100
    WP0101 → WP0011, #WP0011 → WP0101
    WP0101 → WP0011, #WP0011 → WP0110
    WP1001 → WP1000, #WP0011 → WP0111
    WP1001 → WP1000, #WP0011 → WP1000

    WP0110 → WP0001, #WP0100 → WP0100
    WP0110 → WP0011, #WP0100 → WP0101
    WP0110 → WP0011, #WP0100 → WP0110
    WP1001 → WP1000, #WP0100 → WP0111
    WP1001 → WP1000, #WP0100 → WP1000

    WP0111 → WP0011, #WP0101 → WP0101
    WP0111 → WP0011, #WP0101 → WP0110
    WP1001 → WP1000, #WP0101 → WP0111
    WP1001 → WP1000, #WP0101 → WP1000

    WP0111 → WP0011, #WP0110 → WP0110
    WP1001 → WP1000, #WP0110 → WP0111
    WP1001 → WP1000, #WP0110 → WP1000

    WP1001 → WP1000, #WP0111 → WP0111
    WP1001 → WP1000, #WP0111 → WP1000

    WP1001 → WP1000, #WP1000 → WP1000
    ]

    for (a,b) in zip(warlpiri_div_rtb_values,pos_warlpiri)
        @test (a, b) == (b / WRTB, b)
    end

    #test dividing values that round both zero and infinity...  In this  case (2, 0.5]

    warlpiri_rtb_div_values = [
    WP1001 → WP1000, #WP0000 → WP0000  #anything that contains zero will trigger all reals.
    WP1001 → WP1000, #WP0000 → WP0001
    WP1001 → WP1000, #WP0000 → WP0010
    WP1001 → WP1000, #WP0000 → WP0011
    WP1001 → WP1000, #WP0000 → WP0100
    WP1001 → WP1000, #WP0000 → WP0101
    WP1001 → WP1000, #WP0000 → WP0110
    WP1001 → WP1000, #WP0000 → WP0111
    WP1001 → WP1000, #WP0000 → WP1000

    WP1001 → WP1000, #WP0001 → WP0001  (2, 0.5] / (0, 0.5) = allreals
    WP1001 → WP1000, #WP0001 → WP0010
    WP1001 → WP1000, #WP0001 → WP0011
    WP1001 → WP1000, #WP0001 → WP0100
    WP1001 → WP1000, #WP0001 → WP0101
    WP1001 → WP1000, #WP0001 → WP0110
    WP1001 → WP1000, #WP0001 → WP0111
    WP1001 → WP1000, #WP0001 → WP1000

    WP0111 → WP0100, #WP0010 → WP0010  (2, 0.5] / 0.5 = (2, 1]
    WP0111 → WP0100, #WP0010 → WP0011  (2, 0.5] / [0.5, 1) = (2, 1]
    WP0111 → WP0100, #WP0010 → WP0100  (2, 0.5] / [0.5, 1] = (2, 1]
    WP1001 → WP1000, #WP0010 → WP0101  (2, 0.5] / [0.5, 2) = allreals
    WP1001 → WP1000, #WP0010 → WP0110
    WP1001 → WP1000, #WP0010 → WP0111
    WP1001 → WP1000, #WP0010 → WP1000

    WP0111 → WP0011, #WP0011 → WP0011  (2, 0.5] / (0.5, 1) = (2, 1)
    WP0111 → WP0011, #WP0011 → WP0100  (2, 0.5] / (0.5, 1] = (2, 1)
    WP0101 → WP0011, #WP0011 → WP0101  (2, 0.5] / (0.5, 2) = (not 1)
    WP0101 → WP0011, #WP0011 → WP0110  (2, 0.5] / (0.5, 2] = (not 1)
    WP1001 → WP1000, #WP0011 → WP0111
    WP1001 → WP1000, #WP0011 → WP1000

    WP0111 → WP0010, #WP0100 → WP0100  (2, 0.5] / 1        = (2, 0.5]
    WP0101 → WP0010, #WP0100 → WP0101  (2, 0.5] / [1,   2) = (1, 0.5]
    WP0101 → WP0010, #WP0100 → WP0110  (2, 0.5] / [1,   2] = (1, 0.5]
    WP1001 → WP1000, #WP0100 → WP0111  (2, 0.5] / [1, Inf) = allreals
    WP1001 → WP1000, #WP0100 → WP1000

    WP0101 → WP0001, #WP0101 → WP0101  (2, 0.5] / (1,   2) = (1, 0.5)
    WP0101 → WP0001, #WP0101 → WP0110  (2, 0.5] / (1,   2] = (1, 0.5)
    WP1001 → WP1000, #WP0101 → WP0111  (2, 0.5] / (1, Inf) = allreals
    WP1001 → WP1000, #WP0101 → WP1000

    WP0101 → WP0001, #WP0110 → WP0110  (2, 0.5] / 2        = (1, 0.5)
    WP1001 → WP1000, #WP0110 → WP0111  (2, 0.5] / (2, Inf) = allreals
    WP1001 → WP1000, #WP0110 → WP1000

    WP1001 → WP1000, #WP0111 → WP0111
    WP1001 → WP1000, #WP0111 → WP1000

    WP1001 → WP1000, #WP1000 → WP1000
    ]

    for (a,b) in zip(warlpiri_rtb_div_values,pos_warlpiri)
        @test (a, b) == (WRTB / b, b)
    end

    @test (WP1001 → WP1000) == WRTI / WRTB
    @test (WP1001 → WP1000) == WRTZ / WRTB
    @test (WP1001 → WP1000) == WRTB / WRTB

    @test (WP1001 → WP1000) == WRTB / WRTI
    @test (WP1001 → WP1000) == WRTB / WRTZ
    @test (WP1001 → WP1000) == WRTB / WRTB
end

ext: jl

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
export OnlineClustering

mutable struct OnlineClustering <: Clustering
    fields
end

function (olc::OnlineClustering)(feature)
    body
end

max_stars_repo_name: Experica/SpikeSorting

max_issues_repo_name: Experica/SpikeSorting

max_forks_repo_name: Experica/SpikeSorting


text:
#===========================================================m====================
    calibration.jl

    calibration procedure:
    1. multichannel deconvolution
    2. normalize variation between wells in absolute fluorescence

===============================================================================#

import DataStructures.OrderedDict
import StaticArrays: MArray, SArray, SMatrix, SVector
import Memento: debug, error

#===============================================================================
    constants >>
===============================================================================#

## default values
const DEFAULT_CAL_FIRST_WELL        = 1

## preset values
const NORMALIZATION_SCALING_FACTOR  = 3.7           ## used: 9e5, 1e5, 1.2e6, 3.0
const DECONVOLUTION_SCALING_FACTOR  = [1.0, 4.2]    ## used: [1, oneof(1, 2, 3.5, 8, 7, 5.6, 4.2)]
const DECONVOLUTION_BACKUP_K        = K4DCV         ## not used

## convenience array indices
const WATER  = 1
const SIGNAL = 2

#===============================================================================
    function definitions >>
===============================================================================#

"Perform deconvolution between channels and normalize variation between wells."
function calibrate(
    i                       ::Input,
    calibration_data        ::CalibrationData{<: NumberOfChannels, <: Union{Int_T,Float_T}},
    calibration_args        ::CalibrationParameters,
    raw                     ::RawData{<: Union{Int_T,Float_T}},         ## 3D array of raw fluorescence
    data_format             ::DataFormat,               ## array, dict, both
)
    debug(logger, "at calibrate()")

    ## remove MySql dependency
    # calib_info = ensure_ci(db_conn, calib_info)
    # norm_data, norm_well_nums = prep_normalize(db_conn,
    #     calib_info, well_nums_in_req, dye_in, dyes_to_fill)

    ## prepare data for normalization
    (norm_data, norm_wells) = prep_normalize(calibration_data)

    # if length(well_nums_found_in_req) == 0
    #     well_nums_found_in_req = norm_wells
    # end

    ## remove MySql dependency
    #
    # matched_well_idc = find(norm_well_nums) do norm_well_num
    #     norm_well_num in well_nums_found_in_raw
    # end ## do norm_well_num

    ## match (explicit) well numbers in the experimental data
    ## with (presumed complete) set of well numbers in calibration data
    matched_wells = indexin(i.wells, norm_wells)
    matched_exp_well_idc = find(matched_wells)
    matched_calib_well_idc = SVector{length(matched_exp_well_idc)}(
        matched_wells[matched_exp_well_idc])

    #
    ## subtract background
    background_subtracted_data =
        subtract_background(
            raw.data,
            norm_data[:,:,WATER],
            matched_wells,
            matched_exp_well_idc)

    #
    (k_deconv, deconvoluted_data) =
        if  calibration_args.dcv &&
            isa(calibration_data, CalibrationData{DualChannel,<: Union{Int_T,Float_T}})
            ## addition with flexible ratio instead of deconvolution (commented out)
            # k_inv_vec = fill(reshape(DataArray([1, 0, 1, 0]), 2, 2), 16)

            ## removing MySql dependency
            # k4dcv, deconvoluted_data = deconvolute(
            #     1. * background_subtracted_data,
            #     channel_nums, matched_well_idc,
            #     db_conn, calib_info, well_nums_in_req;
            #     out_format="array")
            deconvolute(
                calibration_data,
                calibration_args,
                background_subtracted_data |> cast(bless(Float_T)),
                matched_calib_well_idc,
                norm_wells[matched_calib_well_idc],
                DECONVOLUTION_BACKUP_K,
                DECONVOLUTION_SCALING_FACTOR,
                array)
        else ## !dcv
            DeconvolutionMatrices(calibration_data, calibration_args),
                background_subtracted_data
        end
    #
    calibrated_array =
        map(1:i.num_channels) do channel
            perform_normalize(
                deconvoluted_data[:, :, channel],
                norm_data,
                matched_exp_well_idc,
                matched_calib_well_idc,
                channel,
                calibration_args.subtract_water,
                NORMALIZATION_SCALING_FACTOR)
        end #= do channel =# |>
        splat(tie(3)) ## not gather(tie(3)) because we need 3 dimensions for single channel data
    #
    ## format output
    ## the following line of code needs the keys of calibrated_dict to be in sort order
    calibrated_dict()  =
        OrderedDict(
            map(1:i.num_channels) do channel
                channel => calibrated_array[:, :, channel]
            end)
    calibrated_data =
        if      (data_format == array)  tuple(calibrated_array)
        elseif  (data_format == dict)   tuple(calibrated_dict())
        elseif  (data_format == both)   tuple(calibrated_array, calibrated_dict())
        else                            throw(ArgumentException(
                                            "`data_format` must be array, dict or both"))
        end ## if
    return (
        background_subtracted_data,
        k_deconv,
        deconvoluted_data,
        norm_data,
        norm_wells,
        calibrated_data...)
end ## calibrate()

## called by calibrate()
function subtract_background(
    raw_data                ::Array{<: Union{Int_T,Float_T}, 3},
    water_data              ::SArray{S,<: Union{Int_T,Float_T},2} where {S},
    matched_wells           ::Array{Int_T},
    matched_exp_wells_idc   ::Array{Int_T},
)
    ## vectorized
    # raw_data[:,matched_exp_wells_idc,:] .-
    #     permutedims(cat(3,water_data[matched_wells[matched_exp_wells_idc],:]),[3,1,2])
    ## devectorized
    [   raw_data[ui, wi, ci] - water_data[matched_wells[wi], ci]
        for ui in 1:size(raw_data, 1),
            wi in matched_exp_wells_idc,
            ci in 1:size(raw_data, 3)   ]
end

#==============================================================================#

"Check validity of optical calibration data."
function prep_normalize(calibration_data ::CalibrationData{DualChannel, <: Union{Int_T,Float_T}})
    debug(logger, "at prep_normalize()")
    norm_data = get_norm_data(calibration_data)
    ## index wells in calibration data starting at `DEFAULT_CAL_FIRST_WELL` = 0
    ## issue:
    ## using the current format for the request body there is no well_num information
    ## associated with the calibration data
    num_wells = size(calibration_data.array, 1)
    num_channels = 2
    signal_wells = num_wells |> from(DEFAULT_CAL_FIRST_WELL) |>
        mold(Symbol) |> SVector{num_wells}
    ## check whether signal fluorescence > water fluorescence
    failed = norm_data[:, :, SIGNAL] .<= norm_data[:, :, WATER]
    if any(failed)
        err_msgs = Vector{String}()
        for channel in 1:num_channels
            invalid = find(failed[:, channel])
            if length(invalid) > 0
                failed_wells = join(signal_wells[invalid], ", ")
                push!(err_msgs, "invalid calibration data in channel $channel: " *
                    "fluorescence value of water is greater than or equal to that of dye " *
                    "in the following well(s) - " * failed_wells)
            end ## if invalid
        end ## next channel
        throw(ArgumentError(join(err_msgs, "; ")))
    end ## if any(failed)

    ## issue:
    ## this feature has been temporarily disabled while
    ## removing MySql dependency in get_wva_data because
    ## using the current format for the request body
    ## we cannot subset the calibration data by step_id
    # if length(dyes_to_fill) > 0 ## extrapolate well-to-well variation data for missing channels
    #     println("Preset well-to-well variation data is used to extrapolate " *
    #         "calibration data for missing channels.")
    #     channels_missing = setdiff(channels_in_water, channels_in_signal)
    #     dyes_dyes_to_fill_channels =
    #         map(dye -> DYE2CHST[dye]["channel"], dyes_to_fill) ## DYE2CHST is defined in module scope
    #     check_subset(
    #         Ccsc(channels_missing, "Channels missing well-to-well variation data"),
    #         Ccsc(dyes_dyes_to_fill_channels, "channels corresponding to the dyes " *
    #             "of which well-to-well variation data is needed")
    #     )
    #     # process preset calibration data
    #     preset_step_ids = OrderedDict([
    #         dye => DYE2CHST[dye]["step_id"]
    #         for dye in keys(DYE2CHST)
    #     ])
    #     preset_signal_data_dict, dyes_in_preset =
    #         get_wva_data(PRESET_calib_ids["signal"], preset_step_ids, db_conn, "dye") ## `well_nums` is not passed on
    #     pivot_preset = preset_signal_data_dict[dye_in]
    #     pivot_in = signal_data_dict[DYE2CHST[dye_in]["channel"]]
    #     in2preset = pivot_in ./ pivot_preset
    #     for dye in dyes_to_fill
    #         signal_data_dict[DYE2CHST[dye]["channel"]] = preset_signal_data_dict[dye] .* in2preset
    #     end
    # end ## if

    return (norm_data, signal_wells)
end ## prep_normalize()

function prep_normalize(calibration_data ::CalibrationData{SingleChannel, <: Union{Int_T,Float_T}})
    debug(logger, "at prep_normalize()")
    norm_data = calibration_data.array
    num_wells = size(calibration_data.array, 1)
    signal_wells = num_wells |> from(DEFAULT_CAL_FIRST_WELL) |>
        mold(Symbol) |> SVector{num_wells}
    ## check whether signal fluorescence > water fluorescence
    failed = norm_data[:, 1, SIGNAL] .<= norm_data[:, 1, WATER]
    if any(failed)
        err_msg = "invalid calibration data in channel 1: " *
            "fluorescence value of water is greater than or equal to that of dye " *
            "in the following well(s) - " * join(signal_wells[find(failed)], ", ")
        throw(ArgumentError(err_msg))
    end ## if any(failed)
    return (norm_data, signal_wells)
end ## prep_normalize()

## called by prep_normalize() >>

get_norm_data(data ::CalibrationData{SingleChannel, <: Union{Int_T,Float_T}}) =
    data.array

function get_norm_data(data ::CalibrationData{DualChannel,R}) where {R <: Union{Int_T,Float_T}}
    local (water, dye1, dye2) = 1:3
    SArray{Tuple{size(data.array,1),2,2},R}( ## hcat converts from SArray to Array
        hcat(
            data.array[:, SVector{1}(1), SVector{2}(water, dye1)],    ## channel 1
            data.array[:, SVector{1}(2), SVector{2}(water, dye2)]))   ## channel 2
end

#==============================================================================#

#=
    perform_normalize(): each dye only has data for its target channel;
    calibration/calib/oc - used for `deconvolute` and `normalize`,
    each dye has data for both target and non-target channels.
    Input `fluo` and output: dim2 indexed by well and dim1 indexed by unit,
    which can be cycle (amplification) or temperature point (melt curve).
    Output does not include the automatically created column at index 1
    from rownames of input array as R does
=#

## called from calibrate()
"Normalize variation between wells in absolute fluorescence values."
function perform_normalize(
    fluorescence                    ::Array{<: Union{Int_T,Float_T},2},
    norm_data                       ::SArray{S,<: Union{Int_T,Float_T},3} where {S},
    matched_exp_well_idc            ::Array{Int_T},
    matched_calib_well_idc          ::SVector{L,Int_T} where {L},
    channel                         ::Int_T,
    subtract_water                  ::Bool,
    normalization_scaling_factor    ::Float_T,
)
    debug(logger, "at perform_normalize()")
    #
    ## devectorized code avoids transposing data matrix
    if subtract_water == false
        smw = norm_data[matched_calib_well_idc, channel, SIGNAL]
        return ([
            normalization_scaling_factor * mean(smw) *
                fluorescence[u, matched_exp_well_idc[wi]] / smw[wi]
                    for u  in 1:size(fluorescence, 1),
                        wi in eachindex(matched_exp_well_idc)])
    end
    #
    ## subtract_water == true
    norm_water = norm_data[matched_calib_well_idc, channel, WATER]
    smw = norm_data[matched_calib_well_idc, channel, SIGNAL] .- norm_water
    return ([
        normalization_scaling_factor * mean(smw) *
            (fluorescence[u, matched_exp_well_idc[wi]] - norm_water[wi]) / smw[wi]
                for u  in 1:size(fluorescence, 1),
                    wi in eachindex(matched_exp_well_idc)])
end ## perform_normalize

#==============================================================================#

## called from calibrate()
"Perform multi-channel deconvolution of fluorescence values."
function deconvolute(
    ## ary2dcv dim1 is unit, which can be cycle (amplification), temperature point (melting curve),
    ## or step type (like "water", "channel_1", "channel_2" for calibration experiment);
    ## ary2dcv dim2 must be well, ary2dcv dim3 must be channel
    calibration_data        ::CalibrationData{DualChannel, <: Union{Int_T,Float_T}},
    calibration_args        ::CalibrationParameters,
    ary2dcv                 ::Array{Float_T,3},
    matched_calib_well_idc  ::SVector{L,Int_T} where {L},
    wells                   ::SVector{L,Symbol} where {L},
    k_deconv_backup         ::DeconvolutionMatrices, ## argument not used
    scaling_factor_dcv_vec  ::Vector{Float_T},
    data_format             ::DataFormat, ## array, dict, both

    ## remove MySql dependency
    #
    ## arguments needed if `k` matrix needs to be computed
    ## `db_conn_default` is defined in "__init__.jl"
    # db_conn ::MySQL.MySQLHandle=db_conn_default,
    # calib_info ::Union{Integer,OrderedDict}=calib_info_AIR,
    # well_nums ::AbstractVector=[];
)
    debug(logger, "at deconvolute()")

    ## remove MySql dependency
    # k_deconv = (isa(calib_info, Integer) || begin
    #     step_ids = map(ci_value -> ci_value["step_id"], values(calib_info))
    #     length_step_ids = length(step_ids)
    #     length_step_ids <= 2 || length(unique(step_ids)) < length_step_ids
    # end) ? k_deconv_backup : get_k(db_conn, calib_info, well_nums) ## use default `well_proc` value

    k_deconv = get_k(
        calibration_data,
        calibration_args,
        matched_calib_well_idc,
        wells)
    (num_units, num_channels) = size(ary2dcv,1,3)
    scaled_k_inv_vecs =
        map(eachindex(matched_calib_well_idc)) do wi
            k_deconv.k_inv_vec[wi] .* scaling_factor_dcv_vec
        end
    deconvoluted_array = similar(ary2dcv)
    for ui in 1:num_units, wi in eachindex(matched_calib_well_idc)
        deconvoluted_array[ui, wi, :] =
            scaled_k_inv_vecs[wi] * ary2dcv[ui, wi, :] ## matrix * vector
    end
    deconvoluted_dict() =
        OrderedDict(
            channel => deconvoluted_array[:, :, channel]
            for channel in 1:num_channels)
    ## format output
    deconvoluted_data =
        if      (data_format == array)  tuple(deconvoluted_array)
        elseif  (data_format == dict)   tuple(deconvoluted_dict())
        elseif  (data_format == both)   tuple(deconvoluted_array, deconvoluted_dict())
        else
            throw(ArgumentError("`out_format` must be array, dict or both"))
        end ## if
    return (k_deconv, deconvoluted_data...)
end ## deconvolute()

deconvolute(
    calibration_data        ::CalibrationData{SingleChannel, <: Union{Int_T,Float_T}},
    calibration_args        ::CalibrationParameters,
    ary2dcv                 ::Array{Float_T,3},
    matched_calib_well_idc  ::SVector{L,Int_T} where {L},
    wells                   ::SVector{L,Symbol} where {L},
    k_deconv_backup         ::DeconvolutionMatrices, ## argument not used
    scaling_factor_dcv_vec  ::Vector{Float_T},
    data_format             ::DataFormat, ## array, dict, both
) =
    throw(ArgumentError("cannot perform multi-channel deconvolution " *
        "using single channel calibration data"))

#==============================================================================#

"Calculate the matrix `k` that describes crosstalk between channels."
function get_k(
    calibration_data        ::CalibrationData{DualChannel, <: Union{Int_T,Float_T}},
    calibration_args        ::CalibrationParameters,
    matched_calib_well_idc  ::SVector{L,Int_T} where {L},
    wells                   ::SVector{L,Symbol} where {L};
    save_to                 ::String = "" ## used: "k.jld"

    ## remove MySql dependency
    #
    # db_conn ::MySQL.MySQLHandle,
    ## info on experiment(s) used to calculate matrix K
    ## OrderedDict(
    ##    "water"    =OrderedDict(calibration_id=..., step_id=...),
    ##    "channel_1"=OrderedDict(calibration_id=..., step_id=...),
    ##    "channel_2"=OrderedDict(calibration_id=..., step_id=...))
    # dcv_exp_info ::OrderedDict,
    #
    ## possible  issue:
    ## originally step_ids were provided together with calibration data
    ## I'm not sure that this is a problem now because the calibration data
    ## in the request body is already specific to a single step.
)
    debug(logger, "at get_k()")

    ## remove MySql dependency
    #
    # dcv_exp_info = ensure_ci(db_conn, dcv_exp_info)
    #
    # calib_key_vec = get_ordered_keys(dcv_exp_info)
    # cd_key_vec = calib_key_vec[2:end] # cd = channel of dye, "water" is index 1
    #
    # dcv_data_dict = get_full_calib_data(db_conn, dcv_exp_info, well_nums)
    #
    # water_data, water_well_nums = dcv_data_dict["water"]
    # num_wells = length(water_well_nums)
    #
    ## `dcv_well_nums` is not passed on because it is
    ## expected to be the same as `water_well_nums`,
    ## otherwise error will be raised by `get_full_calib_data`
    # smw_bydye = OrderedDict(map(cd_key_vec) do cd_key
    #    k_data_1dye, dcv_well_nums = dcv_data_dict[cd_key]
    #    return cd_key => k_data_1dye .- water_data
    # end)

    ## water-subtracted calibration data
    smw = SArray{Tuple{2,2,length(matched_calib_well_idc)},Float_T}([
        calibration_data.array[well, channel, dye] - calibration_data.array[well, channel, WATER]
        for channel in 1:2,
            dye     in 2:3,
            well    in matched_calib_well_idc])
    # const smw2 = SArray{Tuple{2,2,length(matched_calib_well_idc)},Float_T}(
    #
    ## check that the water-subtracted signal in the target channel
    ## is greater than that in the non-target channel(s) for each well and each dye
    failed =
        map(DYES) do dye
            target     = dye
            non_target = 3 - target
            smw[target, dye, :] .<= smw[non_target, dye, :]
        end #= next channel =# |>
        gather(hcat)
    if any(failed)
        err_msgs = Vector{String}()
        for dye in DYES
            invalid = find(failed[:, dye])
            if length(invalid) > 0
                target     = dye
                non_target = 3 - target
                wells_invalid = join(wells[invalid], ", ")
                push!(err_msgs,
                    "invalid deconvolution data for the dye targeting channel $target: " *
                    "fluorescence value of non-target channel $non_target " *
                    "is greater than or equal to that of target channel $target " *
                    "in the following well(s) - $wells_invalid")
            end ## if
        end ## next dye
        throw(ArgumentError(join(err_msgs, "; ")))
    end ## if any(failed)
    #
    ## compute inverses and return
    INV_NOTE_PT2 = ": K matrix is singular, using `pinv` instead of `inv` " *
    "to compute inverse matrix of K. Deconvolution result may not be accurate. " *
    "This may be caused by using the same or a similar set of solutions " *
    "in the steps for different dyes."
    (k_s, k_inv_vec, inv_note) =
        calc_kinv(Val{calibration_args.k_method}, smw, DYES, wells)
    k_deconv =
        DeconvolutionMatrices(
            k_s,
            k_inv_vec,
            (length(inv_note) > 0 ? inv_note * INV_NOTE_PT2 : ""))
    (length(save_to) > 0) && save(save_to, "k_deconv", k_deconv)
    return k_deconv
end ## get_k()

get_k(
    calibration_data        ::CalibrationData{SingleChannel, <: Union{Int_T,Float_T}},
    calibration_args        ::CalibrationParameters,
    matched_calib_well_idc  ::SVector{L,Int_T} where {L},
    wells                   ::SVector{L,Symbol} where {L};
    save_to                 ::String = "" ## used: "k.jld"
) =
    throw(ArgumentError("cannot calculate deconvolution matrices " *
        "using single channel calibration data"))

#==============================================================================#

## functions called by get_k() >>

"Calculate a single deconvolution matrix K, averaging across all wells"
function calc_kinv(
    ::Type{Val{well_proc_mean}},
    smw             ::SArray{S,<: Float_T,3} where {S},
    dyes            ::SVector{L,Int_T},
    wells           ::SVector{L,Symbol}
) where {L}
    inv_note = false
    k_s =
        map(eachindex(dyes)) do dye
            SVector{L,Float_T}(sweep(sum)(/)(mean(smw[:,dye,:], 2)))
        end #= next dye =# |>
        gather(hcat)
    k_inv = try
        inv(k_s)
    catch err
        if isa(err, Base.LinAlg.SingularException)
            inv_note = true
            pinv(k_s)
        else
            rethrow(err)
        end ## if isa(err,
    end ## try
    return
        [k_s],
        fill(k_inv, length(wells)),
        inv_note ? "" : "Well mean"
end

"Calculate deconvolution matrix K for each well"
function calc_kinv(
    ::Type{Val{well_proc_vec}},
    smw             ::SArray{S,<: Float_T,3} where {S},
    dyes            ::SVector{L,Int_T},
    wells           ::SVector{M,Symbol} where {M}
) where {L}
    singular_wells = Vector{Symbol}()
    k_s =
        map(eachindex(wells)) do wi
            map(eachindex(dyes)) do dye
                sweep(sum)(/)(smw[:, dye, wi])
            end #= next dye =# |>
            gather(hcat)
        end ## next wi
    k_inv_vec =
        map(eachindex(wells)) do wi
            try
                ## `inv()` is supposedly faster with StaticArray
                inv(k_s[wi])
            catch err
                if  isa(err, Base.LinAlg.SingularException) ||
                    isa(err, Base.LinAlg.LAPACKException)
                        push!(singular_wells, wells[wi])
                        pinv(k_s[wi])
                else
                    rethrow(err)
                end ## if isa(err
            end ## try
        end ## next wi
    inv_note = (length(singular_wells) > 0) ?
        "Well(s) " * string(join(singular_wells, ", ")) : ""
    return k_s, k_inv_vec, inv_note
end

#==============================================================================#

## unused function
# function calib_calib(
#     ## remove MySql dependency
#     #
#     # db_conn_1 ::MySQL.MySQLHandle,
#     # db_conn_2 ::MySQL.MySQLHandle,
#     # calib_info_1 ::OrderedDict,
#     # calib_info_2 ::OrderedDict,
#     # well_nums_1 ::AbstractVector=[],
#     # well_nums_2 ::AbstractVector=[];
#     dye_in      ::Symbol =:FAM,
#     dyes_to_fill ::AbstractVector =[]
# )
#     ## This function is expected to handle situations where `calib_info_1` and `calib_info_2`
#     ## have different combinations of wells, but the number of wells should be the same.
#     if length(well_nums_1) != length(well_nums_2)
#         throw(DimensionMismatch("length(well_nums_1) != length(well_nums_2)"))
#     end

#     ## remove MySql dependency
#     #
#     # calib_dict_1 = get_full_calib_data(db_conn_1, calib_info_1, well_nums_1)
#     # water_well_nums_1 = calib_dict_1["water"][2]
#     #
#     # calib_key_vec_1 = get_ordered_keys(calib_info_1)
#     # cd_key_vec_1 = calib_key_vec_1[2:end] ## cd = channel of dye. "water" is index 1 per original order.
#     # channel_nums_1 = map(cd_key_vec_1) do cd_key
#     #     parse(Int_T, split(cd_key, "_")[2])
#     # end
#
#     const ary2dcv_1 =
#         cat(1,
#             map(values(calib_dict_1)) do value_1
#                 reshape(transpose(fluo_data), 1, size(value_1[1],2,1)...)
#             end...) ## do value_1
#     const (background_subtracted_data_1, k4dcv_2, deconvoluted_data_1,
#         norm_data_2, norm_well_nums_2, calibrated_data_1) =
#         calibrate(
#             ary2dcv_1,
#             db_conn_2,
#             calib_info_2,
#             well_nums_2,
#             well_nums_2,
#             channel_nums_1,
#             true,
#             dye_in,
#             dyes_to_fill;
#             data_format = array)
#     return CalibCalibOutput(
#         ary2dcv_1,
#         background_subtracted_data_1,
#         k4dcv_2,
#         deconvoluted_data_1,
#         norm_data_2,
#         calibrated_data_1)
# end ## calib_calib

## deprecated to remove MySql dependency
#
## get all the data from a calibration experiment,
## including data from all the channels for all the steps
# function get_full_calib_data(
#    db_conn::MySQL.MySQLHandle,
#    calib_info::OrderedDict,
#    well_nums::AbstractVector=[]
#    )
#
#    calib_info = ensure_ci(db_conn, calib_info)
#
#    calib_key_vec = get_ordered_keys(calib_info)
#    cd_key_vec = calib_key_vec[2:end] ## cd = channel of dye. "water" is index 1 per original order.
#    channel_nums = map(cd_key_vec) do cd_key
#        parse(Int_T, split(cd_key, "_")[2])
#    end
#    num_channels = length(channel_nums)
#
#    calib_dict = OrderedDict(map(calib_key_vec) do calib_key
#        exp_id = calib_info[calib_key]["calibration_id"]
#        step_id = calib_info[calib_key]["step_id"]
#        k_qry_2b = "
#            SELECT fluorescence_value, well_num, channel
#                FROM fluorescence_data
#                WHERE
#                    experiment_id = $exp_id AND
#                    step_id = $step_id AND
#                    cycle_num = 1 AND
#                    step_id is not NULL
#                    well_constraint
#                ORDER BY well_num, channel
#        "
#        calib_data_1key, calib_well_nums = get_mysql_data_well(
#            well_nums, k_qry_2b, db_conn, false
#        )
#        if length(well_nums) > 0 && Set(calib_well_nums) != Set(well_nums)
#            error(logger, "experiment $exp_id, step $step_id: " *
#                " calibration data is not found for all the wells requested")
#        end ## if
#        calib_data_1key_chwl = vcat(map(channel_nums) do channel
#            transpose(calib_data_1key[:fluorescence_value][calib_data_1key[:channel] .== channel])
#        end...) ## do channel
#            return an array where rows indexed by channels and columns indexed by wells
#
#        return calib_key => (calib_data_1key_chwl, calib_well_nums)
#    end)
#
#    return calib_dict ## share the same keys as `calib_info`
#
# end ## get_full_calib_data

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#  Functions and methods common to all MixedModel types

"""
    feL(m::MixedModel)

Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerTriangular`
`p × p` matrix.
"""
feL(m::MixedModel) = LowerTriangular(lmm(m).L.data.blocks[end - 1, end - 1])

lmm(m::LinearMixedModel) = m
lmm(m::GeneralizedLinearMixedModel) = m.LMM

"""
    cond(m::MixedModel)

Return the vector of the condition numbers of the blocks of `m.trms`
"""
LinearAlgebra.cond(m::MixedModel) = cond.(reterms(m))

"""
    std{T}(m::MixedModel{T})

Return the estimated standard deviations of the variance components as a `Vector{Vector{T}}`.
"""
function Statistics.std(m::MixedModel{T}) where {T}
    rl =  filter(!isempty, rowlengths.(lmm(m).trms))
    s = sdest(m)
    isfinite(s) ? rmul!(push!(rl, [1.]), s) : rl
end

fixefnames(m::MixedModel) = lmm(m).trms[end - 1].cnames
## methods for generics defined in StatsBase

function StatsBase.coeftable(m::MixedModel)
    co = coef(m)
    se = stderror(m)
    z = co ./ se
    pvalue = ccdf.(Chisq(1), abs2.(z))
    CoefTable(hcat(co, se, z, pvalue), ["Estimate", "Std.Error", "z value", "P(>|z|)"],
        fixefnames(m), 4)
end

"""
    describeblocks(io::IO, m::MixedModel)
    describeblocks(m::MixedModel)

Describe the types and sizes of the blocks in the lower triangle of `m.A` and `m.L`.
"""
function describeblocks(io::IO, m::MixedModel)
    lm = lmm(m)
    A = lm.A
    L = lm.L
    for i in 1 : nblocks(A, 2), j in 1 : i
        println(io, i, ",", j, ": ", typeof(A[Block(i, j)]), " ",
                blocksize(A, (i, j)), " ", typeof(L.data[Block(i, j)]))
    end
end
describeblocks(m::MixedModel) = describeblocks(stdout, m)

"""
    fnames(m::MixedModel)

Return the names of the grouping factors for the random-effects terms.
"""
fnames(m::MixedModel) = map(x -> x.fnm, reterms(m))

function getθ!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    trms = m.trms
    @argcheck(length(v) == sum(nθ, trms), DimensionMismatch)
    offset = 0
    for λ in trms
        nli = nθ(λ)
        getθ!(view(v, offset .+ (1 : nli)), λ)
        offset += nli
    end
    v
end

"""
    getΛ(m::MixedModel)

Return a vector of covariance template matrices for the random effects of `m`
"""
getΛ(m::MixedModel) = getΛ.(reterms(m))

"""
    getθ(m::MixedModel)

Return the current covariance parameter vector.
"""
getθ(m::MixedModel) = getθ(lmm(m).trms)

"""
    grplevels(m::MixedModel)

Return the number of levels in the random-effects terms' grouping factors.
"""
grplevels(m::MixedModel) = nlevs.(reterms(m))

nreterms(m::MixedModel) = sum(t -> isa(t, AbstractFactorReTerm), lmm(m).trms)

reterms(m::MixedModel) = convert(Vector{AbstractFactorReTerm},
                                 filter(t -> isa(t, AbstractFactorReTerm), lmm(m).trms))

"""
    ranef!{T}(v::Vector{Matrix{T}}, m::MixedModel{T}, β, uscale::Bool)

Overwrite `v` with the conditional modes of the random effects for `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise
on the original scale
"""
function ranef!(v::Vector, m::LinearMixedModel{T}, β::AbstractArray{T}, uscale::Bool) where {T}
    Ldat = m.L.data
    @argcheck((k = length(v)) == nreterms(m), DimensionMismatch)
    for j in 1:k
        mulαβ!(vec(copyto!(v[j], Ldat[Block(nblocks(Ldat, 2), j)])), Ldat[Block(k + 1, j)]', β,
            -one(T), one(T))
    end
    for i in k: -1 :1
        Lii = Ldat[Block(i, i)]
        vi = vec(v[i])
        ldiv!(adjoint(isa(Lii, Diagonal) ? Lii : LowerTriangular(Lii)), vi)
        for j in 1:(i - 1)
            mulαβ!(vec(v[j]), Ldat[Block(i, j)]', vi, -one(T), one(T))
        end
    end
    if !uscale
        trms = m.trms
        for (j, vv) in enumerate(v)
            lmul!(trms[j].Λ, vv)
        end
    end
    v
end

ranef!(v::Vector, m::LinearMixedModel, uscale::Bool) = ranef!(v, m, fixef(m), uscale)

"""
    ranef(m::MixedModel; uscale=false, named=true)

Return, as a `Vector{Vector{T}}` (`Vector{NamedVector{T}}` if `named=true`),
the conditional modes of the random effects in model `m`.

If `uscale` is `true` the random effects are on the spherical (i.e. `u`) scale, otherwise on
the original scale.
"""
function ranef(m::MixedModel; uscale=false, named=false)
    LMM = lmm(m)
    T = eltype(LMM.sqrtwts)
    retrms = reterms(LMM)
    v = Matrix{T}[Matrix{T}(undef, vsize(t), nlevs(t)) for t in retrms]
    ranef!(v, LMM, uscale)
    named || return v
    vnmd = map(NamedArray, v)
    for (trm, vnm) in zip(retrms, vnmd)
        setnames!(vnm, trm.cnms, 1)
        setnames!(vnm, string.(levs(trm)), 2)
    end
    vnmd
end

function StatsBase.vcov(m::MixedModel)
    Xtrm = lmm(m).trms[end - 1]
    iperm = invperm(Xtrm.piv)
    p = length(iperm)
    r = Xtrm.rank
    Linv = inv(feL(m))
    permvcov = varest(m) * (Linv'Linv)
    if p == Xtrm.rank
        permvcov[iperm, iperm]
    else
        T = eltype(permvcov)
        covmat = fill(zero(T)/zero(T), (p, p))
        for j in 1:r, i in 1:r
            covmat[i,j] = permvcov[i, j]
        end
        covmat[iperm, iperm]
    end
end

Statistics.cor(m::MixedModel{T}) where {T} = Matrix{T}[stddevcor(t)[2] for t in reterms(m)]

"""
    condVar(m::MixedModel)

Return the conditional variances matrices of the random effects.

The random effects are returned by `ranef` as a vector of length `k`,
where `k` is the number of random effects terms.  The `i`th element
is a matrix of size `vᵢ × ℓᵢ`  where `vᵢ` is the size of the
vector-valued random effects for each of the `ℓᵢ` levels of the grouping
factor.  Technically those values are the modes of the conditional
distribution of the random effects given the observed data.

This function returns an array of `k` three dimensional arrays,
where the `i`th array is of size `vᵢ × vᵢ × ℓᵢ`.  These are the
diagonal blocks from the conditional variance-covariance matrix,

    s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'
"""
function condVar(m::MixedModel)
    lm = lmm(m)
    λ = lm.trms[1]
    L11 = lm.L.data[Block(1, 1)]
    if nreterms(lm) ≠ 1 || !isa(L11, Diagonal{eltype(λ)})
        throw(ArgumentError("code for vector-valued r.e. or more than one term not yet written"))
    end
    ll = λ.Λ[1]
    Ld = L11.diag
    Array{eltype(Ld), 3}[reshape(abs2.(ll ./ Ld) .* varest(m), (1, 1, length(Ld)))]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule atomsk_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("atomsk")
JLLWrappers.@generate_main_file("atomsk", UUID("e89c7bda-2e69-5076-81ab-8e206d1a812b"))
end  # module atomsk_jll


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw2a.py

# import modules
import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    # create the input SpatialReference
    inspatialref = GU.fromEPSG(4269)
    # create the output SpatialReference
    outspatialref = GU.fromEPSG(26912)
    # create the CoordinateTransformation
    coordtrans = GU.createcoordtrans(inspatialref, outspatialref)

    GU.read("tmp/hw2a.shp") do input
        GU.create("tmp/hw2b.shp", "ESRI Shapefile") do output
            inlayer = GU.fetchlayer(input, 0)
            outlayer = GU.createlayer(output, "hw2b", GDAL.wkbPolygon)

            # get the FieldDefn for the county name field
            infeaturedefn = GU.getlayerdefn(inlayer)
            nameindex = GU.getfieldindex(infeaturedefn, "name")
            fielddefn = GU.fetchfielddefn(infeaturedefn, nameindex)

            # add the field to the output shapefile
            GU.createfield(outlayer, fielddefn)
            # get the FeatureDefn for the output shapefile
            outfeaturedefn = GU.getlayerdefn(outlayer)

            # loop through the input features
            for infeature in inlayer
                # get the input geometry
                geom = GU.getgeom(infeature)
                # reproject the geometry
                GU.transform(geom, coordtrans)
                # create a new feature
                GU.createfeature(outfeaturedefn) do outfeature
                    # set the geometry and attribute
                    GU.setgeom(outfeature, geom)
                    inname = GU.fetchfield(infeature, nameindex)
                    GU.setfield(outfeature, 0, inname)
                    # add the feature to the shapefile
                    GU.createfeature(outlayer, outfeature)
                end
            end
        end
    end

    # create the *.prj file
    GU.morphtoesri(outspatialref)
    open("tmp/hw2b.prj", "w") do file
        write(file, GU.toWKT(outspatialref))
    end
end


text:
module DifferentialExpression
export
    de_mannwhitney,
    DESeq2

using GLM: glm, coeftable, LogLink, NegativeBinomial
using HypothesisTests:MannWhitneyUTest
using MultipleTesting
using Statistics: median
using StatsBase: geomean

# scRNASeq Differential Gene Expression with two-sample
# MannWhitneyUTest

function de_mannwhitney(X,Y)
    assert(size(X,1)==size(Y,1))
    n=size(X,1)
    p=zeros(n,1)
    for k in 1:n
        p[k]=pvalue(MannWhitneyUTest(X[k,:],Y[k,:]))
    end
    q=adjust(vec(p), BenjaminiHochberg())
    return p,q
end

add_constant(X) = hcat(ones(size(X, 1)), X)
loge_to_log2(β) = β |> exp |> log2
median_of_ratios_normalization(X) = normalize_geomean(X) |> scale_size_factors
normalize_geomean(X) = X ./ mapslices(geomean, X, dims=1)
scale_size_factors(X) = X .* mapslices(median, X, dims=2)

struct DESeq2_results
    β
    SE
    z
    p
end

function get_DESeq2_results(models)
    # Needs to be fixed to use Log2 in the future
    β  = [coeftable(model).cols[1][2] for model in models]
    SE = [coeftable(model).cols[2][2] for model in models]
    z  = [coeftable(model).cols[3][2] for model in models]
    p  = [coeftable(model).cols[4][2] for model in models]
    DESeq2_results(β, SE, z, p)
end

function DESeq2(X, y)
    X_norm = median_of_ratios_normalization(X)
    models = mapslices(X -> glm(add_constant(X), y, NegativeBinomial(), LogLink()), X_norm, dims=1)
    get_DESeq2_results(models)
end

"""
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8#Sec2

Work in progress
----------------
Some low hanging fruit
1. Convert analysis to log2
2. Offer extra analysis/plots that DESeq2 provides

Challenging work
----------------
1. Ensuring this module behaves exactly like DESeq2
2. Going through the rest of the paper for details
"""

end


text:
include("../test/helper.jl")
include("../test/minitest.jl")
Minitest.do_tests() do
  Minitest.include("test_xpfsa.jl")
end