text:
# legacy content to facilitate transition to AMP

function buildHybridManifoldCallbacks(manif::T) where {T <: Tuple}
  # TODO use multiple dispatch instead -- will be done for second version of system
  addopT = []
  diffopT = []
  getManiMu = []
  getManiLam = []

  for mn in manif
    if mn == :Euclid
      push!(addopT, +)
      push!(diffopT, -)
      push!(getManiMu, KDE.getEuclidMu)
      push!(getManiLam, KDE.getEuclidLambda)
    elseif mn == :Circular
      push!(addopT, addtheta)
      push!(diffopT, difftheta)
      push!(getManiMu, getCircMu)
      push!(getManiLam, getCircLambda)
    else
      error("Unrecognized manifold $(mn)")
    end
  end

  return (addopT...,), (diffopT...,), (getManiMu...,), (getManiLam...,)
end

# """
#     $(SIGNATURES)
#
# Calculate the KDE bandwidths for each dimension independly, as per manifold of each.  Return vector of all dimension bandwidths.
# """
function getKDEManifoldBandwidths(pts::AA,
                                  manif::T1 ) where {AA <: AbstractArray{Float64,2}, T1 <: Tuple}
  #
  ndims = size(pts, 1)
  bws = ones(ndims)

  for i in 1:ndims
    if manif[i] == :Euclid
      bws[i] = getBW( kde!(pts[i,:]) )[1,1]
    elseif manif[i] == :Circular
      bws[i] = getBW( kde!_CircularNaiveCV( pts[i,:] ) )[1,1]
    else
      error("Unrecognized manifold $(manif[i])")
    end
  end

  return bws
end

function ensurePeriodicDomains!( pts::AA, manif::T1 )::Nothing where {AA <: AbstractArray{Float64,2}, T1 <: Tuple}

  i = 0
  for mn in manif
    i += 1
    if manif[i] == :Circular
      pts[i,:] = TUs.wrapRad.(pts[i,:])
    end
  end

  nothing
end

"""
    $(SIGNATURES)

Legacy extension of KDE.kde! function to approximate smooth functions based on samples, using likelihood cross validation for bandwidth selection.  This method allows approximation over hybrid manifolds.
"""
function manikde!(pts::AA2,
                  bws::Vector{Float64},
                  manifolds::T  )::BallTreeDensity where {AA2 <: AbstractArray{Float64,2}, T <: Tuple}
  #
  addopT, diffopT, getManiMu, getManiLam = buildHybridManifoldCallbacks(manifolds)
  KernelDensityEstimate.kde!(pts, bws, addopT, diffopT)
end
function manikde!(pts::AA2,
                  manifolds::T  )::BallTreeDensity where {AA2 <: AbstractArray{Float64,2}, T <: Tuple}
  #
  bws = getKDEManifoldBandwidths(pts, manifolds)
  ensurePeriodicDomains!(pts, manifolds)
  ApproxManifoldProducts.manikde!(pts, bws, manifolds)
end


text:
module MarkovModelAndMMs

include(joinpath(@__DIR__, "MarkovStuff.jl"))
using .MarkovStuff

include(joinpath(@__DIR__, "MMs.jl"))
using .MMs

include(joinpath(@__DIR__, "CostFunctions.jl"))
using .CostFunctions

export MarkovStuff, MMs, CostFunctions

end

ext: jl


text:
module Arrhenius
    using YAML
    using NPZ
    using LinearAlgebra
    using SparseArrays

    include("Constants.jl")
    include("DataStructure.jl")
    include("Solution.jl")
    include("Magic.jl")
    include("Thermo.jl")
    include("Kinetics.jl")
    include("Transport.jl")
end


text:
# sink support
prep!(T, A) = A, 0
prep!(::Type{Union{T, Missing}}, A) where {T} = T[ifelse(ismissing(x), zero(T), x) for x in A]
prep!(::Type{Union{Dates.Date, Missing}}, A) = ODBC.API.SQLDate[ismissing(x) ? ODBC.API.SQLDate() : ODBC.API.SQLDate(x) for x in A], 0
prep!(::Type{Union{Dates.DateTime, Missing}}, A) = ODBC.API.SQLTimestamp[ismissing(x) ? ODBC.API.SQLTimestamp() : ODBC.API.SQLTimestamp(x) for x in A], 0
prep!(::Type{Union{Dec64, Missing}}, A) = Float64[ismissing(x) ? 0.0 : Float64(x) for x in A], 0

getptrlen(x::AbstractString) = pointer(Vector{UInt8}(x)), length(x), UInt8[]
getptrlen(x::WeakRefString{T}) where {T} = convert(Ptr{UInt8}, x.ptr), codeunits2bytes(T, x.len), UInt8[]
getptrlen(x::Missing) = convert(Ptr{UInt8}, C_NULL), 0, UInt8[]
function getptrlen(x::CategoricalArrays.CategoricalValue)
    ref = Vector{UInt8}(String(x))
    return pointer(ref), length(ref), ref
end

prep!(::Type{T}, A) where {T <: AbstractString} = _prep!(T, A)
prep!(::Type{Union{T, Missing}}, A) where {T <: AbstractString} = _prep!(T, A)
prep!(::Type{T}, A) where {T <: CategoricalValue} = _prep!(T, A)
prep!(::Type{Union{T, Missing}}, A) where {T <: CategoricalValue} = _prep!(T, A)

function _prep!(T, column)
    maxlen = maximum(ODBC.clength, column)
    data = zeros(UInt8, maxlen * length(column))
    ind = 1
    for i = 1:length(column)
        ptr, len, ref = getptrlen(column[i])
        unsafe_copyto!(pointer(data, ind), ptr, len)
        ind += maxlen
    end
    return data, maxlen
end

function prep!(column::T, col, columns, indcols) where {T}
    columns[col], maxlen = prep!(eltype(T), column)
    indcols[col] = ODBC.API.SQLLEN[clength(x) for x in column]
    return length(column), maxlen
end

getCtype(::Type{T}) where {T} = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)
getCtype(::Type{Union{T, Missing}}) where {T} = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)
getCtype(::Type{Vector{T}}) where {T} = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)
getCtype(::Type{Vector{Union{T, Missing}}}) where {T} = get(ODBC.API.julia2C, T, ODBC.API.SQL_C_CHAR)

function load!(dsn::DSN, table::AbstractString, itr::T) where {T}
    Tables.istable(T) || throw(ArgumentError("$T doesn't support the required Tables.jl interface"))
    stmt = dsn.stmt_ptr2
    ODBC.execute!(dsn, "select * from $table", stmt)
    df = DataFrame(itr)
    M, N = size(df)
    columns = Vector{Any}(undef, N)
    indcols = Vector{Any}(undef, N)
    for (col, column) in enumerate(Tables.eachcolumn(df))
        rows, len = prep!(column, col, columns, indcols)
        @CHECK stmt API.SQL_HANDLE_STMT API.SQLBindCols(stmt, col, getCtype(eltype(column)), columns[col], M, indcols[col])
    end
    API.SQLSetStmtAttr(stmt, API.SQL_ATTR_ROW_ARRAY_SIZE, M, API.SQL_IS_UINTEGER)
    @CHECK stmt API.SQL_HANDLE_STMT API.SQLBulkOperations(stmt, API.SQL_ADD)
    return
end

# function load(dsn::DSN, table::AbstractString, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T}
#     sink = Data.stream!(T(args...), ODBC.Sink, dsn, table; append=append, transforms=transforms)
#     return Data.close!(sink)
# end
# function load(dsn::DSN, table::AbstractString, source; append::Bool=false, transforms::Dict=Dict{Int,Function}())
#     sink = Data.stream!(source, ODBC.Sink, dsn, table; append=append, transforms=transforms)
#     return Data.close!(sink)
# end

# load(sink::Sink, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(T(args...), sink; append=append, transforms=transforms); return Data.close!(sink))
# load(sink::Sink, source; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))

# function Data.stream!(source, ::Type{Data.Column}, sink::ODBC.Sink, append::Bool=false)
#     Data.types(source) == Data.types(sink) || throw(ArgumentError("schema mismatch: \n$(Data.schema(source))\nvs.\n$(Data.schema(sink))"))
#     rows, cols = size(source)
#     Data.isdone(source, 1, 1) && return sink
#     ODBC.execute!(sink.dsn, "select * from $(sink.table)")
#     stmt = sink.dsn.stmt_ptr
#     types = Data.types(source)
#     columns = Vector{Any}(cols)
#     indcols = Array{Vector{ODBC.API.SQLLEN}}(cols)
#     row = 0
#     # get the column names for a table from the DB to generate the insert into sql statement
#     # might have to try quoting
#     # SQLPrepare (hdlStmt, (SQLTCHAR*)"INSERT INTO customers (CustID, CustName,  Phone_Number) VALUES(?,?,?)", SQL_NTS) ;
#     try
#         # SQLSetConnectAttr(hdlDbc, SQL_ATTR_AUTOCOMMIT, SQL_AUTOCOMMIT_OFF, SQL_NTS)
#         while !Data.isdone(source, row+1, cols+1)
#
#             for col = 1:cols
#                 T = types[col]
#                 # SQLBindParameter(hdlStmt, 1, SQL_PARAM_INPUT, SQL_C_LONG, SQL_INTEGER, 0, 0, (SQLPOINTER)custIDs, sizeof(SQLINTEGER) , NULL);
#                 rows, cT = ODBC.bindcolumn!(source, T, col, columns, indcols)
#                 ret = ODBC.API.SQLBindCols(stmt, col, cT, pointer(columns[col]), sizeof(eltype(columns[col])), indcols[col])
#                 println("$col: $ret")
#             end
#             ODBC.API.SQLSetStmtAttr(stmt, ODBC.API.SQL_ATTR_ROW_ARRAY_SIZE, rows, ODBC.API.SQL_IS_UINTEGER)
#             # SQLSetStmtAttr( hdlStmt, SQL_ATTR_PARAMSET_SIZE, (SQLPOINTER)NUM_ENTRIES, 0 );
#             # ret = SQLExecute(hdlStmt);
#             ODBC.@CHECK stmt ODBC.API.SQL_HANDLE_STMT ODBC.API.SQLBulkOperations(stmt, ODBC.API.SQL_ADD)
#             row += rows
#         end
#         # SQLEndTran(SQL_HANDLE_DBC, hdlDbc, SQL_COMMIT);
#     # finally
#         # SQLSetConnectAttr(hdlDbc, SQL_ATTR_AUTOCOMMIT, SQL_AUTOCOMMIT_ON, SQL_NTS);
#     end
#     Data.setrows!(source, row)
#     return sink
# end


text:
# Example: Use monthly average temperature/salinity fields on poseidon to compute sigma1
using Revise
using ECCOtour, MITgcmTools, MeshArrays

# other useful packages in Reemergence project
# using Statistics, PyPlot, Distributions, FFTW, LinearAlgebra, StatsBase

# get monthly average temperature and salinity
include("examples/read_TS_monavg_poseidon.jl")

p₀ = 1000 # reference pressure = 1000 dbar

# ECCOv4r4 uses approximation for pressure without any horizontal deviations.
# Can precompute pressure for each depth level.
ρ₀ = 1029 # from "data" text file in run directory
g  = 9.81 # from "data"

# depths of MITgcm levels
pathLLC = "../inputs/GRID_LLC90/" # need to input grid earlier: see example
fileZ = "RC" # "R" = radius, "C"= center of grid cell
Z = read_mdsio(pathLLC,fileZ)
Z = vec(Z)

# standard pressures via hydrostatic balance
Pa2dbar = 1/10000
p = -ρ₀ .*g .* Pa2dbar .* Z # 10000 to change Pa to dbar

# mitgcmtools.seawaterdensity is set up for a column.
# random choice of a water column.
ix = 50
iy = 400

# filter out zeroes.
θz = filter(!iszero,θ[ix,iy,:])
Sz = filter(!iszero,S[ix,iy,:])
nz = length(θz)
pz = p[1:nz]
p₀ = 1000

# ρI is in-situ density, ρR is density referenced to p0.
(ρP,ρI,ρR) = SeaWaterDensity(θz,Sz,pz,p₀)
σ₁ = ρR .- 1000

# could add loops for the entire domain.

# to do: plot it


