text:
module Col

import Base: ==

struct GroupedColumn
    columns::Union{Nothing,Vector}
end

Base.hash(colgroup::GroupedColumn, h::UInt64) = hash(colgroup.columns, h)

==(a::GroupedColumn, b::GroupedColumn) = a.columns==b.columns

Base.show(io::IO, gc::GroupedColumn) = print(io, "Column")

index() = GroupedColumn(nothing)

index(xs::T...) where {T <: (Union{Int, Symbol})} = GroupedColumn(collect(Union{Nothing,T}, xs))

struct GroupedColumnValue
    columns::Union{Nothing,Vector}
end

Base.show(io::IO, gc::GroupedColumnValue) = print(io, "Column Value")

value() = GroupedColumnValue(nothing)

value(xs::T...) where {T <: (Union{Int, Symbol})} = GroupedColumnValue(collect(Union{Nothing,T}, xs))

end # module Col

module Row

# represent a row index correspondig to a set of columns
struct GroupedColumnRowIndex
    columns::Union{Nothing,Vector}
end

index() = GroupedColumnRowIndex(nothing)

index(xs::T...) where {T <: (Union{Int, Symbol})} = GroupedColumnRowIndex(collect(Union{Nothing,T}, xs))

end # module Row

# Handle aesthetics aliases and warn about unrecognized aesthetics.
#
# Returns:
#   A new mapping with aliases evaluated and unrecognized aesthetics removed.
#
function cleanmapping(mapping::Dict)
    cleaned = Dict{Symbol, Any}()
    for (key, val) in mapping
        # skip the "order" pesudo-aesthetic, used to order layers
        key == :order && continue

        if haskey(aesthetic_aliases, key)
            # redefining key could have performance impacts
            newkey = aesthetic_aliases[key]
        elseif in(key, valid_aesthetics)
            newkey = key
        else
            @warn "$(string(key)) is not a recognized aesthetic. Ignoring."
            continue
        end

        newval = val == Col.value || val == Col.index || val == Row.index ? val() : val

        cleaned[newkey] = newval
    end
    cleaned
end

# Type to contain fields produced by melting data (and to dispatch on)
struct MeltedData{T}
    data
    melted_data::T
    row_indicators::Array
    col_indicators::Array
    colmap::Dict
end

function meltdata(U::AbstractVector, colgroups_::Vector{Col.GroupedColumn})
    colgroups = Set(colgroups_)

    if length(colgroups) != 1 || first(colgroups).columns!==nothing
        # if every column is of the same length, treat it as a matrix
        if length(Set([length(u for u in U)])) == 1
            return meltdata(cat(U..., dims=2), colgroups_)
        end

        # otherwise it doesn't make much sense
        error("Col.index/Col.value can only be used without arguments when plotting an array of heterogenous arrays")
    end
    colgroup = first(colgroups)
    colmap = Dict{Any, Int}()
    colmap[colgroup] = 1

    V = cat(U..., dims=1)
    col_indicators = Array{Int}(undef, length(V))
    row_indicators = Array{Int}(undef, length(V))
    k = 1
    for i in 1:length(U)
        for j in 1:length(U[i])
            col_indicators[k] = i
            row_indicators[k] = j
            k += 1
        end
    end

    return MeltedData(U, V, row_indicators, col_indicators, colmap)
end

function meltdata(U::AbstractMatrix, colgroups_::Vector{Col.GroupedColumn})
    um, un = size(U)

    colgroups = Set(colgroups_)

    # Figure out the size of the new melted matrix
    allcolumns = BitSet()
    for i in 1:un
        push!(allcolumns, i)
    end

    vm = um
    grouped_columns = BitSet()
    for colgroup in colgroups
        if colgroup.columns===nothing
            vm *= un
            grouped_columns = copy(allcolumns)
        else
            for j in colgroup.columns
                push!(grouped_columns, j)
            end
            vm *= length(colgroup.columns)
        end
    end

    ungrouped_columns = setdiff(allcolumns, grouped_columns)
    vn = length(colgroups) + length(ungrouped_columns)
    V = similar(U, (vm, vn))

    # Indicator columns for each colgroup
    col_indicators = Array{Int}(undef, vm, length(colgroups))
    row_indicators = Array{Int}(undef, vm, length(colgroups))

    colidxs = [colgroup.columns===nothing ? collect(allcolumns) : colgroup.columns
               for colgroup in colgroups]

    vi = 1
    for ui in 1:um, colidx in product(colidxs...)
        # copy grouped columns
        for (vj, uj) in enumerate(colidx)
            V[vi, vj] = U[ui, uj]
            col_indicators[vi, vj] = uj
            row_indicators[vi, vj] = ui
        end

        # copy ungrouped columns
        for (vj, uj) in enumerate(ungrouped_columns)
            V[vi, vj + length(colgroups)] = U[ui, uj]
        end

        vi += 1
    end

    # Map grouped and individual columns in U to columns in V
    colmap = Dict{Any, Int}()
    for (vj, colgroup) in enumerate(colgroups)
        colmap[colgroup] = vj
    end
    for (vj, uj) in enumerate(ungrouped_columns)
        colmap[uj] = vj + length(colgroups)
    end

    return MeltedData(U, V, row_indicators, col_indicators, colmap)
end

# Evaluate one mapping.
evalmapping(source, arg::AbstractArray) = arg
evalmapping(source, arg::Function) = arg
evalmapping(source, arg::Distribution) = arg

evalmapping(source::MeltedData, arg::Integer) = source.melted_data[:,source.colmap[arg]]
evalmapping(source::MeltedData, arg::Col.GroupedColumn) = source.col_indicators[:,source.colmap[arg]]

evalmapping(source::MeltedData, arg::Col.GroupedColumnValue) =
    source.melted_data[:,source.colmap[Col.GroupedColumn(arg.columns)]]
evalmapping(source::MeltedData, arg::Row.GroupedColumnRowIndex) =
    source.row_indicators[:,source.colmap[Col.GroupedColumn(arg.columns)]]
evalmapping(source::MeltedData, arg::Symbol) =
    source.melted_data[:,source.colmap[arg]]
evalmapping(source::MeltedData, arg::AbstractString) =
    source.melted_data[:,source.colmap[Symbol(arg)]]
evalmapping(source::MeltedData, arg::Colon) = source.melted_data

# Share common functionality between general and melted-data-specific methods
function _evalmapping!(mapping::Dict, data_source, data::Data)
    for (k, v) in mapping
        setfield!(data, k, evalmapping(data_source, v))
        data.titles[k] = isa(v, AbstractString) || isa(v, Symbol) ?  string(v) : string(k)
    end

    return data_source
end

# Evalute aesthetic mappings producting a Data instance.
evalmapping!(mapping::Dict, data_source::MeltedData, data::Data) =
    _evalmapping!(mapping, data_source, data)

# Need to do additional work only if data_source is not MeltedData
function evalmapping!(mapping::Dict, data_source, data::Data)
    # Are we doing implicit reshaping?
    colgroups = Col.GroupedColumn[]
    for (k, v) in mapping
        if isa(v, Col.GroupedColumn)
            push!(colgroups, v)
        elseif isa(v, Col.GroupedColumnValue) || isa(v, Row.GroupedColumnRowIndex)
            push!(colgroups, Col.GroupedColumn(v.columns))
        end
    end

    transformed_data_source = if !isempty(colgroups)
        meltdata(data_source, colgroups)
    else
        data_source
    end

    return _evalmapping!(mapping, transformed_data_source, data)
end

function link_dataframes()
    @debug "Loading DataFrames support into Gadfly"
    include("dataframes.jl")
end


text:
export pdfnegbin_truncated

cdfnegbin(k; p, r) = one(k) - cdf(Beta(k+one(k),r),p)
pdfnegbin(k;p,r) = binomial( k+r - one(k), k) * (one(p) - p)^r * (p^k)
meannegbin(p,r) = (p*r) / (1 - p)

function pdfnegbin_truncated(k; p, r, maxwins)
    if k < maxwins
        pdfnegbin(k; p = p, r = r)
    elseif k == maxwins
        pdfnegbin(maxwins; p = p, r = r) + (1 - cdfnegbin(maxwins; p = p, r = r))
    else
        zero(pdfnegbin(k; p = p, r = r))
    end
end

ext: jl

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
" functions for Julia integration with Modeller (http://salilab.org/modeller/)
  adapted from the modeller example scripts
"

module Modeller
export  gen_modeller_script, build_profile, compare, model_single, evaluate_model, align2d, plot_profiles
    using PyCall
    using PyPlot

    # Generator for MODELLER julia scripts.
    # These scripts are based on the basic example scripts from the MODELLER website.
    # Scripts are generated in the current working directory
    # name: The name of the script (minus the extension), e.g. "build_profile"
    function gen_modeller_script(name::AbstractString)
       file = joinpath(dirname(@__FILE__), "..", "MODELLER/modeller-basic-example-julia", "$name.jl")
       if isfile(file)
          cp(file, "./$name.jl")
          println("Generated $name.jl MODELLER script")
       end
    end

    function build_profile(;seq_database_file::AbstractString = "", seq_database_format::AbstractString="PIR", sequence_file::AbstractString = "", sequence_format::AbstractString = "PIR", output_name::AbstractString = "build_profile", output_profile_format::AbstractString="TEXT", output_alignment_format::AbstractString="PIR")
        seq_database_file_out = string(seq_database_file, ".bin")

        pyinitialize()
        @pyimport modeller
        @pyimport _modeller
        mod_lib = _modeller.mod_libdir_get()
        modeller.log[:verbose]()
        env = modeller.environ()

        #-- Prepare the input files

        #-- Read in the sequence database
        sdb = modeller.sequence_db(env)
        sdb[:read](seq_database_file=seq_database_file, seq_database_format=seq_database_format,
                 chains_list="ALL", minmax_db_seq_len=(30, 4000), clean_sequences=true)

        #-- Write the sequence database in binary form
        sdb[:write](seq_database_file=seq_database_file_out, seq_database_format="BINARY",
                  chains_list="ALL")

        #-- Now, read in the binary database
        sdb[:read](seq_database_file=seq_database_file_out, seq_database_format="BINARY",
                 chains_list="ALL")

        #-- Read in the target sequence/alignment
        aln = modeller.alignment(env)
        aln[:append](file=sequence_file, alignment_format=sequence_format, align_codes="ALL")

        #-- Convert the input sequence/alignment into
        #   profile format
        prf = aln[:to_profile]()

        #-- Scan sequence database to pick up homologous sequences
        prf[:build](sdb, matrix_offset=-450, rr_file=string(mod_lib,"/blosum62.sim.mat"),
                  gap_penalties_1d=(-500, -50), n_prof_iterations=1,
                  check_profile=false, max_aln_evalue=0.01)

        #-- Write out the profile in text format
        prf[:write](file=string(output_name, ".prf"), profile_format=output_profile_format)

        #-- Convert the profile back to alignment format
        aln = prf[:to_alignment]()

        #-- Write out the alignment file
        aln[:write](file=string(output_name, ".ali"), alignment_format=output_alignment_format)

    end

    function compare(pdbs)
        @pyimport modeller

        env = modeller.environ()
        aln = modeller.alignment(env)
        for (pdb, chain) in pdbs
            m = modeller.model(env, file=pdb, model_segment=("FIRST:$chain", "LAST:$chain"))
            aln[:append_model](m, atom_files=pdb, align_codes=string(pdb, chain))
        end
        aln[:malign]()
        aln[:malign3d]()
        aln[:compare_structures]()
        aln[:id_table](matrix_file="family.mat")
        env[:dendrogram](matrix_file="family.mat", cluster_cut=-1.0)
    end

    function model_single(alnf::AbstractString, known_structure::AbstractString, seq::AbstractString)

        @pyimport modeller
        @pyimport modeller.automodel as am

        env = modeller.environ()
        a = am.automodel(env, alnfile=alnf,
              knowns=known_structure, sequence=seq,
              assess_methods=(am.assess[:DOPE],
                              #soap_protein_od.Scorer(),
                              am.assess[:GA341]))
        a[:starting_model] = 1
        a[:ending_model] = 5
        a[:make]()
    end

    function evaluate_model(pdbfile::AbstractString, outputfile::AbstractString = "")

        # check for optional argument output and set it to the first part of the name of the PDB file
        if outputfile == ""
            outputfile = string(pdbfile, ".profile")
        end

        @pyimport modeller
        @pyimport modeller.scripts as scripts
        @pyimport _modeller

        mod_lib = _modeller.mod_libdir_get()

        modeller.log[:verbose]()    # request verbose output
        env = modeller.environ()
        env[:libs][:topology][:read](file=string(mod_lib, "/top_heav.lib")) # read topology
        env[:libs][:parameters][:read](file=string(mod_lib, "/par.lib")) # read parameters

        # read model file
        mdl = scripts.complete_pdb(env, pdbfile)

        # Assess with DOPE:
        s = modeller.selection(mdl)   # all atom selection
        s[:assess_dope](output="ENERGY_PROFILE NO_REPORT", file=outputfile,
                  normalize_profile=true, smoothing_window=15)
    end

   function align2d(model_file::AbstractString, model_segment, model_codes::AbstractString, atom_files::AbstractString, sequence_file::AbstractString, sequence_codes::AbstractString, outputname::AbstractString)
       @pyimport modeller

       env = modeller.environ()
       aln = modeller.alignment(env)
       mdl = modeller.model(env, file=model_file, model_segment=model_segment)
       aln[:append_model](mdl, align_codes=model_codes, atom_files=atom_files)
       aln[:append](file=sequence_file, align_codes=sequence_codes)
       aln[:align2d]()
       aln[:write](file=string(outputname, ".ali"), alignment_format="PIR")
       aln[:write](file=string(outputname, ".pap"), alignment_format="PAP")
    end

    function plot_profiles(alignment_file::AbstractString, template_profile::AbstractString, template_sequence::AbstractString, model_profile::AbstractString, model_sequence::AbstractString, plot_file::AbstractString = "dope_profile.png")

        @pyimport modeller

        function r_enumerate(seq)
            """Enumerate a sequence in reverse order"""
            # Note that we don"t use reversed() since Python 2.3 doesn"t have it
            num = pybuiltin(:len)(seq) - 1
            while num >= 0
                produce((num, get(seq, num)))
                num -= 1
            end
        end

        function get_profile(profile_file, seq)
            """Read `profile_file` into a Python array, and add gaps corresponding to
               the alignment sequence `seq`."""
            # Read all non-comment and non-blank lines from the file:
            f = open(profile_file)
            vals = Any[]
            for line in eachline(f)
                if line[1] != '#' && length(line) > 10
                    spl = split(line)
                    push!(vals, float(last(spl)))
                end
            end
            close(f)
            # Insert gaps into the profile corresponding to those in seq:
            taskAtr = @task r_enumerate(seq[:residues])
            for x in taskAtr
                n = x[1]
                res = x[2]
                for gap in range(1, res[:get_leading_gaps]())
                    # Julia arrays start from 1 so we have to add 1 to the position
                    insert!(vals, n+1, nothing)
                end
            end
            # Add a gap at position "0", so that we effectively count from 1:
            insert!(vals, 1, nothing)
            return vals
        end

        e = modeller.environ()
        a = modeller.alignment(e, file=alignment_file)

        template = get_profile(template_profile, get(a, template_sequence))
        model = get_profile(model_profile, get(a, model_sequence))

        # Plot the template and model profiles in the same plot for comparison:
        PyPlot.figure(1, figsize=(10,6))
        PyPlot.xlabel("Alignment position")
        PyPlot.ylabel("DOPE per-residue score")
        PyPlot.plot(model, color="red", linewidth=2, label="Model")
        PyPlot.plot(template, color="green", linewidth=2, label="Template")
        PyPlot.legend()
        PyPlot.savefig(plot_file, dpi=65)
    end
