text:
using MathProgBase
using DataFrames
using Clp
using SparseArrays, LinearAlgebra

import Base.*, Base.-, Base.+, Base./, Base.max

export LinearProgrammingHall, LinearProgrammingShaft, LinearProgrammingRoom, LinearProgrammingHouse
export hallsingle, hall_relabel, hallvalues, hall_duplicate
export shaftsingle, shaftvalues
export roomdiagonal, roomdiagonalintersect, roomsingle, roomchunks, roomempty, roomintersect, room_relabel, room_relabel_parameter, room_duplicate
export varsum, fromindex, fromindexes, discounted
export setobjective!, setconstraint!, setconstraintoffset!, getconstraintoffset, clearconstraint!, setlower!, setupper!, addparameter!, addconstraint!
## Debugging
export getroom, getnonzerorooms
export gethouse, constraining, houseoptimize, summarizeparameters, findinfeasiblepair, varlengths, getconstraintsolution, getparametersolution

# A hallway is a vector of variables
"""
    LinearProgrammingHall

A vector of values, used either to describes the objective function or
the constraint offsets.

# Fields
* `component::Symbol`: The component defining a parameter or variable.
* `name::Symbol`: Either a parameter or variable name.
* `f::Vector{Float64}`: A vector of values, for every entry in the variable.
"""
mutable struct LinearProgrammingHall
    component::Symbol
    name::Symbol
    f::Vector{Float64}
end

"Construct a hall, calling gen with each index."
function hallsingle(model::Model, component::Symbol, name::Symbol, gen::Function, dupover::Vector{Symbol}=Symbol[])
    #println("$(now()) hs $component:$name")
    if isempty(dupover)
        LinearProgrammingHall(component, name, vectorsingle(getdims(model, component, name), gen))
    else
        dupover2 = Bool[dimname in dupover for dimname in getdimnames(model, component, name)]
        LinearProgrammingHall(component, name, vectorsingle(getdims(model, component, name), gen, dupover2))
    end
end

function hallsingle(model::Model, component::Symbol, name::Symbol, constant::Float64)
    LinearProgrammingHall(component, name, vectorsingle(getdims(model, component, name), constant))
end

function hallvalues(model::Model, component::Symbol, name::Symbol, values::Array{Float64})
    #println("$(now()) hv $component:$name")
    gen(inds...) = values[inds...]
    hallsingle(model, component, name, gen)
end

"""Connect a derivative to another component: change the variable component and name to another component."""
function hall_relabel(hall::LinearProgrammingHall, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert hall.name == from "Name mismatch in hall_relabel: $(hall.name) <> $from"

    LinearProgrammingHall(tocomponent, toname, hall.f)
end

"""Connect a derivative to another component while adding an additional dimension"""
function hall_duplicate(hall::LinearProgrammingHall, from::Symbol, tocomponent::Symbol, toname::Symbol, model::Model, newdims::Vector{Bool})
    @assert hall.name == from "Name mismatch in hall_duplicate: $(hall.name) <> $from"

    alldims = getdims(model, tocomponent, toname)

    f = zeros(prod(alldims))
    for kk in 1:prod(alldims[newdims])
        fulliis = expanddims(collect(1:length(hall.f)), alldims, newdims, toindex(kk, alldims[newdims]))
        f[fulliis] = hall.f
    end

    LinearProgrammingHall(tocomponent, toname, f)
end

function -(hall::LinearProgrammingHall)
    LinearProgrammingHall(hall.component, hall.name, -hall.f)
end

function *(hall::LinearProgrammingHall, mm::Number)
    LinearProgrammingHall(hall.component, hall.name, hall.f * mm)
end

function *(mm::Number, hall::LinearProgrammingHall)
    LinearProgrammingHall(hall.component, hall.name, mm * hall.f)
end

function /(hall::LinearProgrammingHall, dd::Number)
    LinearProgrammingHall(hall.component, hall.name, hall.f / dd)
end

"""
    hall1 + hall2

Construct a new LinearProgrammingHall which is the element-wise
sum of the values in hall1 and hall2, which must refer to
the same variable.
"""
function +(hall1::LinearProgrammingHall, hall2::LinearProgrammingHall; skipnamecheck=false)
    if !skipnamecheck
        @assert hall1.name == hall2.name "Hall + Hall name mismatch: $(hall1.name) <> $(hall2.name); use hall_relabel?"
    end
    LinearProgrammingHall(hall1.component, hall1.name, hall1.f + hall2.f)
end

"""
    hall1 - hall2

Construct a new LinearProgrammingHall which is the element-wise
difference in the values in hall1 and hall2, which must refer to
the same variable.
"""
function -(hall1::LinearProgrammingHall, hall2::LinearProgrammingHall; skipnamecheck=false)
    if !skipnamecheck
        @assert hall1.name == hall2.name "Hall - Hall name mismatch: $(hall1.name) <> $(hall2.name); use hall_relabel?"
    end
    LinearProgrammingHall(hall1.component, hall1.name, hall1.f - hall2.f)
end

"""
    max(hall1, hall2)

Construct a new LinearProgrammingHall which is the element-wise
maximum value of each entry in hall1 and hall2, which must refer to
the same variable.
"""
function max(hall1::LinearProgrammingHall, hall2::LinearProgrammingHall; skipnamecheck=false)
    if !skipnamecheck
        @assert hall1.name == hall2.name "Hall - Hall name mismatch: $(hall1.name) <> $(hall2.name); use hall_relabel?"
    end
    LinearProgrammingHall(hall1.component, hall1.name, max(hall1.f, hall2.f))
end

# A shaft is a transpose of a hallway, used for parameters
"""
    LinearProgrammingShaft

A vector of values, used either to describes the objective function or
the constraint offsets.  This acts as the transpose of a Hall,
although this distinction is only important for multiplying a Room by
Shaft, which then returns a Hall.

# Fields
* `component::Symbol`: The component defining a parameter or variable.
* `name::Symbol`: Either a parameter or variable name.
* `x::Vector{Float64}`: A vector of values, for every entry in the variable.
"""
mutable struct LinearProgrammingShaft
    component::Symbol
    name::Symbol
    x::Vector{Float64}
end

function shaftsingle(model::Model, component::Symbol, name::Symbol, gen::Function, dupover::Vector{Symbol}=Symbol[])
    #println("$(now()) ss $component:$name")
    if isempty(dupover)
        LinearProgrammingShaft(component, name, vectorsingle(getdims(model, component, name), gen))
    else
        dupover2 = Bool[dimname in dupover for dimname in getdimnames(model, component, name)]
        LinearProgrammingShaft(component, name, vectorsingle(getdims(model, component, name), gen, dupover2))
    end
end

function shaftsingle(model::Model, component::Symbol, name::Symbol, constant::Float64)
    LinearProgrammingShaft(component, name, vectorsingle(getdims(model, component, name), constant))
end

function shaftvalues(model::Model, component::Symbol, name::Symbol, values::Array{Float64})
    #println("$(now()) sv $component:$name")
    gen(inds...) = values[inds...]
    shaftsingle(model, component, name, gen)
end

function -(shaft::LinearProgrammingShaft)
    LinearProgrammingShaft(shaft.component, shaft.name, -shaft.x)
end

function +(shaft1::LinearProgrammingShaft, shaft2::LinearProgrammingShaft; skipnamecheck=false)
    if !skipnamecheck
        @assert shaft1.name == shaft2.name "Shaft + Shaft name mismatch: $(shaft1.name) <> $(shaft2.name); use shaft_relabel?"
    end
    LinearProgrammingShaft(shaft1.component, shaft1.name, shaft1.x + shaft2.x)
end

"""
    LinearProgrammingRoom

A matrix of values, used to describe a part of the linear programming
problem constraint matrix (a gradient matrix).  The rows of the matrix
correspond to variable indices, and the columns correspond to
parameter indices.  The values are stored as a sparse matrix, since
most parameter indices are assumed to not interact with most variable
indices.

# Fields
* `varcomponent::Symbol`: The component defining the variable.
* `variable::Symbol`: The variable name.
* `paramcomponent::Symbol`: The component defining the parameter.
* `parameter::Symbol`: The parameter name.
* `A::SparseMatrixCSC{Float64, Int64}`: A sparse matrix of values, for every combination of the parameter and variable.
"""
mutable struct LinearProgrammingRoom
    varcomponent::Symbol
    variable::Symbol
    paramcomponent::Symbol
    parameter::Symbol
    A::SparseMatrixCSC{Float64, Int64}
end

## Math to be incorporated below:
# ```math
# \begin{array}{ccc} p_1 & p_2 & \cdots \end{array}
# ```
# ```math
# \begin{array}{c}
#     v_1 \\
#     v_2 \\
#     \vdots
#     \end{array}\left(\begin{array}{ccc}
#         g(1) & 0 & \cdots \\
#         0 & g(2) & \cdots \\
#         \vdots & \vdots & \ddots
#         \end{array}\right)
# ```
"""
    roomdiagonal(model, component, variable, parameter, gen)

Fill in just the diagonal, assuming ``\frac{dv_i}{dp_j} = 0``, if ``i <>
j``.  Requires that the dimensions of `variable` and `parameter` are
the same.  `gen` called for each combination of indices.

# Arguments
* `model::Model`: The model containing `component`.
* `component::Symbol`: The component containing `variable` and `parameter`.
* `variable::Symbol`: The outcome variable, corresponding to matrix rows.
* `parameter::Symbol`: The optimization parameter, corresponding to matrix columns.
* `gen::Function`: The function generating gradient values.
"""
function roomdiagonal(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function, dupover::Vector{Symbol}=Symbol[]; diagdupover::Vector{Symbol}=Symbol[])
    #println("$(now()) rd $component:$variable x $parameter")
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)

    @assert (isempty(dupover) || isempty(diagdupover)) "No current support for both dupover and diagdupover."

    if isempty(dupover) && isempty(diagdupover)
        @assert dimsvar == dimspar "Variable and parameter in roomdiagonal do not have the same dimensions: $component.$variable $dimsvar <> $component.$parameter $dimspar"

        LinearProgrammingRoom(component, variable, component, parameter, matrixdiagonal(dimsvar, gen))
    elseif !isempty(dupover)
        dupoverpar = Bool[dimname in dupover for dimname in getdimnames(model, component, parameter)]
        dupovervar = Bool[dimname in dupover for dimname in getdimnames(model, component, variable)]

        @assert dimsvar[.!dupovervar] == dimspar[.!dupoverpar] "Variable and parameter in roomdiagonal do not have the same dimensions: $component.$variable $dimsvar <> $component.$parameter $dimspar, after dup-ignore $dupover"

        matrix = matrixdiagonal(dimsvar[.!dupovervar], gen)
        matrix2 = matrixduplicate(matrix, dimsvar, dimspar, dupovervar, dupoverpar)

        LinearProgrammingRoom(component, variable, component, parameter, matrix2)
    elseif !isempty(diagdupover)
        @assert dimsvar == dimspar "Variable and parameter in roomdiagonal must be identifical when using diagdupover"
        dupoverparvar = Bool[dimname in diagdupover for dimname in getdimnames(model, component, parameter)]

        matrix = matrixdiagonal(dimsvar[.!dupoverparvar], gen)
        matrix2 = matrixduplicate_diagonal(matrix, dimspar, dupoverparvar)

        LinearProgrammingRoom(component, variable, component, parameter, matrix2)
    end
end