text:
@testset "ritt-wu method" begin
    # Apollonius circle theorem
    hp = @hp begin
        points(A, B, C)
        C := Segment(A, B) ⟂ Segment(A, C)
        M₁ := Midpoint(A, B)
        M₂ := Midpoint(A, C)
        M₃ := Midpoint(B, C)
        H := A ↓ Segment(B, C)
        O := Circle(M₁, M₂, M₃)
    end

    th = @th H ∈ circle(O, M₁)
    p = prove(hp, th)
    @test p.coords[:A] == GTP.Coordinate(0, 0)
    @test p.coords[:A].y == 0
    @test all(isproved(p))

    # test with user given coordinates, example 3 of sec. 6.5 of [CLO15]
    @polyvar u[1:2] x[1:8]

    coords = OrderedDict(
       :A => Coordinate(0, 0),
       :B => Coordinate(u[1], 0),
       :C => Coordinate(0, u[2]),
       :M₁ => Coordinate(x[1], 0),
       :M₂ => Coordinate(0, x[2]),
       :M₃ => Coordinate(x[3], x[4]),
       :H => Coordinate(x[5], x[6]),
       :O => Coordinate(x[7], x[8])
    )
    p = prove(hp, th, RittWuMethod, coords, u, x)
    H_ = [
        2x[1] - u[1],
        2x[2] - u[2],
        2x[3] - u[1],
        2x[4] - u[2],
        -u[2]*x[5] - u[1]*x[6] + u[1]*u[2],
        -u[1]*x[5] + u[2]*x[6],
        x[1]^2 - x[2]^2 - 2*x[1]*x[7] + 2*x[2]*x[8],
        x[1]^2 - 2x[1]*x[7] - x[3]^2 + 2*x[3]*x[7] - x[4]^2 + 2x[4]*x[8],
    ]

    T_ = [
        2x[1] - u[1],
        2x[2] - u[2],
        2x[3] - u[1],
        2x[4] - u[2],
        -(u[1]^2 + u[2]^2) * x[5] + u[1]*u[2]^2,
        -u[1]*x[5] + u[2]*x[6],
        2*(x[1]*x[2] - x[2]*x[3] - x[1]*x[4])*x[7] - x[1]^2*x[2] + x[2]*x[3]^2 + x[1]^2*x[4] - x[2]^2*x[4] + x[2]*x[4]^2,
        x[1]^2 - 2x[1]*x[7] - x[3]^2 + 2*x[3]*x[7] - x[4]^2 + 2x[4]*x[8],
    ]

    @test p.H == H_
    T = map(x -> x ÷ gcd(x.a), p.T)
    @test T == T_

    @test ndg(p) == [-u[1]^2 - u[2]^2, u[2], 4*x[1]*x[2]-4*x[1]*x[4] - 4*x[2]*x[3], 2*x[4]]

    # Properties of parallelogram
    hp = @hp begin
        D := Parallelogram(A, B, C, D)
        O := Segment(A, C) ∩ Segment(B, D)
    end

    th = @th begin
        Segment(A, O) ≅ Segment(O, C)
        Segment(A, B) ∥ Segment(A, C) # this is obviously false
    end

    p = prove(hp, th)
    @test isproved(p) == [true, false]

    thm = Theorem(hp, th)
    @test thm isa Theorem{GTP.NoProof}
    @test !isproved(thm)

    thm = prove(thm)
    @test isproved(thm) == [true, false]

    @test string(thm) == "POINTS:\n------------\nA : free\nB : free\nC : free\nD : dependent by (1)\nO : dependent by (2)\n\nHYPOTHESIS:\n------------\n(1) ABCD parallelogram\n(2) O = AC ∩ BD\n\nTHESIS:\n------------\nAO ≅ OC\nAB ∥ AC\n\nAssigned coordinates:\n---------------------\nA = (0, 0)\nB = (u₁, 0)\nC = (u₂, u₃)\nD = (x₁, x₂)\nO = (x₃, x₄)\n\nGoal 1: success\nGoal 2: fail\n\nNondegeneracy conditions:\n-------------------------\nu₁u₃ + u₂x₂ - u₃x₁ ≠ 0\n-u₁ + x₁ ≠ 0\n"
end


text:
using Test
using DwsimOpt

@test isa(OptProblemDef, Function)