end


text:
immutable Rational{T<:Integer} <: Real
    num::T
    den::T

    function Rational(num::T, den::T)
        if num == 0 && den == 0
            error("invalid rational: 0//0")
        end
        g = den < 0 ? -gcd(den,num) : gcd(den, num)
        new(div(num, g), div(den, g))
    end
end
Rational{T<:Integer}(n::T, d::T) = Rational{T}(n,d)
Rational(n::Integer, d::Integer) = Rational(promote(n,d)...)
Rational(n::Integer) = Rational(n,one(n))

//(n::Integer,  d::Integer ) = Rational(n,d)
//(x::Rational, y::Integer ) = x.num//(x.den*y)
//(x::Integer,  y::Rational) = (x*y.den)//y.num
//(x::Rational, y::Rational) = (x.num*y.den)//(x.den*y.num)
//(x::Complex,  y::Real    ) = complex(real(x)//y,imag(x)//y)
//(x::Real,     y::Complex ) = x*y'//real(y*y')

function //(x::Complex, y::Complex)
    xy = x*y'
    yy = real(y*y')
    complex(real(xy)//yy, imag(xy)//yy)
end

function show(io::IO, x::Rational)
    if isinf(x)
        print(io, x.num > 0 ? "Inf" : "-Inf")
    else
        show(io, num(x)); print(io, "//"); show(io, den(x))
    end
end

convert{T<:Integer}(::Type{Rational{T}}, x::Rational) = Rational(convert(T,x.num),convert(T,x.den))
convert{T<:Integer}(::Type{Rational{T}}, x::Integer) = Rational(convert(T,x), convert(T,1))

convert(::Type{Rational}, x::Rational) = x
convert(::Type{Rational}, x::Integer) = convert(Rational{typeof(x)},x)

convert(::Type{Bool}, x::Rational) = (x!=0) # to resolve ambiguity
convert{T<:Integer}(::Type{T}, x::Rational) = (isinteger(x) ? convert(T, x.num) : throw(InexactError()))

convert(::Type{FloatingPoint}, x::Rational) = float(x.num)/float(x.den)
function convert{T<:FloatingPoint,S}(::Type{T}, x::Rational{S})
    P = promote_type(T,S)
    convert(P,x.num)/convert(P,x.den)
end

function convert{T<:Integer}(::Type{Rational{T}}, x::FloatingPoint)
    r = rationalize(T, x, tol=0)
    x === convert(typeof(x), r) || throw(InexactError())
    r
end
convert(::Type{Rational}, x::Float64) = convert(Rational{Int64}, x)
convert(::Type{Rational}, x::Float32) = convert(Rational{Int}, x)

promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{S}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{Rational{S}}) = Rational{promote_type(T,S)}
promote_rule{T<:Integer,S<:FloatingPoint}(::Type{Rational{T}}, ::Type{S}) = promote_type(T,S)

widen{T}(::Type{Rational{T}}) = Rational{widen(T)}

function rationalize{T<:Integer}(::Type{T}, x::FloatingPoint; tol::Real=eps(x))
    if isnan(x);       return zero(T)//zero(T); end
    if x < typemin(T); return -one(T)//zero(T); end
    if typemax(T) < x; return  one(T)//zero(T); end
    tm = x < 0 ? typemin(T) : typemax(T)
    z = x*tm
    if z <= 0.5 return zero(T)//one(T) end
    if z <= 1.0 return one(T)//tm end
    y = x
    a = d = 1
    b = c = 0
    while true
        f = itrunc(y); y -= f
        p, q = f*a+c, f*b+d
        typemin(T) <= p <= typemax(T) &&
        typemin(T) <= q <= typemax(T) || break
        0 != sign(a)*sign(b) != sign(p)*sign(q) && break
        a, b, c, d = p, q, a, b
        if y == 0 || abs(a/b-x) <= tol
            break
        end
        y = inv(y)
    end
    return convert(T,a)//convert(T,b)
end
rationalize(x::Union(Float64,Float32); tol::Real=eps(x)) = rationalize(Int, x, tol=tol)

num(x::Integer) = x
den(x::Integer) = one(x)
num(x::Rational) = x.num
den(x::Rational) = x.den

sign(x::Rational) = sign(x.num)
signbit(x::Rational) = signbit(x.num)
copysign(x::Rational, y::Real) = copysign(x.num,y) // x.den
copysign(x::Rational, y::Rational) = copysign(x.num,y.num) // x.den

isnan(x::Rational) = false
isinf(x::Rational) = x.den == 0
isfinite(x::Rational) = x.den != 0

typemin{T<:Integer}(::Type{Rational{T}}) = -one(T)//zero(T)
typemax{T<:Integer}(::Type{Rational{T}}) = one(T)//zero(T)

isinteger(x::Rational) = x.den == 1

zero{T}(::Type{Rational{T}}) = zero(T)//one(T)
zero{T}(q::Rational{T}) = zero(T)//one(T)
one{T}(::Type{Rational{T}}) = one(T)//one(T)
one{T}(q::Rational{T}) = one(T)//one(T)

hash(x::Rational) = bitmix(hash(x.num), ~hash(x.den))

-(x::Rational) = (-x.num) // x.den
for op in (:+, :-, :rem, :mod)
    @eval begin
        function ($op)(x::Rational, y::Rational)
            g = gcd(x.den, y.den)
            Rational(($op)(x.num * div(y.den, g), y.num * div(x.den, g)), x.den * div(y.den, g))
        end
    end
end
*(x::Rational, y::Rational) = (x.num*y.num) // (x.den*y.den)
/(x::Rational, y::Rational) = (x.num*y.den) // (x.den*y.num)
/(x::Rational, z::Complex ) = inv(z/x)

==(x::Rational, y::Rational) = (x.den == y.den) & (x.num == y.num)
==(x::Rational, y::Integer ) = (x.den == 1) & (x.num == y)
==(x::Integer , y::Rational) = y == x

# needed to avoid ambiguity between ==(x::Real, z::Complex) and ==(x::Rational, y::Number)
==(z::Complex , x::Rational) = isreal(z) & (real(z) == x)
==(x::Rational, z::Complex ) = isreal(z) & (real(z) == x)

==(x::FloatingPoint, q::Rational) = ispow2(q.den) & (x == q.num/q.den) & (x*q.den == q.num)
==(q::Rational, x::FloatingPoint) = ispow2(q.den) & (x == q.num/q.den) & (x*q.den == q.num)

# TODO: fix inequalities to be in line with equality check
< (x::Rational, y::Rational) = x.den == y.den ? x.num < y.num :
                               widemul(x.num,y.den) < widemul(x.den,y.num)
< (x::Rational, y::Integer ) = x.num < widemul(x.den,y)
< (x::Rational, y::Real    ) = x.num < x.den*y
< (x::Integer , y::Rational) = widemul(x,y.den) < y.num
< (x::Real    , y::Rational) = x*y.den < y.num

<=(x::Rational, y::Rational) = x.den == y.den ? x.num <= y.num :
                               widemul(x.num,y.den) <= widemul(x.den,y.num)
<=(x::Rational, y::Integer ) = x.num <= widemul(x.den,y)
<=(x::Rational, y::Real    ) = x.num <= x.den*y
<=(x::Integer , y::Rational) = widemul(x,y.den) <= y.num
<=(x::Real    , y::Rational) = x*y.den <= y.num

div(x::Rational, y::Rational) = div(x.num*y.den, x.den*y.num)
div(x::Rational, y::Real    ) = div(x.num, x.den*y)
div(x::Real    , y::Rational) = div(x*y.den, y.num)

fld(x::Rational, y::Rational) = fld(x.num*y.den, x.den*y.num)
fld(x::Rational, y::Real    ) = fld(x.num, x.den*y)
fld(x::Real    , y::Rational) = fld(x*y.den, y.num)

itrunc(x::Rational) = div(x.num,x.den)
ifloor(x::Rational) = fld(x.num,x.den)
iceil (x::Rational) = -fld(-x.num,x.den)
iround(x::Rational) = div(x.num*2 + copysign(x.den,x.num), x.den*2)

trunc(x::Rational) = Rational(itrunc(x))
floor(x::Rational) = Rational(ifloor(x))
ceil (x::Rational) = Rational(iceil(x))
round(x::Rational) = Rational(iround(x))

## rational to int coercion ##

for f in (:int8, :int16, :int32, :int64, :int128,
          :uint8, :uint16, :uint32, :uint64, :uint128,
          :signed, :integer, :unsigned, :int, :uint)
    @eval ($f)(x::Rational) = ($f)(iround(x))
end

^(x::Rational, y::Integer) = y < 0 ?
    Rational(x.den^-y, x.num^-y) : Rational(x.num^y, x.den^y)

^(x::Number, y::Rational) = x^(y.num/y.den)
^{T<:FloatingPoint}(x::T, y::Rational) = x^(convert(T,y.num)/y.den)
^{T<:FloatingPoint}(x::Complex{T}, y::Rational) = x^(convert(T,y.num)/y.den)

^{T<:Rational}(z::Complex{T}, n::Bool) = n ? z : one(z) # to resolve ambiguity
^{T<:Rational}(z::Complex{T}, n::Integer) =
    n>=0 ? power_by_squaring(z,n) : power_by_squaring(inv(z),-n)

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using TriMeshRender
using Test

Vertex = TriMeshRender.Vertex

@testset "TriMeshRender.jl" begin
	@test isometric(Vertex(0,0,0)) == Vertex(0,0,0)
	@test isometric(Vertex(0,1,0)) == Vertex(0,1,0)
	@test isometric(Vertex(1,0,0)) ≈ Vertex(cos(deg2rad(30)), sin(deg2rad(30)), 0)
	@test isometric(Vertex(2,0,0)) ≈ Vertex(2cos(deg2rad(30)), sin(2deg2rad(30)), 0)
	@test isometric(Vertex(1,1,0)) ≈ Vertex(cos(deg2rad(30)), 1+sin(deg2rad(30)), 0)
	@test isometric(Vertex(0,0,1)) ≈ Vertex(-cos(deg2rad(30)), sin(2eg2rad(30)), 0)
	@test isometric(Vertex(0,0,2)) ≈ Vertex(-2cos(deg2rad(30)), sin(2deg2rad(30)), 0)
end


text:
# For generating symbols that work as dataframe column names; symbols
# generated by gensym() don't work for this.
global _rvnum = 0

"""
    _make_rvname(name)

Return a unique name for random variable `rv`.
"""
function _make_rvname(name)
    global _rvnum += 1
    return Symbol("$(name)!$(_rvnum)")
end

"""
    _make_dims(args)

Return a vector of dimensions from `args`.
"""
function _make_dims(args)
    dims = Vector{Any}()
    for arg in args
        # println("Arg: $arg")

        if @capture(arg, i_Int)  # scalar (must be integer)
            # println("arg is Int")
            dim = i

        elseif @capture(arg, first_Int:last_)   # last can be an int or 'end', which is converted to 0
            # println("arg is range")
            last = last == :end ? 0 : last
            dim = :($first:$last)

        elseif @capture(arg, first_Int:step_:last_)
            # println("arg is step range")
            last = last == :end ? 0 : last
            dim = :($first:$step:$last)

        elseif @capture(arg, s_Symbol)
            if arg == :(:)
                # println("arg is Colon")
                dim = Colon()
            else
                # println("arg is Symbol")
                dim = :($(QuoteNode(s)))
            end

        elseif @capture(arg, s_String)
            dim = s

        elseif isa(arg, Expr) && arg.head == :tuple  # tuple of Strings/Symbols (@capture didn't work...)
            argtype = typeof(arg.args[1])            # ensure all are same type as first element
            if ! isempty(filter(s -> typeof(s) != argtype, arg.args))
                error("A parameter dimension tuple must be all Strings or all Symbols (got $arg)")
            end
            dim = :(convert(Vector{$argtype}, $(arg.args)))

        else
            error("Unrecognized stochastic parameter specification: $arg")
        end
        push!(dims, dim)
    end
    # println("dims = $dims")
    return dims
end

"""
    defsim(expr::Expr)

Define a Mimi `SimulationDef` with the expressions in `expr`.
"""
macro defsim(expr::Expr)
    let # to make vars local to each macro invocation
        local _rvs        = []
        local _transforms = []
        local _saves      = []
        local _sim_args  = Dict{Symbol, Any}()

        simdatatype = nothing

        # distilled into a function since it's called from two branches below
        function saverv(rvname, distname, distargs)
            expr = :(RandomVariable($(QuoteNode(rvname)), $distname($(distargs...))))
            push!(_rvs, esc(expr))
        end

        @capture(expr, elements__)
        for elt in elements
            # Meta.show_sexpr(elt)
            # println("")
            # e.g.,  rv(name1) = Normal(10, 3)
            if @capture(elt, rv(rvname_) = distname_(distargs__))
                saverv(rvname, distname, distargs)

            elseif @capture(elt, save(vars__))
                for var in vars
                    # println("var: $var")
                    if @capture(var, comp_.datum_)
                        expr = :($(QuoteNode(comp)), $(QuoteNode(datum)))
                        push!(_saves, esc(expr))
                    else
                        error("Save arg spec must be of the form comp_name.datum_name; got ($var)")
                    end
                end

            # handle vector of distributions
            elseif @capture(elt, extvar_ = [items__])
                for pair in items
                    if (@capture(pair, [dims__] => distname_(distargs__)) ||
                        @capture(pair,     dim_ => distname_(distargs__)))

                        dims = _make_dims(dims === nothing ? [dim] : dims)

                        rvname = _make_rvname(extvar)
                        saverv(rvname, distname, distargs)

                        if @capture(extvar, compname_.name_)
                            expr = :(TransformSpec($(QuoteNode(compname)), $(QuoteNode(name)), :(=), $(QuoteNode(rvname)), [$(dims...)]))
                        else
                            expr = :(TransformSpec($(QuoteNode(extvar)), :(=), $(QuoteNode(rvname)), [$(dims...)]))
                        end
                        push!(_transforms, esc(expr))
                    end
                end

            # e.g., ext_var5[2010:2050, :] *= name2
            elseif (@capture(elt, extvar_  = rvname_Symbol) ||
                    @capture(elt, extvar_ += rvname_Symbol) ||
                    @capture(elt, extvar_ *= rvname_Symbol) ||
                    @capture(elt, extvar_  = distname_(distargs__)) ||
                    @capture(elt, extvar_ += distname_(distargs__)) ||
                    @capture(elt, extvar_ *= distname_(distargs__)))

                # For "anonymous" RVs, e.g., ext_var2[2010:2100, :] *= Uniform(0.8, 1.2), we
                # gensym a name based on the model parameter name and process it as a named RV.
                if rvname === nothing
                    param_name = @capture(extvar, name_[args__]) ? name : extvar
                    rvname = _make_rvname(param_name)
                    saverv(rvname, distname, distargs)
                end

                op = elt.head
                # println("Ref:  $name, $args")
                # Meta.show_sexpr(extvar)
                # println("")
                # if extvar.head == :ref, extvar.args must be one of the following,
                # where the extvar could be paramname or compname.paramname:
                    # - a scalar value, e.g., name[2050] => (:ref, :name, 2050)
                    #   convert to tuple of dimension specifiers (:name, 2050)
                    # - a slice expression, e.g., name[2010:2050] => (:ref, :name, (:(:), 2010, 2050))
                    #   convert to (:name, 2010:2050) [convert it to actual UnitRange instance]
                    # - a tuple of symbols, e.g., name[(US, CHI)] => (:ref, :name, (:tuple, :US, :CHI))
                    #   convert to (:name, (:US, :CHI))
                    # - combinations of these, e.g., name[2010:2050, (US, CHI)] => (:ref, :name, (:(:), 2010, 2050), (:tuple, :US, :CHI))
                    #   convert to (:name, 2010:2050, (:US, :CHI))

                if @capture(extvar, compname_.name_[args__])
                    dims = _make_dims(args)
                    expr = :(TransformSpec($(QuoteNode(compname)), $(QuoteNode(name)), $(QuoteNode(op)), $(QuoteNode(rvname)), [$(dims...)]))
                elseif @capture(extvar, compname_.name_)
                    expr = :(TransformSpec($(QuoteNode(compname)), $(QuoteNode(name)), $(QuoteNode(op)), $(QuoteNode(rvname))))

                elseif @capture(extvar, name_[args__])
                    dims = _make_dims(args)
                    expr = :(TransformSpec($(QuoteNode(name)), $(QuoteNode(op)), $(QuoteNode(rvname)), [$(dims...)]))
                else
                    expr = :(TransformSpec($(QuoteNode(extvar)), $(QuoteNode(op)), $(QuoteNode(rvname))))
                end
                push!(_transforms, esc(expr))

            elseif @capture(elt, sampling(simdatatype_, simargs__))
                # General method for setting an SA method's required parameters
                for arg in simargs
                    if !@capture(arg, name_=value_)
                        error("simargs must be given as keyword arguments (got $simargs)")
                    end
                    _sim_args[name] = __module__.eval(value)
                end
            else
                error("Unrecognized expression '$elt' in @defsim")
            end
        end

        # set default
        simdatatype = (simdatatype == nothing) ? nameof(MCSData) : simdatatype

        # call constructor on given args
        data = :($simdatatype(; $(_sim_args)...))

        return :(SimulationDef{$simdatatype}(
                    [$(_rvs...)],
                    TransformSpec[$(_transforms...)],
                    Tuple{Symbol, Symbol}[$(_saves...)],
                    $data))
    end
end

#
# Simulation Definition update methods
#

function _update_nt_type!(sim_def::SimulationDef{T}) where T <: AbstractSimulationData
    names = (keys(sim_def.rvdict)...,)
    types = [eltype(fld) for fld in values(sim_def.rvdict)]
    sim_def.nt_type = NamedTuple{names, Tuple{types...}}
    nothing
end

"""
    delete_RV!(sim_def::SimulationDef, name::Symbol)

Delete the random variable with name `name` from the Simulation definition `sim_def`.
Transformations using this RV are deleted, and the Simulation's
NamedTuple type is updated to reflect the dropped RV.
"""
function delete_RV!(sim_def::SimulationDef, name::Symbol)
    if !haskey(sim_def.rvdict, name)
        @warn("Simulation def does not have RV :$name. Nothing being deleted.")
    else
        delete_transform!(sim_def, name)
        delete!(sim_def.rvdict, name)
        _update_nt_type!(sim_def)
    end
end

"""
    add_RV!(sim_def::SimulationDef, rv::RandomVariable)

Add random variable definition `rv` to Simulation definition `sim_def`. The
Simulation's NamedTuple type is updated to include the RV.
"""
function add_RV!(sim_def::SimulationDef, rv::RandomVariable)
    name = rv.name
    haskey(sim_def.rvdict, name) && error("Simulation def already has RV :$name. Use replace_RV! to replace it.")
    sim_def.rvdict[name] = rv
    _update_nt_type!(sim_def)
end

"""
    add_RV!(sim_def::SimulationDef, name::Symbol, dist::Distribution)

Add a random variable with name `name` and distribution `dist` to Simulation definition
`sim_def`. The Simulation's NamedTuple type is updated to include the RV.
"""
add_RV!(sim_def::SimulationDef, name::Symbol, dist::Distribution) = add_RV!(sim_def, RandomVariable(name, dist))

"""
    replace_RV!(sim_def::SimulationDef, rv::RandomVariable)

Replace the RV with the given `rv`s name in the Simulation definition `sim_def` with
`rv` and update the Simulation's NamedTuple type accordingly.
"""
function replace_RV!(sim_def::SimulationDef, rv::RandomVariable)
    name = rv.name
    haskey(sim_def.rvdict, name) || error("Simulation def does not have has RV :$name. Use add_RV! to add it.")
    sim_def.rvdict[rv.name] = rv
    _update_nt_type!(sim_def)
end

"""
    replace_RV!(sim_def::SimulationDef, name::Symbol, dist::Distribution)

Replace the rv with name `name` in the Simulation definition `sim_def` with a new RV
with `name` and distribution `dist`. Update the Simulation's NamedTuple
type accordingly.
"""
replace_RV!(sim_def::SimulationDef, name::Symbol, dist::Distribution) = replace_RV!(sim_def, RandomVariable(name, dist))

"""
    delete_transform!(sim_def::SimulationDef, name::Symbol)

Delete all data transformations--i.e., replacement, addition or multiplication
of original data values with values drawn from the RV named `name`. Update the
Simulation definition's NamedTuple type accordingly.
"""
function delete_transform!(sim_def::SimulationDef, name::Symbol)
    pos = findall(t -> t.rvname == name, sim_def.translist)
    if isempty(pos)
        @warn("Simulation def does not have and transformations using RV :$name. Nothing being deleted.")
    else
        deleteat!(sim_def.translist, pos)
        _update_nt_type!(sim_def)
    end
end

"""
    add_transform!(sim_def::SimulationDef, t::TransformSpec)

Add the data transformation `t` to the Simulation definition `sim_def`, and update the
Simulation's NamedTuple type. The TransformSpec `t` must refer to an
existing RV.
"""
function add_transform!(sim_def::SimulationDef, t::TransformSpec)
    push!(sim_def.translist, t)
    _update_nt_type!(sim_def)
end

"""
    add_transform!(sim_def::SimulationDef, paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T

Create a new TransformSpec based on `paramname`, `op`, `rvname` and `dims` to the
Simulation definition `sim_def`, and update the Simulation's NamedTuple type.
The symbol `rvname` must refer to an existing random variable, and `paramname`
must refer to an existing shared model parameter that can be accessed by that
name. Use the  signature that includes `compname` if your `paramname`
is an unshared model parameter specific to a component. If `dims` are
provided, these must be legal subscripts of `paramname`. Op must be one of :+=, :*=,
or :(=).
"""
function add_transform!(sim_def::SimulationDef, paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T
    add_transform!(sim_def, TransformSpec(nothing, paramname, op, rvname, dims))
end

"""
    add_transform!(sim_def::SimulationDef, compname_paramname::Tuple{Symbol, Symbol}, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T

Create a new TransformSpec based on `compname`, `paramname`, `op`, `rvname` and `dims` to the
Simulation definition `sim_def`, and update the Simulation's NamedTuple type. The symbol
`rvname` must refer to an existing RV, and `compname` and `paramname` must refer to
an existing parameter for a given component, parameter pair. If `dims` are provided,
these must be legal subscripts of `paramname`. Op must be one of :+=, :*=, or :(=).
"""
function add_transform!(sim_def::SimulationDef, compname::Symbol, paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T
    add_transform!(sim_def, TransformSpec(compname, paramname, op, rvname, dims))
end

"""
    delete_save!(sim_def::SimulationDef, key::Tuple{Symbol, Symbol})

Delete from Simulation definition `sim_def` a "save" instruction for the given `key`, comprising
component `comp_name` and parameter or variable `datum_name`. This result will no
longer be saved to a CSV file at the end of the simulation.
"""
function delete_save!(sim_def::SimulationDef, key::Tuple{Symbol, Symbol})
    pos = findall(isequal(key), sim_def.savelist)
    isempty(pos) ? @warn("Simulation def doesn't have $key in its save list. Nothing being deleted.") : deleteat!(sim_def.savelist, pos)
    _update_nt_type!(sim_def)
end

"""
    delete_save!(sim_def::SimulationDef, comp_name::Symbol, datum_name::Symbol)

Delete from Simulation definition `sim_def` a "save" instruction for component `comp_name` and parameter
or variable `datum_name`. This result will no longer be saved to a CSV file at the end
of the simulation.
"""
delete_save!(sim_def::SimulationDef, comp_name::Symbol, datum_name::Symbol) = delete_save!(sim_def, (comp_name, datum_name))

"""
    add_save!(sim_def::SimulationDef, key::Tuple{Symbol, Symbol})

Add to Simulation definition `sim_def` a "save" instruction for the given `key`, comprising
component `comp_name` and parameter or variable `datum_name`. This result will
be saved to a CSV file at the end of the simulation.
"""
function add_save!(sim_def::SimulationDef, key::Tuple{Symbol, Symbol})
    pos = findall(isequal(key), sim_def.savelist)
    !isempty(pos) ? @warn("Simulation def already has $key in its save list. Nothing being added.") : push!(sim_def.savelist, key)
end

"""
    add_save!(sim_def::SimulationDef, comp_name::Symbol, datum_name::Symbol)

Add to Simulation definition`sim_def` a "save" instruction for component `comp_name` and parameter or
variable `datum_name`. This result will be saved to a CSV file at the end of the simulation.
"""
add_save!(sim_def::SimulationDef, comp_name::Symbol, datum_name::Symbol) = add_save!(sim_def, (comp_name, datum_name))

"""
    get_simdef_rvnames(sim_def::SimulationDef, name::Union{String, Symbol})

A helper function to support finding the keys in a Simulation Definition `sim_def`'s
rvdict that contain a given `name`. This can be particularly helpful if the random
variable was set via the shortcut syntax ie. my_parameter = Normal(0,1) and thus the
`sim_def` has automatically created a key in the format `:my_parameter!xx`. In this
case this function is useful to get the key and carry out modification functions
such as `replace_rv!`.
"""
function get_simdef_rvnames(sim_def::SimulationDef, name::Union{String, Symbol})
    names = String.([keys(sim_def.rvdict)...])
    matches = Symbol.(filter(x -> occursin(String(name), x), names))
end


text:
"""
Collection of the methods used for generating the cuts and plates
that make up the PP-G2KP model.
"""
module Enumeration

# TODO: Consider removing flag ignore_2th_dim. Reasoning: nobody does that,
# even Furini only consider the pieces that fit a plate before discretizing it,
# even if ignoring the 2th dimension would save a lot of effort in the
# preprocessing nobody believes would be a valid trade-off to do so (the
# preprocessing phase participation on the total time is negligible), and it
# makes the code far more complex (because I have to have code that respects
# and ignores a basic assumption).

using DocStringExtensions
using UnPack: @unpack

using ...Utilities: SortedLinkedLW, switched_dims

export gen_cuts, ByproductPPG2KP
export should_extract_piece_from_plate

"""
    discretize(d::[D], l::[S], w::[S], L::S, W::S; kwargs...) :: [S]

!!! **Internal use.**

Given the pieces demand, length, and width (`d`, `l`, `w`), as well the plate
dimensions (`L`, `W`), return a vector of all linear combinations (constrained
by demand) of the length of pieces that fit the plate (the returned vector is
sorted by increasing value and no value is greater than `L`). To discretize the
width just swap `l` with `w` and `L` with `W`.

# Keyword arguments

All keyword arguments are of type `Bool` and have `false` as default.

* `only_single_pieces`: Instead of discretizing the length, return the
  piece lengths that exist only as the length of a piece and not as the
  combination of two or more smaller piece lengths.
* `ignore_W`: Do not use `W` to exclude pieces that do not fit the plate.
* `ignore_d`: Do not use `d`emand information to discretize (consider
  an unlimited amount of each piece type available).
"""
function discretize(
	d :: Vector{D}, l :: Vector{S}, w :: Vector{S}, L :: S, W :: S;
	only_single_pieces :: Bool = false, ignore_W :: Bool = false,
	ignore_d :: Bool = false
) where {D, S}
	# If two pieces have the same dimension they should be merged into a single
	# piece (summing their demands) for performance reasons.
	#@assert l == unique(l)
	# Each piece in 1:N has a demand and a dimension.
	@assert length(d) == length(l)
	N = length(d)
	# Remove items that cannot fit into L anyway.
	d_, l_, w_ = similar.((d, l, w))
	j = 0
	for i in 1:N
		if l[i] <= L && (ignore_W || w[i] <= W)
			j += 1
			d_[j] = d[i]
			l_[j] = l[i]
			w_[j] = w[i]
		end
	end
	resize!(d_, j)
	resize!(l_, j)
	resize!(w_, j)
	d, l, w = d_, l_, w_
	N = j
	iszero(N) && return S[]
	# marks: for each unit of the plate dimension, if there can be a cut there
	# (considering the pieces available) or not.
	marks = fill(false, L)
	# Mark the cuts of the first piece.
	!only_single_pieces && (marks[l[1]] = true)
	y = l[1] # y: used to iterate capacity, inherited from knapsack papers
	for _ = 2:(ignore_d ? L ÷ l[1] : min(d[1], L ÷ l[1]))
		marks[y += l[1]] = true
	end
	# If the lengths of the single pieces are not marked, the pairs are needed.
	if only_single_pieces
		for pii = 2:N, pij = 1:(pii - 1) # PIece `j` (as opposed to `i`)
			lij = l[pii] + l[pij]
			lij <= L && (marks[lij] = true)
		end
	end
	# Mark the cuts of all other pieces.
	for pii = 2:N # PIece Index
		li = l[pii] # length of i
		di = d[pii] # demand of i
		for y = (L - li):-1:1
			if marks[y]
				yrli = y + li # yrli: y + repeated lengths of i
				marks[yrli] = true
				for r = 2:(ignore_d ? L : di)
					yrli += li
					yrli > L && break
					marks[yrli] = true
				end
			end
		end
		# Mark cuts using just multiple copies of the same piece type (could be
		# done using a dummy cut at index zero with value true, but the array has
		# no index zero).
		y = li
		!only_single_pieces && (marks[y] = true)
		for _ = 2:(ignore_d ? L ÷ li : min(di, L ÷ li))
			marks[y += li] = true
		end
	end

	cuts = Vector{S}()
	sizehint!(cuts, L)
	if only_single_pieces
		# If we just want to check if the piece lengths cannot be obtained from
		# linear combinations of other pieces, we just iterate the unique piece
		# lengths in increasing order, and if they are not marked, then there is
		# no linear combination of other pieces that give the same length.
		usl = unique!(sort(l))
		for y in usl
			!marks[y] && push!(cuts, y)
		end
	else
		for (position, is_marked) in enumerate(marks)
			is_marked && push!(cuts, position)
		end
	end
	cuts
end

"""
    should_extract_piece_from_plate(pii, L, W, sllw, symm = 0x03) :: Bool

!!! **Internal use.**

Considering the enhanced PP-G2KP model, it is necessary to known if there
will be an 'extraction' variable representing the extraction of piece
with index `pii` from a plate with length `L` and width `L`.

A piece should be extracted from a plate if:
1. The piece fits inside the plate.
2. It is not possible to extract the piece in consideration together with any
   other piece (even other copy of the same piece) from the plate in
   consideration (for the cut orientations allowed by `symm` which by default
   are both vertical and horizontal).

# Arguments

1. `pii::D`: The index of the piece.
2. `L::S`: The length of the plate.
3. `W::S`: The width of the plate.
4. `sllw::SortedLinkedLW{D, S}`: The SortedLinkedLW struct for the pieces.
5. `symm::UInt8 = 0x03`: which cut orientations are allowed: `0x01` means
   'allow only horizontal cuts', `0x02` means 'allow only vertical cuts',
   and `0x03` means 'allow both kinds of cuts'.
"""
function should_extract_piece_from_plate(
	pii :: D, L :: S, W :: S, sllw :: SortedLinkedLW{D, S}, symm :: UInt8 = 0x03
) :: Bool where {D, S}
	li = sllw.l[pii]
	wi = sllw.w[pii]
	# If the piece does not fit the plate, then it is not even possible to
	# extract the piece from the plate.
	(li > L || wi > W) && return false

	n = length(sllw.l)

	# If we cut the plate by the length of the piece (a horizontal cut), does
	# this leaves space for another piece? If it does then we postpone to extract
	# the piece from a child of the current plate. Note that we can only do that
	# if the plate is marked to allow horizontal cuts or any cuts.
	# (The 'f' in the variables is from 'frontal'.)
	if symm == 1 || symm == 3
		fl = L - li
		fw = W

		for sli = 1:n # for the pieces in increase-or-stay order of their length
			# if the current piece does not fit the length none after it will fit
			sllw.sl[sli] > fl && break
			# if a piece fits the frontal child, then postpone extraction
			sllw.w[sllw.sli2pii[sli]] <= fw && return false
		end
	end

	if symm == 2 || symm == 3
		# The same as above but for the width of the piece (a vertical cut).
		# (The 'l' in the variables is from 'lateral'.)
		ll = L
		lw = W - wi

		# the same comments of the last loop are valid for this one
		for swi = 1:n
			sllw.sw[swi] > lw && break
			sllw.l[sllw.swi2pii[swi]] <= ll && return false
		end
	end

	# If no piece is small enough to make the plate-parameter be cut in a length
	# (or width) that leaves a second child large enough for extracting the
	# piece-parameter, then we cannot postpone and should allow extracting the
	# piece-parameter from the plate-parameter.
	return true
end

"""
    ub_num_pieces_fit(d, l, w, a, L, W, A, cutoff; ignore_2th_dim = false)

!!! **Internal use.**

Given a piece set (defined by `d`emand, `l`ength, `w`idth, and `a`rea),
and a plate (defined by `L`ength, `W`idth, and `A`rea), returns the `min`
between `cutoff` and an upper bound on the number of pieces that fit inside
the plate. If `ignore_2th_dim` is true, do not use `w` and `W` to filter
pieces that do not even fit the plate when alone.

In other words, `ub_num_pieces_fit(..., 6)` will return `6` if it is
possible (but not guaranteed) for `6` or more pieces to fit in the plate,
and will return `x` if `x < 6` is the greatest number of pieces that have
yet some possibility of being packed together; `x + 1` is already
guaranteed to be impossible.
"""
function ub_num_pieces_fit(
	d :: Vector{D},
	l :: Vector{S},
	w :: Vector{S},
	a :: Vector{P},
	L :: S,
	W :: S,
	A :: P,
	cutoff :: D;
	ignore_2th_dim = false
) :: D where {D, S, P}
	idxs = collect(1:length(a))
	sort!(idxs, by = i -> a[i])
	ub = zero(D)
	for i in idxs
		if l[i] <= L && (ignore_2th_dim || w[i] <= W)
			for _ = 1:d[i]
				A < a[i] && return ub
				A -= a[i]
				ub += 1
				ub > cutoff && return cutoff
			end
		end
	end
	return ub
end

"""
    no_chance_to_fit_6_piece(d, l, w, a, L, W, A; ignore_2th_dim = false)

Returns true if it is demonstrably impossible to fit six pieces into the plate;
and false if there is some possibility (but no guarantee) of fitting six pieces
in the plate.

If `ignore_2th_dim` is true, `w` and `W` are ignored by the algorithm.
"""
function no_chance_to_fit_6_piece(
	d :: Vector{D},
	l :: Vector{S},
	w :: Vector{S},
	a :: Vector{P},
	L :: S,
	W :: S,
	A :: P;
	ignore_2th_dim = false
) :: Bool where {D, S, P}
	ub = ub_num_pieces_fit(
		d, l, w, a, L, W, A, 6; ignore_2th_dim = ignore_2th_dim
	)
	return ub < 6
end

"""
    reduce2fit_usl(sl, sli2pii, w, L, W) :: [S]

!!! **Internal use.**

Given a plate (`L`, `W`), and the pieces sorted by length (`sl`, `sli2pii`,
`w`), return an increasing vector of unique lengths that pertain to a piece
that fits the plate. If two or more pieces share length and fit the plate, only
one copy of that length is included in the returned vector.

# Arguments

1. `sl::Vector{S}`: The piece lengths sorted by increase-or-stay order.
2. `sli2pii::Vector{D}`: If `sli2pii[i] == j` then `sl[i]` and `w[j]`
   correspond to the same item.
3. `w::Vector{S}`: The piece widths (in 'original' order).
4. `L::S`: The length of the plate.
5. `W::S`: The width of the plate.

"""
function reduce2fit_usl(
	sl :: Vector{S},
	sli2pii :: Vector{D},
	w :: Vector{S},
	L :: S,
	W :: S
) :: Vector{S} where {D, S, P}
	xs = empty(sl)
	for (i, x) in enumerate(sl)
		x > L && break # stop before the first lenght greater than L
		!isempty(xs) && last(xs) == x && continue # do not add duplicate lenghts
		w[sli2pii[i]] <= W && push!(xs, x) # add only if fits BOTH dimensions
	end

	xs
end

"""
    fits_at_least_one(sllw, L, W) :: Bool

!!! **Internal use.**

True if at least one piece in `sllw` fits inside a `L`x`W` plate;
false otherwise.
"""
function fits_at_least_one(
	sllw :: SortedLinkedLW{D, S}, L :: S, W :: S
) :: Bool where {D, S}
	for (i, v) in enumerate(sllw.sl)
		v > L && return false
		sllw.w[sllw.sli2pii[i]] <= W && return true
	end

	return false
end

"""
    filter_symm_pos(disc, dim, disc_copy = copy(disc)) -> disc_copy

!!! **Internal use.**

Implement the symmetry-breaking described in Section 2.1 of
10.1287/ijoc.2016.0710 (just after equation 10). Given a list of discretized
lengths (or widths) `disc` and the plate length (or width) `dim`, return a
copy of `disc` without all positions `i` in which ```disc[i] > div(dim/2)
\\land (\\exists. disc[j] == dim - disc[i])```.
"""
function filter_symm_pos(
	disc :: Vector{S},
	dim :: S,
	disc_copy = copy(disc) :: Vector{S}
) :: Vector{S} where {S}
	half_dim = dim ÷ 2
	for xy in disc
		(xy > half_dim || (xy == dim - xy)) && break
		symm_pos = dim - xy
		symm_idx = searchsortedlast(disc, symm_pos)
		!iszero(symm_idx) && disc[symm_idx] == symm_pos && (disc_copy[symm_idx] = 0)
	end
	return filter!(!iszero, disc_copy)
end

# From Furini2016 supplement: "if one (or more) of the flags of plate j has
# value -1, do not perform trim cuts on j in the flag orientation."
function filter_redundant_cuts!(
	# A list of cuts (parent plate, first child, second child). Only the
	# second child may be waste, and if it is, the index of the second child
	# is zero. All cuts share the same orientation (vertical or horizontal).
	nnn :: Vector{NTuple{3, P}},
	# Self Flag: corresponds to 'sh' or 'sv'. The orientation must match 'nnn'.
	sf :: Vector{Int8},
	# Father Flag: corresponds to 'fh' or 'fv'. Analogue to Self Flag.
	ff :: Vector{Int8}
) :: Vector{NTuple{3, P}} where {P}
	return filter!(nnn) do cut
		pp, _, sc = cut # parent plate, first child, second child
		# A cut should be kept if it is not a trim cut (!iszero(sc)) or,
		# in the case it is a trim cut, if none of its flags indicate a trim
		# cut was already made to obtain pp (check sf), or to obtain every possible
		# father plate of pp (check ff).
		!iszero(sc) || (sf[pp] > -1 && ff[pp] > -1)
	end
end

#=
Old documentation of the gen_cuts method. It was superseeded by the
ByproductPPG2KP documentation, but it is more complete than it and so
was kept here for now.

`pli_lwb`: List of all generated plates. The triple is: plate length, plate
	width, and plate bound (how many of the plate fit inside the original plate).
	The index/identifier of a plate is the same as its position in this list.
`hcuts`: If `(pp, cp1, cp2)` is present in this list, then the model has an
   horizontal cut that divides parent plate `pp` into child plates `cp1` and
   `cp2`. The three values are indexes in the list of all generated plates
   (except `cp2` may be zero, i.e., a trim/waste).
`vcuts`: The same as the second returned value, but for the vertical cuts.
`np`: A list of plate-piece pairs in which the plate may be sold as the
   respective piece. If `faithful2furini2016` is `true`, then this list
   has always the same size than the number of piece types and each
   plate that may be sold as a piece has the exact dimensions of the respective
   piece. If `faithful2furini2016` is `false` then this is the list of
   'extraction variables' and the plates may have dimensions slight
   different from their respective pieces, and each plate may be sold as
   one of many distinct pieces, and each piece may be extracted from many
   different plates.
=#
# TODO: Optional. Consider if changing the order of the variables can impact
# the solver performance. The only problem with this is that current
# ByproductPPG2KP forces all horizontal cuts to come before all vertical
# cuts. Changing the order of the constraints, however, needs the complete
# recomputation of ByproductPPG2KP.cuts and ByproductPPG2KP.np.
# Some criteria to be considered for the np variables (mostly):
# how much is the absolute profit obtained; how much relative area is wasted;
# how much is the profit of the a squared unit of the plate used; or group
# them by pii just for making it easier to the demand constraint.
"""
Collection of internal datastructures used to create the PP-G2KP model and
needed to assemble the solution given by the values of the model variables.

$TYPEDFIELDS

"""
struct ByproductPPG2KP{D, S, P}
	# pp (parent plate), fc (first child), sc (second child)
	"If `(pp, fc, sc)` is in `cuts`, then `pp` may be cut into `fc` and `sc`."
	cuts :: Vector{NTuple{3, P}}
	"If hybridize-with-restricted, has the extracted piece for each `cuts` elem."
	cut_extraction :: Vector{D}
	"Every cut before this position is horizontal, the rest are vertical."
	first_vertical_cut_idx :: P
	"If `(n, p)` is in `np` then plate `n` may be sold as piece `p`."
	np :: Vector{Tuple{P, D}} # TODO: rename np to pe (piece extraction)?
	"Indexed by a plate index, values are the plate length, width, and bound."
	pli_lwb :: Vector{Tuple{S, S, P}}
	"The demand of the pieces."
	d :: Vector{D}
	"The length of the pieces."
	l :: Vector{S}
	"The width of the pieces."
	w :: Vector{S}
	#p :: Vector{P} # for now, it is not necessary
	"The length of the original plate."
	L :: S
	"The width of the original plate."
	W :: S
	"If all plates were rotated so length is never greater than width."
	mirror_plates :: Bool
end
Base.broadcastable(x :: ByproductPPG2KP{D, S, P}) where {D, S, P} = Ref(x)

# `_piece_only_positions` is similar to `guarantee_discretization!` but:
# 1. It is only necessary if hybridize-with-restricted is enabled.
# 2. It is eager instead of lazy (i.e., it computes the whole data
#    structure instead of filling the gaps as the enumeration goes).
#    The reasoning is that: (i) different from the discretization,
#    the number of distinct 'piece-only positions' is bounded by the
#    number of pieces (linear); (ii) it is almost certain that we will
#    need all of them.
# The return is a `Vector{Vector{S}}`, lets call it `pols`.
# `pows` outer Vector is indexed by the width available and gives the
# inner vector that is all the lengths that can only be obtained by
# single pieces (not by any combination) considering the length available.
# To obtain its counterpart, `pows`, indexed by the width available and
# returning the lengths only obtainable by single pieces, just switc
# `L` and `W` as well all fields of SortedLinkedLW for the call.
function _piece_only_positions(
	L :: S, W :: S, sllw :: SortedLinkedLW{D, S}, d; ignore_d = false
) :: Vector{Vector{S}} where {D, S}
	@unpack l, w = sllw
	pols = Vector{Vector{S}}(undef, W)
	empty_list = S[]
	# If less width is available than any piece the only valid value
	# is an empty list.
	for x in 1:(first(sllw.sw) - 1)
		pols[x] = empty_list
	end

	for (swi, available_w) in pairs(sllw.sw)
		not_is_last = swi < lastindex(sllw.sw)
		if not_is_last && sllw.sw[swi + 1] == available_w
			continue # i.e., only use the last in a block of repeated values
		end
		pol = discretize(
			d, l, w, L, available_w; only_single_pieces = true,
			ignore_d = ignore_d
		)
		before_next_w = not_is_last ? sllw.sw[swi + 1] - one(S) : W
		for x in available_w:before_next_w
			pols[x] = pol
		end
	end

	return pols
end

# Same as _piece_only_positions, but does the work for both dimensions
# instead of just the length dimension.
# The keyword arguments are passed to `_piece_only_positions` which
# pass them to `discretize`.
function piece_only_positions(
	L :: S, W :: S, sllw :: SortedLinkedLW{D, S}, d; kwargs...
) :: NTuple{2, Vector{Vector{S}}} where {D, S}
	pols = _piece_only_positions(L, W, sllw, d; kwargs...)
	rev_sllw = switched_dims(sllw)
	pows = _piece_only_positions(W, L, rev_sllw, d; kwargs...)

	return pols, pows
end

function guarantee_discretization!(
	dls, dws, pll :: S, plw :: S, L :: S, W :: S, d, l, w,
	round2disc, ignore_2th_dim, ignore_d, level :: UInt8 = UInt8(1)
) :: Tuple{S, S} where {S}
	# If our understanding of the discretization is correct, it is impossible
	# for this method to recurse more than one time in a row.
	@assert level == 1 || level == 2
	(iszero(pll) || iszero(plw)) && return zero(S), zero(S)

	# First of all, let us get the two discretizations, memoized or not.
	dls_plw = if isassigned(dls, plw)
		dls[plw]
	else
		discretize(d, l, w, L, plw; ignore_W = ignore_2th_dim, ignore_d = ignore_d)
	end
	dws_pll = if isassigned(dws, pll)
		dws[pll]
	else
		discretize(d, w, l, W, pll; ignore_W = ignore_2th_dim, ignore_d = ignore_d)
	end

	# The code for fully initializing `dls[plw]` needs partially
	# initializing `dws[L]`, if we also wanted to fully initialize
	# `dws[L]` then we need to save the information that it is not
	# fully initialized (to remember initializing it after).
	dws_L_full_init_pending = pll == L && !isassigned(dws, L)

	# The following two `if`s do essentially the same thing but one is for `dlw`
	# and the other for `dlw`. They find a range of the second (non-discretized)
	# dimension for which the discretization (of the first dimension) is valid.
	# Then assign the same discretization to all positions in the range, possibly
	# saving the effort for new calls (an eager memoization).
	if !isassigned(dls, plw)
		dws_L_full_init_pending && (dws[L] = dws_pll)
		idx_lte_plw_in_dws_L = searchsortedlast(dws[L], plw)
		if iszero(idx_lte_plw_in_dws_L)
			# `searchsortedlast` returns zero if the first value in `dws[L]` is
			# already greater than `plw`, i.e., there is no piece or piece set
			# that fit in so little width (the minimum width, or first
			# discretized position, is after `plw`).
			@assert isempty(dls_plw) # can only be an empty discretization set
			# Assign the empty set to these positions so we to not lose time
			# recomputing them. Those are the first spaces no piece fits.
			!isempty(dws[L]) && (dls[1:(first(dws[L]) - 1)] .= (dls_plw,))
		elseif idx_lte_plw_in_dws_L == length(dws[L])
			# If `idx_lte_plw_in_dws_L` is the last discretization in `dws[L]`
			# this means that from there to `W` all positions can only pack the
			# same pieces/'piece sets'.
			first_valid_width = dws[L][idx_lte_plw_in_dws_L]
			dls[first_valid_width:W] .= (dls_plw,)
		else
			# If `idx_lte_plw_in_dws_L` is nor zero, nor the last discretization,
			# then `dls_plw` is valid for the range `dws[L][idx_lte_plw_in_dws_L]`
			# to `dws[L][idx_lte_plw_in_dws_L + 1] - 1` (i.e., one width unit
			# before the next discretized width position).
			first_valid_width = dws[L][idx_lte_plw_in_dws_L]
			last_valid_width = dws[L][idx_lte_plw_in_dws_L + 1] - 1
			dls[first_valid_width:last_valid_width] .= (dls_plw,)
		end
	end
	if !isassigned(dws, pll) || dws_L_full_init_pending
		@assert plw != W || isassigned(dls, W)
		idx_lte_pll_in_dls_W = searchsortedlast(dls[W], pll)
		# See comments of the similar `if` block above.
		if iszero(idx_lte_pll_in_dls_W)
			@assert isempty(dws_pll)
			!isempty(dls[W]) && (dws[1:(first(dls[W]) - 1)] .= (dws_pll,))
		elseif idx_lte_pll_in_dls_W == length(dls[W])
			first_valid_length = dls[W][idx_lte_pll_in_dls_W]
			dws[first_valid_length:L] .= (dws_pll,)
		else
			first_valid_length = dls[W][idx_lte_pll_in_dls_W]
			last_valid_length = dls[W][idx_lte_pll_in_dls_W + 1] - 1
			dws[first_valid_length:last_valid_length] .= (dws_pll,)
		end
	end

	if round2disc
		r2d_pll_idx = searchsortedlast(dls[plw], pll)
		r2d_plw_idx = searchsortedlast(dws[pll], plw)
		r2d_pll = iszero(r2d_pll_idx) ? zero(S) : dls[plw][r2d_pll_idx]
		r2d_plw = iszero(r2d_plw_idx) ? zero(S) : dws[pll][r2d_plw_idx]
		if iszero(r2d_pll) || iszero(r2d_plw)
			r2d_pll = r2d_plw = zero(S)
		end
		#println("Before recursion")
		#@show plw, pll, r2d_plw, r2d_pll
		r2d_not_guaranteed = (!iszero(r2d_pll) && !isassigned(dws, r2d_pll)) ||
			(!iszero(r2d_plw) && !isassigned(dls, r2d_plw))
		if r2d_not_guaranteed
			r2d_pll, r2d_plw = guarantee_discretization!(
				dls, dws, r2d_pll, r2d_plw, L, W, d, l, w,
				round2disc, ignore_2th_dim, ignore_d, level + one(level)
			)
		end
		@assert iszero(r2d_pll) || isassigned(dws, r2d_pll)
		@assert iszero(r2d_plw) || isassigned(dls, r2d_plw)
		return r2d_pll, r2d_plw
	else
		return pll, plw
	end
end

# mirror plate
_mp(::Val{true}, l, w) = l > w ? (w, l) : (l, w)
_mp(::Val{false}, l, w) = l, w

macro mp(v, l, w)
	v_, l_, w_ = esc(v), esc(l), esc(w)#gensym(v), gensym(l), gensym(w)# v, l, w #
	return :(_mp($v_, $l_, $w_)...)
end

function do_cut!(
	is_vertical :: Bool, # If the cut is vertical (i.e., on a width position)
	n :: P, # The last plate id, incremented with each plate created and returned.
	m :: Val{M}, # If the plates are mirrored (i.e., length always smaller).
	fcl :: S, # The length of the first child (if round2disc then it may shrink).
	fcw :: S, # The width of the first child (if round2disc then it may shrink).
	scl :: S, # The length of the second child (if round2disc then it may shrink).
	scw :: S, # The width of the second child (if round2disc then it may shrink).
	next_idx :: P, # The index of next and sfhv for the current parent plate.
	dls :: Vector{Vector{S}}, # Discretized lengths set for each width.
	dws :: Vector{Vector{S}}, # Discretized widths set for each length.
	pols :: Vector{Vector{S}}, # Piece-Only Length Sets (by available width).
	pows :: Vector{Vector{S}}, # Piece-Only Widths Sets (by available length).
	L :: S, # The original plate length.
	W :: S, # The original plate width.
	d :: Vector{D}, # The pieces demand (needed for guaranteeing discretization).
	l :: Vector{S}, # The pieces length.
	w :: Vector{S}, # The pieces width.
	sllw :: SortedLinkedLW{D, S}, # Used for `fits_at_least_one`
	sfhv :: NTuple{4, Vector{Int8}}, # Used for the Redundant-Cut optimization.
	next :: Vector{Tuple{S, S, P}}, # New plates are to be pushed here.
	plis :: Array{P, 2}, # plis[pll, plw] == pli (if plate already known)
	nnn :: Vector{NTuple{3, P}}, # Either vnnn or hnnn.
	dce :: Vector{D}, # Either vdce or hdce.
	round2disc :: Bool, # Passed to guarantee_discretization!
	ignore_2th_dim :: Bool, # Passed to guarantee_discretization!
	ignore_d :: Bool, # Passed to guarantee_discretization!
	faithful2furini2016 :: Bool,
	hybridize_with_restricted :: Bool,
	aggressive_hybridization :: Bool,
	no_redundant_cut :: Bool,
) :: Tuple{P, P, P, P} where {D, S, P, M}
	aggressive_hybridization_extra_vars = zero(P)
	# fcl: first child length
	# fcw: first child width
	fcl, fcw = guarantee_discretization!(
		dls, dws, fcl, fcw, L, W, d, l, w,
		round2disc, ignore_2th_dim, ignore_d
	)
	# scl: second child length
	# scw: second child width
	scl, scw = guarantee_discretization!(
		dls, dws, scl, scw, L, W, d, l, w,
		round2disc, ignore_2th_dim, ignore_d
	)

	# dcpiis: Double-Cut PIece IndexeS. Has a single zero if no hybridization
	#   has happened, has one or more piece ids if hybridizations have happened.
	# fcdpars: if multiple hybridizations happen, then multiple `fcdpar` exist,
	#   so after the `if hybridize_with_restricted` block below, we loop over
	#   `fcdpars`. If there was no hybridization, then the only `fcdpar` is
	#   pushed into the vector; if there is one or more hybridizations, they
	#   are pushed into the `fcdpar`.
	dcpiis = D[]
	fcdpars = S[]
	# pop -- Piece-Only Position (always perpendicular to cut)
	# fcdper -- First Child Dimension PERpendicular to cut
	# fcdpar -- First Child Dimension PARallel to cut
	# pii2dper -- PIece Index TO Dimension PERpendicular to cut
	# pii2dpar -- PIece Index TO Dimension PARallel to cut
	# sddper -- Sorted Discretized Dimensions PERpendicular to cut
	# sddpar -- Sorted Discretized Dimensions PARallel to cut
	# sdd2pii -- Sorted Discretized Dimensions TO PIece Index (perpendicular)
	# First rename the dimensions to be cut-oriented.
	fcdper, fcdpar = is_vertical ? (fcw, fcl) : (fcl, fcw)
	if hybridize_with_restricted
		# Now rename all the other relevant data structures to be cut-oriented.
		pop, pii2dper, pii2dpar, sddper, sddpar, sdd2pii = if is_vertical
			(ignore_2th_dim ? pows[L] : pows[fcl]), sllw.w, sllw.l,
				sllw.sw, sllw.sl, sllw.swi2pii
		else
			(ignore_2th_dim ? pols[W] : pols[fcw]), sllw.l, sllw.w,
				sllw.sl, sllw.sw, sllw.sli2pii
		end
		# If `fcdper` is in `pop`, then we have one or more hybridizations,
		# otherwise we have no hybridization.
		fcdper_pop_idx = searchsorted(pop, fcdper)
		if isempty(fcdper_pop_idx) # no hybridization
			push!(dcpiis, 0)
			push!(fcdpars, fcdpar)
		else # one or more hybridizations
			spiis = searchsorted(sddper, fcdper)
			# If there is a single hybridization possibility, then we always
			# hybridize. However, if there are multiple possibilities, we only
			# hybridize if aggressive_hybridization is enabled, to avoid
			# increasing the model size too much.
			if isone(length(spiis)) || aggressive_hybridization
				for spii in spiis
					pii = sdd2pii[spii]
					if pii2dpar[pii] <= fcdpar
						push!(dcpiis, convert(D, pii))
						push!(fcdpars, fcdpar - pii2dpar[pii])
					end
				end
				aggressive_hybridization_extra_vars += length(spiis) - 1
			else # multiple hybridization candidates but !aggressive_hybridization
				# Abstain from any hybridization to keep optimality and model size.
				push!(dcpiis, 0)
				push!(fcdpars, fcdpar)
			end
			@assert !isempty(dcpiis)
			@assert !isempty(fcdpars)
		end
	else
		# If !hybridize_with_restricted, we fill dcpiis and fcdpars with
		# sentinel values, that are the correct if no hybridization has ocurred.
		push!(dcpiis, 0)
		push!(fcdpars, fcdpar)
	end # if hybridize_with_restricted

	qt_first_born_kills :: P = zero(P)
	for (dcpii, fcdpar) in zip(dcpiis, fcdpars)
		fcl, fcw = is_vertical ? (fcdpar, fcdper) : (fcdper, fcdpar)

		# If `dcpii` is not zero, then this is an hybridization, and the plate
		# may be further reduced by round2disc.
		if !iszero(dcpii) && round2disc
			fcl, fcw = guarantee_discretization!(
				dls, dws, fcl, fcw, L, W, d, l, w,
				round2disc, ignore_2th_dim, ignore_d
			)
		end

		n, killed_first_born = do_cut!(is_vertical,
			n, m, fcl, fcw, scl, scw, next_idx, dls, dws,
			sllw, sfhv, next, plis, nnn, faithful2furini2016,
			hybridize_with_restricted, no_redundant_cut,
		)
		qt_first_born_kills += killed_first_born

		# Finally, if hybridize_with_restricted is enabled, and the current cut
		# was transformed by this optimization in a double cut, we need to save
		# the extracted piece in dce (otherwise zero is our sentinel).
		hybridize_with_restricted && push!(dce, dcpii)
	end

	qt_hybridizations :: P = zero(P)
	if !iszero(first(dcpiis))
		qt_hybridizations = convert(P, length(dcpiis))
	end

	return (n, qt_hybridizations, qt_first_born_kills,
		aggressive_hybridization_extra_vars)
end

function do_cut!(
	is_vertical :: Bool, # If the cut is vertical (i.e., on a width position)
	n :: P, # The last plate id, incremented with each plate created and returned.
	m :: Val{M}, # If the plates are mirrored (i.e., length always smaller).
	fcl :: S, # The length of the first child (if round2disc then it may shrink).
	fcw :: S, # The width of the first child (if round2disc then it may shrink).
	scl :: S, # The length of the second child (if round2disc then it may shrink).
	scw :: S, # The width of the second child (if round2disc then it may shrink).
	next_idx :: P, # The index of next and sfhv for the current parent plate.
	dls :: Vector{Vector{S}}, # Discretized lengths set for each width.
	dws :: Vector{Vector{S}}, # Discretized widths set for each length.
	sllw :: SortedLinkedLW{D, S}, # Used for `fits_at_least_one`
	sfhv :: NTuple{4, Vector{Int8}}, # Used for the Redundant-Cut optimization.
	next :: Vector{Tuple{S, S, P}}, # New plates are to be pushed here.
	plis :: Array{P, 2}, # plis[pll, plw] == pli (if plate already known)
	nnn :: Vector{NTuple{3, P}}, # Either vnnn or hnnn.
	faithful2furini2016 :: Bool,
	hybridize_with_restricted :: Bool,
	no_redundant_cut :: Bool,
) :: Tuple{P, Bool} where {D, S, P, M}
	# Put some info readily available.
	pll, plw, pli = next[next_idx] # PLate Length, Width, and Index
	if !no_redundant_cut
		sh, sv, fh, fv = sfhv
		sh_j, sv_j, fh_j, fv_j = getindex.(sfhv, next_idx)
	end
	# If the preprocessing is faithful2furini2016, the plates are cut until
	# they have the same size as pieces and, consequently, there exist the
	# concept of trim cut (i.e., a cut in which the second child is waste).
	trim_cut = faithful2furini2016 && !fits_at_least_one(sllw, scl, scw)
	# The trim_cut flag is used below outside of 'if faithful2furini2016'
	# because a true value implicates that 'faithful2furini2016 == true'.

	# For now, the only reason for a first child to be unable to pack
	# any piece is the hybridize_with_restricted optimization.
	killed_first_born = false
	if hybridize_with_restricted && !fits_at_least_one(sllw, fcl, fcw)
		killed_first_born = true
	elseif iszero(plis[@mp m fcl fcw]) # If the first child does not yet exist.
		push!(next, (fcl, fcw, n += 1)) # Create the first child.
		if !no_redundant_cut
			# For now, we do not know if hybridize_with_restricted is
			# entirely compatible with Redundant-Cut, so we only consider
			# a trim_cut if it does not clash with hybridize_with_restricted.
			if trim_cut && !hybridize_with_restricted
				# From Furini2016 supplement: "if a NEW plate j_1 ∈ J is obtained
				# from j through a trim cut with orientation v:"
				if is_vertical
					fchild_sfhv = (0, -1, sh_j, sv_j)
				else
					fchild_sfhv = (-1, 0, sh_j, sv_j)
				end
			else
				# If a plate (NEW or existing) j_1 is obtained from j without a trim
				# cut: set all flags of j 1 to 1.
				fchild_sfhv = (1, 1, 1, 1)
			end
			push!.(sfhv, fchild_sfhv)
		end
		plis[@mp m fcl fcw] = n # Mark plate existence.
	elseif !no_redundant_cut
		# If plate already exists, and we are implementing Redundant-Cut.
		# Here we do not care for hybridize_with_restricted because
		# setting positions to positive (i.e., `1`) is disabling
		# Redundant-Cut optimizations not enabling them.
		if trim_cut
			# From Furini2016 supplement: "if an existing plate j_1 ∈ J is
			# obtained from j through a trim cut:"
			sh_j > -1 && (fh[plis[@mp m fcl fcw]] = 1)
			sv_j > -1 && (fv[plis[@mp m fcl fcw]] = 1)
		else
			# From Furini2016 supplement: "if a plate (new or EXISTING) j_1 is
			# obtained from j without a trim cut: set all flags of j_1 to 1."
			setindex!.(sfhv, (1, 1, 1, 1), plis[@mp m fcl fcw])
		end
	end

	# If the second child is not waste (!trim_cut), nor was already
	# generated, then it is a new plate obtained without a trim cut.
	if !trim_cut && iszero(plis[@mp m scl scw])
		# Assert meaning: if the second child is not waste, then it must
		# have an associated size.
		@assert !iszero(scw)
		# Save the plate to be processed later, and mark its existence.
		push!(next, (scl, scw, n += 1))
		plis[@mp m scl scw] = n
		# From Furini2016 supplement: "if a plate (NEW or existing) j 1 is
		# obtained from j without a trim cut: set all flags of j_1 to 1."
		!no_redundant_cut && push!.(sfhv, (1, 1, 1, 1))
	end
	# Add the cut to the cut list. Check if the second child plate is waste
	# before trying to get its plate index.
	push!(nnn, (pli,
		killed_first_born ? 0 : plis[@mp m fcl fcw],
		trim_cut ? 0 : plis[@mp m scl scw]
	))

	return n, killed_first_born
end

"""
    reflect!(l, L)

!!! **Internal use.**

Change `l` to the sorted union of `[li ∈ l ∧ li ≤ L÷2 | li]` and `[li ∈ l ∧ li
> L÷2 ∧ li < L | L - li]`; assumes `l` is in increase-or-stay order.

The new value of `l` has an important property that: any pair of plates
generated by a cut using the original `l` will also be generated by some cut in
the modified set, even with the modified set only having cuts in the first half
of the plate. Note any values equal to or above L are discarded.
"""
function reflect!(l, L)
	mid_L = L ÷ 2
	qt_elems = 0
	for (i, v) in pairs(l)
		v >= L && break
		qt_elems += 1
		v > mid_L && (l[i] = L - v)
	end
	resize!(l, qt_elems)
	return unique!(sort!(l))
end

# !!! **Internal use.**
#
#     _pop_smaller_than!(l, min_l)
#
# Assumes a increase-or-stay ordered `Vector` `l` and delete values smaller
# than `min_l`. Returns `l`.
function _pop_smaller_than!(l, min_l)
	# This may seem terrible performance-wise, but the way Julia stores
	# `Vectors` allows `popfirst!` to be on par with `pop!`.
	while !isempty(l) && first(l) < min_l
		popfirst!(l)
	end
	return l
end

"""
    gen_cuts(::Type{P}, d, sllw, L, W; [kwargs])

!!! **Internal use.**

The main method responsible for generating all the cuts used to create both the
original PP-G2KP model as its enhanced version (in this case, it generates the
piece extractions too).

# Arguments

## Positional and required arguments

1. `::Type{P}`: The type of the piece profits (and area).
2. `d::Vector{D}`: The pieces demand.
3. `sllw::SortedLinkedLW{D, S}`: The SortedLinkedLW containing the length and
   width of the pieces.
4. `L::S`: The plate length.
5. `W::S`: The plate width.

## Keyword arguments

All keyword arguments are of type `Bool`.

* `ignore_2th_dim`: Default: false. Ignore the dimension not being discretized
  during discretization.
* `ignore_d`: Default: false. Ignore the demand information during
  discretization.
* `round2disc`: Default: true. Round the size of the second child of a
  cut to a discretized position.
* `faithful2furini2016`: Default: false. Tries to be the most faithful possible
  to the description in 10.1287/ijoc.2016.0710.
* `no_redundant_cut`: Default: false. Disables the Redundant-Cut reduction
  described in 10.1287/ijoc.2016.0710.
* `no_cut_position`: Default: false. Disables the Cut-Position reduction
  described in 10.1287/ijoc.2016.0710.
* `no_furini_symmbreak`: Default: false. Ignored if `faithful2furini2016`
  is `false`. Disables the symmetry-breaking used in 10.1287/ijoc.2016.0710
  (consequently all discretized positions are used, none is removed).
"""
function gen_cuts(
	::Type{P}, d :: Vector{D}, sllw :: SortedLinkedLW{D, S}, L :: S, W :: S,
	m :: Val{M} = Val(false);
	ignore_2th_dim = false, ignore_d = false, round2disc = true,
	hybridize_with_restricted = false,
	aggressive_hybridization = false,
	faithful2furini2016 = false,
	no_redundant_cut = false, no_cut_position = false,
	no_furini_symmbreak = false, quiet = false, verbose = false
) :: ByproductPPG2KP{D, S, P} where {D, S, P, M}
	mirror_plates = M # the Bool inside the Val
	# If both verbose and quiet are passed, then quiet wins.
	verbose = verbose & !quiet
	!quiet && faithful2furini2016 && round2disc && @warn(
		"Enabling both faithful2furini2016 and round2disc is allowed, but you" *
		" are not being entirely faithful to Furini2016 if you do so."
	)
	!faithful2furini2016 && no_redundant_cut && !quiet && @warn(
		"The Redundant-Cut is only used when faithful2furini2016" *
		" is enabled. As flag faithful2furini2016 is not enabled, flag" *
		" no-redundant-cut has no effect."
	)
	if !no_redundant_cut && mirror_plates
		!quiet && @warn(
			"Flag mirror-plates has disabled the Redundant-Cut reduction."
		)
		no_redundant_cut = true
	end

	l = sllw.l
	w = sllw.w
	@assert length(d) == length(l)
	@assert length(d) == length(w)
	# unique versions
	usl = unique(sllw.sl)
	usw = unique(sllw.sw)
	a = l .* w
	A = L * W
	max_piece_type = convert(D, length(l))
	#max_num_plates = convert(P, L) * convert(P, W)
	min_pil = minimum(l)
	min_piw = minimum(w)

	# If (n1, n2, n3) is in nnn, then n1 may be partitioned in n2 and n3.
	# Not every possible partition is present. If the model is trying to be
	# the most faithful to Furini, then redundant-cuts is implemented. If the
	# model is not restricted to Furini's description, then it may cut even more.
	hnnn = Vector{NTuple{3, P}}()
	vnnn = Vector{NTuple{3, P}}()
	# dce (double cut extractions): if !hybridize_with_restricted these are kept
	# as empty Vectors; otherwise, some cuts are double cuts, that generate
	# up to two child plates as usual, but also a piece extraction. If vdce[i]
	# is zero, then hnnn[i] is not a double cut, otherwise hnnn[i] is a double
	# cut and hdci[i] indicate the piece extracted. Analogue for vnnn and vdce.
	hdce = Vector{D}()
	vdce = Vector{D}()
	# The np vector has the pairs plate-piece for which the respective plate is
	# allowed to be sold as the respective piece.
	# If the model is faithful to Furini, it just has a single plate for each
	# piece type and the plate has the exact size of the piece.
	# Otherwise, if the piece fits the plate size and the plate does not allow
	# any other piece by its side (i.e., there is no second piece that can be
	# placed in the same plate vertically or horizontally) then (plate, piece) is
	# in np.
	np = Vector{Tuple{P, D}}()
	# The list of plates attributes: plate length, plate width, and plate bound.
	# The plate index is the same as the index in pli_lwb.
	pli_lwb = Vector{Tuple{S, S, P}}()
	# plis: matrix of the plate dimensions in which zero means "never seen that
	# plate before" and nonzero means "this nonzero number is the plate index".
	if mirror_plates
		plis = zeros(P, max(L, W), max(L, W))
	else
		plis = zeros(P, L, W)
	end
	plis[@mp m L W] = one(P)
	# next: plates already indexed but not yet processed, starts with (L, W, 1).
	# Storing the index as the third value is not necessary (as it could be
	# queried from plis) but this is probably more efficient this way.
	next = Vector{Tuple{S, S, P}}()
	next_idx = one(P)
	#sizehint!(next, max_num_plates)
	push!(next, (L, W, one(P)))
	sfhv = (Vector{Int8}(), Vector{Int8}(), Vector{Int8}(), Vector{Int8}())
	if !no_redundant_cut
		# If the preprocessing is faithful2furini2016 and Redundant-Cut is
		# enabled, then four auxiliary trim cut flag vectors are necessary.
		sh, sv, fh, fv = sfhv
		push!.(sfhv, (1, 1, 1, 1))
	end
	# n: The amount of plates (the index of the highest plate type).
	n = one(P) # there is already the original plate
	# Memoized discretizations. The discretized lengths for every plate width.
	# The discretized widths for every plate length.
	dls = Vector{Vector{S}}(undef, W)
	dws = Vector{Vector{S}}(undef, L)
	# piece only length/width -- for use with hybridize_with_restricted,
	# has the length and widths that only appear on pieces and not on
	# piece combinations. We start simple, without the optimization below.
	# TODO: `pol` and `pow` could be a set of discretizations like `dls`/`dws`
	# because in smaller plates some piece length/width may stop being
	# attainable by combinations of smaller plates (e.g., a set of pieces
	# matched another piece length, but in the smaller plate some piece of
	# the set does not fit because of the width, while the piece that had
	# the length matched yet fits at that width).
	if hybridize_with_restricted
		pols, pows = piece_only_positions(L, W, sllw, d; ignore_d = ignore_d)
	else
		# If `hybridize_with_restricted` these will never be used, but
		# it is practical if they exist and are the right type to be passed
		# as arguments (and type-stable).
		pols = pows = Vector{S}[]
	end

	# There are two discretizations for each dimension: dl1/dw1 and dl2/dw2.
	# They are only different if the Cut-Position reduction is being used. In
	# this case, dl1/dw1, which is used for cutting the plate, may be a
	# restricted cut set (i.e., it does not have piece size combinations, just
	# single piece sizes). The dl2/dw2 is always the complete discretization,
	# and may be needed (even if Cut-Position is used) for reducing the size of
	# the second child to the last discretized position (if round2disc is
	# enabled).
	guarantee_discretization!(
		dls, dws, L, W, L, W, d, l, w, round2disc, ignore_2th_dim, ignore_d
	)
	hybridizations = 0
	first_borns_killed_by_hybridization = 0
	aggressive_hybridization_extra_vars = 0
	while next_idx <= length(next)
		#@show next_idx
		#@show length(next)
		@assert (no_redundant_cut || [length(next)] == unique(length.(sfhv)))

		pll, plw, pli = next[next_idx] # PLate Length, Width, and Index
		if !faithful2furini2016
			for pii in 1:max_piece_type # pii: PIece Index
				if should_extract_piece_from_plate(pii, pll, plw, sllw)
					push!(np, (pli, pii))
				end
			end
		end

		plb = (L ÷ pll) * (W ÷ plw) # PLate Bound
		# It is not necessary to store the plate id in pli_lwb because they are added
		# in order (a queue is used), so the array index is the plate index.
		push!(pli_lwb, (pll, plw, plb))
		if ignore_2th_dim
			dl1 = dl2 = dls[W]
			dw1 = dw2 = dws[L]
		else
			#@show plw
			@assert isassigned(dls, plw)
			#@show pll
			@assert isassigned(dws, pll)
			#=
			pll, plw = guarantee_discretization!(
				dls, dws, pll, plw, L, W, d, l, w, round2disc, ignore_2th_dim, ignore_d
			)
			=#
			dl1 = dl2 = dls[plw]
			dw1 = dw2 = dws[pll]
		end

		if !no_cut_position
			# In the code below, we had a long-standing bug in the interaction
			# between the Cut-Position dominance and the enhanced model. The bug only
			# occured in very specific circunstances, i.e., ALL optimal solutions
			# must depend on (i) the existence of an intermediary plate that does not
			# fit 6 pieces and must also depend on (ii) the existence of a cut over
			# such plate that is not the size of any piece smaller than half the
			# plate.  For example, in the specific case we found this bug we had a
			# 2x12 plate that fitted at most 5 pieces (and, in fact, in the studied
			# optimal packing it had 5 pieces inside), and the first cut to be done
			# had to divide it from 2x12 into 2x4 and 2x8. The Cut-Position changed
			# dw1 to [1, 2, 3, 8, 12], the original discretization had a cut at 4
			# because 1 + 3, but as Cut-Position only allows restricted cuts and
			# there is no piece of width 4. In the faithful/original model this is
			# not a problem, because every position in the array is considered, but
			# the same is not true for enhanced: only positions before the middle of
			# the plate are cut. Disabling Cut-Position for the enhanced formulation
			# is an obvious fix, but there are many cases in which the Cut-Position
			# and the enhanced formulation do not conflict, and it is in our interest
			# to reduce every cut possible while guaraneeing optimality. The solution
			# found was to "reflect" the restricted cut set to the first half of the
			# plate; so, if some piece larger-than-half-plate needs to be cut from a
			# do-not-fit-6-piece plate, then there will be a cut in the first half
			# that leaves a second child of the exact length (or width) of the piece
			# we need to obtain. It is also possible to compare the enhanced cut set
			# without Cut-Position with this restricted-"reflected" one, and chose
			# the best (i.e., the smallest). The advantage of the
			# restricted-"reflected" approach is to keep the assumption the cuts will
			# only be done at the first half of the plate (from the enhanced
			# formulation), while univocally having the same benefit than the one
			# given by the original restricted approach (in fact, given the removal
			# of perfect symmetries used by the faithful formulation, things like
			# the vertical cuts 4 and 8 in a plate of width 12 were already being
			# conflated into a single cut).
			if no_chance_to_fit_6_piece(
				d, l, w, a, pll, plw, pll * plw; ignore_2th_dim = ignore_2th_dim
			)
				if ignore_2th_dim
					# This `deepcopy` is a ridiculous solution, but `ignore_2th_dim' is
					# never used seriously and I refuse to create yet another variable
					# just for it.
					dl1 = _pop_smaller_than!(reflect!(deepcopy(usl), pll), min_pil)
					dw1 = _pop_smaller_than!(reflect!(deepcopy(usw), plw), min_piw)
				else
					dl1 = reduce2fit_usl(sllw.sl, sllw.sli2pii, w, pll, plw)
					@assert dl1 ⊆ dl2
					dw1 = reduce2fit_usl(sllw.sw, sllw.swi2pii, l, plw, pll)
					@assert dw1 ⊆ dw2
					# Here is the workaround we talk about above.
					# NOTE: if dl1/dw1 stops being generated for each new plate
					# we will need to `deepcopy` the variables.
					# NOTE 2: the `_pop_smaller_than` is necessary for not breaking
					# our assumption of "no cut position exists before the
					# length (width) of the piece with smallest length (width) that
					# fits the current plate". If this assumption is broken the
					# round2disc/plate-normalization breaks because it normalizes
					# child plates from such cuts to 0x0 plates (rightfully so,
					# a 0x0 plate is the smallest plate that fits an empty set of
					# pieces). These cuts are never needed in the enhanced (even
					# if we are only allowed to use restricted cuts) because
					# in such cases there is an extraction variable to extract the
					# larger-than-half-plate piece that, when reflected, generated
					# such these thin cuts (that broken our assumptions).
					if !faithful2furini2016
						if !isempty(dl1)
							plate_min_pil = first(dl1)
							dl1 = _pop_smaller_than!(reflect!(dl1, pll), plate_min_pil)
						end
						if !isempty(dw1)
							plate_min_piw = first(dw1)
							dw1 = _pop_smaller_than!(reflect!(dw1, plw), plate_min_piw)
						end
					end
				end
				#= Debug. Showing something is not an error.
				if dl1 != dl2[1:length(dl1)] || dw1 != dw2[1:length(dw1)]
					fpr_count += 1
					@show pll
					@show plw
					if dl1 != dl2[1:length(dl1)]
						@show dl1
						@show dl2
					end
					if dw1 != dw2[1:length(dw1)]
						@show dw1
						@show dw2
					end
				end
				=#
			end
		end

		# If faithful2furini2016 flag is enabled, then the whole discretization is
		# iterated, except by the cuts on the second half of the plate that have an
		# exact symmetry on the fist half (same distance from plate midpoint). The
		# calls below create a new vector without such symmetric points.
		if faithful2furini2016 && !no_furini_symmbreak
			dl1 = filter_symm_pos(dl1, pll)
			dw1 = filter_symm_pos(dw1, plw)
		end

		for x in dw1
			# See comment above filter_symm_pos calls.
			if faithful2furini2016
				x >= plw && break
			else
				x > plw ÷ 2 && break
			end

			@assert pll >= min_pil
			@assert x >= min_piw
			@assert faithful2furini2016 || plw - x >= min_piw

			fcl, fcw, scl, scw = pll, x, pll, plw - x
			n, qt_hybridizations, qt_first_born_kills, agg_hyb_extra = do_cut!(true,
				n, m, fcl, fcw, scl, scw, next_idx, dls, dws, pols, pows,
				L, W, d, l, w, sllw, sfhv, next, plis,
				vnnn, # Either vnnn or hnnn, i.e., the correct orientation of the two
				vdce, # Either vdce or hdce, i.e., the correct orientation of the two
				round2disc, ignore_2th_dim, ignore_d,
				faithful2furini2016, hybridize_with_restricted,
				aggressive_hybridization, no_redundant_cut
			)
			hybridizations += qt_hybridizations
			first_borns_killed_by_hybridization += qt_first_born_kills
			aggressive_hybridization_extra_vars += agg_hyb_extra
		end

		for y in dl1
			# See comment above filter_symm_pos calls.
			if faithful2furini2016
				y >= pll && break
			else
				y > pll ÷ 2 && break
			end

			@assert plw >= min_piw
			@assert y >= min_pil
			@assert faithful2furini2016 || pll - y >= min_pil

			fcl, fcw, scl, scw = y, plw, pll -y, plw
			n, qt_hybridizations, qt_first_born_kills, agg_hyb_extra = do_cut!(false,
				n, m, fcl, fcw, scl, scw, next_idx, dls, dws, pols, pows,
				L, W, d, l, w, sllw, sfhv, next, plis,
				hnnn, # Either vnnn or hnnn, i.e., the correct orientation of the two
				hdce, # Either vdce or hdce, i.e., the correct orientation of the two
				round2disc, ignore_2th_dim, ignore_d,
				faithful2furini2016, hybridize_with_restricted,
				aggressive_hybridization, no_redundant_cut
			)

			hybridizations += qt_hybridizations
			first_borns_killed_by_hybridization += qt_first_born_kills
			aggressive_hybridization_extra_vars += agg_hyb_extra
		end

		# Goes to the next unprocessed plate.
		next_idx += 1
	end

	if faithful2furini2016
		# Apply Redundant-Cut.
		if !no_redundant_cut
			# red_cut_flags = collect(zip(sh, sv, fh, fv))
			# println("red_cut_flags begin")
			# show(stdout, "text/plain", red_cut_flags)
			# println()
			# println("red_cut_flags end")
			qt_cuts_before_redundant_cut = length(hnnn) + length(vnnn)
			filter_redundant_cuts!(hnnn, sh, fh)
			filter_redundant_cuts!(vnnn, sv, fv)
			qt_cuts_after_redundant_cut = length(hnnn) + length(vnnn)
			qt_cuts_removed_by_redundant_cut = qt_cuts_before_redundant_cut -
				qt_cuts_after_redundant_cut
			if verbose
				print(
					"""
					qt_cuts_before_redundant_cut = $qt_cuts_before_redundant_cut
					qt_cuts_after_redundant_cut = $qt_cuts_after_redundant_cut
					qt_cuts_removed_by_redundant_cut = $qt_cuts_removed_by_redundant_cut
					"""
				)
			end
		end

		# Create conversion table of the index of the plate with same size as
		# a piece to the corresponding piece index.
		for pii in 1:max_piece_type # pii: PIece Index
			# Assert meaning: for every distinct piece dimensions there should be
			# a plate with the exact same dimensions.
			#@assert !iszero(get(plis, (l[pii], w[pii]), 0))
			#push!(np, (plis[(l[pii], w[pii])], pii))
			@assert !iszero(plis[@mp m l[pii] w[pii]])
			push!(np, (plis[@mp m l[pii] w[pii]], pii))
		end
	end

	# The two commented 'for' loops below print all variable indexes. They were
	# used to check the variable generated against the Martin's code. They were
	# kept commented given the necessity arises again, but may be removed in the
	# future (or a flag may be created to output them).
	#=
	for hcut in hnnn
		pp, fc, _ = hcut
		ppl, ppw, _ = pli_lwb[pp]
		q = pli_lwb[fc][1]
		println("$(ppl) $(ppw) 1 $(q)")
	end
	for vcut in vnnn
		pp, fc, _ = vcut
		ppl, ppw, _ = pli_lwb[pp]
		q = pli_lwb[fc][2]
		println("$(ppl) $(ppw) 0 $(q)")
	end
	=#

	if verbose
		@show hybridizations
		@show first_borns_killed_by_hybridization
		@show aggressive_hybridization_extra_vars
	end

	first_vertical_cut_idx = length(hnnn) + 1
	return ByproductPPG2KP(
		append!(hnnn, vnnn), append!(hdce, vdce),
		first_vertical_cut_idx,
		np, pli_lwb, deepcopy(d), deepcopy(l), deepcopy(w),
		deepcopy(L), deepcopy(W), mirror_plates
	)
end

end # module


text:
# Autogenerated wrapper script for Kokkos_jll for x86_64-linux-gnu-cxx11
export libkokkoscontainers, libkokkoscore

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Kokkos")
JLLWrappers.@declare_library_product(libkokkoscontainers, "libkokkoscontainers.so.3.5")
JLLWrappers.@declare_library_product(libkokkoscore, "libkokkoscore.so.3.5")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libkokkoscontainers,
        "lib/libkokkoscontainers.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libkokkoscore,
        "lib/libkokkoscore.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MultiAssayExperiments
using Documenter
import SummarizedExperiments
import DataStructures
import DataFrames

makedocs(
    sitename="MultiAssayExperiments.jl",
    modules = [MultiAssayExperiments],
    pages=[
        "Home" => "index.md"
    ]
)

deploydocs(;
    repo="github.com/LTLA/MultiAssayExperiments.jl",
)


text:
# useful guide here http://math.mit.edu/~stevenj/18.303/fall13/lecture-10.html

diff1(n, d=1.) = inv(d) .* [ [1.0 zeros(1,n-1)]; diagm(1=>ones(n-1)) - diagm(0=>ones(n)) ]
project_on_diff1(n) = [[1 zeros(1,n-1)]; diagm(1=>0.5.*ones(n-1)) + diagm(0=>vcat(0.5.*ones(n-1), [1.0]))]
project_on_Dx(nz,nx) = kron(project_on_diff1(nx), diagm(0=>ones(nz)))
project_on_Dz(nz,nx) = kron(diagm(0=>ones(nx)), project_on_diff1(nz))
diff1x(nz,nx,dz=1.,dx=1.) = kron(diff1(nx,dx), diagm(0=>ones(nz)))
diff1z(nz,nx,dz=1.,dx=1.) = kron(diagm(0=>ones(nx)), diff1(nz,dz))

#=
x=repeat([1,2,3], 1,4)
xx=project_on_Dz(3,4)*vec(x)

function Laplacian(Nx, Ny)

	dy=inv(Ny-1);
	dx=inv(Nx-1);
   Dx = diff1(Nx) / dx
   Dy = diff1(Ny) / dy
   Ax = Dx' * Dx
   Ay = Dy' * Dy
   return kron(diagm(0=>ones(Ny)), Ax) + kron(Ay, diagm(0=>ones(Nx)))
end

n=5
A=Array(Param(n,n,ones(n*n)).A)[1:n*n,1:n*n]
B=Laplacian(n,n)

x=zeros(n,n)
x[div(n,2)+1,div(n,2)+1]=1.0
x=vec(x)

xx1=reshape(A*x,n,n)
xx2=reshape(B*x,n,n)
=#

"""
* `mpars` : medium parameters for forward or inverse mode
"""
function updateA!(pa::Param, mpars; A=pa.A)
	println("ffffwewe")
	nz=pa.nz; nx=pa.nx; dz=pa.dz; dx=pa.dx
	nznx=nz*nx
	dx=pa.dx
	dz=pa.dz

	Ax=diff1x(nz,nx,dz,dx)' * diagm(0=>project_on_Dx(nz,nx)*vec(mpars))  * diff1x(nz,nx,dz,dx)
	Az=diff1z(nz,nx,dz,dx)' * diagm(0=>project_on_Dz(nz,nx)*vec(mpars))  * diff1z(nz,nx,dz,dx)

	copyto!(A,Ax.+Az)

	# update At
	for i in 1:nznx+1
		for j in 1:nznx
			pa.At[j,i]=A[i,j]
		end
	end
	return pa
end

"""
transpose of dAdx
"""
function dAdmpars(dx,dz,nx,nz,T)

	ddz = inv(dz*dz); # for simplicity below
	ddx = inv(dx*dx);
	X=[spzeros(T,(nz*nx)+1,nz*nx) for i in 1:nz*nx]  # note the size, it has appended zeros

	m1=randn(nz*nx)
	m2=zeros(nz*nx)

	for i in 1:nz*nx
		copyto!(m2,m1)
		m2[i]+=0.1
		Ax1=diff1x(nz,nx,dz,dx)' * diagm(0=>project_on_Dx(nz,nx)*m1)  * diff1x(nz,nx,dz,dx)
		Az1=diff1z(nz,nx,dz,dx)' * diagm(0=>project_on_Dz(nz,nx)*m1)  * diff1z(nz,nx,dz,dx)

		Ax2=diff1x(nz,nx,dz,dx)' * diagm(0=>project_on_Dx(nz,nx)*m2)  * diff1x(nz,nx,dz,dx)
		Az2=diff1z(nz,nx,dz,dx)' * diagm(0=>project_on_Dz(nz,nx)*m2)  * diff1z(nz,nx,dz,dx)
		copyto!(X[i],Ax1+Az1-Ax2-Az2)
		rmul!(X[i],10)
	end

	 return X
end


text:
using Pkg

Pkg.activate(".")

using Documenter, SCARGC
using DocumenterTools: Themes

Themes.compile(joinpath(@__DIR__,"src/assets/scargc-light.css"), joinpath(@__DIR__,"src/assets/themes/documenter-light.css"))

makedocs(
        format = Documenter.HTML(
            assets=[asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css)]
        ),
        sitename="SCARGC.jl",
        authors = "Gabriel Marinho",
        pages = [
            "Home" => "index.md",
            "Functions" => "functions.md"
        ]
    )

deploydocs(
    repo = "github.com/MarinhoGabriel/SCARGC.jl.git",
)


text:
struct ExposureResult
    val::Dict{RiskFactor, Float64}

    function ExposureResult(dict::Dict{RiskFactor, Float64})
        # creates a copy insted of reusing the instance
        dict_copy = Dict{RiskFactor, Float64}()

        for (k, v) in dict
            dict_copy[k] = v
        end

        return new(dict_copy)
    end
end

function ExposureResult()
    return ExposureResult(Dict{RiskFactor, Float64}())
end

Base.:(==)(ex1::ExposureResult, ex2::ExposureResult) = ex1.val == ex2.val
Base.hash(ex::ExposureResult) = hash(length(ex.val))

#
# ExposureResult API
#

function Base.getindex(r::ExposureResult, k::RiskFactor) :: Float64
    return r.val[k]
end

function Base.haskey(r::ExposureResult, k::RiskFactor) :: Bool
    return haskey(r.val, k)
end

function Base.setindex!(r::ExposureResult, v::Float64, k::RiskFactor)
    r.val[k] = v
    nothing
end

function Base.setindex!(r::ExposureResult, v::Integer, k::RiskFactor)
    r.val[k] = Float64(v)
    nothing
end

Base.length(r::ExposureResult) = length(r.val)
Base.isempty(r::ExposureResult) = isempty(r.val)
Base.keys(r::ExposureResult) = keys(r.val)
Base.values(r::ExposureResult) = values(r.val)
Base.iterate(r::ExposureResult) = iterate(r.val)
Base.iterate(r::ExposureResult, state) = iterate(r.val, state)
Base.delete!(r::ExposureResult, k::RiskFactor) = delete!(r.val, k)
Base.get(r::ExposureResult, k::RiskFactor, default::Float64) = get(r.val, k, default)

function Base.:*(s::Real, expo::ExposureResult) :: ExposureResult
    result = ExposureResult()

    if s ≈ 0.0
        return result
    end

    for (k, v) in expo
        result[k] = s * v
    end

    return result
end

function Base.:*(expo::ExposureResult, s::Real) :: ExposureResult
    return s * expo
end

function Base.:+(ex1::ExposureResult, ex2::ExposureResult) :: ExposureResult
    result = ExposureResult()

    for (k, v) in ex1
        # copies ex1 to result
        result[k] = v
    end

    for (k, v) in ex2
        # increment result by ex2
        result[k] = v + get(result, k, 0.0)
    end

    _clean_zero_exposures!(result)
    return result
end

function _clean_zero_exposures!(expo::ExposureResult)
    for (k, v) in expo
        if v ≈ 0.0
            delete!(expo, k)
        end
    end

    nothing
end