function roomdiagonal(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, constant::Float64, dupover::Vector{Symbol}=Symbol[]; diagdupover::Vector{Symbol}=Symbol[])
    #println("$(now()) rd $component:$variable x $parameter")
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)

    @assert (isempty(dupover) || isempty(diagdupover)) "No current support for both dupover and diagdupover."

    if isempty(dupover) && isempty(diagdupover)
        @assert dimsvar == dimspar "Variable and parameter in roomdiagonal do not have the same dimensions: $component.$variable $dimsvar <> $component.$parameter $dimspar"

        LinearProgrammingRoom(component, variable, component, parameter, matrixdiagonal(dimsvar, constant))
    elseif !isempty(dupover)
        dupoverpar = Bool[dimname in dupover for dimname in getdimnames(model, component, parameter)]
        dupovervar = Bool[dimname in dupover for dimname in getdimnames(model, component, variable)]

        @assert dimsvar[.!dupovervar] == dimspar[.!dupoverpar] "Variable and parameter in roomdiagonal do not have the same dimensions: $component.$variable $dimsvar <> $component.$parameter $dimspar, after dup-ignore $dupover"

        matrix = matrixdiagonal(dimsvar[.!dupovervar], constant)
        matrix2 = matrixduplicate(matrix, dimsvar, dimspar, dupovervar, dupoverpar)
        LinearProgrammingRoom(component, variable, component, parameter, matrix2)
    elseif !isempty(diagdupover)
        @assert dimsvar == dimspar "Variable and parameter in roomdiagonal must be identifical when using diagdupover"
        dupoverparvar = Bool[dimname in diagdupover for dimname in getdimnames(model, component, parameter)]

        matrix = matrixdiagonal(dimsvar[.!dupoverparvar], constant)
        matrix2 = matrixduplicate_diagonal(matrix, dimspar, dupoverparvar)

        LinearProgrammingRoom(component, variable, component, parameter, matrix2)
    end
end

"""
gen returns a full matrix for all shared variables, in variable dimension order.
"""
function roomdiagonalintersect(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function)
    #println("$(now()) rd $component:$variable x $parameter")
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)
    dimnamesvar = getdimnames(model, component, variable)
    dimnamespar = getdimnames(model, component, parameter)

    A = matrixdiagonalintersect(dimsvar, dimspar, dimnamesvar, dimnamespar, gen)
    LinearProgrammingRoom(component, variable, component, parameter, A)
end

"""
Fill in every element
"""
function roomsingle(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function, vardupover::Vector{Symbol}=Symbol[], pardupover::Vector{Symbol}=Symbol[])
    #println("$(now()) rs $component:$variable x $parameter")
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)

    if isempty(vardupover) && isempty(pardupover)
        LinearProgrammingRoom(component, variable, component, parameter, matrixsingle(dimsvar, dimspar, gen))
    else
        vardupover2 = Bool[dimname in vardupover for dimname in getdimnames(model, component, variable)]
        pardupover2 = Bool[dimname in pardupover for dimname in getdimnames(model, component, parameter)]
        LinearProgrammingRoom(component, variable, component, parameter, matrixsingle(dimsvar, dimspar, gen, vardupover2, pardupover2))
    end
end

"""
Call with matrices for each combination of chunk variables
    Chunked variables must be consecutive at the end.
"""
function roomchunks(model::Model, component::Symbol, variable::Symbol, parameter::Symbol, gen::Function, varchunk::Vector{Symbol}, parchunk::Vector{Symbol})
    #println("$(now()) rc $component:$variable x $parameter")
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)

    dimvarnames = getdimnames(model, component, variable)
    dimparnames = getdimnames(model, component, parameter)

    varchunk2 = 0
    while varchunk2 < length(dimsvar)
        if !(dimvarnames[end - varchunk2] in varchunk)
            break
        end
        varchunk2 += 1
    end

    parchunk2 = 0
    while parchunk2 < length(dimspar)
        if !(dimparnames[end - parchunk2] in parchunk)
            break
        end
        parchunk2 += 1
    end

    LinearProgrammingRoom(component, variable, component, parameter, matrixchunks(dimsvar, dimspar, gen, varchunk2, parchunk2))
end

"""
    roomintersect(model, component, variable1, variable2, gen)

Generate a room at the intersection of two variables.  Variable1 is
the constraint, and variable2 has the optimization parameter.

Call `gen` for each shared index, passing an array to be filled for unshared indexes.
This version assumes both variables come from the same component.

See `matrixintersect` for the matrix generation logic.
"""
function roomintersect(model::Model, component::Symbol, variable1::Symbol, variable2::Symbol, gen::Function, dupover1::Vector{Symbol}=Symbol[], dupover2::Vector{Symbol}=Symbol[])
    #println("$(now()) ri $component:$variable1 x $variable2")
    dims1 = getdims(model, component, variable1)
    dims2 = getdims(model, component, variable2)
    dimnames1 = getdimnames(model, component, variable1)
    dimnames2 = getdimnames(model, component, variable2)

    if isempty(dupover1) && isempty(dupover2)
        LinearProgrammingRoom(component, variable1, component, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen))
    else
        dupover12 = Bool[dimname in dupover1 for dimname in getdimnames(model, component, variable1)]
        dupover22 = Bool[dimname in dupover2 for dimname in getdimnames(model, component, variable2)]
        LinearProgrammingRoom(component, variable1, component, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen, dupover12, dupover22))
    end
end

"""
    roomintersect(model, component, variable1, component2, variable2, gen)

Generate a room at the intersection of two variables.

Call `gen` for each shared index, passing an array to be filled for unshared indexes.
This version allows the variables to come from different component.

See `matrixintersect` for the matrix generation logic.
"""
function roomintersect(model::Model, component1::Symbol, variable1::Symbol, component2::Symbol, variable2::Symbol, gen::Function, dupover1::Vector{Symbol}=Symbol[], dupover2::Vector{Symbol}=Symbol[])
    #println("$(now()) ri $component1:$variable1 x $component2:$variable2")
    dims1 = getdims(model, component1, variable1)
    dims2 = getdims(model, component2, variable2)
    dimnames1 = getdimnames(model, component1, variable1)
    dimnames2 = getdimnames(model, component2, variable2)

    if isempty(dupover1) && isempty(dupover2)
        LinearProgrammingRoom(component1, variable1, component2, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen))
    else
        dupover12 = Bool[dimname in dupover1 for dimname in getdimnames(model, component1, variable1)]
        dupover22 = Bool[dimname in dupover2 for dimname in getdimnames(model, component2, variable2)]
        LinearProgrammingRoom(component1, variable1, component2, variable2, matrixintersect(dims1, dims2, dimnames1, dimnames2, gen, dupover12, dupover22))
    end
end

"""
Create a room to be filled in later
"""
function roomempty(model::Model, component::Symbol, variable::Symbol, parameter::Symbol)
    dimsvar = getdims(model, component, variable)
    dimspar = getdims(model, component, parameter)
    LinearProgrammingRoom(component, variable, component, parameter, matrixempty(dimsvar, dimspar))
end