text:
#=
This file defines some simple formatting functions that may be used
in a [`Field`](#) by passing it as a `formatfn`.
=#

"""
    string_format(str)

Ensure that a string is shown with double quotes.

```julia-repl
julia> string_format("hi")
"\"hi\""
"""
string_format(x) = sprint(io -> show(io, string(x)))

code_format(x) = Markdown.parse("`$x`")

md_format(x) = Markdown.parse(x)

type_format(::T) where T = "$(T.name.name)"

@testset "Formatting functions" begin
    @testset "string_format" begin
        @test string_format("hi") == "\"hi\""

    end
    @testset "code_format" begin
        @test code_format("hi") == Markdown.parse("`hi`")
    end
end


text:
# This file is not included by any module and is not intended to exist in
# the long-term. It is just a home for julia code that was developed and
# needs to be renamed, moved to a new module, or maybe deleted.

function build(
	model, ::Type{P}, d :: Vector{D}, l :: Vector{S}, w :: Vector{S},
	L :: S, W :: S
) where {D, S, P}
	@assert length(l) == length(w) && length(w) == length(d)
	a = l .* w # area of the piece types

	partitions(P, d, l, w, L, W)
	@variables model begin
		x[1:N, 1:T], Bin  # true if piece type is placed at node, false otherwise
		v[1:N, 1:DL], Bin # true if node is a vertical cut, false otherwise
		h[1:N, 1:DW], Bin # true if node is a horizontal cut, false otherwise
		l_[1:N] >= 0 # length available to node N
		w_[1:N] >= 0 # width available to node N
	end

	@objective(model, Max, sum(p[j] * x[i, j] for i = 1:N, j = 1:T))

	# ub0: trivial area upper bound
	@constraint(model, sum(a[j] * x[i, j] for i = 1:N, j = 1:T) <= L*W)

	# c1: Each node is either: 1) a placed piece (leaf); 2) a horizontal cut
	# (branch); 3) a vertical cut (branch).
	for i = 1:N
		@constraint(model,
			sum(x[i, j] for j = 1:T) + sum(v[i, j] for j = 1:DL) +
			sum(h[i, j] for j = 1:DW) <= 1
		)
	end
	# c2: Each non-root node can only be a branch or leaf if the node parent
	# is a branch.
	for i = 2:N
		@constraint(model,
			sum(x[i, j] for j = 1:T) + sum(v[i, j] for j = 1:DL) +
			sum(h[i, j] for j = 1:DW) <= sum(v[div(i, 2), j] for j = 1:DL) +
			sum(h[div(i, 2), j] for j = 1:DW)
		)
	end

	# c3 and c4: The root node is limited by the size of the original plate.
	fix(l_[1], L; force = true)
	fix(w_[1], W; force = true)

	# c5: There is a limit on the number of pieces available for each piece type.
	for j = 1:T
		@constraint(model,
			sum(x[i, j] for i = 1:N) <= d[j]
		)
	end

	# c6 and c7: The cuts reduce the space available for their children.
	# c9 and c10: the length (width) of a plate is always at max the length
	# NOTE: it is not possible to tighten the big-M using dl[1] and dw[1] because
	# they can be children that yet keep one dimension at the same size as root.
	for i = 2:N
		i2 = div(i, 2)
		if iseven(i) # the first child needs a big-M to restrict its dimensions
			@constraints model begin
				l_[i] <= sum(dl[j] * v[i2, j] for j = 1:DL) + L * sum(h[i2, j] for j = 1:DW)
				w_[i] <= sum(dw[j] * h[i2, j] for j = 1:DW) + W * sum(v[i2, j] for j = 1:DL)
				# these are c9 and c10, they are implicit for the isodd(i) case
				l_[i] <= l_[i2]
				w_[i] <= w_[i2]
			end
		else
			# the second child does not need a big-M to restrict its dimensions,
			# nor does it need c9 and c10
			@constraints model begin
				l_[i] <= l_[i2] - sum(dl[j] * v[i2, j] for j = 1:DL)
				w_[i] <= w_[i2] - sum(dw[j] * h[i2, j] for j = 1:DW)
			end
		end
	end
	# c11 and c12: the placed pieces respect length (width) of the plate where
	# they are placed.
	for i = 1:N
		@constraints model begin
			sum(l[j] * x[i, j] for j = 1:T) <= l_[i]
			sum(w[j] * x[i, j] for j = 1:T) <= w_[i]
		end
	end

	model
end # build

# HIGH LEVEL EXPLANATION OF THE MODEL
#
# Variables:
#
# `pieces_sold[pii]`: Integer. The number of pieces `pii` sold. Is the minimum
#   between the demand of the piece and the amount of plates generated and not
#   used that have exactly the same size as the piece.
# `cuts_made[n1, n2, n3]`: Integer. The number of subplates of type `n1` that
#   are cut into subplates `n2` and `n3` (horizontal and vertical cuts are
#   together for now). As this is a symmetry-breaking model, the plate types
#   are not only each distinct `l` and `w` but each different `l`, `w`, and
#   `symm` (that marks if the plate can be cut only horizontally, only
#   vertically, or both ways).
#
# Objective function:
#
# Maximize the profit of the pieces sold.
#   sum(p[pii] * pieces_sold[pii])
#
# Constraints:
#
# There is exactly one of the original plate, which may be used for cutting
# or extracting a piece.
#   sum(pieces_sold[plates exactly the size of the original plate]) +
#   sum(cuts_made[1, _,  _]) <= 1
# The number of subplates available depends on the number of plates that have
# it as children.
#   sum(cuts_made[n1>1, _, _]) <= sum(cuts_made[_, n2, n3])
#     where n2 == n1 or n3 == n1, doubling cuts_made[_, n2, n3] if n2 == n3
# The number of pieces sold is bounded both by the demand of the piece type and
# the the number of unused plates with the same size as the piece.
#   sum(pieces_sold[pii]) <= d[pii]
#   sum(pieces_sold[pii]) <= cuts_made[_, n2, n3] - cuts_made[n2 or n3, _, _]
#     where n2 or n3 has the same size as pii, fixing for when n2 == n3
#
# Unnecessary constraints:
#
# The amount of times a plate type may be cut is bounded by how many of them
# could fit the original plate. Note that we ignore the symmetry tag here
# and group all the plates with the same `l` and `w` but distinct symmetry tag.
#   sum(cuts_made[plates sharing `l` and `w`, _, _]) <= (L ÷ l) * (W ÷ w)
function build_model_with_symmbreak(
	model, d :: Vector{D}, p :: Vector{P}, l :: Vector{S}, w :: Vector{S},
	L :: S, W :: S; only_binary = false, use_c25 = false,
	ignore_2th_dim = false, ignore_d = false, round2disc = true
) where {D, S, P}
	@assert length(d) == length(l) && length(l) == length(w)
	num_piece_types = convert(D, length(d))

	sllw = SortedLinkedLW(D, l, w)
	pli2lwsb, hcuts, vcuts, pii2plis, pli2piis, same_size_plis =
		gen_cuts_sb(P, d, sllw, L, W; ignore_2th_dim = ignore_2th_dim,
		ignore_d = ignore_d,
		round2disc = round2disc
	)
	num_plate_types = length(pli2lwsb)
	hvcuts = vcat(hcuts, vcuts)

	# The vectors below all have the same length as the number of plate types (to
	# allow indexing by plate type). The value of a position is a vector of
	# arbitrary length and irrelevant index, the values of this inner vector are
	# cut indexes. Such cut indexes are related to the plate that is the index of
	# the outer vector.
	# any CHILD plate to respective CUT indexes
	child2cut = [Vector{P}() for _ = 1:num_plate_types]
	# any PARENT plate to respective CUT indexes
	parent2cut = [Vector{P}() for _ = 1:num_plate_types]

	# Initialize all inverse indexes.
	for i in eachindex(hvcuts)
		parent, fchild, schild = hvcuts[i]
		push!(parent2cut[parent], i)
		push!(child2cut[fchild], i)
		!iszero(schild) && push!(child2cut[schild], i)
	end

	# If all pieces have demand one, a binary variable will suffice to make the
	# connection between a piece type and the plate it is extracted from.
	naturally_only_binary = all(di -> di <= 1, d)
	if naturally_only_binary || only_binary
		# only_binary is equal to naturally_only_binary for now, but the idea is
		# that only_binary will expand the number of binary variables to account
		# for piece solds that can repeat (i.e., have demand more than one).
		# NOTE that using only_binary with this model will restrict much more
		# than using only_binary with the model without symmetry, unless the demand
		# of all pieces is naturally_only_binary the model will give much worse
		# results.
		@variable(model, pieces_sold[1:num_piece_types], Bin)
	else
		@variable(model,
			0 <= pieces_sold[pii = 1:num_piece_types] <= d[pii],
		Int)
	end

	if only_binary
		@variable(model, cuts_made[1:length(hvcuts)], Bin)
	else
		@variable(model, cuts_made[1:length(hvcuts)] >= 0, Int)
	end

	# The objective function maximizes the profit of the pieces sold.
	@objective(model, Max,
		sum(p[pii] * sum(pieces_sold[pii]) for pii = 1:num_piece_types)
	)

	# c1: There is just one of the original plate, and so it can be only used to
	# extract a single piece (this is rare, because there would need to be a
	# piece the exact same size as the original plate) xor make a single cut that
	# would make two new subplates available.
	@constraint(model,
		sum(pieces_sold[pli2piis[1]]) + sum(cuts_made[parent2cut[1]]) <= 1
	)

	# c2: for each subplate type that is not the original plate, such subplate
	# type may be cut at most the number of times it was a child of another cut.
	for pli in 2:num_plate_types
		@constraints model begin
			sum(cuts_made[parent2cut[pli]]) <= sum(cuts_made[child2cut[pli]])
		end
	end

	if use_c25
		@error "not sure if correctly implemented, check before"
		# TODO: check if the below is corret. What seems to be wrong is that ssplis
		# is a vector of plate indexes, while cuts_made should be indexed by
		# vectors of cut indexes.

		# c2.5: The amount of each subplate type generated by cuts (and used either
		# as a piece or as a intermediary plate) is bounded by the amount that can be
		# cut from the original plate.
		for ssplis in same_size_plis
			@assert !iszero(length(ssplis))
			@assert isone(length(unique!(map(i -> pli2lwsb[i][1], ssplis))))
			@assert isone(length(unique!(map(i -> pli2lwsb[i][2], ssplis))))
			@constraint(model,
				sum(cuts_made[ssplis]) <= pli2lwsb[ssplis[1]][4]
			)
		end
	end

	# c3: finally, for each piece type, the amount of pieces sold of that type is
	# at most the number of plates with the piece exact size that were not cut to
	# make smaller plates.
	for pii in 1:num_piece_types
		@constraint(model, pieces_sold[pii] <= sum(cuts_made[vcat(child2cut[pii2plis[pii]]...)]) - sum(cuts_made[vcat(parent2cut[pii2plis[pii]]...)]))
	end

	model, hvcuts, pli2lwsb, pii2plis, pli2piis
end # build_model_with_symmbreak

#=
# CODE that lived inside the script method, and printed the solutions
# from the models. Need to be broken in many model-specific methods
# that take a solved model, and maybe the model building output, and
# create a common and structured representation of the solution.
		if p_args["flow-model"]
			ps = value.(m[:edge][pii] for pii = 1:length(d) if is_valid(m, m[:edge][pii]))
			ps_nz = [iv for iv in enumerate(ps) if iv[2] > 0.001]
			@show ps_nz
		else
			if p_args["break-hvcut-symmetry"]
				ps = m[:pieces_sold]
				cm = m[:cuts_made]
				ps_nz = [(i, value(ps[i])) for i = 1:length(ps) if is_valid(m, ps[i]) && value(ps[i]) > 0.001]
				cm_nz = [(i, value(cm[i])) for i = 1:length(cm) if is_valid(m, cm[i]) && value(cm[i]) > 0.001]
				@show ps_nz
				@show cm_nz
				println("(piece length, piece width) => (piece index, amount in solution, profit of single piece, total profit contributed in solution) ")
				foreach(ps_nz) do e
					pii, v = e
					println((l[pii], w[pii]) => (pii, v, p[pii], v * p[pii]))
				end
				println("(parent plate length, parent plate width) => ((first child plate length, first child plate width), (second child plate length, second child plate width))")
				foreach(cm_nz) do e
					i, _ = e
					parent, fchild, schild = hvcuts[i]
					if iszero(schild)
						println((pli2lwsb[parent][1], pli2lwsb[parent][2]) => ((pli2lwsb[fchild][1], pli2lwsb[fchild][2]), (0, 0)))
					else
						println((pli2lwsb[parent][1], pli2lwsb[parent][2]) => ((pli2lwsb[fchild][1], pli2lwsb[fchild][2]), (pli2lwsb[schild][1], pli2lwsb[schild][2])))
					end
				end
			else # same as: if !p_args["break-hvcut-symmetry"]
				p_args["relax2lp"] && @warn "relax2lp flag used, be careful regarding solution"
				ps_nz = Vector{Tuple{Int, Float64}}()
				cm_nz = Vector{Tuple{Int, Float64}}()
				for i = 1:length(ps)
					if is_valid(m, ps[i]) && value(ps[i]) > 0.0001
						push!(ps_nz, (i, value(ps[i])))
					end
				end
				for i = 1:length(cm)
					if is_valid(m, cm[i]) && value(cm[i]) > 0.0001
						push!(cm_nz, (i, value(cm[i])))
					end
				end
				#ps_nz = [(i, value(ps[i])) for i = 1:length(ps) if (is_valid(m, ps[i]) && value(ps[i]) > 0.001)]
				#cm_nz = [(i, value(cm[i])) for i = 1:length(cm) if (is_valid(m, cm[i]) && value(cm[i]) > 0.001)]
				@show ps_nz
				@show cm_nz
				println("(plate length, plate width) => (number of times this extraction happened, piece length, piece width)")
				foreach(ps_nz) do e
					i, v = e
					pli, pii = np[i]
					println((pli_lwb[pli][1], pli_lwb[pli][2]) => (v, l[pii], w[pii]))
				end
				println("(parent plate length, parent plate width) => (number of times this cut happened, (first child plate length, first child plate width), (second child plate length, second child plate width))")
				foreach(cm_nz) do e
					i, v = e
					parent, fchild, schild = hvcuts[i]
		@assert !iszero(parent)
		@assert !iszero(fchild)
		if iszero(schild)
			@assert p_args["faithful2furini2016"]
			if pli_lwb[parent][1] == pli_lwb[fchild][1]
							println((pli_lwb[parent][1], pli_lwb[parent][2]) => (v, (pli_lwb[fchild][1], pli_lwb[fchild][2]), (pli_lwb[parent][1], pli_lwb[parent][2] - pli_lwb[fchild][2])))
			else
				@assert pli_lwb[parent][2] == pli_lwb[fchild][2]
							println((pli_lwb[parent][1], pli_lwb[parent][2]) => (v, (pli_lwb[fchild][1], pli_lwb[fchild][2]), (pli_lwb[parent][1] - pli_lwb[fchild][1], pli_lwb[parent][2])))
			end
		else
						println((pli_lwb[parent][1], pli_lwb[parent][2]) => (v, (pli_lwb[fchild][1], pli_lwb[fchild][2]), (pli_lwb[schild][1], pli_lwb[schild][2])))
		end
				end
			end
		end
	end
=#

# A situation that occurred with some LP solving methods was the following:
# The method below intends to fix an unproductive behavior that happened with
# some LP methods. The method would give a very small positive reduced profit
# for a few variables, and the _iterative_pricing would keep executing because
# of it. The values were so small and the variables so few that they could not
# affect the objective value of the LP significantly, but triggered a
# resolving anyway.
function _is_significant(LP, cuts, plate_cons)
	delta_to_next_int = (one(LP) + floor(LP)) - LP
	improvement = sum(_reduced_profit.(cuts, (plate_cons,)))
	@show improvement
	@show delta_to_next_int
	return improvement >= delta_to_next_int
end

#################### START OF THE OLD WARM_START FILE ####################
"""
    _search_approx_cut(pp, fcl, fcw, max_diff, approx_l, nnn, pli_lwb) :: P

!!! **Internal use.**

"""
function _search_approx_cut(
	pp :: P, # parent plate
	fcl :: S, # first child length
	fcw :: S, # first child width
	max_diff :: S, # how much smaller the approx dimension may be
	approx_l :: Bool, # if true, length is approx, otherwise, w is approx
	nnn :: Vector{NTuple{3, P}},
	pli_lwb :: Vector{Tuple{S, S, P}},
) :: P where {D, S, P}
	for i in one(P):convert(P, length(nnn))
		(pp_, fc, _) = nnn[i]
		pp_ != pp && continue
		if approx_l
			pli_lwb[fc][1] <= fcl && pli_lwb[fc][1] >= (fcl - max_diff) &&
			pli_lwb[fc][2] == fcw && return i
		else
			pli_lwb[fc][2] <= fcw && pli_lwb[fc][2] >= (fcw - max_diff) &&
			pli_lwb[fc][1] == fcl && return i
		end
	end
	@show pp
	@show pli_lwb[pp]
	@show fcl
	@show fcw
	@show max_diff
	@assert false # this should not be reachable
end

"""
    _search_cut_or_symmetry(pp, fcl, fcw, nnn, pli_lwb)

!!! **Internal use.**

"""
function _search_cut_or_symmetry(
	pp :: P, # parent plate
	fcl :: S, # first child length
	fcw :: S, # first child width
	nnn :: Vector{NTuple{3, P}},
	pli_lwb :: Vector{Tuple{S, S, P}},
) :: P where {D, S, P}
	pll, plw, _ = pli_lwb[pp]
	@assert fcl == pll || fcw == plw
	for i in one(P):convert(P, length(nnn))
		(pp_, fc, _) = nnn[i]
		pp_ != pp && continue
		pli_lwb[fc][1] == fcl && pli_lwb[fc][2] == fcw && return i
	end
	@assert false # this should not be reachable
end

"""
    _search_cut(pp, fcl, fcw, nnn, pli_lwb)

!!! **Internal use.**

"""
function _search_cut(
	pp :: P, # parent plate
	fcl :: S, # first child length
	fcw :: S, # first child width
	nnn :: Vector{NTuple{3, P}},
	pli_lwb :: Vector{Tuple{S, S, P}},
) :: P where {D, S, P}
	for i in one(P):convert(P, length(nnn))
		(pp_, fc, _) = nnn[i]
		pp_ != pp && continue
		pli_lwb[fc][1] == fcl && pli_lwb[fc][2] == fcw && return i
	end
	@assert false # this should not be reachable
end

# TODO: Check why this method is used if a structure like SortedLinkedLW
# would answer this more efficiently and be aware of the type used for the
# index.
"""
    _min_l_fitting_piece(l, w, L, W)

!!! **Internal use.**

Given a plate `L`x`W` and two pieces dimensions list (`l`, `w`),
return the index of the piece of smallest length that fits the
plate (the width dimension may preclude this from being just
the piece of smallest length).
"""
function _min_l_fitting_piece(l, w, L, W)
	@assert length(l) == length(w)
	n = length(l)
	min_i = zero(n)
	min_l = zero(eltype(l))
	for i = one(n):n
		(l[i] > L || w[i] > W) && continue
		if iszero(min_i) || l[i] < min_l
			min_i = i
			min_l = l[i]
		end
	end
	min_i
end

# Warm start faithful2furini.
function warm_start_f2f(
	model, l, w, L, W,
	pat :: Vector{Vector{D}},
	bp :: ByproductPPG2KP{D, S, P}
	#round2disc wait to see if this is needed
	# which other model building options will need to be passed to this?
) where {D, S, P}
	# TODO: implement
end

# TODO: check if the piece will be immediatelly extracted or does it need
# an extra cut to make a plate small enough to extract the piece?
# TODO: if the piece needs an extra trim, the things are even worse than
# I thought at first, there is no guarantee the trim will be obtained
# in a single cut. It may be necessary to make many fake trim cuts to
# get a plate of the right size to make an extraction. It is just much
# simpler to add fake extraction variables for every plate that would need
# fake trims than to add all cuts needed to arrive at it. On the other side,
# this is something interesting to comment on the paper, the lack of guarantee
# on the number of cuts needed to extract the pieces from the plates, a
# 'downside' that seems not to be a problem for the solver, of have common
# worst-cases.
# TODO: comment all code this method is specially hard to follow
# TODO: the warm-start for the flag faithful2furini enabled and disabled will
# need to be different? the rules for which plates exist and which do not are
# different from one to another.
# NOTE: this method only work for simple patterns in which:
# (1) the cuts are two-staged (i.e., the pattern is justa a vector of vector);
# (1.5) by consequence, the cuts are restricted;
# (2) the first stage cuts vertically (width strips);
# (3) the first piece of each strip gives the width of the strip;
# If you need to warm-start with a more complex pattern, create another
# method with the same name, and another type for parameter `pat`.
"""
    warm_start(model, l, w, L, W, pat, pli_lwb, nnn, np; faithful2furini = false)

!!! **Internal use.**
"""
function warm_start(
	model, l, w, L, W,
	pat :: Vector{Vector{D}},
	pli_lwb :: Vector{Tuple{S, S, P}},
	nnn :: Vector{NTuple{3, P}},
	np :: Vector{Tuple{P, D}};
	faithful2furini2016 = false
	#round2disc wait to see if this is needed
	# which other model building options will need to be passed to this?
) where {D, S, P}
	@assert !faithful2furini2016
	# the initial residual plate is L, W
	# visit the outer vector in reverse
	# if the current head of stripe is smaller than half residual plate
	# then search for a vertical cut on the residual plate, with the right width
	#   for the first child and enable it, change the residual plate to
	#   be the second child
	# else assert this is the last stripe, just use the remaining plate (second
	#   child of the last cut, or the whole root if there is just one stripe)
	# after finishing the strip processing, for each plate that is a stripe:
	#   do the same as the first stage, but for the subplate and the opposite cut
	#   orientation (as the inner vectors are not sorted by length, we need to
	#   sort them ourselves);
	# finally, for every subplate that will become a piece, connect it to a piece
	#   for faithful2furini2016 we need to trim the plate and have it with exact
	#   plate size; for !faithful2furini2016 we just link directly to np
	rw = W # remaining width, initialized with root plate width
	rpli = 1 # NOTE: the root plate is guaranteed to have index 1
	cut_var_vals = Dict{P, Int}() # the amount of times each cut was made
	first_stage_plates = Vector{P}() # the plate index of all stripes
	rpat = reverse(pat)
	final = false
	for stripe in rpat
		@assert !final
		@assert !isempty(stripe)
		ws = w[first(stripe)]
		if ws <= div(rw, 2)
			cut = _search_cut(rpli, L, ws, nnn, pli_lwb)
			cut_var_vals[cut] = 1 + get(cut_var_vals, cut, 0)
			_, fc, rpli = nnn[cut]

			rw -= ws
			push!(first_stage_plates, fc)
			fcl, fcw = pli_lwb[fc]
			println("0\ti\t$(fcl)\t$(fcw)")
		else
			# If the stripe width is already more than half residual plate, then it
			# is the last stripe (the stripes are ordered in increase-or-keep order,
			# and cannot exist a larger width stripe in the remaining space).
			push!(first_stage_plates, rpli)
			scl, scw = pli_lwb[rpli]
			println("0\tl\t$(scl)\t$(scw)")
			final = true
		end
	end
	# TODO: both loops need to be by index, and stop before the last element,
	# the last element need to be checked against the min_piece_l/w and the
	# remaining space in that dimension, there may be needed to make a cut
	# that creates an unused first child and a used second child.
	# The elements before the last do not have this problem because the
	# ordering guarantee that they are smaller than half plate and therefore
	# a cut for them exist.
	# TODO: Consider: would it be simpler to just create missing variables?
	# this would work for any model type, and would make this code unaffedcted
	# by new variable reductions; the increase in the number of variables
	# would be very small.

	println("start of second stage")
	pli2pii_qt = Dict{Tuple{P, D}, Int}()
	for i in 1:length(first_stage_plates)
		# For each stripe, reset the remaining info.
		rpli = first_stage_plates[i] # get the index of the plate/stripe
		rl = pli_lwb[rpli][1] # get the length of the plate/stripe
		ws = pli_lwb[rpli][2] # stripe width (do not diminish)

		pli2pii = Vector{Tuple{P, D}}()
		for piece in sort(rpat[i], by = pii -> l[pii])
			lp = l[piece]
			wp = w[piece]
			if lp <= div(rl, 2)
				cut = _search_cut(rpli, lp, ws, nnn, pli_lwb)
				cut_var_vals[cut] = 1 + get!(cut_var_vals, cut, 0)
				_, fc, rpli = nnn[cut]

				fcl, fcw, _ = pli_lwb[fc]
				println("$(i)\ti\t$(fcl)\t$(fcw)")
				println("$(piece)\tp\t$(lp)\t$(wp)")

				min_w = _min_l_fitting_piece(w, l, fcw, fcl)
				if faithful2furini2016
				elseif fcw >= wp + min_w
					# The piece may have a small width relative to the plate width
					# (i.e., the piece of smallest width could fit there) if this is
					# the case, we cannot extract it directly from tne plate.
					# There are two options: (1) if the piece is smaller than half
					# the plate, it can be cut as first child; (2) otherwise,
					# we need to make an extra cut (a fake trim) in which
					# the second child will be a plate "not too big" to allow a
					# direct extraction of the piece from it.
					if wp <= div(fcw, 2)
						# Cut the fc again, at the piece width this time, use the
						# first child to extract the piece.
						cut2 = _search_cut(fc, lp, wp, nnn, pli_lwb)
						cut_var_vals[cut2] = 1 + get!(cut_var_vals, cut2, 0)
						_, fc2, _ = nnn[cut2]

						fc2l, fc2w, _ = pli_lwb[fc2]
						println("$(i)\ti\t$(scl)\t$(scw)")
						println("$(piece)\tp\t$(lp)\t$(wp)")

						pli2pii_qt[(fc2, piece)] = 1 + get(pli2pii_qt, (fc2, piece), 0)
					else
						trimw = fcw - wp
						# There is no guarantee of existence of a cut in which the
						# second child will have the exact size of the piece, however
						# there is guarantee that a cut between the exact size and
						# the exact size - (min_w - 1) exists. Note that cuts smaller
						# than this lower bound do not interest us, as they would have
						# the same problem again (the second child would need a fake
						# trim cut again).
						cut2 = _search_approx_cut(
							rpli, lp, trimw, min_w - 1, false, nnn, pli_lwb
						)
						cut_var_vals[cut2] = 1 + get!(cut_var_vals, cut2, 0)
						_, _, sc = nnn[cut2]

						scl, scw, _ = pli_lwb[sc]
						println("$(i)\ti\t$(scl)\t$(scw)")
						println("$(piece)\tp\t$(lp)\t$(wp)")

						pli2pii_qt[(sc, piece)] = 1 + get(pli2pii_qt, (sc, piece), 0)
					end
				else
					pli2pii_qt[(fc, piece)] = 1 + get(pli2pii_qt, (fc, piece), 0)
				end

				rl -= lp
			else # if it is the last piece of a strip
				pll, plw, _ = pli_lwb[rpli]
				@assert plw == ws
				min_l_pii = _min_l_fitting_piece(l, w, pll, plw)
				min_l = iszero(min_l_pii) ? pll : l[min_l_pii]
				if pll >= lp + min_l
					triml = pll - lp
					cut = _search_approx_cut(
						rpli, triml, ws, min_l - 1, true, nnn, pli_lwb
					)
					cut_var_vals[cut] = 1 + get!(cut_var_vals, cut, 0)
					_, _, sc = nnn[cut]

					scl, scw, _ = pli_lwb[sc]
					println("$(i)\tl\t$(scl)\t$(scw)")
					println("$(piece)\tp\t$(lp)\t$(wp)")
					rpli = sc
				end
				pll, plw, _ = pli_lwb[rpli]
				@assert plw == ws
				min_w_pii = _min_l_fitting_piece(w, l, plw, pll)
				min_w = iszero(min_w_pii) ? plw : l[min_w_pii]
				if plw >= wp + min_w
					trimw = plw - wp
					cut = _search_approx_cut(
						rpli, pll, trimw, min_w - 1, false, nnn, pli_lwb
					)
					cut_var_vals[cut] = 1 + get!(cut_var_vals, cut, 0)
					_, _, sc = nnn[cut]

					scl, scw, _ = pli_lwb[sc]
					println("$(i)\tl\t$(scl)\t$(scw)")
					println("$(piece)\tp\t$(lp)\t$(wp)")
					rpli = sc
				end

				pli2pii_qt[(rpli, piece)] = 1 + get(pli2pii_qt, (rpli, piece), 0)
			end
		end
	end
	for (var_index, var_value) in cut_var_vals
		set_start_value(model[:cuts_made][var_index], var_value)
	end
	if faithful2furini2016
		# HERE WE NEED TO DO THE TRIM CUT IF FAITHFUL2FURINI IS ENABLED
	else
		for (pair, value) in pli2pii_qt
			@show pli_lwb[pair[1]]
			@show l[pair[2]], w[pair[2]]
			picut_idx = findfirst(isequal(pair), np)
			@show picut_idx
			set_start_value(model[:picuts][picut_idx], value)
		end
	end

	model
end

# TODO: Consider if this should be enabled back again considering that, for
# now, the idea of a symmetry-breaking version of PP-G2KP was abandoned.
# TODO: internal change. The symmetry tag dimension of plis should be the
# first one, given how columnwise memory layout works.
#=
function gen_cuts_sb(
	::Type{P}, d :: Vector{D}, sllw :: SortedLinkedLW{D, S}, L :: S, W :: S;
	ignore_2th_dim = false, ignore_d = false, round2disc = true
) where {D, S, P}
	(ignore_d || ignore_2th_dim || round2disc) && @error "ignore_2th_dimm, ignore_d, and round2disc are not yet implemented for gen_cuts_sb, also, first improve the performance by memoizing the discretizations"
	l = sllw.l
	w = sllw.w
	@assert length(d) == length(l)
	@assert length(d) == length(w)
	only_single_l = discretize(
		d, l, w, L, W; only_single_pieces = true, ignore_W = ignore_2th_dim,
		ignore_d = ignore_d
	)
	only_single_w = discretize(
		d, w, l, W, L; only_single_pieces = true, ignore_W = ignore_2th_dim,
		ignore_d = ignore_d
	)
	max_piece_type = convert(D, length(l))
	max_num_plates = convert(P, L) * convert(P, W) * 3
	min_pil = minimum(l)
	min_piw = minimum(w)
	# If (n1, n2, n3) is in hcut (vcuts), then n1 may be partitioned in n2 and n3
	# by the means of a horizontal (vertical) cut. Not every possible partition
	# is present, just the ones which each child plate can fit at least one
	# piece, and the ones where one child plate has the same dimension as a piece
	# and the other is the dummy waste plate.
	hcuts = Vector{NTuple{3, P}}()
	vcuts = Vector{NTuple{3, P}}()
	# If the plate has the exact same size that a piece, then it is added to
	# the vector at pii2plis[piece_type].
	pii2plis = [Vector{P}() for _ = 1:max_piece_type]
	# If the piece has the exact same size that a plate, then it is added to
	# the vector at pli2piis[plate_type].
	pli2piis = Vector{D}[]
	# The list of plates attributes: plate length, plate width, plate symmetry,
	# and plate bound. The plate index is the same as the index in pli_lwsb.
	pli2lwsb = Vector{Tuple{S, S, UInt8, 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".
	plis = zeros(P, L, W, 3)
	plis[L, W, 3] = one(P)
	# next: plates already indexed but not yet processed, starts with (L, W,
	# 3, 1). The size of the original plate, the code for being able to cut
	# in any direction (same as third dimension of plis), and the plate index.
	# 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, UInt8, P}}()
	next_idx = one(P)
	#sizehint!(next, max_num_plates)
	push!(next, (L, W, 3, one(P)))
	# n: The amount of plates (the index of the highest plate type).
	n = one(P) # there is already the original plate
	piece_cuts_avoided_by_symmb = 0
	while next_idx <= length(next)
		# PLate Length, Width, Symmetry, and Index
		pll, plw, pls :: UInt8, pli = next[next_idx]
		next_idx += 1
		plb = (L ÷ pll) * (W ÷ plw) # PLate Bound
		# It is not necessary to store the plate id in pli2lwsb because they are
		# added in order (a queue is used), so the array index is the plate index.
		push!(pli2lwsb, (pll, plw, pls, plb))
		# the pli2piis vector grows with the number of plates processed
		push!(pli2piis, D[])
		for pii in 1:max_piece_type # pii: PIece Index
			pil, piw = l[pii], w[pii] # PIece Length and Width (homophones, I know)
			if pil == pll && piw == plw
				# If the current plate is exactly the size of a plate, add it to
				# the list of plate codes that correspond to a piece.
				push!(pii2plis[pii], pli)
				push!(pli2piis[pli], pii)
			elseif should_extract_piece_from_plate(pii, pll, plw, sllw, pls)
				# If there is no way to cut the current piece from the current plate
				# except if we cut after the half of the plate... (i.e., the most
				# unbalanced cut possible, that is, making a cut with the minimal
				# length or width among the pieces, removes the possibility of getting
				# such piece from the larger child plate)
				if pil < pll && (pls == 1 || pls == 3)
					# If the current plate is allowed to be cut horizontally (i.e.,
					# reducing the length) and we need to reduce the length to make the
					# plate into the piece, then we do so.
					if iszero(plis[pil, plw, 2])
						push!(next, (pil, plw, 2, n += 1))
						plis[pil, plw, 2] = n
					end
					push!(hcuts, (pli, plis[pil, plw, 2], 0))
				elseif piw < plw && (pls == 2 || pls == 3)
					# If the current plate is allowed to be cut vertically (i.e.,
					# reducing the width) and we need to reduce the width to make the
					# plate into the piece, then we do so.
					if iszero(plis[pll, piw, 1])
						push!(next, (pll, piw, 1, n += 1))
						plis[pll, piw, 1] = n
					end
					push!(vcuts, (pli, plis[pll, piw, 1], 0))
				else
					piece_cuts_avoided_by_symmb += 1
				end
			end
		end
		# TODO: after is working, refactor the code below to use a single loop
		# using a vector saved in a variable, or an inner function
		sorted_in(a, v) = searchsortedfirst(a, v) <= length(a)
		if pls == 1 && sorted_in(only_single_w, plw)
			indexes = filter(i -> w[i] == plw && l[i] < pll, collect(1:max_piece_type))
			#@assert !isempty(indexes)
			for y in l[indexes]
				if iszero(plis[y, plw, 2])
					push!(next, (y, plw, 2, n += 1))
					plis[y, plw, 2] = n
				end
				if iszero(plis[pll - y, plw, 3])
					push!(next, (pll - y, plw, 3, n += 1))
					plis[pll - y, plw, 3] = n
				end
				push!(hcuts, (pli, plis[y, plw, 2], plis[pll - y, plw, 3]))
			end
		elseif pls == 1 || pls == 3
			for y in discretize(
				d, l, w, pll ÷ 2, plw; ignore_2th_dim = ignore_2th_dim
			)
				#y > pll ÷ 2 && break
				@assert plw >= min_piw
				@assert y >= min_pil
				@assert pll - y >= min_pil
				if iszero(plis[y, plw, 2])
					push!(next, (y, plw, 2, n += 1))
					plis[y, plw, 2] = n
				end
				if iszero(plis[pll - y, plw, 3])
					push!(next, (pll - y, plw, 3, n += 1))
					plis[pll - y, plw, 3] = n
				end
				push!(hcuts, (pli, plis[y, plw, 2], plis[pll - y, plw, 3]))
			end
		end
		if pls == 2 && sorted_in(only_single_l, pll)
			indexes = filter(i -> l[i] == pll && w[i] < plw, collect(1:max_piece_type))
			#@assert !isempty(indexes)
			for x in w[indexes]
				if iszero(plis[pll, x, 1])
					push!(next, (pll, x, 1, n += 1))
					plis[pll, x, 1] = n
				end
				if iszero(plis[pll, plw - x, 3])
					push!(next, (pll, plw - x, 3, n += 1))
					plis[pll, plw - x, 3] = n
				end
				push!(vcuts, (pli, plis[pll, x, 1], plis[pll, plw - x, 3]))
			end
		elseif pls == 2 || pls == 3
			for x in discretize(
				d, w, l, plw ÷ 2, pll; ignore_2th_dim = ignore_2th_dim
			)
				#x > plw ÷ 2 && break
				@assert pll >= min_pil
				@assert x >= min_piw
				@assert plw - x >= min_piw
				if iszero(plis[pll, x, 1])
					push!(next, (pll, x, 1, n += 1))
					plis[pll, x, 1] = n
				end
				if iszero(plis[pll, plw - x, 3])
					push!(next, (pll, plw - x, 3, n += 1))
					plis[pll, plw - x, 3] = n
				end
				push!(vcuts, (pli, plis[pll, x, 1], plis[pll, plw - x, 3]))
			end
		end
	end
	#@show piece_cuts_avoided_by_symmb

	# For each possible subplatex size, same_size_plis groups the plate indexes
	# that refer to the plates with the same size, and distinct symmetry tags.
	same_size_plis = Vector{Vector{P}}()
	# The length of same_size_plis may be as low as one third of length(pli2lwsb)
	# and as high as length(pli2lwsb).
	sizehint!(same_size_plis, length(pli2lwsb))
	for j = 1:W, i = 1:L # column-major for performance
		temp = Vector{P}()
		!iszero(plis[i, j, 1]) && push!(temp, plis[i, j, 1])
		!iszero(plis[i, j, 2]) && push!(temp, plis[i, j, 2])
		!iszero(plis[i, j, 3]) && push!(temp, plis[i, j, 3])
		!iszero(length(temp)) && (push!(same_size_plis, temp); temp = Vector{P}())
	end
	@assert last(same_size_plis) == P[1]
	pop!(same_size_plis)

	return pli2lwsb, hcuts, vcuts, pii2plis, pli2piis, same_size_plis
end
=#

# CODE USED TO SERIALIZE AND UNSERIALIZE THE CONSTRAINT DUALS
# USED IN FILE src/ppg2kp/pricing.jl and METHOD _restricted_final_pricing!
#=
	@assert has_duals(model)
	plate_cons = model[:plate_cons]
	#plate_cons_duals = Serialization.deserialize("./plate_cons_duals_gcut9_simplex")
	plate_cons_duals = dual.(plate_cons)
	#println("plate_cons_duals stats")
	#vector_summary(plate_cons_duals)
	#Serialization.serialize("./plate_cons_duals_gcut9_simplex", plate_cons_duals)
	# This is a little costy, and not necessary anymore.
	#println("stats on reduced profit values of restricted final pricing")
	#rps = _reduced_profit.(rc_cuts, (plate_cons,))
	rps = _reduced_profit.(rc_cuts, (plate_cons_duals,))
=#


text:
################################################################################
# PrimalDualTraj
################################################################################

mutable struct PrimalDualTraj{KN,n,m,T,SVd}
    probsize::ProblemSize
    pr::Traj{n,m,T,KN} # Primal trajectory
    du::Vector{SVd} # Dual trajectory
end

function PrimalDualTraj(probsize::ProblemSize, dt::T; f=rand, amplitude=1e-8) where {T}
    N = probsize.N
    n = probsize.n
    m = probsize.m
    p = probsize.p
    pr = Traj(n,m,dt,N)
    du = [[amplitude*f(SVector{n,T}) for k=1:N-1] for i=1:p]
    for k = 1:N
        pr[k].z = amplitude*f(n+m)
    end
    TYPE = (eltype(pr),n,m,T,eltype(du))
    return PrimalDualTraj{TYPE...}(probsize, pr, du)
end

################################################################################
# Methods
################################################################################

function init_traj!(pdtraj::PrimalDualTraj{KN,n,m,T,SVd}; x0=1e-8*rand(SVector{n,T}),
    s::Int=2^10, f=rand, amplitude=1e-8) where {KN,n,m,T,SVd}
    N = pdtraj.probsize.N
    p = pdtraj.probsize.p

    for k = 1:N
        pdtraj.pr[k].z = (k+s<=N) ? pdtraj.pr[k+s].z : amplitude*f(SVector{n+m,T})
    end
    for i = 1:p
        for k = 1:N-1
            pdtraj.du[i][k] = (k+s<=N-1) ? pdtraj.du[i][k+s] : amplitude*f(SVector{n,T})
        end
    end
    RobotDynamics.set_state!(pdtraj.pr[1], x0)
    return nothing
end

function set_traj!(core::NewtonCore, Δpdtraj::PrimalDualTraj{KN,n,m,T,SVd},
    Δtraj::AbstractVector) where {KN,n,m,T,SVd}
    N = core.probsize.N
    p = core.probsize.p
    pu = core.probsize.pu
    # Primals
    indu = zeros(Int,m)
	stamp = HStamp()
    for k = 1:N-1
        # States
		stampify!(stamp, :x, 1, k+1)
        ind = horizontal_idx(core, stamp)
        RobotDynamics.set_state!(Δpdtraj.pr[k+1], Δtraj[ind])
        # Controls
        for i = 1:p
			stampify!(stamp, :u, i, k)
            indu[pu[i]] = horizontal_idx(core, stamp)
        end
        RobotDynamics.set_control!(Δpdtraj.pr[k], Δtraj[indu])
    end
    # Duals
    for i = 1:p
        for k = 1:N-1
			stampify!(stamp, :λ, i, k)
            ind = horizontal_idx(core, stamp)
            Δpdtraj.du[i][k] = Δtraj[ind]
        end
    end
    return nothing
end

function get_traj!(core::NewtonCore, Δpdtraj::PrimalDualTraj{KN,n,m,T,SVd},
    Δtraj::AbstractVector) where {KN,n,m,T,SVd}
    N = core.probsize.N
    p = core.probsize.p
    pu = core.probsize.pu
    # Primals
	indu = zeros(Int,m)
	stamp = HStamp()
	for k = 1:N-1
		# States
		stampify!(stamp, :x, 1, k+1)
		ind = horizontal_idx(core, stamp)
		Δtraj[ind] = state(Δpdtraj.pr[k+1])
        # Controls
        for i = 1:p
			stampify!(stamp, :u, i, k)
            indu[pu[i]] = horizontal_idx(core, stamp)
        end
        Δtraj[indu] = control(Δpdtraj.pr[k])
    end
    # Duals
    for i = 1:p
        for k = 1:N-1
			stampify!(stamp, :λ, i, k)
            ind = horizontal_idx(core, stamp)
            Δtraj[ind] = Δpdtraj.du[i][k]
        end
    end
    return nothing
end

function update_traj!(target::PrimalDualTraj{KN,n,m,T,SVd},
    source::PrimalDualTraj{KN,n,m,T,SVd}, α::T,
    Δ::PrimalDualTraj{KN,n,m,T,SVd}) where {KN,n,m,T,SVd,SVx}
    N = target.probsize.N
    p = target.probsize.p
    # Primals
    for k = 1:N-1
        # States
        RobotDynamics.set_state!(target.pr[k+1], state(source.pr[k+1]) + α*state(Δ.pr[k+1]))
        # Controls
        RobotDynamics.set_control!(target.pr[k], control(source.pr[k]) + α*control(Δ.pr[k]))
    end
    # Duals
    for i = 1:p
        for k = 1:N-1
            target.du[i][k] = source.du[i][k] + α*Δ.du[i][k]
        end
    end
    return nothing
end

function Δ_step(pdtraj::PrimalDualTraj, α::T) where {T}
	s = 0.0
	N = pdtraj.probsize.N
	n = pdtraj.probsize.n
	m = pdtraj.probsize.m
	p = pdtraj.probsize.p
	for k = 1:N-1
		s += norm(state(pdtraj.pr[k+1]), 1) # xk+1
		s += norm(control(pdtraj.pr[k]), 1) # uk
		for i = 1:p
			# s += norm(pdtraj.du[i][k], 1) # λik
		end
	end
	s *= α
	# s /= (N-1)*(n+m+n*p)
	s /= (N-1)*(n+m)
	return s
end

function reset_duals!(pdtraj::PrimalDualTraj) where {T}
	N = pdtraj.probsize.N
	p = pdtraj.probsize.p
	for k = 1:N-1
		for i = 1:p
			pdtraj.du[i][k] *= 0.0
		end
	end
	return nothing
end


text:
immutable FDist <: ContinuousUnivariateDistribution
    ndf::Float64
    ddf::Float64
    function FDist(d1::Real, d2::Real)
    	d1 > zero(d1) && d2 > zero(d2) ||
    	    error("Numerator and denominator degrees of freedom must be positive")
  	new(float64(d1), float64(d2))
    end
end

@_jl_dist_2p FDist f

@continuous_distr_support FDist 0.0 Inf

mean(d::FDist) = 2.0 < d.ddf ? d.ddf / (d.ddf - 2.0) : NaN

median(d::FDist) = quantile(d, 0.5)

mode(d::FDist) = d.ndf <= 2. ? 0.0 : (d.ndf - 2.) / d.ndf * d.ddf / (d.ddf + 2.)

var(d::FDist) = d.ddf > 4. ?  2.0 * d.ddf^2 *
		       (d.ndf + d.ddf - 2.0) /
		       (d.ndf * (d.ddf - 2.0)^2 * (d.ddf - 4.0)) : NaN

skewness(d::FDist) = d.ddf > 6. ?  (2.d.ndf + d.ddf - 2.) * sqrt(8.(d.ddf - 4.)) /
           ((d.ddf - 6.) * sqrt(d.ndf * (d.ndf + d.ddf - 2.))) : NaN

function kurtosis(d::FDist)
    d.ddf <= 8. && return NaN
    a = d.ndf * (5. * d.ddf - 22.) * (d.ndf + d.ddf - 2.) +
        (d.ddf - 4.) * (d.ddf - 2.)^2
    b = d.ndf * (d.ddf - 6.) * (d.ddf - 8.) * (d.ndf + d.ddf - 2.)
    12. * a / b
end

entropy(d::FDist) = (log(d.ddf) -log(d.ndf)
                     +lgamma(0.5*d.ndf) +lgamma(0.5*d.ddf) -lgamma(0.5*(d.ndf+d.ddf))
                     +(1.0-0.5*d.ndf)*digamma(0.5*d.ndf) +(-1.0-0.5*d.ddf)*digamma(0.5*d.ddf)
                     +0.5*(d.ndf+d.ddf)*digamma(0.5*(d.ndf+d.ddf)))


text:
struct BlockIndices
    i1::UnitRange{Int}
    i2::UnitRange{Int}
    isdiag::Bool
end

# Remove isdiag field from comparison for dict
Base.isequal(b1::BlockIndices, b2::BlockIndices) = b1.i1 == b2.i1 && b1.i2 == b2.i2
Base.hash(b::BlockIndices, h::UInt) = hash((b.i1, b.i2), h)

"""
    SparseBlocks

A type for speeding up expressions like the following:

    A[rows, cols] .= data
    A[rows, cols] .+= data

where `A` is a `SparseMatrixCSC` and `rows` and `cols` are `UnitRange`s. The
`data` is usually a `Matrix` but will work with any `AbstractArray` or
collection that supports broadcasting.

The key assumption is that the sparsity in `A` is known and fixed, specified by
a list of "blocks," or a chunk of consecutive rows and columns, addressable by
two `UnitRange` objects.

Usage consists of two phases: 1) initialization, where the sparsity structure is
specified, and 2) usage, where the precomputed data in this struct is used to
speed up indexing operations.

## Initialization
To specify the sparsity structure, use [`addblock!`](@ref). For example, we could
do the following:

    sb = SparseBlocks(10,10)  # blocks for a 10x10 sparse matrix
    addblock!(sb, 1:2, 1:2)   # nonzeros block in top left corner
    addblock!(sb, 1:3, 8:10)  # nonzeros block in top right corner

After adding all of our blocks, we initialize the sparsity structure via
[`initialize!`](@ref):

    A = spzeros(10,10)
    initialize!(sb, A)

After calling `initialize!`, the matrix `A` will have the correct number of nonzero
entries. The speed of the indexing operations is given by using cached conversions
from the entries of each block to the storage location in the nonzeros vector of `A`,
so it's critical that the sparsity structure of `A` does not change!

## Usage
Once initialized, we can use the cached information in this object to speed up
writing to our sparse matrix `A`. Basic usage is as simple as wrapping our indices
with an indexing call to our `SparseBlocks` object:

    A[sb[1:2, 1:2]] .= randn(2,2)

Note that these indices must exactly match one of the argument pairs passed to
[`addblock`](@ref). For best performance, you can cache the [`SparseBlockIndex`](@ref)
returned from this object:

    inds = sb[1:2, 1:2]  # cache this somewhere for future indexing calls
    A[inds] .= randn(2,2)

All of these methods avoid dynamic memory allocations and are about 5-20x faster
indexing directly into the sparse matrix.
"""
struct SparseBlocks
    n::Int
    m::Int
    inds::Dict{BlockIndices,Matrix{Int}}
    SparseBlocks(n,m) = new(n,m,Dict{BlockIndices,Matrix{Int}}())
end

"""
    addblock!(blocks, rows, cols)

Add a nonzero block for the sparse matrix, spanning the block defined by
   consecutive `rows` and `cols` ranges. All calls to `addblock!` should be
   made before calling [`initialize!`](@ref).
"""
function addblock!(blocks::SparseBlocks, i1::UnitRange, i2::UnitRange, isdiag::Bool=false)
    if i1.start < 1 || i1.stop > blocks.n || i2.start < 1 || i2.stop > blocks.m
        throw(BoundsError(spzeros(blocks.n, blocks.m), (i1, i2)))
    end
    block = BlockIndices(i1, i2, isdiag)
    if isdiag
        blocks.inds[block] = zeros(Int, min(length(i1), length(i2)), 1)
    else
        blocks.inds[block] = zeros(Int, length(i1), length(i2))
    end
    blocks
end

"""
    initialize!(blocks::SparseBlocks, A::SparseMatrixCSC)

Initialize the sparsity structure specified by `blocks` (via previous calls
to [`addblock!`](@ref)). Writes epsilon values to the nonzero blocks in `A`
to initialize the data storage for the sparse matrix. The indices into this
storage is computed and stored in `blocks`.
"""
function initialize!(blocks::SparseBlocks, A::SparseMatrixCSC{T}) where T
    n,m = blocks.n, blocks.m
    if size(A) != (n,m)
        throw(DimensionMismatch("Dimension of sparse matrix doesn't match expected size of ($n,$m)."))
    end
    @assert size(A) == (blocks.n, blocks.m)
    for block in keys(blocks.inds)
        if block.isdiag
            for (i,j) in zip(block.i1, block.i2)
                A[i,j] = eps(T)
            end
        else
            A[block.i1, block.i2] .= eps(T)
        end
    end
    inds = blocks.inds
    for (block,inds) in pairs(blocks.inds)
        rows = block.i1
        cols = block.i2
        n = length(rows)
        if block.isdiag
            for (j,(r,c)) in enumerate(zip(rows, cols))
                nzind = getnzind(A, r, c)
                if !isempty(nzind)
                    inds[j] = nzind[1]
                end
            end
        else
            for (j,col) in enumerate(cols)
                nzind = getnzind(A, rows[1], col)
                if !isempty(nzind)
                    inds[:,j] .= nzind[1] - 1 .+ (1:n)
                end
            end
        end
    end
    blocks
end

# Get the index into the nonzero vector for the given row and column
function getnzind(A,row,col)
    rows = view(A.rowval, nzrange(A,col))
    istart = searchsorted(rows, row) .+ A.colptr[col] .- 1
    return istart
end

"""
    SparseBlockIndex

A custom index for sparse arrays with a fixed sparsity structure, especially for
    those whose nonzeros appear as dense blocks. This is created by indexing
    into a [`SparseBlocks`] object, using indices to one of the cached blocks.

This object can then be used as a index for those same indices, creating a
    [`SparseBlockView`](@ref) object that efficiently maps cartesian indices
    into the nonzeros vector of the sparse matrix.
"""
struct SparseBlockIndex
    block::BlockIndices
    nzinds::Matrix{Int}
end

function Base.getindex(blocks::SparseBlocks, i1::UnitRange, i2::UnitRange)
    block0 = BlockIndices(i1, i2, false)
    block = getkey(blocks.inds, block0, nothing)
    if isnothing(block)
        throw(KeyError((i1, i2)))
    end
    SparseBlockIndex(block, blocks.inds[block])
end

"""
    SparseBlockView

A custom view into a sparse matrix with known sparsity, and whose nonzero entries
appears as dense blocks. Can be used a normal array, where operations on the elements
write and read directly from the nonzeros vector of the sparse matrix.
"""
struct SparseBlockView{Tv,Ti} <: AbstractMatrix{Tv}
    data::SparseMatrixCSC{Tv,Ti}
    block::SparseBlockIndex
end

# Needed to support expressions like A[B] .= data
Broadcast.dotview(A::AbstractSparseMatrix, block::SparseBlockIndex) = SparseBlockView(A, block)

# Define this type when indexing into a SparseMatrixCSC with a SparseBlockIndex
Base.getindex(A::AbstractSparseMatrix, block::SparseBlockIndex) = SparseBlockView(A, block)

# Array interface
Base.size(B::SparseBlockView) = size(B.block.nzinds)
Base.getindex(B::SparseBlockView, index) = getindex(B.data.nzval, B.block.nzinds[index])
Base.getindex(B::SparseBlockView, I::Vararg{Int,2}) = getindex(B.data.nzval, B.block.nzinds[I...])
Base.setindex!(B::SparseBlockView, v, index) = setindex!(B.data.nzval, v, B.block.nzinds[index])
Base.setindex!(B::SparseBlockView, v, I::Vararg{Int,2}) = setindex!(B.data.nzval, v, B.block.nzinds[I...])
Base.IndexStyle(::SparseBlockView) = IndexLinear()

# Broadcasting interface
struct SparseViewStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:SparseBlockView}) = SparseViewStyle()

# Always use the new broadcasting style over others
Base.BroadcastStyle(::SparseViewStyle, ::Broadcast.BroadcastStyle) = SparseViewStyle()

# Handle generic expressions like B .= f.(arg)
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{Broadcast.DefaultMatrixStyle})
    for i in eachindex(bc)
        B[i] = bc[i]
    end
    B