"""Connect a gradient to another component: change the variable component and name to another component."""
function room_relabel(room::LinearProgrammingRoom, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert room.variable == from "Variable name mismatch in room_relabel: $(room.variable) <> $from"

    LinearProgrammingRoom(tocomponent, toname, room.paramcomponent, room.parameter, room.A)
end

"""Connect a gradient to another component: change the variable component and name to another component."""
function room_relabel_parameter(room::LinearProgrammingRoom, from::Symbol, tocomponent::Symbol, toname::Symbol)
    @assert room.parameter == from "Parameter name mismatch in room_relabel_parameter: $(room.variable) <> $from"

    LinearProgrammingRoom(room.varcomponent, room.variable, tocomponent, toname, room.A)
end

"""Connect a derivative to another component while adding an additional dimension"""
function room_duplicate(room::LinearProgrammingRoom, tovariable::Symbol, toparameter::Symbol, model::Model)
    parolddimnames = getdimnames(model, room.paramcomponent, room.parameter)
    parnewdimnames = getdimnames(model, room.paramcomponent, toparameter)
    @assert issubset(parolddimnames, parnewdimnames)
    varolddimnames = getdimnames(model, room.varcomponent, room.variable)
    varnewdimnames = getdimnames(model, room.varcomponent, tovariable)
    @assert issubset(varolddimnames, varnewdimnames)

    rowdupover = Bool[]
    jj = 1
    for ii in 1:length(varnewdimnames)
        if jj <= length(varolddimnames) && varnewdimnames[ii] == varolddimnames[jj]
            jj += 1
            push!(rowdupover, false)
        else
            push!(rowdupover, true)
        end
    end
    @assert jj == length(varolddimnames) + 1 "Got $jj after matching $varolddimnames in $varnewdimnames for $rowdupover"

    coldupover = Bool[]
    jj = 1
    for ii in 1:length(parnewdimnames)
        if jj <= length(parolddimnames) && parnewdimnames[ii] == parolddimnames[jj]
            jj += 1
            push!(coldupover, false)
        else
            push!(coldupover, true)
        end
    end
    @assert jj == length(parolddimnames) + 1 "Got $jj after matching $parolddimnames in $parnewdimnames for $coldupover"

    A = matrixduplicate(room.A, getdims(model, room.varcomponent, tovariable), getdims(model, room.paramcomponent, toparameter), rowdupover, coldupover)
    LinearProgrammingRoom(room.varcomponent, tovariable, room.paramcomponent, toparameter, A)
end

"""
    Construct a new LinearProgrammingHall, with each entry as a weighted sum of the partials for a given variable.
        This is equivalent to translating a constraint by a functional relationship.
"""
function *(hall::LinearProgrammingHall, room::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert hall.name == room.variable "Hall * Room name mismatch: $(hall.name) <> $(room.variable); use room_relabel?"
    end
    LinearProgrammingHall(room.paramcomponent, room.parameter, vec(transpose(hall.f) * room.A))
end

# Returns a hall, since it then contains variable/constraint information
function *(room::LinearProgrammingRoom, shaft::LinearProgrammingShaft; skipnamecheck=false)
    if !skipnamecheck
        @assert shaft.name == room.parameter "Room * Shaft name mismatch: $(room.parameter) <> $(hall.name); use room_relabel?"
    end
    LinearProgrammingHall(room.varcomponent, room.variable, room.A * shaft.x)
end

function *(room1::LinearProgrammingRoom, room2::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert room1.parameter == room2.variable "Room * Room name mismatch: $(room1.parameter) <> $(room2.variable); use room_relabel?"
    end
    LinearProgrammingRoom(room1.varcomponent, room1.variable, room2.paramcomponent, room2.parameter, room1.A * room2.A)
end

function +(room1::LinearProgrammingRoom, room2::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert room1.parameter == room2.parameter "Room + Room parameter name mismatch: $(room1.parameter) <> $(room2.parameter); use room_relabel?"
        @assert room1.variable == room2.variable "Room + Room variable name mismatch: $(room1.variable) <> $(room2.variable); use room_relabel?"
    end
    if size(room1.A)[1] == 0 || size(room1.A)[2] == 0
        LinearProgrammingRoom(room1.varcomponent, room1.variable, room1.paramcomponent, room1.parameter, spzeros(size(room1.A)...))
    else
        LinearProgrammingRoom(room1.varcomponent, room1.variable, room1.paramcomponent, room1.parameter, room1.A + room2.A)
    end
end

function -(room1::LinearProgrammingRoom, room2::LinearProgrammingRoom; skipnamecheck=false)
    if !skipnamecheck
        @assert room1.parameter == room2.parameter "Room - Room parameter name mismatch: $(room1.parameter) <> $(room2.parameter); use room_relabel?"
        @assert room1.variable == room2.variable "Room - Room variable name mismatch: $(room1.variable) <> $(room2.variable); use room_relabel?"
    end
    if size(room1.A)[1] == 0 || size(room1.A)[2] == 0
        LinearProgrammingRoom(room1.varcomponent, room1.variable, room1.paramcomponent, room1.parameter, spzeros(size(room1.A)...))
    else
        LinearProgrammingRoom(room1.varcomponent, room1.variable, room1.paramcomponent, room1.parameter, room1.A - room2.A)
    end
end

function -(room::LinearProgrammingRoom)
    LinearProgrammingRoom(room.varcomponent, room.variable, room.paramcomponent, room.parameter, -room.A)
end

"""
Construct a hall from a room by summing over rows.
        This is equivalent to a constraint for the sum of the variables in the original room.
"""
function varsum(room::LinearProgrammingRoom)
    LinearProgrammingHall(room.paramcomponent, room.parameter, vec(sum(room.A, dims=1)))
end

function varsum(room::LinearProgrammingRoom, axis::Int64, model::Model, newvar::Symbol)
    vardims = getdims(model, room.varcomponent, room.variable)
    iis, jjs, vvs = findnz(room.A)

    indexes = toindexes(iis, vardims)
    newiis = fromindexes(indexes[:, 1:length(vardims) .!= axis], vardims[1:length(vardims) .!= axis])

    useiis = []
    usejjs = []
    usevvs = []
    for newii in unique(newiis)
        withinii = newiis .== newii
        for newjj in unique(jjs[withinii])
            withinij = jjs[withinii] .== newjj
            push!(useiis, newii)
            push!(usejjs, newjj)
            push!(usevvs, sum(vvs[withinii][withinij]))
        end
    end

    A = sparse(useiis, usejjs, usevvs, prod(vardims[1:length(vardims) .!= axis]), size(room.A, 2))
    LinearProgrammingRoom(room.varcomponent, newvar, room.paramcomponent, room.parameter, A)
end

"""
Discount economic costs in room.  Only rows are discounted (i.e., doesn't matter when action made, just when cost incurred).
"""
function discounted(model::Model, room::LinearProgrammingRoom, rate::Float64)
    dims = getdims(model, room.varcomponent, room.variable)

    dupover2 = Bool[dimname != :time for dimname in getdimnames(model, room.varcomponent, room.variable)]
    outers, inners = interpretdupover(dupover2)
    dupouter = prod(dims[outers])
    dupinner = prod(dims[inners])

    timevalues = collect(1:getindexcount(model, :time))
    alltime = repeat(timevalues, inner=[dupinner], outer=[dupouter])

    discount = exp.(-rate * (alltime - 1))
    # discounted = room.A .* discount
    iis, jjs, vvs = findnz(room.A)
    vvs2 = vvs .* discount[iis]
    discounted = sparse(iis, jjs, vvs2, size(room.A, 1), size(room.A, 2))

    LinearProgrammingRoom(room.varcomponent, room.variable, room.paramcomponent, room.parameter, discounted)
end

## Math to be incorporated below:
# The linear programming that is solved is always:
# ```math
# \max f' x
# ```
# ```math
# A x \le b
# ```
# ```math
# x_{lower} \le x \le x_{upper}
# ```
"""
    LinearProgrammingHouse

The full description of a linear programming problem, including all
its variables, parameters, the constraints, and the objective.

# Fields
* `model::Model`: The model containing all components
* `paramcomps::Vector{Symbol}`: The components defining each variable.
* `parameters::Vector{Symbol}`: The names of each variable.
* `constcomps::Vector{Symbol}`: The components defining each parameter.
* `constraints::Vector{Symbol}`: The names of each parameter.
* `namedictionary::Dict{Symbol, Symbol}`: The names used in `parameters` and `constraints` must be unique, but can refer to the same parameter or variable by including an entry in this disctionary mapping the unique name to the true variable name.
* `lowers::Vector{Float64}`: The lower bound for each parameter.
* `uppers::Vector{Float64}`: The upper bound for each parameter.
* `f::Vector{Float64}`: The derivative of the objective function for each parameter.
* `A::SparseMatrixCSC{Float64, Int64}`: Each row describes the derivatives of a given row for each parameter.
* `b::Vector{Float64}`: The maximum value for ``A x`` for parameter values ``x``.
"""
mutable struct LinearProgrammingHouse
    model::Model
    paramcomps::Vector{Symbol}
    parameters::Vector{Symbol}
    constcomps::Vector{Symbol}
    constraints::Vector{Symbol}
    namedictionary::Dict{Symbol, Symbol}
    lowers::Vector{Float64}
    uppers::Vector{Float64}
    f::Vector{Float64} # length of parameters
    A::SparseMatrixCSC{Float64, Int64} # constraint rows, parameter columns
    b::Vector{Float64} # length of constraints
end

function LinearProgrammingHouse(model::Model, paramcomps::Vector{Symbol}, parameters::Vector{Symbol}, constcomps::Vector{Symbol}, constraints::Vector{Symbol}, namedictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    paramlen = sum(varlengths(model, paramcomps, parameters, namedictionary))
    variablelen = sum(varlengths(model, constcomps, constraints, namedictionary))
    A = spzeros(variablelen, paramlen)
    f = zeros(paramlen)
    b = zeros(variablelen)
    LinearProgrammingHouse(model, paramcomps, parameters, constcomps, constraints, namedictionary, zeros(length(f)), Inf * ones(length(f)), f, A, b)
end

function addparameter!(house::LinearProgrammingHouse, component::Symbol, parameter::Symbol)
    paramlen = prod(getdims(house.model, component, parameter))

    append!(house.paramcomps, [component])
    append!(house.parameters, [parameter])
    append!(house.lowers, zeros(paramlen))
    append!(house.uppers, Inf * ones(paramlen))
    append!(house.f, zeros(paramlen))
    house.A = [house.A spzeros(length(house.b), paramlen)]
end

function addconstraint!(house::LinearProgrammingHouse, component::Symbol, constraint::Symbol, variable::Symbol)
    varlen = prod(getdims(house.model, component, variable))

    append!(house.constcomps, [component])
    append!(house.constraints, [constraint])
    house.namedictionary[constraint] = variable
    append!(house.b, zeros(varlen))
    house.A = [house.A; spzeros(varlen, length(house.f))]
end

function setobjective!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) .& (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.f[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.f[before+ii] = hall.f[ii]
    end
end

function setconstraint!(house::LinearProgrammingHouse, room::LinearProgrammingRoom)
    kk = findfirst((house.paramcomps .== room.paramcomponent) .& (house.parameters .== room.parameter))
    @assert kk > 0 "$(room.paramcomponent).$(room.parameter) not a known parameter"
    ll = findfirst((house.constcomps .== room.varcomponent) .& (house.constraints .== room.variable))
    @assert ll > 0 "$(room.varcomponent).$(room.variable) not a known variable"

    paramspans = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)
    @assert size(room.A, 1) == constspans[ll] "Length of variable $(room.variable) unexpected: $(size(room.A, 1)) given <> $(constspans[ll])"
    @assert size(room.A, 2) == paramspans[kk] "Length of parameter $(room.parameter) unexpected: $(size(room.A, 2)) given <> $(paramspans[kk])"

    if size(room.A, 2) == 0 # expressions below mishandle this case
        return
    end

    parambefore = sum(paramspans[1:kk-1])
    constbefore = sum(constspans[1:ll-1])

    house.A[constbefore+1:constbefore+constspans[ll], parambefore+1:parambefore+paramspans[kk]] = room.A
end

function getroom(house::LinearProgrammingHouse, varcomponent::Symbol, variable::Symbol, paramcomponent::Symbol, parameter::Symbol)
    @assert parameter in house.parameters "$parameter not a known parameter"
    @assert variable in house.constraints "$variable not a known variable"

    # Determine where to index into vector
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)
    rowkk = findfirst((house.constcomps .== varcomponent) .& (house.constraints .== variable))
    rowbefore = sum(constspans[1:rowkk-1])

    paramspans = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    colkk = findfirst((house.paramcomps .== paramcomponent) .& (house.parameters .== parameter))
    colbefore = sum(paramspans[1:colkk-1])

    # Collect the offset
    subA = house.A[(rowbefore + 1):(rowbefore + constspans[rowkk]), (colbefore + 1):(colbefore + paramspans[colkk])]

    LinearProgrammingRoom(varcomponent, variable, paramcomponent, parameter, subA)
end

"""Return a list of all rooms (parameter+component and variable+component) that are related by the derivative matrix.
"""
function getnonzerorooms(house::LinearProgrammingHouse)
    df = DataFrame(varcomponent=Symbol[], variable=Symbol[], paramcomponent=Symbol[], parameter=Symbol[])
    for kk in 1:length(house.paramcomps)
        for ll in 1:length(house.constcomps)
            room = getroom(house, house.constcomps[ll], house.constraints[ll], house.paramcomps[kk], house.parameters[kk])
            if !isempty(nonzeros(room.A))
                push!(df, [house.constcomps[ll], house.constraints[ll], house.paramcomps[kk], house.parameters[kk]])
            end
        end
    end

    df
end

"""
Set offset values from a `LinearProgrammingHall`.  See the other `setconstraintoffset!`
"""
function setconstraintoffset!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    setconstraintoffset!(house, hall.component, hall.name, hall.f)
end

"""
    setconstraintoffset!(house, component, variable, f)

Set offset values within a `LinearProgrammingHouse`.

# Arguments
* `house::LinearProgrammingHouse`: The house to set values within.
* `component::Symbol`: The component for the constraint variable.
* `variable::Symbol`: The variable for the constraint.
* `f::Vector{Float64}`: The values, with all dimensions collapsed into a single vector.
"""
function setconstraintoffset!(house::LinearProgrammingHouse, component::Symbol, variable::Symbol, f::Vector{Float64})
    @assert variable in house.constraints "$(variable) not a known variable"
    kk = findfirst((house.constcomps .== component) .& (house.constraints .== variable))
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)
    @assert length(f) == constspans[kk] "Length of parameter $(variable) unexpected: $(length(f)) <> $(constspans[kk])"
    before = sum(constspans[1:kk-1])
    for ii in 1:constspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.b[before+ii] = f[ii]
    end
end

"""
    getconstraintoffset(house, component, variable, reshp)

Return the values for a constraint, optionally reshaped to the original dimensions.

# Arguments
* `house::LinearProgrammingHouse`: The house from which to get the values.
* `component::Symbol`: The component for the constraint variable.
* `variable::Symbol`: The variable for the constraint.
* `reshp::Bool`: Should it be reshaped to the original variable dimensions? (default: false)
"""
function getconstraintoffset(house::LinearProgrammingHouse, component::Symbol, variable::Symbol; reshp::Bool=false)
    @assert variable in house.constraints "$variable not a known variable"

    # Determine where to index into vector
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)
    kk = findfirst((house.constcomps .== component) .& (house.constraints .== variable))
    before = sum(constspans[1:kk-1])

    # Collect the offset
    constraintoffset = house.b[(before + 1):(before + constspans[kk])]

    if reshp
        # Reshape to original dimensions
        reshape(constraintoffset, getdims(house.model, component, variable)...)
    else
        constraintoffset
    end
end

function clearconstraint!(house::LinearProgrammingHouse, component::Symbol, variable::Symbol)
    @assert variable in house.constraints "$(variable) not a known variable"

    ll = findfirst((house.constcomps .== component) .& (house.constraints .== variable))
    constspans = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)

    constbefore = sum(constspans[1:ll-1])

    house.A[constbefore+1:constbefore+constspans[ll], :] = 0
    house.b[constbefore+1:constbefore+constspans[ll]] = 0
end

function setlower!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) .& (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.lowers[before+ii] = hall.f[ii]
    end
end

function setupper!(house::LinearProgrammingHouse, hall::LinearProgrammingHall)
    @assert hall.name in house.parameters "$(hall.name) not a known parameter"
    kk = findfirst((house.paramcomps .== hall.component) .& (house.parameters .== hall.name))
    paramspans = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    @assert length(hall.f) == paramspans[kk] "Length of parameter $(hall.name) unexpected: $(length(hall.f)) <> $(paramspans[kk])"
    before = sum(paramspans[1:kk-1])
    for ii in 1:paramspans[kk]
        #@assert house.b[before+ii] == 0 "Overwrite existing gradient in setobjective!"
        house.uppers[before+ii] = hall.f[ii]
    end
end

function gethouse(house::LinearProgrammingHouse, rr::Int64, cc::Int64)
    # Determine the row and column names
    varii = findlast(cumsum(varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)) .< rr) + 1
    parii = findlast(cumsum(varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)) .< cc) + 1
    rrrelative = rr - sum(varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)[1:varii-1])
    ccrelative = cc - sum(varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)[1:parii-1])
    vardims = getdims(house.model, house.constcomps[varii], get(house.namedictionary, house.constraints[varii], house.constraints[varii]))
    pardims = getdims(house.model, house.paramcomps[parii], house.parameters[parii])

    println("$(house.constcomps[varii]).$(house.constraints[varii])$(toindex(rrrelative, vardims)), $(house.paramcomps[parii]).$(house.parameters[parii])$(toindex(ccrelative, pardims)) = $(house.A[rr, cc])")
end

function constraining(house::LinearProgrammingHouse, solution::Vector{Float64}; subset=[])
    # Determine which constraint (if any) is stopping an increase or decrease of each
    if subset == []
        df = DataFrame(solution=solution)
        df[:component] = :na
        df[:parameter] = :na
    else
        df = DataFrame(solution=solution[subset])
    end
    df[:abovefail] = ""
    df[:belowfail] = ""

    # Produce names for all constraints
    varlens = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)
    names = ["" for ii in 1:sum(varlens)]
    for kk in 1:length(house.constcomps)
        ii0 = sum(varlens[1:kk-1])
        for ii in 1:varlens[kk]
            names[ii0 + ii] = "$(house.constraints[kk]).$ii"
        end
    end

    varlens = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    baseconsts = house.A * solution

    println("Ignore:")
    println(join(names[findall(baseconsts .> house.b)], ", "))
    ignore = baseconsts .> house.b

    if subset == []
        for kk in 1:length(house.paramcomps)
            println(kk / length(house.paramcomps))
            ii0 = sum(varlens[1:kk-1])
            for ii in 1:varlens[kk]
                df[ii0 + ii, :component] = house.paramcomps[kk]
                df[ii0 + ii, :parameter] = house.parameters[kk]

                newconst = baseconsts + house.A[:, ii0 + ii] * 1e-6 + (house.A[:, ii0 + ii] .> 0) * 1e-6
                df[ii0 + ii, :abovefail] = join(names[findall((newconst .> house.b) .& .!ignore)], ", ")

                newconst = baseconsts - house.A[:, ii0 + ii] * 1e-6 - (house.A[:, ii0 + ii] .> 0) * 1e-6
                df[ii0 + ii, :belowfail] = join(names[findall((newconst .> house.b) .& .!ignore)], ", ")
            end
        end
    else
        for ii in 1:length(subset)
            println(ii / length(subset))

            newconst = baseconsts + house.A[:, subset[ii]] * 1e-6 + 1e-6
            df[ii, :abovefail] = join(names[findall((newconst .> house.b) .& !ignore)], ", ")

            newconst = baseconsts - house.A[:, subset[ii]] * 1e-6 - 1e-6
            df[ii, :belowfail] = join(names[findall((newconst .> house.b) .& !ignore)], ", ")
        end
    end

    df
end

houseoptimize(house::LinearProgrammingHouse) = linprog(-house.f, house.A, '<', house.b, house.lowers, house.uppers, ClpSolver())
houseoptimize(house::LinearProgrammingHouse, solver) = linprog(-house.f, house.A, '<', house.b, house.lowers, house.uppers, solver)
houseoptimize(house::LinearProgrammingHouse, solver, subset::Vector{Int64}) = linprog(-house.f, house.A[subset, :], '<', house.b[subset], house.lowers, house.uppers, solver)

"""
    findinfeasiblepair(house, solver)

Finds a range within the matrix for which the results become minimally
infeasible.  In other words, suppose that the full linear programming
matrix is ``A``.  It returns ``i``, ``j``, such that ``A[1:i, :]`` is
infeasible, but ``A[1:i-1, :]`` is not, and ``A[j:end, :]`` is infeasible
but ``A[j-1:end, :]`` is not.

# Arguments
* `house::LinearProgrammingHouse`: An infeasible LinearProgrammingHouse.
* `solver`: A solver object which finds the infeasibility.
"""
function findinfeasiblepair(house::LinearProgrammingHouse, solver)
    top = findinfeasiblepairhelper(house, solver, "top", 2, length(house.b))
    bottom = findinfeasiblepairhelper(house, solver, "bottom", 1, length(house.b) - 1)

    [top, bottom]
end

"""
Look for the row that makes the problem feasible
"""
function findinfeasiblepairhelper(house::LinearProgrammingHouse, solver, direction, searchtop, searchbottom)
    searchpoint = round(Int, (searchtop + searchbottom) / 2 + .01) # Always round up
    println([direction, searchtop, searchbottom])
    isfeasible = checkfeasibility(house, solver, direction, searchpoint)

    if !isfeasible
        # Look for a shorter span
        if searchpoint == searchbottom
            if direction == "bottom" || checkfeasibility(house, solver, direction, searchtop)
                searchbottom
            else
                searchtop
            end
        elseif direction == "top"
            findinfeasiblepairhelper(house, solver, direction, searchtop, searchpoint)
        else
            findinfeasiblepairhelper(house, solver, direction, searchpoint, searchbottom)
        end
    else
        # Look for a longer span
        if searchpoint == searchbottom
            searchtop
        elseif (direction == "top")
            findinfeasiblepairhelper(house, solver, direction, searchpoint, searchbottom)
        else
            findinfeasiblepairhelper(house, solver, direction, searchtop, searchpoint)
        end
    end
end

function checkfeasibility(house::LinearProgrammingHouse, solver, direction, searchpoint)
    if (direction == "top")
        sol = linprog(-house.f, house.A[1:searchpoint, :], '<', house.b[1:searchpoint], house.lowers, house.uppers, solver)
    else
        sol = linprog(-house.f, house.A[searchpoint:end, :], '<', house.b[searchpoint:end], house.lowers, house.uppers, solver)
    end

    sol.status != :Infeasible
end

function summarizeparameters(house::LinearProgrammingHouse, solution::Vector{Float64})
    # Look at parameter values
    varlens = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)
    for ii in 1:length(house.parameters)
        println(house.parameters[ii])
        index1 = sum(varlens[1:ii-1]) + 1
        index2 = sum(varlens[1:ii])

        values = solution[index1:index2]

        if (sum(values .!= 0) == 0)
            println("All zero.")
        else
            println(values[1:min(100, index2 - index1 + 1)])
            println("Sum: $(sum(values))")
        end
    end
end

"""
    getparametersolution(house, solution, parameter)

Return the array of solution values for a given parameter.
"""
function getparametersolution(house::LinearProgrammingHouse, solution::Vector{Float64}, parameter::Symbol)
    @assert parameter in house.parameters "Unknown parameter."

    varlens = varlengths(house.model, house.paramcomps, house.parameters, house.namedictionary)

    ii = findall(house.parameters .== parameter)[1]

    index1 = sum(varlens[1:ii-1]) + 1
    index2 = sum(varlens[1:ii])

    solution[index1:index2]
end

"""
    getparametersolution(house, sol, constraint)

Return the array of solution values for a given constraint.
"""
function getconstraintsolution(house, sol, constraint)
    constvalues = house.A * sol.sol

    varlens = varlengths(house.model, house.constcomps, house.constraints, house.namedictionary)

    ii = findall(house.constraints .== constraint)[1]

    index1 = sum(varlens[1:ii-1]) + 1
    index2 = sum(varlens[1:ii])

    constvalues[index1:index2]
end

## Helpers