end

# Handle generic expressions like B .= f.(B, g.(args...))
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{SparseViewStyle})
    for i in eachindex(bc)
        B[i] = bc.f(B[i], bc[i])
    end
    B
end

# Specialization of B .= data::AbstractMatrix
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{Broadcast.DefaultMatrixStyle,Nothing,typeof(identity),<:Tuple{<:AbstractMatrix}})
    data = bc.args[1]
    for i in eachindex(B)
        B[i] = data[i]
    end
    B
end

# Specialization of B .= f.(B, data::AbstractMatrix)
#   which includes B .+= data::AbstractMatrix and variants
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{SparseViewStyle,Nothing,<:Any,<:Tuple{<:SparseBlockView,<:AbstractArray}})
    data = bc.args[2]
    for i in eachindex(B)
        B[i] = bc.f(B[i], data[i])
    end
    B
end

# Specialization of B .= D::Diagonal
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{<:LinearAlgebra.StructuredMatrixStyle{<:Diagonal}, Nothing, typeof(identity), <:Tuple{<:Diagonal}})
    D = bc.args[1]
    if B.block.block.isdiag
        for i in eachindex(B)
            B[i] = D.diag[i]
        end
    else
        for i = 1:minimum(size(B))
            B[i,i] = D.diag[i]
        end
    end
end

function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{SparseViewStyle,Nothing,<:Any,<:Tuple{<:SparseBlockView,<:Diagonal}})
    D = bc.args[2]
    f = bc.f
    if B.block.block.isdiag
        for i in eachindex(B)
            B[i] = f(B[i], D.diag[i])
        end
    else
        for i = 1:minimum(size(B))
            B[i,i] = f(B[i], D.diag[i])
        end
    end
end

# Specialize B' .= bc
function Broadcast.materialize!(B::Adjoint{T, <:SparseBlockView{T}}, bc::Broadcast.Broadcasted{Broadcast.DefaultMatrixStyle,Nothing,typeof(identity),<:Tuple{<:AbstractMatrix}}) where T
    data = bc.args[1]
    for i in eachindex(bc)
        it = CartesianIndex(i[2], i[1])
        B.parent[it] = data[i]
    end
    B
end

# Specialize B = UpperTriangular(data)
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{<:LinearAlgebra.StructuredMatrixStyle{<:UpperTriangular}})
    for r = 1:size(B,1)
        for c = r:size(B,2)
            i = CartesianIndex(r,c)
            B[r,c] = bc[i]
        end
    end
end

# Specialize B = LowerTriangular(data)
function Broadcast.materialize!(B::SparseBlockView, bc::Broadcast.Broadcasted{<:LinearAlgebra.StructuredMatrixStyle{<:LowerTriangular}})
    for c = 1:size(B,2)
        for r = c:size(B,1)
            i = CartesianIndex(r,c)
            B[r,c] = bc[i]
        end
    end