function rangeof(m::Model, name, components, names, vardictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    varlens = varlengths(m, components, names, vardictionary)
    kk = findfirst(name .== names)
    sum(varlens[1:kk-1])+1:sum(varlens[1:kk])
end

"""
Translate an offset value (+1) to an index vector.
This faster if dims is a vector.  Use ind2sub if dims is tuple.
"""
function toindex(ii::Int64, dims::Vector{Int64})
    indexes = Vector{Int64}(undef, length(dims))
    offset = ii - 1
    for dd in 1:length(dims)
        indexes[dd] = offset % dims[dd] + 1
        offset = div(offset, dims[dd])
    end

    return indexes
end

"""
As `toindex`, but for many iis.  Ruturns N x D
"""
function toindexes(iis::Vector{Int64}, dims::Vector{Int64})
    indexes = Matrix{Int64}(undef, length(iis), length(dims))
    offsets = iis - 1
    for dd in 1:length(dims)
        indexes[:, dd] = offsets .% dims[dd] + 1
        offsets = div.(offsets, dims[dd])
    end

    return indexes
end

"""
Translate an index vector to an offset (+1).
This seems to be faster than sub2ind irrespective of the tupleness of dims.
"""
function fromindex(index::Vector{Int64}, dims::Vector{Int64})
    offset = index[end]
    for ii in length(dims)-1:-1:1
        offset = (offset - 1) * dims[ii] + index[ii]
    end

    return offset
end

"""
As `fromindex` but for many indexes (N x D)
"""
function fromindexes(indexes::Matrix{Int64}, dims::Vector{Int64})
    offsets = indexes[:, end]
    for ii in length(dims)-1:-1:1
        offsets = (offsets - 1) * dims[ii] + indexes[:, ii]
    end

    return offsets
end

"""
Insert an extra dimension somewhere, adjusting offset (+1) values accordingly
"""
function insertdim(iis::Vector{Int64}, dims::Vector{Int64}, insertat::Int64, dimvalue::Int64, dimsize::Int64)
    @assert insertat >= 1 && insertat <= length(dims) + 1
    @assert dimvalue >= 1 && dimvalue <= dimsize
    if insertat == 1
        (iis - 1) * dimsize + dimvalue
    elseif insertat == length(dims) + 1
        (dimvalue - 1) * prod(dims) + iis
    else
        proddimsleft = prod(dims[1:insertat-1])
        (div.(iis .- 1, proddimsleft) * dimsize .+ dimvalue .- 1) * proddimsleft .+ (iis .- 1) .% proddimsleft .+ 1
    end
end

"""
Insert as many dims as needed to expand to the alldims
"""
function expanddims(iis::Vector{Int64}, alldims::Vector{Int64}, newdims::AbstractVector{Bool}, newdimvalues::Vector{Int64})
    insertat = findfirst(newdims)
    iis2 = insertdim(iis, alldims[.!newdims], insertat, newdimvalues[1], alldims[insertat])
    if length(newdimvalues) == 1
        return iis2
    else
        newdims[insertat] = false
        return expanddims(iis2, alldims, newdims, newdimvalues[2:end])
    end
end

"Return the total span occupied by each variable or parameter."
function varlengths(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, vardictionary::Dict{Symbol, Symbol}=Dict{Symbol, Symbol}())
    Int64[prod(getdims(model, components[ii], get(vardictionary, names[ii], names[ii]))) for ii in 1:length(components)]
end

## Matrix methods

"Splits dupover into inner and outer dimensions"
function interpretdupover(dupover::Vector{Bool})
    ii = 1

    # Initial trues
    inners = Int64[]
    if dupover[1]
        while ii < length(dupover)
            push!(inners, ii)
            ii += 1
            if !dupover[ii]
                break
            end
        end
    end
    if dupover[ii] # fell off the end
        error("Cannot duplicate all dimensions.")
    end

    # Middle falses
    while ii < length(dupover)
        ii += 1
        if dupover[ii]
            break
        end
    end
    if !dupover[ii]
        return Int64[], inners
    end

    # Final trues
    outers = Int64[]
    while ii <= length(dupover)
        if !dupover[ii]
            error("Cannot support middle dimension duplication.")
        end

        push!(outers, ii)
        ii += 1
    end

    return outers, inners
end

"Construct a vector of values corresponding to entries in a matrix with the given dimensions, calling gen for each element."
function vectorsingle(dims::Vector{Int64}, gen::Function)
    f = Vector{Float64}(undef, prod(dims))
    for ii in 1:length(f)
        f[ii] = gen(toindex(ii, dims)...)
    end

    f
end

function vectorsingle(dims::Vector{Int64}, constant::Float64)
    return constant * ones(prod(dims))
end

function vectorsingle(dims::Vector{Int64}, gen::Function, dupover::Vector{Bool})
    outers, inners = interpretdupover(dupover)
    dupouter = prod(dims[outers])
    dupinner = prod(dims[inners])

    f = Vector{Float64}(undef, prod(dims[.!dupover]))
    for ii in 1:length(f)
        f[ii] = gen(toindex(ii, dims[.!dupover])...)
    end

    repeat(f, inner=[dupinner], outer=[dupouter])
end

## Example to be included below:
# # Examples
# ```jldoctest
# using OptiMimi

# dims = [3, 2]
# A = OptiMimi.matrixdiagonal(dims, (ii, jj) -> 1)
# sum(A)

# # output
# 6.0
# ```
"""
    matrixdiagonal(dims, gen)
"""
function matrixdiagonal(dims::Vector{Int64}, gen::Function)
    dimlen = prod(dims)
    A = spzeros(dimlen, dimlen)
    index = ones(Int64, length(dims))
    for ii in 1:dimlen
        A[ii, ii] = gen(index...)

        index[1] += 1
        for kk in 1:(length(dims) - 1)
            if index[kk] <= dims[kk]
                break
            end
            index[kk] = 1
            index[kk+1] += 1
        end
        #A[ii, ii] = gen(toindex(ii, dims)...) # TOO SLOW
    end

    A
end

function matrixdiagonal(dims::Vector{Int64}, constant::Float64)
    dimlen = prod(dims)
    sparse(1.0I, dimlen, dimlen) * constant
end

function matrixdiagonal(dims::Vector{Int64}, gen::Function, dupover::Vector{Bool})
    diag = vectorsingle(dims, gen, dupover)

    dimlen = prod(dims)
    A = spzeros(dimlen, dimlen)
    for ii in 1:dimlen
        if diag[ii] != 0
            A[ii, ii] = diag[ii]
        end
    end

    A
end

"""
Identifies all the common dimensions, calling the generator function
with an intersection matrix once for each value of the shared
dimensions (filling out their diagonal).

gen is called with each combination of unshared variables, and
returns a full matrix for all shared variables, in variable dimension order.
"""
function matrixdiagonalintersect(rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdimnames::Vector{Symbol}, coldimnames::Vector{Symbol}, gen::Function)
    common = intersect(rowdimnames, coldimnames)

    commonrow = Bool[dimname in common for dimname in rowdimnames]
    commoncol = Bool[dimname in common for dimname in coldimnames]

    remainrowdims = rowdims[.!commonrow]
    remaincoldims = coldims[.!commoncol]

    fullrowsub = repmat([0], length(rowdims))
    fullcolsub = repmat([0], length(coldims))

    alliis = []
    alljjs = []
    allvvs = []
    for ii in 1:prod(remainrowdims)
        for jj in 1:prod(remaincoldims)
            remainrowsub = toindex(ii, remainrowdims)
            remaincolsub = toindex(jj, remaincoldims)

            Apart = gen(remainrowsub..., remaincolsub...)
            kks = collect(1:prod(rowdims[commonrow]))

            if (sum(.!commonrow) == 0)
                append!(alliis, kks)
            else
                append!(alliis, expanddims(kks, rowdims, .!commonrow, remainrowsub))
            end
            if (sum(.!commoncol) == 0)
                append!(alljjs, kks)
            else
                append!(alljjs, expanddims(kks, coldims, .!commoncol, remaincolsub))
            end
            append!(allvvs, vec(Apart))
        end
    end

    sparse(alliis, alljjs, allvvs, prod(rowdims), prod(coldims))
end

"""
    matrixsingle(vardims, pardims, gen)

Call the generate function for every element.  `gen` is given the
dimensions for the rows/variable followed by the dimensions for the
columns/parameters, all in order.  It should return a Float64.
"""
function matrixsingle(vardims::Vector{Int64}, pardims::Vector{Int64}, gen::Function)
    vardimlen = prod(vardims)
    pardimlen = prod(pardims)
    A = spzeros(vardimlen, pardimlen)
    for ii in 1:vardimlen
        for jj in 1:pardimlen
            A[ii, jj] = gen(toindex(ii, vardims)..., toindex(jj, pardims)...)
        end
    end

    A
end

function matrixsingle(vardims::Vector{Int64}, pardims::Vector{Int64}, gen::Function, vardupover::Vector{Bool}, pardupover::Vector{Bool})
    outers, inners = interpretdupover(vardupover)
    vardupouter = prod(vardims[outers])
    vardupinner = prod(vardims[inners])

    outers, inners = interpretdupover(pardupover)
    pardupouter = prod(pardims[outers])
    pardupinner = prod(pardims[inners])

    vardimlen = prod(vardims[.!vardupover])
    pardimlen = prod(pardims[.!pardupover])
    A = zeros(vardimlen, pardimlen) # Not sparse
    for ii in 1:vardimlen
        for jj in 1:pardimlen
            A[ii, jj] = gen(toindex(ii, vardims[.!vardupover])..., toindex(jj, pardims[.!pardupover])...)
        end
    end

    sparse(repeat(A, inner=[vardupinner, pardupinner], outer=[vardupouter, pardupouter]))
end

"""
matrixchunks(rowdims, coldims, gen, rowchunks, colchunks)
Call the `gen` function with all combinations of chunked indices.
    rowchunks and colchunks is number of row and col vars at end to chunk.
"""
function matrixchunks(rowdims::Vector{Int64}, coldims::Vector{Int64}, gen::Function, rowchunk::Int64, colchunk::Int64)
    alliis = Int64[]
    alljjs = Int64[]
    allvvs = Float64[]
    #A = spzeros(prod(rowdims), prod(coldims))

    chunkrowdims = rowdims[end-rowchunk+1:end]
    chunkrowdimlen = prod(chunkrowdims)
    chunkcoldims = coldims[end-colchunk+1:end]
    chunkcoldimlen = prod(chunkcoldims)
    for ii in 1:chunkrowdimlen
        rowindex = toindex(ii, chunkrowdims)
        for jj in 1:chunkcoldimlen
            colindex = toindex(jj, chunkcoldims)

            topii = fromindex([ones(Int64, length(rowdims) - rowchunk); rowindex], rowdims)
            bottomii = fromindex([rowdims[1:length(rowdims) - rowchunk]; rowindex], rowdims)
            leftii = fromindex([ones(Int64, length(coldims) - colchunk); colindex], coldims)
            rightii = fromindex([coldims[1:length(coldims) - colchunk]; colindex], coldims)

            subA = gen(rowindex..., colindex...)
            @assert size(subA)[1] == bottomii - topii + 1 "Got matrix of size $(size(subA)) for rows $topii - $bottomii"
            @assert size(subA)[2] == rightii - leftii + 1 "Got matrix of size $(size(subA)) for columns $leftii - $rightii"

            iis, jjs, vvs = findnz(subA)
            append!(alliis, iis .+ (topii - 1))
            append!(alljjs, jjs .+ (leftii - 1))
            append!(allvvs, vvs)
            #A[topii:bottomii, leftii:rightii] = gen(rowindex..., colindex...)
        end
    end

    sparse(alliis, alljjs, allvvs, prod(rowdims), prod(coldims))
    #A
end

"""
    matrixintersect(rowdims, coldims, rowdimnames, coldimnames, gen)

Call the `gen` function with all combinations of the shared indices.
Shared dimensions must come in the same order and at the end
of the dimensions lists.

# Arguments
* `rowdims::Vector{Int64}` and `coldims::Vector{Int64}`: the number of elements in each dimension of the row and column variables.
* `rowdimnames::Vector{Symbol}` and `coldimnames::Vector{Symbol}`: The names of the dimensions; note that the last 1 or more dimensions should be the same in both of these lists.
* `gen::Function`: a function of arguments (`A`, `indices`...), where `A` is a matrix of the un-shared indices and `indices` are multiple indexes describing all shared indices.

# Examples
```jldoctest
using OptiMimi

rowdims = [3, 2]
coldims = [4, 2]
rowdimnames = [:region, :time]
coldimnames = [:person, :time]
function gen(A, tt)
    @assert tt == 1 || tt == 2
    @assert size(A) == (3, 4)
    A[2, 3] = tt
end
A = OptiMimi.matrixintersect(rowdims, coldims, rowdimnames, coldimnames, gen)
sum(A)

# output
3.0
```
"""
function matrixintersect(rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdimnames::Vector{Symbol}, coldimnames::Vector{Symbol}, gen::Function)
    A = spzeros(prod(rowdims), prod(coldims))

    # Determine shared dimensions: counting from 0 = end
    numshared = 0
    while numshared < min(length(rowdims), length(coldims))
        if rowdims[end - numshared] != coldims[end - numshared]
            break
        end
        if rowdimnames[end - numshared] != coldimnames[end - numshared]
            break
        end
        numshared += 1
    end

    sharedims = rowdims[end - numshared + 1:end]
    sharedimslen = prod(sharedims)
    for ii in 1:sharedimslen
        index = toindex(ii, sharedims)
        topii = fromindex([ones(Int64, length(rowdims) - numshared); index], rowdims)
        bottomii = fromindex([rowdims[1:length(rowdims) - numshared]; index], rowdims)
        leftii = fromindex([ones(Int64, length(coldims) - numshared); index], coldims)
        rightii = fromindex([coldims[1:length(coldims) - numshared]; index], coldims)
        subA = view(A, topii:bottomii, leftii:rightii)
        gen(subA, index...)
    end

    A
end

function matrixintersect(rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdimnames::Vector{Symbol}, coldimnames::Vector{Symbol}, gen::Function, rowdupover::Vector{Bool}, coldupover::Vector{Bool})
    # Figure out how many of dupovers are shared
    rowdupovershared = zeros(Bool, length(rowdupover))
    coldupovershared = zeros(Bool, length(coldupover))
    for ii in 0:(min(length(rowdims), length(coldims))-1)
        if rowdupover[end-ii] && coldupover[end-ii] && rowdimnames[end-ii] == coldimnames[end-ii]
            rowdupovershared[end-ii] = true
            coldupovershared[end-ii] = true
        else
            break
        end
    end

    # Generate, without any dups
    A = matrixintersect(rowdims[.!rowdupover], coldims[.!coldupover], rowdimnames[.!rowdupover], coldimnames[.!coldupover], gen)

    A2 = matrixduplicate(A, rowdims[.!rowdupovershared], coldims[.!coldupovershared], rowdupover[.!rowdupovershared], coldupover[.!coldupovershared])

    alliis, alljjs, allvvs = findnz(A2)

    # Create the shared-dupped portion

    shareddupouter = prod(rowdims[rowdupovershared])
    rownotshared = prod(rowdims[.!rowdupovershared])
    colnotshared = prod(coldims[.!coldupovershared])

    allvvs2 = repeat(allvvs, outer=[shareddupouter])
    alliis2 = zeros(Int64, length(alliis) * shareddupouter)
    alljjs2 = zeros(Int64, length(alljjs) * shareddupouter)

    for kk in 1:shareddupouter
        alliis2[(kk - 1) * length(alliis) .+ (1:length(alliis))] = (kk - 1) * rownotshared .+ alliis
        alljjs2[(kk - 1) * length(alljjs) .+ (1:length(alljjs))] = (kk - 1) * colnotshared .+ alljjs
    end

    sparse(alliis2, alljjs2, allvvs2, prod(rowdims), prod(coldims))
end

"""
Create an empty matrix for the given variables
"""
function matrixempty(vardims::Vector{Int64}, pardims::Vector{Int64})
    vardimlen = prod(vardims)
    pardimlen = prod(pardims)
    spzeros(vardimlen, pardimlen)
end

"""
General duplication of a known matrix over additional dimensions.
Falls back on matrixduplicate_general if not all inner-most and outer-most dupovers
"""
function matrixduplicate(A::SparseMatrixCSC{Float64, Int64}, rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdupover::Vector{Bool}, coldupover::Vector{Bool})
    if sum(abs.(diff([true; rowdupover]))) <= 2 && sum(abs.(diff([true; coldupover]))) <= 2
        matrixduplicate_extremes(A, rowdims, coldims, rowdupover, coldupover)
    else
        matrixduplicate_general(A, rowdims, coldims, rowdupover, coldupover)
    end
end

"""
Duplication of a known matrix over inner-most and outer-most dimensions
"""
function matrixduplicate_extremes(A::SparseMatrixCSC{Float64, Int64}, rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdupover::Vector{Bool}, coldupover::Vector{Bool})
    outers, inners = interpretdupover(rowdupover)
    rowdupouter = prod(rowdims[outers])
    rowdupinner = prod(rowdims[inners])

    outers, inners = interpretdupover(coldupover)
    coldupouter = prod(coldims[outers])
    coldupinner = prod(coldims[inners])

    # Create fully-dupped portion
    iis, jjs, vvs = findnz(A)

    kkdupnum = rowdupouter*rowdupinner*coldupouter*coldupinner

    allvvs = repeat(vvs, inner=[kkdupnum])
    # Need to fill these next two in
    alliis = zeros(Int64, length(iis) * kkdupnum)
    alljjs = zeros(Int64, length(jjs) * kkdupnum)

    nrows = size(A)[1]
    ncols = size(A)[2]

    for kk in 1:length(iis)
        ## order is vec([(slow, fast) for fast in 1:N, slow in 1:M])
        rowduped = vec([(rowoo - 1) * nrows * rowdupinner + (iis[kk] - 1) * rowdupinner + rowii for rowii in 1:rowdupinner, rowoo in 1:rowdupouter])
        colduped = vec([(coloo - 1) * ncols * coldupinner + (jjs[kk] - 1) * coldupinner + colii for colii in 1:coldupinner, coloo in 1:coldupouter])
        alliis[(kk - 1) * kkdupnum .+ (1:kkdupnum)] = repeat(rowduped, outer=[coldupouter*coldupinner])
        alljjs[(kk - 1) * kkdupnum .+ (1:kkdupnum)] = repeat(colduped, inner=[rowdupouter*rowdupinner])
    end

    sparse(alliis, alljjs, allvvs, prod(rowdims), prod(coldims))
end

"""
General duplication of a known matrix over additional dimensions.
This is not constrained to have dupover only on inner-most and outer-most dimensions.
"""
function matrixduplicate_general(origA::SparseMatrixCSC{Float64, Int64}, rowdims::Vector{Int64}, coldims::Vector{Int64}, rowdupover::Vector{Bool}, coldupover::Vector{Bool})
    iis, jjs, vvs = findnz(origA)

    alliis = Int64[]
    alljjs = Int64[]
    allvvs = Float64[]

    if !any(rowdupover)
        newdimsizes = coldims[coldupover]
        for kk in 1:prod(coldims[coldupover])
            fulljjs = expanddims(jjs, coldims, coldupover, toindex(kk, newdimsizes))
            append!(alliis, iis)
            append!(alljjs, fulljjs)
            append!(allvvs, vvs)
        end
    elseif !any(coldupover)
        newdimsizes = rowdims[rowdupover]
        for kk in 1:prod(rowdims[rowdupover])
            fulliis = expanddims(iis, rowdims, rowdupover, toindex(kk, newdimsizes))
            append!(alliis, fulliis)
            append!(alljjs, jjs)
            append!(allvvs, vvs)
        end
    else
        newrowsizes = rowdims[rowdupover]
        newcolsizes = coldims[coldupover]
        for rowkk in 1:prod(rowdims[rowdupover])
            for colkk in 1:prod(coldims[coldupover])
                fulliis = expanddims(iis, rowdims, rowdupover, toindex(rowkk, newrowsizes))
                fulljjs = expanddims(jjs, coldims, coldupover, toindex(colkk, newcolsizes))
                append!(alliis, fulliis)
                append!(alljjs, fulljjs)
                append!(allvvs, vvs)
            end
        end
    end

    sparse(alliis, alljjs, allvvs, prod(rowdims), prod(coldims))
end

"""
General duplication of a known matrix in a block diagonal over additional dimensions.
All off-diagonal entries are block 0s.
"""
function matrixduplicate_diagonal(origA::SparseMatrixCSC{Float64, Int64}, rowcoldims::Vector{Int64}, rowcoldupover::Vector{Bool})
    iis, jjs, vvs = findnz(origA)

    alliis = Int64[]
    alljjs = Int64[]
    allvvs = Float64[]

    newrowcolsizes = rowcoldims[rowcoldupover]
    for rowcolkk in 1:prod(rowcoldims[rowcoldupover])
        fulliis = expanddims(iis, rowcoldims, rowcoldupover, toindex(rowcolkk, newrowcolsizes))
        fulljjs = expanddims(jjs, rowcoldims, rowcoldupover, toindex(rowcolkk, newrowcolsizes))
        append!(alliis, fulliis)
        append!(alljjs, fulljjs)
        append!(allvvs, vvs)
    end

    sparse(alliis, alljjs, allvvs, prod(rowcoldims), prod(rowcoldims))
end


text:
@doc raw"""
    quarnetGoFtest!(net::HybridNetwork, df::DataFrame, optbl::Bool; quartetstat=:LRT, correction=:simulation, seed=1234, nsim=1000, verbose=false, keepfiles=false)
    quarnetGoFtest!(net::HybridNetwork, dcf::DataCF,   optbl::Bool; kwargs...)

Goodness-of-fit test for the adequacy of the multispecies network coalescent,
to see if a given population or species network explains the
quartet concordance factor data adequately.
The network needs to be of level 1 at most (trees fullfil this condition),
and have branch lengths in coalescent units.
The test assumes a multinomial distribution for the observed quartet
concordance factors (CF), such that information on the number of genes for each
four-taxon set (`ngenes` field) must be present.

For each four-taxon set, an outlier p-value is obtained by comparing a
test statistic (-2log likelihood ratio by default) to a chi-square distribution
with 2 degrees of freedom (see below for other options).

The four-taxon sets are then categorized as outliers or not, according to their
outlier p-values (outlier if p<0.05).
Finally, a one-sided goodness-of-fit test is performed on the frequency of
outlier 4-taxon sets to see if there are more outliers than expected.
The z-value for this test corresponds to the null hypothesis that
5% outlier p-values are < 0.05 (versus more than 5%):
``z = \frac{proportion_of_outliers - 0.05}{\sqrt(0.05 * 0.95/number_4_taxon_sets)}``
This z-value corresponds to a test that assumes independent outlier p-values
across 4-taxon sets: it makes no correction for dependence.

To correct for dependence with `correction=:simulation`, the distribution of
z-values is obtained by simulating gene trees under the coalescent along the
network (after branch length optimization if `optbl=true`) using `hybrid-Lambda`.
The z-score is calculated on each simulated data set.
Under independence, these z-scores have mean 0 and variance 1.
Under dependence, these z-scores still have mean 0, but an inflated variance.
This variance σ² is estimated from the simulations, and the corrected p-value
is obtained by comparing the original z value to N(0,σ).
When `correction=:none`, σ is taken to be 1 (independence): *not* recommended!

- The first version takes a `DataFrame` where each row corresponds to a given
  four-taxon set. The data frame is modified by having an additional another column
  containing the p-values corresponding to each four-taxon set.
- The second version takes a `DataCF` object and modifies it by updating
  the expected concordance factors stored in that object.

# arguments

- `optbl`: when `false`, `net.loglik` is updated but branch lengths are
  taken as is. When `true`, a copy of `net` is used to conduce the test,
  with optimized branch lengths (in coalescent units) and updated loglik.
  This network is returned.

# keyword arguments

- `quartetstat`: the test statistic used to obtain an outlier p-value for
  each four-taxon set, which is then compared to a chi-squared distribution
  with 2 degrees of freedom to get a p-value.
  The default is `:LRT` for the likelihood ratio:
  ``2n_{genes} \sum_{j=1}^3 {\hat p}_j (\log{\hat p}_j - \log p_j)`` where ``p_j``
  is the quartet CF expected from the network, and ``{\hat p}_j`` is the
  quartet CF observed in the data.
  Alternatives are `:Qlog` for the Qlog statistics (Lorenzen, 1995):
  ``2n_{genes} \sum_{j=1}^3 \frac{({\hat p}_j - p_j)^2}{p_j (\log{\hat p}_j - \log p_j)}``
  and `:pearson` for Pearon's chi-squared statistic, which behaves poorly when
  the expected count is low (e.g. less than 5):
  ``n_{genes} \sum_{j=1}^3 \frac{({\hat p}_j - p_j)^2 }{p_j}``
- `correction=:simulation` to correct for dependence across 4-taxon.
  Use `:none` to turn off simulations (and the correction for dependence).
- `seed=1234`: master seed to control the seeds for gene tree simulations.
- `nsim=1000`: number of simulated data sets. Each data set is simulated to have the
  median number of genes that each 4-taxon sets has data for.
- `verbose=false`: output from `hybrid-Lambda` and other output is suppressed.
  Turn to `true` to diagnose potential issues.
- `keepfiles=false`: files generated by `hybrid-Lambda` are stored in a temporary
  folder, whose name starts with `jl_` and placed in the current directory.
  Each of the 1000 simulations creates a file, that is deleted when is no longer
  necessary. Turn `keepfiles=true` to keep these gene tree files.

# output

1. p-value of the overall goodness-of-fit test (corrected for dependence if requested)
2. test statistic (z value), uncorrected
3. estimated σ for the test statistic, used for the correction (1.0 if no correction)
4. a vector of outlier p-values, one for each four-taxon set
5. network (first and second versions):
   `net` with loglik field updated if `optbl` is false;
    copy of `net` with optimized branch lengths and loglik if `optbl` is true
6. in case `correction = :simulation`, vector of simulated z values
   (`nothing` if `correction = :none`). These z-values could be used to
   calculate an empirical p-value (instead of the p-value in #1), as the
   proportion of simulated z-values that are ⩾ the observed z-value (in #2).

See also: [`ticr!`](@ref).

# references

- Lorenzen (1995).
  A new family of goodness-of-fit statistics for discrete multivariate data.
  Statistics & Probability Letters, 25(4):301-307.
  doi: [10.1016/0167-7152(94)00234-8](https://doi.org/10.1016/0167-7152(94)00234-8)
"""
function quarnetGoFtest!(net::HybridNetwork,  df::DataFrame, optbl::Bool; kwargs...)
    d = readTableCF(df);
    res = quarnetGoFtest!(net, d, optbl; kwargs...);
    df[!,:p_value] .= res[4] # order in "res": overallpval, uncorrected z-value, sigma, pval, ...
    return res
end

function quarnetGoFtest!(net::HybridNetwork, dcf::DataCF, optbl::Bool;
                         quartetstat::Symbol=:LRT, correction::Symbol=:simulation,
                         seed=1234::Int, nsim=1000::Int, verbose=false::Bool, keepfiles=false::Bool)
    correction in [:simulation, :none] || error("correction ($correction) must be one of :none or :simulation")
    quartetstat in [:LRT, :Qlog, :pearson] || error("$quartetstat is not a valid quartetstat option")
    if optbl
        net = topologyMaxQPseudolik!(net,dcf);
    else
        topologyQPseudolik!(net,dcf);
    end
    outlierp_fun! = ( quartetstat ==  :LRT ? multinom_lrt! :
                     (quartetstat == :Qlog ? multinom_qlog! : multinom_pearson!))
    gof_zval, outlierpvals = quarnetGoFtest(dcf.quartet, outlierp_fun!)
    sig = 1.0 # 1 if independent: correction for dependence among 4-taxon set outlier pvalues
    sim_zvals = nothing
    if correction == :simulation
        # expCF in dcf correspond to net: assumption made below
        sig, sim_zvals = quarnetGoFtest_simulation(net, dcf, outlierp_fun!, seed, nsim, verbose, keepfiles)
    end
    overallpval = normccdf(gof_zval/sig) # one-sided: P(Z > z)
    return (overallpval, gof_zval, sig, outlierpvals, net, sim_zvals)
end

"""
    quarnetGoFtest(quartet::Vector{Quartet}, outlierp_fun!::Function)
    quarnetGoFtest(outlier_pvalues::AbstractVector)

Calculate an outlier p-value for each `quartet` according to function `outlierp_fun!`
(or take outlier-values as input: second version) and calculate the z-value
to test the null hypothesis that 5% of the p-values are < 0.05, versus
the one-sided alternative of more outliers than expected.

See [`quarnetGoFtest!`](@ref) for more details.

Output:
- z-value
- outlier p-values (first version only)
"""
function quarnetGoFtest(quartet::Vector{Quartet}, outlierp_fun!::Function)
   for q in quartet
        q.ngenes > 0 || error("quartet $q does not have info on number of genes")
    end
    pval = fill(-1.0, length(quartet))
    outlierp_fun!(pval, quartet) # calculate outlier p-values
    gof_zval = quarnetGoFtest(pval)
    return (gof_zval, pval)
end

function quarnetGoFtest(pval::AbstractVector)
    # uncorrected z-value: from one-sided test after binning p-values into [0,.05) and [.05,1.0]
    nsmall = count(p -> p < 0.05, pval)
    ntot = count(p -> isfinite(p), pval)
    ntot == length(pval) || @warn "$(length(pval) - ntot) outlier p-value(s) is/are Inf of NaN..."
    gof_zval = (nsmall/ntot - 0.05)/sqrt(0.0475/ntot) # 0.0475 = 0.05 * (1-0.05)
    return gof_zval
end

"""
    quarnetGoFtest_simulation(net::HybridNetwork, dcf::DataCF, outlierp_fun!::Function,
                              seed::Int, nsim::Int, verbose::Bool, keepfiles::Bool)

Run external program `hybrid-Lambda` (downloaded at installation, with
path stored in `QuartetNetworkGoodnessFit.hybridlambda`) to simulate
gene trees under the multispecies coalescent model along network `net`.
The quartet concordance factors (CFs) from these simulated gene trees are used
as input to `outlierp_fun!` to categorize each 4-taxon set as an outlier
(p-value < 0.05) or not.
For each simulated data set, a goodness-of-fit z-value is calculated by
comparing the proportion of outlier 4-taxon sets to 0.05. The standard deviation
of these z-values (assuming a mean of 0), and the z-values themselves are returned.

Used by [`quarnetGoFtest!`](@ref).

**Warning**: The quartet CFs expected from `net` are assumed to be stored in
`dcf.quartet[i].qnet.expCF`. This is *not* checked.
"""
function quarnetGoFtest_simulation(net::HybridNetwork, dcf::DataCF, outlierp_fun!::Function,
        seed::Int, nsim::Int, verbose::Bool, keepfiles::Bool)
    ngenes = ceil(Int, median(q.ngenes for q in dcf.quartet))
    # how to handle case when different #genes across quartets? need taxon set for each gene, not just ngenes.
    seed!(seed) # master seed
    hlseed = rand(1:10_000_000_000, nsim) # 1 seed for each simulation
    netHL = hybridlambdaformat(net)
    netHLquoted = "'$netHL'"
    nq = length(dcf.quartet)
    pval = fill(-1.0, nq) # to be re-used across simulations, but NOT shared between processes
    sim_zval = SharedArray{Float64}(nsim) # to be shared between processes
    # expected CFs: in dcf.quartet[i].qnet.expCF[j] for 4-taxon set i and resolution j
    # BUT countquartetsintrees might list 4-taxon sets in a different order, and might list
    # the 4 taxa within a set in a different order -> need to re-order resolutions within a set
    # extra problem: hybrid-Lambda changes tip names: adds suffix "_1" to all of them (if 1 individual/tip)
    expCF, taxa = expectedCF_ordered(dcf, net, "_1") # 'taxa' gives the map i -> taxon[i]
    hyblamdir = mktempdir(pwd()) # temporary directory to store hybrid-Lambda output files

    @sync @distributed for irep in 1:nsim
        verbose && @info "starting replicate $irep"
        gt = joinpath(hyblamdir, "genetrees_rep$irep") # root name given to hybrid-Lambda
        gtcu = gt * "_coal_unit"                 # name of output file created by hybrid-Lambda
        hlcommand = `$hybridlambda -spcu $netHLquoted -num $ngenes -seed $(hlseed[irep]) -o $gt`
        hlout = ( verbose ? stdout : devnull )
        run(pipeline(hlcommand; stderr = hlout));
        # problem: hybrid-Lambda replaces individual names like "s5" by "s5_1", since 1 individual / tip
        ## solution below: platform specific. Instead, deal with suffix in expectedCF_ordered
        ## run(`sed -i "" 's/_1//g' $gtcu`); # changes individual names like "s5_1" back into "s5"
        treelist = readMultiTopology(gtcu);
        length(treelist) == ngenes || @warn "unexpected number of gene trees, file $gtcu" # sanity check
        obsCF, t = countquartetsintrees(treelist; showprogressbar=verbose)
        irep == 1 && (taxa == t || error("different order of taxa used by countquartetsintrees"))
        outlierp_fun!(pval, (q.data for q in obsCF), expCF) # calculate outlier p-values
        sim_zval[irep] = quarnetGoFtest(pval)
        keepfiles || rm(gtcu) # delete file for replicate irep
    end
    keepfiles || rm(hyblamdir) # delete (empty) directory that had output of hybrid-Lambda
    sigma = sqrt(sum(sim_zval.^2)/nsim) # estimated sigma, assuming mean=0
    return sigma, sim_zval
end

"""
    expectedCF_ordered(dcf::DataCF, net::HybridNetwork, suffix=""::AbstractString)

Expected quartet concordance factors in `dcf`, but ordered as they would be if
output by `PhyloNetworks.countquartetsintrees`.
Output:
- 2-dimentional `SharedArray` (number of 4-taxon sets x 3).
  `dcf.quartet[i].qnet.expCF[j]` for 4-taxon set `i` and resolution `j`
  is stored in row `qi` and column `k` if `qi` is the rank of 4-taxon set `i`
  (see `PhyloNetworks.quartetrank`). This rank depends on how taxa are ordered.
- vector of taxon names, whose order matters. These are tip labels in `net` with
  suffix `suffix` added, then ordered alphabetically, or
  numerically if taxon names can be parsed as integers.
"""
function expectedCF_ordered(dcf::DataCF, net::HybridNetwork, suffix=""::AbstractString)
    nq = length(dcf.quartet)
    expCF = SharedArray{Float64, 2}(nq,3)
    # careful ordering: ["t","t_0"] ordered differently after we add suffix "_1"
    taxa = PhyloNetworks.sort_stringasinteger!(tipLabels(net) .* suffix) # same as done in countquartetsintrees
    nsuff = length(suffix)
    taxonnumber = Dict(chop(taxa[i];tail=nsuff) => i for i in eachindex(taxa))
    ntax = length(taxa)
    nCk = PhyloNetworks.nchoose1234(ntax) # matrix used to ranks 4-taxon sets
    nq == nCk[ntax+1,4] || error("dcf is assumed to contain ALL $(nCk[ntax+1,4]) four-taxon sets, but contains $nq only.")
    ptype = Vector{Int8}(undef,4) # permutation vector to sort the 4 taxa
    resperm = Vector{Int8}(undef,3) # permutation to sort the 3 resolutions accordingly
    for q in dcf.quartet
        tn = map(i->taxonnumber[i], q.taxon)
        sortperm!(ptype, tn)
        qi = PhyloNetworks.quartetrank(tn[ptype]..., nCk)
        # next: find the corresponding permutation of the 3 CFs.
        # this permutation is invariant to permuting ptype like: [a b c d] -> [c d a b] or [b a d c]
        # modify ptype to have 1 first, then only 6 permutations, corresponding to 6 permutations of 3 resolutions
        ptype .-= 0x01
        if ptype[1]==0x00
            resperm[1] = ptype[2]
            resperm[2] = ptype[3]
            resperm[3] = ptype[4]
        elseif ptype[2]==0x00 # consider ptype[[2,1,4,3]] to put 0(=1-1) first
            resperm[1] = ptype[1]
            resperm[2] = ptype[4]
            resperm[3] = ptype[3]
        elseif ptype[3]==0x00 # consider ptype[[3,4,1,2]] to put 0(=1-1) first
            resperm[1] = ptype[4]
            resperm[2] = ptype[1]
            resperm[3] = ptype[2]
        else # consider ptype[[4,3,2,1]] to put 0(=1-1) first instead of last
            resperm[1] = ptype[3]
            resperm[2] = ptype[2]
            resperm[3] = ptype[1]
        end
        expCF[qi,:] = q.qnet.expCF[resperm]
    end
    return expCF, taxa
end

"""
    multinom_pearson!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    multinom_pearson!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})

Calculate outlier p-values (one per four-taxon set)
using Pearson's chi-squared statistic under a multinomial distribution
for the observed concordance factors.
"""
function multinom_pearson!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    nq = length(quartet)
    for i in 1:nq
        qt = quartet[i]
        phat = qt.obsCF
        p = qt.qnet.expCF
        ngenes = qt.ngenes
        ipstat = ngenes * sum((phat .- p).^2 ./ p)
        pval[i] = chisqccdf(2, ipstat)
    end
    return nothing
end
function multinom_pearson!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})
    for (i,o) in enumerate(obsCF) # o = [observedCF, ngenes]
        p = expCF[i,:]
        ipstat = o[4] * sum((o[1:3] .- p).^2 ./ p)
        pval[i] = chisqccdf(2, ipstat)
    end
end

"""
    multinom_qlog!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    multinom_qlog!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})

Calculate outlier p-values (one per four-taxon set)
using the Qlog statistic (Lorenzen, 1995),
under a multinomial distribution for the observed concordance factors.
"""
function multinom_qlog!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    nq = length(quartet)
    for i in 1:nq
        qt = quartet[i]
        phat = qt.obsCF
        p = qt.qnet.expCF
        ngenes = qt.ngenes
        mysum = 0.0
        for j in 1:3
            if phat[j] > eps(Float64) && !isapprox(phat[j],p[j]; atol=1e-20)
                mysum += (phat[j]-p[j])^2/(p[j]*(log(phat[j])-log(p[j])))
            end
        end
        iqstat = 2 * ngenes * mysum
        pval[i] = chisqccdf(2,iqstat)
    end
    return nothing
end
function multinom_qlog!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})
    for (i,o) in enumerate(obsCF) # o = [observedCF, ngenes]
        p = expCF[i,:]
        mysum = 0.0
        for j in 1:3
            if o[j] > eps(Float64) && !isapprox(o[j],p[j]; atol=1e-20)
                mysum += (o[j]-p[j])^2/(p[j]*(log(o[j])-log(p[j])))
            end
        end
        iqstat = 2 * o[4] * mysum
        pval[i] = chisqccdf(2,iqstat)
    end