end


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using Plots
# setup
finep = 4
coarsep = 1
numbercomponents = 1
dimension = 1
mesh = Mesh1D(1.0)

ctofbasis = TensorH1LagrangePProlongationBasis(coarsep+1, finep+1, numbercomponents, dimension)

# diffusion operators
finediffusion = GalleryOperator("diffusion", finep+1, finep+1, mesh)
coarsediffusion = GalleryOperator("diffusion", coarsep+1, finep+1, mesh)

# Chebyshev smoother
nos = Schwarz(finediffusion)

# p-multigrid preconditioner
multigrid = PMultigrid(finediffusion, coarsediffusion, nos, [ctofbasis])

# full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# setup
ω = [2.0/3.0]
eigenvalues = zeros(numbersteps, finep)

# compute
for i in 1:numbersteps
    θ = [θ_range[i]]
    if abs(θ[1]) >  π/512
        A = computesymbols(multigrid, ω, [1,1], θ)
        #currenteigenvalues = [real(val) for val in eigvals(I-A)]
        currenteigenvalues = [real(val) for val in eigvals(A)]
        eigenvalues[i, :] = currenteigenvalues
    end
end

# plot
xrange = θ_range/π
plot(
    xrange,
    xlabel="θ/π",
    xtickfont=font(12, "Courier"),
    eigenvalues,
    ytickfont=font(12, "Courier"),
    ylabel="λ",
    linewidth=3,
    legend=:none,
    title="Spectrum of NOS Symbol",
    palette=palette(:tab10)
)
ylims!(min(0.0, eigenvalues...) * 1.1, max(eigenvalues...) * 1.1)
savefig("pmg_nos_spectrum_4_1")


text:
# TODO: add finalizers for all disposable objects
# TODO: ensure all require _poll methods are invoked in fixed intervals

mutable struct KafkaClient
    conf::Dict{Any, Any}
    typ::Cint
    rk::Ptr{Cvoid}
end

# # TODO: save task and stop it in finalizer
# @async begin
#     while true
#         println("calling kafka_poll")
#         kafka_poll(rk, 1000)
#     end
# end

function KafkaClient(typ::Integer, conf::Dict=Dict())
    c_conf = kafka_conf_new()
    for (k, v) in conf
        kafka_conf_set(c_conf, string(k), string(v))
    end
    rk = kafka_new(c_conf, Cint(typ))
    client = KafkaClient(conf, typ, rk)
    # seems like `kafka_destroy` also destroys its config, so we don't attempt it twice
    #finalizer(client, client -> kafka_destroy(rk))
    return client
end

Base.show(io::IO, kc::KafkaClient) = print(io, "KafkaClient($(kc.typ))")

mutable struct KafkaTopic
    conf::Dict{Any, Any}
    topic::String
    rkt::Ptr{Cvoid}
end

function KafkaTopic(kc::KafkaClient, topic::String, conf::Dict=Dict())
    c_conf = kafka_topic_conf_new()
    for (k, v) in conf
        kafka_topic_conf_set(c_conf, string(k), string(v))
    end
    rkt = kafka_topic_new(kc.rk, "mytopic", c_conf)
    topic = KafkaTopic(conf, topic, rkt)
    #finalizer(topic, topic -> (kafka_topic_destroy(rkt)))
    return topic
end

Base.show(io::IO, kt::KafkaTopic) = print(io, "KafkaTopic($(kt.topic))")

## KafkaProducer

mutable struct KafkaProducer
    client::KafkaClient
    topics::Dict{String, KafkaTopic}
end

function KafkaProducer(conf::Dict)
    kc = KafkaClient(KAFKA_TYPE_PRODUCER, conf)
    return KafkaProducer(kc, Dict())
end

function KafkaProducer(bootstrap_servers::String, conf::Dict=Dict())
    conf["bootstrap.servers"] = bootstrap_servers
    return KafkaProducer(conf)
end

function Base.show(io::IO, p::KafkaProducer)
    bootstrap_servers = p.client.conf["bootstrap.servers"]
    print(io, "KafkaProducer($bootstrap_servers)")
end

function produce1(kt::KafkaTopic, partition::Integer, key, payload)
    #produce(kt.rkt, partition, key, payload)
    produce(kt.rkt, partition, Vector{UInt8}(key), Vector{UInt8}(payload))
end

function produce1(p::KafkaProducer, topic::String, partition::Integer, key, payload)
    if !haskey(p.topics, topic)
        p.topics[topic] = KafkaTopic(p.client, topic, Dict())
    end
    produce1(p.topics[topic], partition, key, payload)
end

## consumer

struct Message{K,P}
    err::Int
    topic::KafkaTopic
    partition::Int32
    key::Union{K, Cvoid}
    payload::Union{P, Cvoid}
    offset::Int64
end

function Message{K,P}(c_msg::CKafkaMessage) where {K,P}
    topic = KafkaTopic(Dict(), "<unkown>", c_msg.rkt)
    if c_msg.err == 0
        key = convert(K, unsafe_load_array(c_msg.key, c_msg.key_len))
        payload = convert(P, unsafe_load_array(c_msg.payload, c_msg.len))
        return Message(Int(c_msg.err), topic, c_msg.partition, key, payload, c_msg.offset)
    else
        return Message{K,P}(Int(c_msg.err), topic, c_msg.partition, nothing, nothing, c_msg.offset)
    end
end

Base.show(io::IO, msg::Message) = print(io, "Message($(msg.key): $(msg.payload))")

mutable struct PartitionList
    rkparlist::Ptr{Cvoid}
end

function PartitionList()
    ptr = kafka_topic_partition_list_new()
    parlist = PartitionList(ptr)
    #finalizer(parlist, parlist -> kafka_topic_partition_list_destroy(ptr))
    return parlist
end

mutable struct KafkaConsumer
    client::KafkaClient
    parlist::PartitionList
end

function KafkaConsumer(conf::Dict)
    @assert haskey(conf, "bootstrap.servers") "`bootstrap.servers` should be specified in conf"
    @assert haskey(conf, "group.id") "`group.id` should be specified in conf"
    client = KafkaClient(KAFKA_TYPE_CONSUMER, conf)
    parlist = PartitionList()
    return KafkaConsumer(client, parlist)