end

"""
    multinom_lrt!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    multinom_lrt!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})

Calculate outlier p-values (one per four-taxon set)
using the likelihood ratio test under a multinomial distribution
for the observed concordance factors.
"""
function multinom_lrt!(pval::AbstractVector{Float64}, quartet::Vector{Quartet})
    nq = length(quartet)
    for i in 1:nq
        qt = quartet[i]
        phat = qt.obsCF
        p = qt.qnet.expCF
        ngenes = qt.ngenes
        mysum = 0.0
        for j in 1:3
            if phat[j] > eps(Float64)
                mysum += phat[j]*log(phat[j]/p[j])
            end
        end
        igstat = 2 * ngenes * mysum
        pval[i] = chisqccdf(2,igstat)
    end
    return nothing
end
function multinom_lrt!(pval::AbstractVector{Float64}, obsCF, expCF::AbstractMatrix{Float64})
    for (i,o) in enumerate(obsCF) # o = [observedCF, ngenes]
        p = expCF[i,:]
        mysum = 0.0
        for j in 1:3
            if o[j] > eps(Float64)
                mysum += o[j]*log(o[j]/p[j])
            end
        end
        igstat = 2 * o[4] * mysum
        pval[i] = chisqccdf(2,igstat)
    end
    return nothing
end

ext: jl


text:
import ThreeJS

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    ThreeJS.outerdiv() <<
    (ThreeJS.initscene() <<
        [
            ThreeJS.mesh(0.0, 0.0, 0.0) <<
            [
                #Any accessible URL can be provided to the texture attribute to render it as texture
                ThreeJS.box(2.0, 2.0, 2.0),
                ThreeJS.material(Dict(:kind=>"texture",:texture=>"assets/crate.png"))
            ],
            ThreeJS.mesh(3.0, 3.0, 3.0) <<
            [
                ThreeJS.plane(2.0, 2.0),
                ThreeJS.material(Dict(:kind=>"texture",:texture=>"assets/doge.png"))
            ],
            ThreeJS.pointlight(10.0, 10.0, 10.0),
            ThreeJS.camera(0.0, 0.0, 20.0)
        ]
    )