end

function KafkaConsumer(bootstrap_servers::String, group_id::String, conf::Dict=Dict())
    conf["bootstrap.servers"] = bootstrap_servers
    conf["group.id"] = group_id
    return KafkaConsumer(conf)
end

function Base.show(io::IO, c::KafkaConsumer)
    group_id = c.client.conf["group.id"]
    bootstrap_servers = c.client.conf["bootstrap.servers"]
    print(io, "KafkaConsumer($group_id @ $bootstrap_servers)")
end

function subscribe(c::KafkaConsumer, topic_partitions::Vector{Tuple{String, Int}})
    rkparlist = c.parlist.rkparlist
    for (topic, partition) in topic_partitions
        kafka_topic_partition_list_add(rkparlist, topic, partition)
    end
    kafka_subscribe(c.client.rk, rkparlist)
end

function poll(::Type{K}, ::Type{P}, c::KafkaConsumer, timeout::Int=1000) where {K,P}
    c_msg_ptr = kafka_consumer_poll(c.client.rk, timeout)
    if c_msg_ptr != nothing
        c_msg = unsafe_load(c_msg_ptr)
        msg = Message{K,P}(c_msg)
        kafka_message_destroy(c_msg_ptr)
        return msg
    else
        return nothing
    end
end

poll(c::KafkaConsumer, timeout::Int=1000) = poll(Vector{UInt8}, Vector{UInt8}, c, timeout)


text:
# (coeff, featurenum, lag) for each factor and possibly an intercept.
"""
A lagged factor linear regression model for time series prediction. Each factor has
a coefficient, uses a given feature of the data and with a given lag back
in time. By setting a minlag we can ensure that the model can be used that
many time steps in advance, since no factor with less than `minlag` can be
used. The `maxlag` specifies how far back in time we are allowed to look,
since very large lag values can limit the amount of prediction that can be done
in short time series.
"""
mutable struct LaggedFactorLITRM <: AbstractLinearITRM
    nfactors::Int
    nfeatures::Int
    minlag::Int
    maxlag::Int
    maxcoef::Float64
    mincoef::Float64
    hasintercept::Bool
    floorcoefs::Bool
    ndigits::Int # When printing
    params::Union{Nothing,Vector{Float64}}
    featurenames::Vector{String}
    LaggedFactorLITRM(nfactors::Int, nfeatures::Int;
        minlag::Int = 3, maxlag::Int = 10,
        maxcoef = 1e3, mincoef = -maxcoef,
        hasintercept = true, floorcoefs = false, ndigits = 2,
        params = nothing,
        featurenames = String[]) = begin
        new(nfactors, nfeatures, minlag, maxlag, maxcoef, mincoef,
            hasintercept, floorcoefs, ndigits, params, featurenames)
    end
end

function isvalidparams(m::LaggedFactorLITRM, ps::Vector{Float64})
    lagvals = lags(m, ps)
    minl, maxl = extrema(lagvals)
    (minlag(m) <= minl) && (maxl <= maxlag(m))
end

actualminlag(m::LaggedFactorLITRM, ps::Vector{Float64} = params(m)) =
    minimum(lags(m, ps))

actualmaxlag(m::LaggedFactorLITRM, ps::Vector{Float64} = params(m)) =
    maximum(lags(m, ps))

lag(m::LaggedFactorLITRM, fi::Int, ps::Vector{Float64} = params(m)) =
    floor(Int, ps[fi*3])
lag(m::LaggedFactorLITRM, lagval::Float64) = floor(Int, lagval)
lags(m::LaggedFactorLITRM, ps::Vector{Float64} = params(m)) =
    Int[lag(m, ps[fi*3]) for fi in 1:nfactors(m)]

horizon(m::LaggedFactorLITRM) = actualminlag(m)

const AlmostRoundUpTo1 = 0.9999999

function parambounds(m::LaggedFactorLITRM)
    bounds = Tuple{Float64, Float64}[]
    for _ in 1:nfactors(m)
        push!(bounds, (m.mincoef, m.maxcoef))
        # Since we will floor(Int, v) the following we add almost 1.0 to:
        push!(bounds, (1.0, float(m.nfeatures) + AlmostRoundUpTo1))
        push!(bounds, (float(m.minlag), float(m.maxlag) + AlmostRoundUpTo1))
    end
    hasintercept(m) && push!(bounds, (m.mincoef, m.maxcoef))
    bounds
end

function predictforrow(m::LaggedFactorLITRM, X::AbstractMatrix{N}, row::Int,
    ps::Vector{Float64} = params(m)) where {N<:Number}
    prediction = hasintercept(m) ? ps[end] : 0.0
    for fi in 1:nfactors(m)
        coef, featureidx, lag = getfactor(m, fi, ps)
        prediction += (coef * X[row-lag, featureidx])
    end
    return prediction
end

function getfactor(m::LaggedFactorLITRM, fi::Int, ps::Vector{Float64} = params(m))
    i = 1 + (fi-1)*3
    c, f, l = ps[i:(i+2)]
    coef = m.floorcoefs ? floor(Int, c) : c
    return coef, floor(Int, f), floor(Int, l)
end

modelname(m::LaggedFactorLITRM) =
    "LaggedFactor LITRM (Linear Interpretable Timeseries Regression Model)"

function writefactor(io::IO, m::LaggedFactorLITRM, fi::Int)
    c, f, l = getfactor(m, fi)
    coef = r(m, c)
    if coef != 0
        op = (coef > 0.0) ? "+" : "-"
        write(io, " $op $(abs(coef)) * $(featurename(m, f))[T-$(l)]")
    end
end


text:
using Flux

# test model with linear coefficient
struct TestModel
    coeff
end
TestModel(coeff::Number) = TestModel([coeff])
Flux.trainable(m::TestModel) = (m.coeff,)
(m::TestModel)(x) = x .* m.coeff
Flux.@functor TestModel

# test data

function testbatch(batchsize, coeff)
    xs = rand(1.:100., batchsize)
    return (xs, xs .* coeff)
end

function testbatches(n::Int, coeff, batchsize = 8)
    (testbatch(batchsize, coeff) for _ ∈ 1:n)
end

function testlearner(args...; nbatches = 16, coeff = 3, batchsize = 8, kwargs...)
    model = TestModel(rand())
    data = collect(testbatches(nbatches, coeff, batchsize))
    opt = Descent(0.001)
    Learner(
        model,
        (data, data),
        opt,
        Flux.mae,
        args...;
        usedefaultcallbacks = false,
        kwargs...)
end


text:
# PROX_L1    The proximal operator of the l1 norm.
#
#   prox_l1(v,lambda) is the proximal operator of the l1 norm
#   with parameter lambda.
function x = prox_l1(v, lambda)
    x = max(0, v - lambda) - max(0, -v - lambda)
    return x
end


text:
using Sigma
using Distributions

import Sigma: Exponential, uniform, mvexponential

# Data
λreal = 1.5
n = 10
data = rand(Exponential(λreal),n)

λ = uniform(0,2)
x = mvexponential(λ, n)
observations = x == data
posterior_samples = rand(λ, observations, 10)


text:
# Copyright (c) 2016, Jernej Kovacic
# Licensed under the MIT "Expat" License.
#
# See LICENSE.md for more details.

"""
Required precision for iterative algorithms,
it depends on type 't' that must be derived
from AbstractFloat.
"""
function SPECFUN_ITER_TOL(t)

    if ( Float64 == t )
        retVal = 1e-12
    elseif ( Float32 == t )
        retVal = 1f-6
    else   # t == Float16
        retVal = eps(Float16)
    end

    return retVal
end

"""
Maximum allowed number of iterations for
iterative algorithms.
"""
const SPECFUN_MAXITER = 10000


text:
# findnewton.jl
#
#   Faster search for sorted array
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 09.12.2016
# Last Edit: 18.04.2019

" findnewton{T<:Any}(x,target::Array{T,1};max_iterations::Int,verbose::Bool)

In a sorted data set, search iteratively for the target value/s. Uses Newtonian method, look in the center for a match. If value is greater look halfway between end and current value. If less than look between current value and previous low value. Range will be halved each time until the same index or the sentinal value (iteration max) have been reached.

---

#### Keywords:\n
* max_iterations::Int = Maximum number of iterations to attempt, 1000 is default
* verbose::Bool = Display imformation as the function runs, FALSE is default\n\n

---

\n
#### Related Functions:\n
* searchsorted
* searchsortedfirst
* searchsortedlast"
function findnewton(x,target::Array{T,1};max_iterations::Int=1000,verbose::Bool=false) where T <: Any
	answer = Array{Int}(undef,length(target))
	answer[1] = findnewton(x,target[1];max_iterations=max_iterations,verbose=verbose)
	lastgood = Int(1)
	temp = []
	bad_vals = false
	for i=2:1:length(target)
		temp = findnewton(x,target[i];max_iterations=max_iterations,verbose=verbose)
		if temp <= answer[lastgood] + 1
			# If the answer is the same or just one more than the previous value
			# then a gap has been found, mark it for later deletion
			answer[i] = -1
			bad_vals = true
		else
			answer[i] = temp
			lastgood = i
		end
	end

	# Remove Bad Values
	if bad_vals
		answer = answer[findall(answer .!= -1)]
	end

	return answer
end

function findnewton(x,target;max_iterations::Int=1000,verbose::Bool=false)
	if verbose
		println("Min = " * string(findmin(x)[1]))
		println("Max = " * string(findmax(x)[1]))
		println("Target = " * string(target) * "\n")
	end

	############################
	##  Initialize Variables  ##
	############################
	lower = 1
	upper = length(x)
	test = false
	sentinal = 1
	answer = []

	########################
	##  Newtonian Search  ##
	########################
	while !test && sentinal < max_iterations
		mid = Int(lower + floor((upper-lower)/2))
		verbose ? println("  Iteration (" * string(sentinal) * ") = " * string(x[mid])) : []
		if x[mid] > target
			# Target is higher than guess
			if upper == mid
				test = true
				answer = mid
			else
				upper = mid
			end
		elseif x[mid] < target
			# Target is lower than guess
			if lower == mid
				test = true
				answer = mid
			else
				lower = mid
			end
		else
			# Exact Match!
			test = true
			answer = mid
		end

		sentinal += 1
	end

	# Sentinal Catch Case
	if sentinal >= max_iterations
		warn("Sentinal value reached, returned value may not be the best answer.")
		answer = mid
	end

	if verbose
		println("\nIndex = " * string(Int(answer)))
		println("Answer = " * string(x[answer]) * "\n")
	end
	return Int(answer)
end


text:
using Test
using CompScienceMeshes

h = 2π / 51

for T in [Float32, Float64]
    Γ = meshcircle(T(1.0), T(h),2) ## creating the mesh first
    points=vertices(Γ)     ## get the vertices
    num_points=length(points)
    segments=cellarray(Γ)  ## get the faces
    Γ2 = barycentric_refinement(Γ) ## create the refinment of that one
    points2=vertices(Γ2)     ## get the vertices
    segments2=cellarray(Γ2)  ## get the faces

    @test points[1] == points2[1]               # test if the first point in the original and refined mesh are the same
    @test points[num_points]==points2[num_points]# test if the last point in the original has the same index in the refined one
    @test segments[1,1] ==segments2[1,1]        # test if the face are starting with the same point in test
    @test segments[1,2] ==segments2[2,2]        # test if the original second point in the first cell got separated by a single point

    # test barycentric refinement of surface mesh
    local m = meshrectangle(T(1.0), T(1.0), T(0.25), 3)
    f = barycentric_refinement(m)

    @test CompScienceMeshes.refines(f,m)
    @test numcells(f) == 6*numcells(m)

    local m1 = skeleton(m,1)
    f1 = skeleton(f,1)
    @test numcells(f1) == 2*numcells(m1) + 6*numcells(m)

    m0 = skeleton(m,0)
    @test numvertices(f) == numcells(m0) + numcells(m1) + numcells(m)

    ## test bisecting referinment of surfacic meshes
    b = bisecting_refinement(m)
    @test numcells(b) == 4*numcells(m)

    b1 = skeleton(b,1)
    @test numcells(b1) == 2*numcells(m1) + 3*numcells(m)
    @test numvertices(b) == numcells(m0) + numcells(m1)
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
module SIIPExamples

export print_struct
export print_tree

using Literate
using JSON3
using AbstractTrees
using InteractiveUtils
# lightweight type tree printing
AbstractTrees.children(x::Type) = subtypes(x)

repo_directory = dirname(joinpath(@__DIR__))

"""
`print_struct()`

Prints the definition of a struct.
"""
function print_struct(type)
    mutable = type.mutable ? "mutable" : ""
    println("$mutable struct $type")
    for (fn, ft) in zip(fieldnames(type), fieldtypes(type))
        println("    $fn::$ft")
    end
    println("end")
end

function read_json(filename)
    return open(filename) do io
        JSON3.read(io, Dict)
    end
end

abstract type AbstractOS end
abstract type Unix <: AbstractOS end
abstract type BSD <: Unix end

abstract type Windows <: AbstractOS end
abstract type MacOS <: BSD end
abstract type Linux <: BSD end

if Sys.iswindows()
    const os = Windows
elseif Sys.isapple()
    const os = MacOS
else
    const os = Linux
end

abstract type TestData end

"""
Download Data from `branch="master"` name into a "data" folder in given argument path.
Skip the actual download if the folder already exists and force=false.
Defaults to the root of the PowerSystems package.

Returns the downloaded folder name.
"""
function download(
    repo::AbstractString,
    folder::AbstractString = abspath(joinpath(@__DIR__, "..")),
    branch::String = "master",
    force::Bool = false,
)
    if Sys.iswindows()
        DATA_URL = "$repo/archive/$branch.zip"
    else
        DATA_URL = "$repo/archive/$branch.tar.gz"
    end
    directory = abspath(normpath(folder))
    reponame = splitpath(repo)[end]
    data = joinpath(directory, "$reponame-$branch")
    if !isdir(data) || force
        @info "Downloading $DATA_URL"
        tempfilename = Base.download(DATA_URL)
        mkpath(directory)
        @info "Extracting data to $data"
        unzip(os, tempfilename, directory)
        # mv(joinpath(directory, "$reponame-$branch"), data, force = true)
    end

    return data
end

function unzip(::Type{<:BSD}, filename, directory)
    @assert success(`tar -xvf $filename -C $directory`) "Unable to extract $filename to $directory"
end

function unzip(::Type{Windows}, filename, directory)
    path_7z = if Base.VERSION < v"0.7-"
        "$JULIA_HOME/7z"
    else
        sep = Sys.iswindows() ? ";" : ":"
        withenv(
            "PATH" => string(
                joinpath(Sys.BINDIR, "..", "libexec"),
                sep,
                Sys.BINDIR,
                sep,
                ENV["PATH"],
            ),
        ) do
            Sys.which("7z")
        end
    end
    @assert success(`$path_7z x $filename -y -o$directory`) "Unable to extract $filename to $directory"
end

const empty_nb = "{\n \"cells\": [\n  {\n   \"outputs\": [],\n   \"cell_type\": \"markdown\",\n   \"source\": [\n    \"---\\n\",\n    \"\\n\",\n    \"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*\"\n   ],\n   \"metadata\": {}\n  }\n ],\n \"nbformat_minor\": 3,\n \"metadata\": {\n  \"language_info\": {\n   \"file_extension\": \".jl\",\n   \"mimetype\": \"application/julia\",\n   \"name\": \"julia\",\n   \"version\": \"1.3.1\"\n  },\n  \"kernelspec\": {\n   \"name\": \"julia-1.3\",\n   \"display_name\": \"Julia 1.3.1\",\n   \"language\": \"julia\"\n  }\n },\n \"nbformat\": 4\n}\n"
const empty_script = "# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl\n\n"

function rm_if_empty(filepath::String)
    s = open(filepath) do file
        read(file, String)
    end
    if s == empty_nb || s == empty_script
        @info "Removing empty file: $filepath"
        rm(filepath)
    end
end

function postprocess_notebook(nb)
    txt = read(nb, String)
    open(nb, "w") do f
        write(
            f,
            replace(
                txt,
                r"\"outputs\":\ \[\]\,\n\ \ \ \"cell_type\":\ \"markdown\"" =>
                    "\"cell_type\": \"markdown\"",
            ),
        )
    end
end

"""
`literate_file(folder::AbstractString, file::AbstractString)`

Checks if the file has been modified since the last Weave and updates the notebook and tests accordingly.
* `folder` = Name of the folder the tutorial is in
* `file` = Name of the tutorial file
* `force` = foce weave irrespective of file changes
"""
function literate_file(folder, file; force = false, kwargs...)
    filename = splitext(file)[1]
    srcpath = joinpath(repo_directory, "script", folder, file)
    testpath = joinpath(repo_directory, "test", folder)
    notebookpath = joinpath(repo_directory, "notebook", folder)
    notebookfilepath = joinpath(notebookpath, join([filename, ".ipynb"]))
    configpath = joinpath(repo_directory, "script", folder, filename * "_config.json")

    config = get(kwargs, :config, Dict())
    if isfile(configpath)
        @info "found config file for $filename"
        config = read_json(configpath)
    end

    literate = get(config, "literate", true)

    if literate
        make_test = get(config, "test", true)
        make_notebook = get(config, "notebook", true)

        if make_test &&
           (mtime(srcpath) > mtime(testpath) || mtime(testpath) == 0.0 || force)
            @warn "Updating tests for $filename."
            fn = Literate.script(srcpath, testpath; config = config, kwargs...)
            rm_if_empty(fn)
        else
            @warn "Skipping tests for $filename."
        end
        if make_notebook && (
            mtime(srcpath) > mtime(notebookfilepath) ||
            mtime(notebookfilepath) == 0.0 ||
            force
        )
            @warn "Converting $filename to Jupyter Notebook."
            fn = Literate.notebook(srcpath, notebookpath; config = config, kwargs...)
            postprocess_notebook(fn)
            rm_if_empty(fn)
        else
            @warn "Skipping Jupyter Notebook for $filename."
        end
    else
        @warn "Skipping literate for $filename per config"
    end
end

"""
`literate_folder(folder::AbstractString)`

Checks the files present in the specified folder for modifications and updates the notebook and tests accordingly.

* `folder` = Name of the folder to check
"""
function literate_folder(folder; force = false, kwargs...)
    for file in readdir(joinpath(repo_directory, "script", folder))
        if splitext(file)[end] == ".jl"
            println("")
            println("Building $(joinpath(folder, file))")
            try
                literate_file(folder, file, force = force; kwargs...)
            catch
                @error "failed to build $(joinpath(folder, file))"
            end
            println("")
        end
    end
end

"""
`literate_all()`

Checks every tutorial for modifications and updates the notebook accordingly.
"""
function literate_all(; force = false, kwargs...)
    for folder in readdir(joinpath(repo_directory, "script"))
        literate_folder(folder; force = force, kwargs...)
    end
end

end # module


text:
"""
    DiffusionNormal{Tμ,TΣ,Te,TB,TD,TM<:Manifolds.AbstractGroupManifold} <:
    Manifolds.MPointDistribution{TM}

# Constructor

    DiffusionNormal(
        M::Manifolds.AbstractGroupManifold,
        μ,
        Σ,
        e = Manifolds.identity(M, μ),
        basis = Manifolds.DefaultOrthogonalBasis(),
        direction = Manifolds.LeftAction()
    )

Construct a `DiffusionNormal` distribution on a group manifold ``M``, with points
represented with the same type as the identity element ``e``, which may be user-specified.
The mode of the distribution is at ``μ``, and ``Σ`` is a diffusion matrix applied to the
coordinates of the Lie algebra ``𝔪 = Tₑ M`` using the orthogonal `basis`.
"""
struct DiffusionNormal{Tμ,TΣ,Te,TB,TD,TM<:Manifolds.AbstractGroupManifold} <:
       Manifolds.MPointDistribution{TM}
    manifold::TM
    μ::Tμ
    Σ::TΣ
    e::Te
    basis::TB
    direction::TD

    function DiffusionNormal(M, μ, Σ, e, basis, direction)
        pdΣ = _pdmat(M, Σ)
        return new{
            typeof(μ),typeof(pdΣ),typeof(e),typeof(basis),typeof(direction),typeof(M)
        }(
            M, μ, pdΣ, e, basis, direction
        )
    end
end

function DiffusionNormal(
    M,
    μ,
    Σ;
    e=identity(M, μ),
    basis=Manifolds.DefaultOrthogonalBasis(),
    direction=Manifolds.LeftAction(),
)
    return DiffusionNormal(M, μ, Σ, e, basis, direction)
end

const SO3DiffusionNormal{Tμ,TΣ,Te,TB,TD} = DiffusionNormal{
    Tμ,TΣ,Te,TB,TD,SpecialOrthogonal{3}
}
const SE3DiffusionNormal{Tμ,TΣ,Te,TB,TD} = DiffusionNormal{
    Tμ,TΣ,Te,TB,TD,SpecialEuclidean{3}
}

Base.eltype(::Type{<:DiffusionNormal{Tμ,TΣ,Te}}) where {Tμ,TΣ,Te} = Te

Distributions.insupport(d::DiffusionNormal, p) = Manifolds.is_manifold_point(d.manifold, p)

Distributions.mode(d::DiffusionNormal) = d.μ

# TODO: overload convolve for Dirac and DiffusionNormal

function inversion(d::DiffusionNormal)
    M = d.manifold
    μ = inv(M, d.μ)
    direction = Manifolds.switch_direction(d.direction)
    return DiffusionNormal(M, μ, d.Σ, d.e, d.basis, direction)
end

function Distributions._rand!(rng::AbstractRNG, d::DiffusionNormal, q)
    # sample using Euler-Maruyama scheme
    # TODO: use logdetjac of group_exp to tune n so that
    n = 100
    M = d.manifold
    B = d.basis
    e = d.e
    dliealg = Distributions.MvNormal(d.Σ / n)
    Xᵛ = rand(rng, dliealg)
    X = Manifolds.get_vector(M, e, Xᵛ, B)
    p = Manifolds.group_exp!(M, q, X)
    for i in 1:(n - 1)
        rand!(rng, dliealg, Xᵛ)
        Manifolds.get_vector!(M, X, e, Xᵛ, B)
        Manifolds.group_exp!(M, p, X)
        Manifolds.compose!(M, q, q, p)
    end
    Manifolds.translate!(M, q, d.μ, q, d.direction)
    return q
end

function Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:DiffusionNormal})
    d = s.self
    q = Manifolds.allocate_result(d.manifold, rand, d.e)
    return Random.rand!(rng, d, q)
end


text:
using PyPlot
function tplot(p, t)
    # Plot triangular mesh with nodes `p` and triangles `t`
    tris = convert(Array{Int64}, hcat(t...)')
    tripcolor(first.(p), last.(p), tris .- 1, 0*tris[:,1],
              cmap="Set3", edgecolors="k", linewidth=1)
    axis("equal")
    return
end

p = [[0,0], [1,0], [0,1], [1,2], [2,1]]
t = [[1,2,5], [1,5,3], [3,5,4]]
tplot(p,t)

using Pkg; Pkg.add("Triangle") # In case not already installed

using Triangle

function delaunay(p)
    # Delaunay triangulation `t` of array of nodes `p`
    vertex_map = Array{Int64,1}(1:size(p,1))
    t = Triangle.basic_triangulation(convert(Array{Float64}, hcat(p...)'), vertex_map)
end

t = delaunay(p)
tplot(p,t)

p = [ randn(2) for i = 1:50 ]
t = delaunay(p)
tplot(p,t)


text:
# simple kernel_function

function kappa(x1,x2,y1,y2,rho)
    d = sqrt((x1-y1)^2+(x2-y2)^2);
    return (1+d/rho+d^2/(3*rho^2))*exp(-d/rho)
end

function Delta_x_kappa(x1,x2,y1,y2,rho)
    d = sqrt((x1-y1)^2+(x2-y2)^2);
    return (d^2-2*rho*(rho+d))/(3*rho^4)*exp(-d/rho)
end

function Delta_x_y_kappa(x1,x2,y1,y2,rho)
    d = sqrt((x1-y1)^2+(x2-y2)^2);
    return (8*rho^2+d^2-7*rho*d)/(3*rho^6)*exp(-d/rho)
end


text:
"""
    InitialValueProblem{S <: AbstractSystem, XT} <: AbstractSystem

Parametric composite type for initial value problems. It is parameterized in the
system's type and the initial state's type

### Fields

- `s`  -- system
- `x0` -- initial state

### Examples

The linear system ``x' = -x`` with initial condition ``x₀ = [-1/2, 1/2]``:

```jldoctest
julia> s = LinearContinuousSystem([-1.0 0.0; 0.0 -1.0]);

julia> x₀ = [-1/2, 1/2];

julia> p = InitialValueProblem(s, x₀);

julia> initial_state(p) # same as p.x0
2-element Array{Float64,1}:
 -0.5
  0.5

julia> statedim(p)
2

julia> inputdim(p)
0
```
"""
struct InitialValueProblem{S <: AbstractSystem, XT} <: AbstractSystem
    s::S
    x0::XT
end

statedim(ivp::InitialValueProblem) = statedim(ivp.s)
stateset(ivp::InitialValueProblem) = stateset(ivp.s)
inputdim(ivp::InitialValueProblem) = inputdim(ivp.s)
inputset(ivp::InitialValueProblem) = inputset(ivp.s)
islinear(ivp::InitialValueProblem) = islinear(ivp.s)
isaffine(ivp::InitialValueProblem) = isaffine(ivp.s)
ispolynomial(ivp::InitialValueProblem) = ispolynomial(ivp.s)
state_matrix(ivp::InitialValueProblem) = state_matrix(ivp.s)
input_matrix(ivp::InitialValueProblem) = input_matrix(ivp.s)
noise_matrix(ivp::InitialValueProblem) = noise_matrix(ivp.s)
affine_term(ivp::InitialValueProblem) = affine_term(ivp.s)

"""
    initial_state(ivp::InitialValueProblem)

Return the initial state of an initial-value problem.

### Input

- `ivp` -- initial-value problem

### Output

The initial state of an initial-value problem.
"""
initial_state(ivp::InitialValueProblem) = ivp.x0

"""
    system(ivp::InitialValueProblem)

Return the system wrapped by an initial-value problem.

### Input

- `ivp` -- initial-value problem

### Output

The system of the given initial-value problem.
"""
system(ivp::InitialValueProblem) = ivp.s

"""
    IVP

`IVP` is an alias for `InitialValueProblem`.
"""
const IVP = InitialValueProblem