end


text:
"""
## loc_to_glob!

This subroutine transforms the local end reactions and
moments into the global system (3-d).

### Function
```julia
loc_to_glob!(loc, glob, gamma, coord)
```

### Arguments
```julia
* loc::Vector{Float64}       : Local force and momemts (Updated)
* glob::Vector{Float64}      : Globale forces and moments
* gamma::Float64             : Element orientation angle (3D)
* coord::Matrix{Float64}     : Nodal coordinates
```
"""
function loc_to_glob!(loc::Vector{Float64}, glob::Vector{Float64},
  gamma::Float64, coord::LinearAlgebra.Adjoint{Float64,Array{Float64,2}})
  #
  # This subroutine transforms the local end reactions and
  # moments into the global system (3-d).
  #
  # Updates glob[:]
  #
  t = zeros(12, 12)
  r0 = zeros(3, 3)

  x1 = coord[1,1]
  y1 = coord[1,2]
  z1 = coord[1,3]
  x2 = coord[2,1]
  y2 = coord[2,2]
  z2 = coord[2,3]
  xl = x2 - x1
  yl = y2 - y1
  zl = z2 - z1
  ell = sqrt(xl*xl+yl*yl+zl*zl)
  gamrad = gamma*pi/180.0
  cg = cos(gamrad)
  sg = sin(gamrad)
  den = ell*sqrt(xl*xl+zl*zl)
  if den !== 0.0
   r0[1,1] = xl/ell
   r0[2,1] = yl/ell
   r0[3,1] = zl/ell
   r0[1,2] = (-xl*yl*cg-ell*zl*sg)/den
   r0[2,2] = den*cg/(ell*ell)
   r0[3,2] = (-yl*zl*cg+ell*xl*sg)/den
   r0[1,3] = (xl*yl*sg-ell*zl*cg)/den
   r0[2,3] = -den*sg/(ell*ell)
   r0[3,3] = (yl*zl*sg+ell*xl*cg)/den
  else
   r0[1,1] = 0.0
   r0[3,1] = 0.0
   r0[2,2] = 0.0
   r0[2,3] = 0.0
   r0[2,1] = 1.0
   r0[1,2] = -cg
   r0[3,3] = cg
   r0[3,2] = sg
   r0[1,3] = sg
  end
  for i in 1:3
   for j in 1:3
     x = r0[i,j]
     for k in 0:3:9
       t[i+k, j+k] = x
     end
   end
  end
  for i in 1:12
   sum = 0.0
   for j in 1:12
     sum = sum+t[i,j]*loc[j]
   end
   glob[i] = sum
  end
  glob
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Primitive{Dim,T}

We say that a geometry is a primitive when it can be expressed as a single
entity with no parts (a.k.a. atomic). For example, a sphere is a primitive
described in terms of a mathematical expression involving a metric and a radius.
See https://en.wikipedia.org/wiki/Geometric_primitive.
"""
abstract type Primitive{Dim,T} <: Geometry{Dim,T} end

include("primitives/box.jl")
include("primitives/ball.jl")
include("primitives/sphere.jl")
include("primitives/cylinder.jl")


text:
# convert used as a generalized reinterpret

@inline Base.convert(::Type{Unsigned}, ::Type{Float16}) = UInt16
@inline Base.convert(::Type{Unsigned}, ::Type{Float32}) = UInt32
@inline Base.convert(::Type{Unsigned}, ::Type{Float64}) = UInt64

@inline Base.convert(::Type{SysFloat}, ::Type{UInt16}) = Float16
@inline Base.convert(::Type{SysFloat}, ::Type{UInt32}) = Float32
@inline Base.convert(::Type{SysFloat}, ::Type{UInt64}) = Float64

@inline Base.convert(::Type{Unsigned}, x::Float16) = reinterpret(Unsigned, x)
@inline Base.convert(::Type{Unsigned}, x::Float32) = reinterpret(Unsigned, x)
@inline Base.convert(::Type{Unsigned}, x::Float64) = reinterpret(Unsigned, x)

@inline Base.convert(::Type{SysFloat}, x::UInt16) = reinterpret(Float16, x)
@inline Base.convert(::Type{SysFloat}, x::UInt32) = reinterpret(Float32, x)
@inline Base.convert(::Type{SysFloat}, x::UInt64) = reinterpret(Float64, x)


text:
#!/usr/bin/env julia

insert(LOAD_PATH, 1, "../src/")
load("compose.jl")

const phi = (1 + sqrt(5)) / 2

function golden_rect(n::Int)
    if n == 0; return empty_canvas; end
    c = canvas() << (rectangle() << fill(LCHab(90, 80, 70 - 15n)))
    c | (canvas(0, -1/phi, 1h, 1/phi, Rotation(pi/2, 0, 1)) | golden_rect(n - 1))
end

@upon SVG("golden_rect.svg", phi * 3inch, 3inch) begin
    draw(golden_rect(10) << (fill(nothing) | stroke("white") | linewidth(0.2mm)))
end


text:
@testset "gaussian" begin
    seed!(42)
    @testset "Likelihood" begin
        σ² = 1e-3
        l = GaussianLikelihood(σ²)
        f = 0.5
        y = 0.2
        @test AGP.noise(l) == σ²
        @test l(y, f) == pdf(Normal(f, sqrt(σ²)), y)
        @test loglikelihood(l, y, f) == logpdf(Normal(f, sqrt(σ²)), y)
        @test repr(l) == "Gaussian likelihood (σ² = $(AGP.noise(l)))"
        @test AGP.AugmentedKL(l, []) == 0.0
        @test AGP.num_latent(l) == 1
    end
    N, d = 20, 2
    k = SqExponentialKernel() ∘ ScaleTransform(10.0)
    σ = 0.1
    X, f = generate_f(N, d, k)
    y = f + σ * randn(N)
    floattypes = [Float64]
    @testset "GP" begin
        for floattype in floattypes
            model = GP(X, y, k; opt_noise = true, verbose = 0)
            @test eltype(model) == floattype
            @test AGP.likelihood(model) isa GaussianLikelihood
            @test AGP.inference(model) isa Analytic
            @test AGP.getf(model) isa AGP.LatentGP
            @test AGP.output(model) isa AbstractVector
            @test AGP.input(model) isa AbstractVector
            @test AGP.nLatent(model) == 1
            L = AGP.objective(model)
            @test_nowarn train!(model, 10)
            @test L < AGP.objective(model)
            @test testconv(model, "Regression", X, f, y)
            @test all(proba_y(model, X)[2] .> 0)
        end
    end
    @testset "VGP" begin
        @test_throws ErrorException VGP(
            X,
            y,
            k,
            GaussianLikelihood(),
            AnalyticVI(),
        )
        @test_throws ErrorException VGP(
            X,
            y,
            k,
            GaussianLikelihood(),
            QuadratureVI(),
        )
        @test_throws ErrorException VGP(
            X,
            y,
            k,
            GaussianLikelihood(),
            MCIntegrationVI(),
        )
    end
    @testset "SVGP" begin
        for floattype in floattypes
            @testset "AnalyticVI" begin
                model = SVGP(
                    X,
                    y,
                    k,
                    GaussianLikelihood(),
                    AnalyticVI(),
                    10,
                    optimiser = false,
                    verbose = 0,
                )
                @test eltype(model) == floattype
                @test AGP.likelihood(model) isa GaussianLikelihood
                @test AGP.inference(model) isa AnalyticVI
                @test AGP.getf(model) isa NTuple{1, AGP.SparseVarLatent}
                @test AGP.output(model) isa AbstractVector
                @test AGP.input(model) isa AbstractVector
                @test AGP.nLatent(model) == 1
                @test AGP.likelihood(model) isa GaussianLikelihood
                model_opt = SVGP(
                    X,
                    y,
                    k,
                    GaussianLikelihood(opt_noise = true),
                    AnalyticVI(),
                    10,
                    optimiser = true,
                    Zoptimiser = true,
                    verbose = 0,
                )
                tests(model, model_opt, X, f, y, "Regression")
            end
            @test_throws ErrorException SVGP(
                X,
                y,
                k,
                GaussianLikelihood(),
                QuadratureVI(),
                20,
            )
            @test_throws ErrorException SVGP(
                X,
                y,
                k,
                GaussianLikelihood(),
                MCIntegrationVI(),
                20,
            )
        end
    end
    @testset "MCGP" begin
        @test_throws ErrorException MCGP(
            X,
            y,
            k,
            GaussianLikelihood(),
            GibbsSampling(),
        )
        @test_throws ErrorException MCGP(
            X,
            y,
            k,
            GaussianLikelihood(),
            HMCSampling(),
        )
    end
end


text:
"""
    init_model_triangular_dm_c3!(
        J :: Vector{Vector{Float64}},
        l :: Lattice
        ) :: Nothing

Init spin model with out-of-plane Dzyaloshinskii-Moriya (DM) interactions on the triangular lattice.
The DM interaction is chosen such that it switches sign between associated neighboring bonds, thus breaking the sixfold rotation symmetry of the triangular lattice down to threefold.
Here, J[n] = [Jxx, Jzz, Jxy] (n <= 3), with
* `Jxx` : n-th nearest neighbor coupling between the xx and yy spin components
* `Jzz` : n-th nearest neighbor coupling between the zz spin components
* `Jxy` : n-th nearest neighbor coupling between the xy spin components (aka DM interaction)
"""
function init_model_triangular_dm_c3!(
    J :: Vector{Vector{Float64}},
    l :: Lattice
    ) :: Nothing

    # sanity checks
    @assert l.name == "triangular" "Model requires triangular lattice."
    @assert length(J) <= 3 "Model initialization only works up to third-nearest neighbors."

    for n in eachindex(J)
        @assert length(J[n]) == 3 "Jxx, Jzz and Jxy all need to be specified for J[$(n)]."
    end

    # increase test set according to J
    max_nbs = get_nbs(length(J), l.sites[1], l.sites)
    metric  = maximum(Int64[get_metric(l.sites[1], l.sites[i], l.uc) for i in max_nbs])
    grow_test_sites!(l, metric)

    # save reference vectors to ensure uniform initialization
    ref_vecs = Vector{Float64}[]

    # iterate over sites and add couplings to lattice bonds
    for i in eachindex(l.sites)
        for n in eachindex(J)
            # find n-th nearest neighbors
            nbs = get_nbs(n, l.sites[i], l.sites)

            # determine couplings
            Jxx, Jzz, Jxy = J[n][1], J[n][2], J[n][3]

            # determine (and save) reference vector
            if length(ref_vecs) < n
                push!(ref_vecs, l.sites[nbs[1]].vec .- l.sites[i].vec)
            end

            ref_vec = ref_vecs[n]

            # iterate over neighbors and overwrite bond matrices
            for j in nbs
                vec = l.sites[j].vec .- l.sites[i].vec
                p   = round(dot(ref_vec, vec) / (norm(ref_vec) * norm(vec)), digits = 8)

                # disentanle bond dependence of DM coupling, preserving C3 but not C6 symmetry
                if acos(p) % (2.0 * pi / 3.0) < 1e-8
                    add_bond!(+Jxy, l.bonds[i, j], 1, 2)
                    add_bond!(-Jxy, l.bonds[i, j], 2, 1)
                else
                    add_bond!(-Jxy, l.bonds[i, j], 1, 2)
                    add_bond!(+Jxy, l.bonds[i, j], 2, 1)
                end

                # add other couplings to bond
                add_bond!(Jxx, l.bonds[i, j], 1, 1)
                add_bond!(Jxx, l.bonds[i, j], 2, 2)
                add_bond!(Jzz, l.bonds[i, j], 3, 3)
            end
        end
    end

    return nothing
end