text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Base: @propagate_inbounds

# SIMD-compatible KBN-summation

@inline function kbn_add(a::NTuple{2, Real}, b::Real)
    s = a[1] + b
    c = a[2] + ifelse(
        abs(a[1]) >= abs(b),
        (a[1] - s) + b,
        (b - s) + a[1]
    )
    (s, c)
end

@inline function kbn_add(a::NTuple{2, Real}, b::NTuple{2, Real})
    s = a[1] + b[1]
    c = ifelse(
        abs(a[1]) > abs(b[1]),
        (((a[1] - s) + b[1]) + b[2]) + a[2],
        (((b[1] - s) + a[1]) + a[2]) + b[2]
    )
    (s, c)
end

@doc doc"""
    OnlineMvMean{T<:AbstractFloat} <: AbstractVector{T}

*BAT-internal, not part of stable public API.*

Multivariate mean implemented via Kahan-Babuška-Neumaier summation.
"""
mutable struct OnlineMvMean{T<:AbstractFloat} <: AbstractVector{T}
    m::Int
    sum_w::DoubleFloat{T}
    S::Vector{T}
    C::Vector{T}

    OnlineMvMean{T}(m::Integer) where {T<:AbstractFloat} =
        new{T}(m, zero(DoubleFloat{T}), zeros(T, m), zeros(T, m))
end

OnlineMvMean(m::Integer) = OnlineMvMean{Float64}(m::Integer)

@inline Base.size(omn::OnlineMvMean) = size(omn.S)

@propagate_inbounds function Base.getindex(omn::OnlineMvMean{T}, idxs::Integer...) where {T}
    T((DoubleFloat{T}(omn.S[idxs...]) + DoubleFloat{T}(omn.C[idxs...])) / omn.sum_w)
end

function Base.empty!(target::OnlineMvMean{T}) where T
    target.sum_w = zero(DoubleFloat{T})
    fill!(target.S, zero(T))
    fill!(target.C, zero(T))

    target
end

function Base.merge!(target::OnlineMvMean, others::OnlineMvMean...)
    for x in others
        target.m != x.m && throw(ArgumentError("can't merge OnlineMvMean instances with different size"))
    end
    for x in others
        target.sum_w += x.sum_w
        target_S = target.S; target_C = target.C
        x_S = x.S; x_C = x.C
        @assert eachindex(target_S) == eachindex(x.S) == eachindex(target_C) == eachindex(x.C)  # TODO: Use exception instead of assert
        @inbounds @simd for i in eachindex(target_S)
            target_S[i], target_C[i] = kbn_add(
                (target_S[i], target_C[i]),
                (x.S[i], x.C[i])
            )
        end
    end
    target
end

@inline function Base.push!(
    omn::OnlineMvMean{T}, data::AbstractVector{<:Real}, weight::Real = one(T)
) where {T}
    # Workaround for lack of promotion between, e.g., Float32 and DoubleFloat{Float64}
    weight_conv = T(weight)

    m = omn.m
    S = omn.S
    C = omn.C

    idxs = axes(data, 1)
    @assert idxs == axes(S, 1) == axes(C, 1)  # TODO: Use exception instead of assert

    omn.sum_w += weight_conv

    @inbounds @simd for i in idxs
        x = weight_conv * data[i]
        S[i], C[i] = kbn_add((S[i], C[i]), x)
    end

    omn
end

@doc doc"""
    OnlineMvCov{T<:AbstractFloat,W} <: AbstractMatrix{T}

*BAT-internal, not part of stable public API.*

Implementation based on variance calculation Algorithms of Welford and West.

`W` must either be `Weights` (no bias correction) or one of `AnalyticWeights`,
`FrequencyWeights` or `ProbabilityWeights` to specify the desired bias
correction method.
"""
mutable struct OnlineMvCov{T<:AbstractFloat,W} <: AbstractMatrix{T}
    m::Int
    n::Int64
    sum_w::DoubleFloat{T}
    sum_w2::DoubleFloat{T}
    Mean_X::Vector{T}
    New_Mean_X::Vector{T}
    S::Matrix{T}

    OnlineMvCov{T,W}(m::Integer) where {T<:AbstractFloat,W} =
        new{T,W}(
            m, zero(Int64), zero(DoubleFloat{T}), zero(DoubleFloat{T}),
            zeros(T, m), zeros(T, m), zeros(T, m, m)
        )
end

OnlineMvCov(m::Integer) = OnlineMvCov{Float64, ProbabilityWeights}(m::Integer)

@inline Base.size(ocv::OnlineMvCov) = size(ocv.S)

@propagate_inbounds function Base.getindex(ocv::OnlineMvCov{T, Weights}, idxs::Integer...) where {T}
    sum_w = ocv.sum_w
    ifelse(
        sum_w > zero(T),
        T(ocv.S[idxs...] / sum_w),
        T(NaN)
    )
end

@propagate_inbounds function Base.getindex(ocv::OnlineMvCov{T, AnalyticWeights}, idxs::Integer...) where {T}
    sum_w = ocv.sum_w
    sum_w2 = ocv.sum_w2
    d = sum_w - sum_w2 / sum_w
    ifelse(
        sum_w > zero(T) && d > zero(T),
        T(ocv.S[idxs...] / d),
        T(NaN)
    )
end

@propagate_inbounds function Base.getindex(ocv::OnlineMvCov{T, FrequencyWeights}, idxs::Integer...) where {T}
    sum_w = ocv.sum_w
    ifelse(
        sum_w > one(T),
        T(ocv.S[idxs...] / (sum_w - one(T))),
        T(NaN)
    )
end

@propagate_inbounds function Base.getindex(ocv::OnlineMvCov{T, ProbabilityWeights}, idxs::Integer...) where {T}
    n = ocv.n
    sum_w = ocv.sum_w
    ifelse(
        n > one(T) && sum_w > zero(T),
        T(ocv.S[idxs...] * n / ((n - one(T)) * sum_w)),
        T(NaN)
    )
end

function Base.empty!(target::OnlineMvCov{T,W}) where {T,W}
    target.n = zero(Int64)
    target.sum_w = zero(DoubleFloat{T})
    target.sum_w2 = zero(DoubleFloat{T})
    fill!(target.Mean_X, zero(T))
    fill!(target.New_Mean_X, zero(T))
    fill!(target.S, zero(T))

    target
end

function Base.merge!(target::OnlineMvCov{T,W}, others::OnlineMvCov{T,W}...) where {T,W}
    # TODO: Improve implementation

    for x in others
        target.m != x.m && throw(ArgumentError("can't merge OnlineMvCov instances with different size"))
    end
    m = target.m
    for x in others
        idxs = Base.OneTo(m)
        scl = (target.sum_w * x.sum_w) / (target.sum_w + x.sum_w)
        for j in idxs
            for i in idxs
                target.S[i, j] = scl *
                    (target.Mean_X[i] - x.Mean_X[i]) *
                    (target.Mean_X[j] - x.Mean_X[j]) +
                    target.S[i, j] + x.S[i, j]
            end
        end

        d = inv(T(target.sum_w + x.sum_w))
        for i in idxs
            target.Mean_X[i] = (target.Mean_X[i] * target.sum_w + x.Mean_X[i] * x.sum_w) * DoubleFloat(d)
            target.New_Mean_X[i] = target.Mean_X[i]
        end

        target.n += x.n
        target.sum_w += x.sum_w
        target.sum_w2 += x.sum_w2
    end
    target
end

@inline function Base.push!(
    ocv::OnlineMvCov{T,W}, data::AbstractVector{<:Real}, weight::Real = one(T)
) where {T,W}
    # Ignore zero weights (can't be handled)
    if weight ≈ 0
        return ocv
    end

    # Workaround for lack of promotion between, e.g., Float32 and DoubleFloat{Float64}
    weight_conv = T(weight)

    m = ocv.m
    n = ocv.n
    sum_w = ocv.sum_w
    sum_w2 = ocv.sum_w2
    Mean_X = ocv.Mean_X
    New_Mean_X = ocv.New_Mean_X
    S = ocv.S

    idxs = axes(data, 1)
    @assert idxs == axes(Mean_X, 1) == axes(New_Mean_X, 1) == axes(S, 1) == axes(S, 2)  # TODO: Use exception instead of assert

    n += one(n)
    sum_w += weight_conv
    sum_w2 += weight_conv^2

    weight_over_sum_w = T(weight_conv / sum_w)

    @inbounds @simd for i in idxs
        x = data[i]
        mean_X = Mean_X[i]
        New_Mean_X[i] = muladd(x - mean_X, weight_over_sum_w, mean_X)
    end

    @inbounds for j in idxs
        new_dx_j = data[j] - New_Mean_X[j]

        j_offs = (LinearIndices(S))[last(idxs), j] - last(idxs)
        #@assert (LinearIndices((S))[last(idxs), j] == last(idxs) + j_offs  # TODO: Use exception instead of assert
        @simd for i in idxs
            dx_i = data[i] - Mean_X[i]
            S[i + j_offs] = muladd(dx_i, weight_conv * new_dx_j, S[i + j_offs])
        end
    end

    @inbounds @simd for i in idxs
        Mean_X[i] = New_Mean_X[i]
    end

    ocv.n = n
    ocv.sum_w = sum_w
    ocv.sum_w2 = sum_w2

    ocv
end

@doc doc"""
    BasicMvStatistics{T<:Real,W}

*BAT-internal, not part of stable public API.*

`W` must either be `Weights` (no bias correction) or one of `AnalyticWeights`,
`FrequencyWeights` or `ProbabilityWeights` to specify the desired bias
correction method.
"""
mutable struct BasicMvStatistics{T<:Real,W}
    m::Int
    mean::OnlineMvMean{T}
    cov::OnlineMvCov{T,W}
    maximum::Vector{T}
    minimum::Vector{T}

    BasicMvStatistics{T,W}(
        mean::OnlineMvMean{T},
        cov::OnlineMvCov{T,W},
        maximum::Vector{T},
        minimum::Vector{T}
    ) where {T<:Real,W} = new(mean, cov, maximum, minimum)

    BasicMvStatistics{T,W}(m::Integer) where {T<:Real,W} =
        new(m, OnlineMvMean{T}(m), OnlineMvCov{T,W}(m), fill(typemin(T), m), fill(typemax(T), m))
end

function Base.empty!(target::BasicMvStatistics{T}) where T
    m = target.m
    empty!(target.mean)
    empty!(target.cov)
    fill!(target.maximum, typemin(T))
    fill!(target.minimum, typemax(T))

    target
end

function Base.merge!(target::BasicMvStatistics, others::BasicMvStatistics...)
    for x in others
        merge!(target.mean, x.mean)
        merge!(target.cov, x.cov)
        target.maximum = max.(target.maximum, x.maximum)
        target.minimum = min.(target.minimum, x.minimum)
    end
    target
end

function Base.push!(
    stats::BasicMvStatistics{T,W}, data::AbstractVector{<:Real}, weight::Real = one(T)
) where {T,W}
    push!(stats.mean, data, weight)
    push!(stats.cov, data, weight)

    max_v = stats.maximum
    min_v = stats.minimum

    idxs = axes(data, 1)
    @assert idxs == axes(max_v, 1) == axes(min_v, 1)  # TODO: Use exception instead of assert

    @inbounds @simd for i in idxs
        x = data[i]
        max_v[i] = max(max_v[i], x)
        min_v[i] = min(min_v[i], x)
    end

    stats
end

const OnlineMvStatistic = Union{OnlineMvMean, OnlineMvCov, BasicMvStatistics}

function Base.append!(target::OnlineMvStatistic, data::VectorOfSimilarVectors)
    @uviews data push!.(Ref(target), data)
    target
end

function Base.append!(target::OnlineMvStatistic, data::VectorOfSimilarVectors, weights::AbstractVector)
    @uviews data weights push!.(Ref(target), data, weights)
    target
end

Base.merge(x::S, others::S...) where {S <: OnlineMvStatistic} = merge!(deepcopy(x), others...)


text:
@kwdef mutable struct SDivHandler <: AbstractHandler
    id = 0
end

@register(SDivHandler, /, begin
    [SBit, SBit] => [SBit]
    [SBit, Vector{SBit}] => [Vector{SBit}]
    [Vector{SBit}, SBit] => [Vector{SBit}]
    [Vector{SBit}, Vector{SBit}] => [Vector{SBit}]
    [SBit, Matrix{SBit}] => [Matrix{SBit}]
    [Matrix{SBit}, SBit] => [Matrix{SBit}]
    [Matrix{SBit}, Matrix{SBit}] => [Matrix{SBit}]
end)

function (handler::SDivHandler)(netlist::Netlist,
                                inputs::Vector{Variable},
                                outputs::Vector{Variable})
    # update netlist with inputs
    setsigned!(netlist, getname(inputs[1]), true)
    setsigned!(netlist, getname(inputs[2]), true)

    # compute output size
    lname, rname, outsize = handlebroadcast(inputs[1].name, inputs[2].name,
                                            getsize(netlist, getname(inputs[1])),
                                            getsize(netlist, getname(inputs[2])))

    # update netlist with output
    setsigned!(netlist, getname(outputs[1]), true)

    # add internal nets to netlist
    update!(netlist, Net(name = "div$(handler.id)_pp", size = outsize))
    update!(netlist, Net(name = "div$(handler.id)_mp", size = outsize))

    outstring = """
        $stdcomment
        // BEGIN div$(handler.id)
        stoch_div_mat #(
                .COUNTER_SIZE(8),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(handler.id)_pp (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_p"))),
                .B($(rname("_p"))),
                .Y(div$(handler.id)_pp)
            );
        stoch_div_mat #(
                .COUNTER_SIZE(8),
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(handler.id)_mp (
                .CLK(CLK),
                .nRST(nRST),
                .A($(lname("_m"))),
                .B($(rname("_p"))),
                .Y(div$(handler.id)_mp)
            );
        stoch_sat_sub_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(handler.id)_p (
                .CLK(CLK),
                .nRST(nRST),
                .A(div$(handler.id)_pp),
                .B(div$(handler.id)_mp),
                .Y($(outputs[1].name)_p)
            );
        stoch_sat_sub_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) div$(handler.id)_m (
                .CLK(CLK),
                .nRST(nRST),
                .A(div$(handler.id)_mp),
                .B(div$(handler.id)_pp),
                .Y($(outputs[1].name)_m)
            );
        // END div$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

ext: jl


text:
__precompile__()

"""
    TextWrap

This module provides the function [`wrap`](@ref) which parses an input text
and reorganizes its white space so that it can be printed with a fixed screen width, optionally indenting it.
It also provides the two convenience functions [`print_wrapped`](@ref) and [`println_wrapped`](@ref).
"""
module TextWrap

using Compat

if VERSION < v"0.7.0-DEV.2949"
    # horrible, but we don't need to reimplement all of the logic here
    # TODO: remove it if Compat implements this...
    Base.nextind(s, i, j) = (@assert i == 0; chr2ind(s, j))
end

export
    wrap,
    print_wrapped,
    println_wrapped

# A regex to match any sequence of spaces except non-breakable spaces (\xA0)
const spaceregex = r"((?!\xA0)\s)+"

function _expand_tabs(text::AbstractString, i0::Int)
    out_buf = IOBuffer()
    i = i0 % 8
    for c in text
        if c == '\t'
            while i < 8
                print(out_buf, ' ')
                i += 1
            end
            i = 0
        elseif c == '\n' || c == '\r'
            print(out_buf, c)
            i = 0
        else
            print(out_buf, c)
            i = (i + 1) % 8
        end
    end
    return String(take!(out_buf))
end

function _check_width(width::Integer)
    if width <= 0
        error("invalid width $width (must be > 0)")
    end
    return true
end
function _check_indent(indent::Integer, width::Integer)
    if indent < 0 || indent >= width
        error("invalid intent $indent (must be an integer between 0 and width-1, or an AbstractString)")
    end
    return true
end
function _check_indent(indent::AbstractString, width::Integer)
    if length(indent) >= width
        error("invalid intent (must be shorter than width-1)")
    end
end

function _put_chunks(chunk::AbstractString, out_str,
                    cln, cll, bol, soh,
                    width, initial_indent, subsequent_indent,
                    break_on_hyphens, break_long_words)

    # This function just performs breaks-on-hyphens and passes
    # individual chunks to _put_chunk

    _hyphen_re = r"^(\w+-)*?\w+[^0-9\W]+-(?=\w+[^0-9\W])"

    while break_on_hyphens
        m = match(_hyphen_re, chunk)
        if m === nothing
            break
        end
        c = m.match
        cln, cll, bol, lcise = _put_chunk(c, out_str,
                    cln, cll, bol, soh,
                    width, initial_indent, subsequent_indent,
                    break_long_words)
        soh = ""
        chunk = chunk[m.offset+lastindex(c):end]
    end

    cln, cll, bol, lcise = _put_chunk(chunk, out_str,
                cln, cll, bol, soh,
                width, initial_indent, subsequent_indent,
                break_long_words)
    return cln, cll, bol, lcise
end

function _put_chunk(chunk::AbstractString, out_str,
                    cln, cll, bol, soh,
                    width, initial_indent, subsequent_indent,
                    break_long_words)

    # Writes a chunk to out_str, based on the current position
    # as encoded in (cln, cll, bol) = (current_line_number,
    # current_line_length, beginning_of_line), and returns the
    # updated position (plus a flag to signal that an end-of-sentence
    # was detected).
    # The argument soh (=space_on_hold) is the spacing which should
    # go in front of chunk, and it may or may not be printed.
    # The rest are options.

    liindent = length(initial_indent)
    lsindent = length(subsequent_indent)
    lchunk = length(chunk)
    lsoh = length(soh)

    if cll + lsoh > width
        soh = ""
        lsoh = 0
        print(out_str, "\n")
        cln += 1
        cll = 0
        bol = true
    end

    if bol
        if cln == 1
            indent = initial_indent
            lindent = liindent
        else
            indent = subsequent_indent
            lindent = lsindent
        end

        print(out_str, indent)
        cll += lindent

        soh = ""
        lsoh = 0
    end

    # is there enough room for the chunk? or is this the
    # beginning of the text and we cannot break words?
    if cll + lsoh + lchunk <= width ||
            (cln == 1 && bol && !break_long_words)
        print(out_str, soh, chunk)
        cll += lchunk + lsoh
        bol = false
    # does the chunk fit into next line? or are we
    # forced to put it there?
    elseif lchunk <= width - lsindent || !break_long_words
        print(out_str, bol ? "" : "\n", subsequent_indent, chunk)
        cll = lsindent + lchunk
        cln += 1
        bol = false
    # break it until it fits
    else
        while cll + lsoh + lchunk > width
            if width - cll - lsoh > 0
                print(out_str, soh, chunk[1:nextind(chunk, 0, width-cll-lsoh)], "\n",
                        subsequent_indent)
                chunk = chunk[nextind(chunk, 0, width-cll-lsoh+1):end]
                lchunk = length(chunk)
            else
                print(out_str, "\n", subsequent_indent)
            end
            cll = lsindent
            cln += 1
            soh = ""
            lsoh = 0
        end
        print(out_str, chunk)
        cll += lchunk
        bol = false
    end

    # detect end-of-sentences
    _sentence_end_re = r"\w([\.\!\?…]|\.\.\.)[\"\'´„]?\Z"
    lcise = occursin(_sentence_end_re, chunk)

    return cln, cll, bol, lcise
end

"""
    wrap(string; keywords...)

Parses `string` and returns a new string in which newlines are inserted as appropriate in order
for each line to fit within a specified width.

The behaviour can be controlled via optional keyword arguments:

* `width` (deafult=`70`): the maximum width of the wrapped text, including indentation.
* `initial_indent` (default=`""`): indentation of the first line. This can
   be any string (shorter than `width`), or it can be an integer number (lower than `width`).
* `subsequent_indent` (default=`""`): indentation of all lines except the first. Works the same as
   `initial_indent`.
* `break_on_hyphens` (default=`true`): this flag determines whether words can be broken on hyphens,
  e.g. whether "high-precision" can be split into "high-" and "precision".
* `break_long_words` (default=`true`): this flag determines what to do when a word is too long to
  fit in any line. If `true`, the word will be broken, otherwise it will go beyond the desired text
  width.
* `replace_whitespace` (default=`true`): if this flag is `true`, all whitespace characters in the
  original text (including newlines) will be replaced by spaces.
* `expand_tabs` (default=`true`): if this flag is `true`, tabs will be expanded in-place into spaces.
  The expansion happens before whitespace replacement.
* `fix_sentence_endings` (default=`false`): if this flag is `true`, the wrapper will try to recognize
  sentence endings in the middle of a paragraph and put two spaces before the next sentence in case
  only one is present.
"""
function wrap(text::AbstractString;
              width::Int = 70,
              initial_indent::Union{Integer,AbstractString} = "",
              subsequent_indent::Union{Integer,AbstractString} = "",
              expand_tabs::Bool = true,
              replace_whitespace::Bool = true,
              fix_sentence_endings::Bool = false,
              break_long_words::Bool = true,
              break_on_hyphens::Bool = true)

    # Reformat the single paragraph in `text` so it fits in lines of
    # no more than `width` columns, and return an AbstractString.

    # Sanity checks
    _check_width(width)
    _check_indent(initial_indent, width)
    _check_indent(subsequent_indent, width)

    if isa(initial_indent, Integer)
        initial_indent = " "^initial_indent
    end
    if isa(subsequent_indent, Integer)
        subsequent_indent = " "^subsequent_indent
    end

    # State variables initialization
    cln = 1 # current line number
    cll = 0 # current line length
    bol = true # beginning of line
    lcise = false # last chunk is sentence ending
    soh = "" # space on hold

    # We iterate over the text, looking for whitespace
    # where to split.
    i = firstindex(text)
    l = lastindex(text)
    out_str = IOBuffer()

    wsrng = findnext(spaceregex, text, i)
    j, k = wsrng ≢ nothing ?
        (first(wsrng), nextind(text, last(wsrng))) :
        (0, -1)
    while 0 < j <= l
        if i < k
            if i < j
                # This is non-whitespace. We write it out according
                # to the current cursor position and the leading space.
                chunk = text[i:prevind(text,j)]

                cln, cll, bol, lcise = _put_chunks(chunk, out_str,
                            cln, cll, bol, soh,
                            width, initial_indent, subsequent_indent,
                            break_on_hyphens, break_long_words)
            end
            i = k
        end

        # This is whitespace. We mangle it (expand tabs, fix
        # sentence endings, replace it with single spaces) and
        # then we keep it on hold.
        soh = text[j:prevind(text,k)]
        if expand_tabs && occursin(r"\t", soh)
            soh = _expand_tabs(soh, cll)
        end
        if fix_sentence_endings && lcise && soh == " "
            soh = "  "
        end
        if replace_whitespace
            soh = " "^length(soh)
        end

        # Continue the search

        k <= j && (k = nextind(text,j))
        wsrng = findnext(spaceregex, text, k)
        j, k = wsrng ≢ nothing ?
            (first(wsrng), nextind(text, last(wsrng))) :
            (0, -1)
    end
    if i ≤ ncodeunits(text)
        # Some non-whitespace is left at the end.
        chunk = text[i:end]
        cln, cll, bol = _put_chunks(chunk, out_str,
                    cln, cll, bol, soh,
                    width, initial_indent, subsequent_indent,
                    break_on_hyphens, break_long_words)
    end
    return String(take!(out_str))
end

# print functions signature:
#   first arg: IO
#   last arg: Options
#   inbetween: anything printable
#
#   all arguments are optional
#
function _print_wrapped(newline::Bool, args...; kwargs...)
    if !isempty(args) && isa(args[1], IO)
        io = args[1]
        args = args[2:end]
    else
        io = stdout
    end

    if !isempty(args)
        ws = wrap(string(args...); kwargs...)
    else
        ws = ""
    end

    if newline
        println(io, ws)
    else
        print(io, ws)
    end
end

"""
    print_wrapped([io,] text...; keywords...)

This is just like the standard `print` function (it prints multiple arguments and accepts
an optional `IO` first argument), except that it wraps the result, and accepts keyword
arguments to pass to [`wrap`](@ref).
"""
print_wrapped(args...; kwargs...) = _print_wrapped(false, args...; kwargs...)

"""
    println_wrapped([io,] text...; keywords...)

Like [`print_wrapped`](@ref), but adds a newline at the end.
"""
println_wrapped(args...; kwargs...) = _print_wrapped(true, args...; kwargs...)

end # module TextWrap


text:
function helper_fn(z, attempt)
	single_data_name = z
	proc_id = myid();

	helper_start_time = time()

	train_transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_transaction_counts.csv") ) )
	train_assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_assortments.csv") ) )

	test_transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_R2_temporally_split/week3_to_week6/", single_data_name, "_transaction_counts.csv") ) )
	test_assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_R2_temporally_split/week3_to_week6/", single_data_name, "_assortments.csv") ) )

	num_train_transactions = sum(train_transaction_counts)
	num_test_transactions = sum(test_transaction_counts)

	num_train_assortments = size(train_assortments, 1)
	num_test_assortments = size(test_assortments, 1)

	N = 10

	time_limit = 1200

	perTol = 1e-6

	## RANKING-BASED MODEL ##
	CS_size_vec = [2,3,4,5,6,7,8,9]

	for CS_size in CS_size_vec

		ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking_CS(N, train_assortments, train_transaction_counts, time_limit, perTol, CS_size)
		ranking_K = length(ranking_lambda)

		test_ranking_predict, test_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, test_assortments)
		ranking_KL, ranking_noinf_count, ranking_KL_avg = tcm_evaluateKL(N, test_ranking_predict, test_transaction_counts)

		avg_CS_size, weighted_avg_CS_size, max_CS_size = tcm_rankingStats(N, ranking_lambda, orderings)

		row_string = string(single_data_name,  ",",
								num_train_assortments, ",",
								num_test_assortments, ",",
								num_train_transactions, ",",
								num_test_transactions, ",",
								proc_id, ",",
								attempt, ",",
								"ranking", ",",
								CS_size, ",",
								ranking_loglik, ",",
								ranking_K, ",",
								perTol, ",",
								ranking_KL, ",",
								ranking_noinf_count, ",",
								ranking_KL_avg, ",",
								ranking_elapsed_time, ",",
								avg_CS_size, ",",
								weighted_avg_CS_size, ",",
								max_CS_size)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(ranking_outcsvhandle, row_string);
			flush(ranking_outcsvhandle);
		end
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
include("main.jl")

REWRITE_DATA = true

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# Test dependency on variance estimation confidence
casename = "out_of_sample_lindist"
case_path = "results/$(casename)"

if REWRITE_DATA && isdir(case_path)
    rm(case_path, recursive=true)
end

if !(isdir(case_path)) || REWRITE_DATA
    # Only run if there are no saved results
    mkdir(case_path)

    sample_size = 750
    # Create global set of true observations to make tests comparable
    OBSERVATIONS = create_data(ERROR_VARIANCES_TRUE, sample_size)

    objective, bus_results, line_results = Solve_LinDist(BUSES, LINES, GENERATORS)
    det_alpha = vcat([1], zeros(length(BUSES)-1))
    # In Deterministic Case all imbalance is controlled by the substation
    bus_results = hcat(bus_results, DataFrame(α=det_alpha))
    ldtest = test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC", print_on=false)

    ldtest_dict = DataFrame()

    ldtest_dict = vcat(ldtest_dict, ldtest)

    writetable("$(case_path)/ldtest_dict.csv", ldtest_dict)

else
    # If there is saved data load it
    info("Data already exists at $case_path")
    cctest_dict = readtable("$(case_path)/cctest_dict.csv")
    ddtest_dict = readtable("$(case_path)/ddtest_dict.csv")
end


text:
using Test

@testset "GaN-SNAP-LAMMPS" begin

    params = get_conf_params("../examples/GaN-SNAP-LAMMPS/")

    # Get DFT data
    dft_train_data, dft_val_data = generate_data("dft", params)

    # Get reference data
    ref_train_data, ref_val_data = generate_data("ref", params)

    # Get potential learning problem (e.g. A β = b)
    snap = learning_problem(dft_train_data, ref_train_data, params)

    # Solve potential learning problem (e.g. β = A \ b)
    learn(snap, params)

    # Validate potentials, forces, and stresses
    m = error_metrics(snap, dft_val_data - ref_val_data, params)

    @test m["energy"]["max_rel_error"] < 0.1

end


text:
function groundstate(H::LocalHamiltonian, Ψ₀::FiniteCMPS;
                        optalg::OptimKit.OptimizationAlgorithm =
                            LBFGS(30; verbosity = 2, gradtol = 1e-3),
                        boundaries = :fixed,
                        left_boundary = boundaries,
                        right_boundary = boundaries,
                        optimnodes = nodes(Ψ₀.Q),
                        Kmax = 60,
                        tol = eps(),
                        callback = (numiter, E, normgrad, Ψ, ρL, ρR, e) -> nothing)

    if !(left_boundary ∈ (:free, :fixed) && right_boundary ∈ (:free, :fixed))
         throw(ArgumentError("unknown boundary conditions, expected values are `:free` or `:fixed`"))
    end

    envnodes = nodes(Ψ₀.Q)
    optimindices::Vector{Int} = indexin(optimnodes, envnodes)
    @assert !any(isnothing, optimindices)
    @assert optimindices[1] == 1
    @assert optimindices[end] == length(envnodes)

    function fg(x)
        Ψ, HL, HR, ρL, ρR, hL, hR, E, e = x
        g = gradient(H, (Ψ, ρL, ρR), HL, HR; gradindices = optimindices,
                    left_boundary = left_boundary, right_boundary = right_boundary)
        return E, g
    end
    inner(x, g1, g2) = 2*real(dot(g1, g2))
    function retract(x, g, α)
        Ψ, = x
        Qold, Rsold, vLold, vRold = Ψ
        dQ, dRs, dvL, dvR = g
        Q = Qold + α * dQ
        Rs = Rsold .+ α .* dRs
        if left_boundary == :free
            ndvL = norm(dvL)
            sα, cα = sincos(α*ndvL)
            vL = normalize!(axpby!(sα/ndvL, dvL, cα, copy(vLold)))
            dvL = axpby!(-sα*ndvL, vLold, cα, copy(dvL))
        else
            vL = vLold
        end
        if right_boundary == :free
            ndvR = norm(dvR)
            sα, cα = sincos(α*ndvR)
            vR = normalize!(axpby!(sα/ndvR, dvR, cα, copy(vRold)))
            dvR = axpby!(-sα*ndvR, vRold, cα, copy(dvR))
        else
            vR = vRold
        end
        Ψ = FiniteCMPS(Q, Rs, vL, vR)
        HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR =
            environments!(H, Ψ; Kmax = Kmax, tol = tol)
        infoHL || @warn "HL not converged"
        infoHR || @warn "HR not converged"
        infoρL || @warn "ρL not converged"
        infoρR || @warn "ρR not converged"
        g = FiniteCMPSTangent(Ψ, dQ, dRs, dvL, dvR, optimindices)
        return (Ψ, HL, HR, ρL, ρR, hL, hR, E, e), g
    end
    function transport!(ξ, x, η, α, x′)
        Ψ, = x
        Ψ′, = x′
        Qold, Rolds, vLold, vRold = Ψ
        dQ, dRs, dvL, dvR = ξ
        _, _, dvL′, dvR′ = η
        if left_boundary == :free
            ndvL′ = norm(dvL′)
            sα, cα = sincos(α*ndvL′)
            d1 = dot(vLold, dvL)
            d2 = dot(dvL′, dvL)/ndvL′
            dvL = axpy!((cα-1)*d1 - sα*d2, vLold, dvL)
            dvL = axpy!(((cα-1)*d2 + sα*d1)/ndvL′, dvL′, dvL)
        end
        if right_boundary == :free
            ndvR′ = norm(dvR′)
            sα, cα = sincos(α*ndvR′)
            d1 = dot(vRold, dvR)
            d2 = dot(dvR′, dvR)/ndvR′
            dvR = axpy!((cα-1)*d1 - sα*d2, vRold, dvR)
            dvR = axpy!(((cα-1)*d2 + sα*d1)/ndvR′, dvR′, dvR)
        end
        return FiniteCMPSTangent(Ψ′, dQ, dRs, dvL, dvR, ξ.indices)
    end
    function finalize!(x, f, g, numiter)
        Ψ, HL, HR, ρL, ρR, hL, hR, E, e = x
        normgrad = sqrt(inner(x, g, g))
        callback(numiter, E, normgrad, Ψ, ρL, ρR, e)
        return x, f, g
    end

    x = let (Q₀, R₀s, vL₀, vR₀) = Ψ₀
        Q = PiecewiseLinear(collect(Q₀.nodes), copy.(Q₀.values))
        Rs = map(R₀s) do R₀
            PiecewiseLinear(collect(R₀.nodes), copy.(R₀.values))
        end
        vL = vL₀/norm(vL₀)
        vR = vR₀/norm(vR₀)
        Ψ = FiniteCMPS(Q, Rs, vL, vR)
        HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR =
            environments!(H, Ψ; Kmax = Kmax, tol = tol)
        (Ψ, HL, HR, ρL, ρR, hL, hR, E, e)
    end
    x, _, grad, numfg, history = optimize(fg, x, optalg; retract = retract, inner = inner,
                                            transport! = transport!, finalize! = finalize!)

    return let (Ψ, HL, HR, ρL, ρR, hL, hR, E, e) = x
        normgrad = sqrt(inner(x, grad, grad))
        Ψ, ρL, ρR, E, e, normgrad, numfg, history
    end
end

function groundstate2(H::LocalHamiltonian, Ψ₀::FiniteCMPS;
                        optalg::OptimKit.OptimizationAlgorithm =
                            LBFGS(30; verbosity = 2, gradtol = 1e-3),
                        boundaries = :fixed,
                        left_boundary = boundaries,
                        right_boundary = boundaries,
                        optimnodes = nodes(Ψ₀.Q),
                        Kmax = 60,
                        tol = eps(),
                        callback = (numiter, E, normgrad, Ψ, ρL, ρR, e) -> nothing)

    if !(left_boundary ∈ (:free, :fixed) && right_boundary ∈ (:free, :fixed))
         throw(ArgumentError("unknown boundary conditions, expected values are `:free` or `:fixed`"))
    end

    envnodes = nodes(Ψ₀.Q)
    optimindices::Vector{Int} = indexin(optimnodes, envnodes)
    @assert !any(isnothing, optimindices)
    @assert optimindices[1] == 1
    @assert optimindices[end] == length(envnodes)

    function fg(x)
        Ψ, HL, HR, ρL, ρR, hL, hR, E, e = x
        g = gradient(H, (Ψ, ρL, ρR), HL, HR; gradindices = optimindices,
                    left_boundary = left_boundary, right_boundary = right_boundary)
        return E, g
    end
    inner(x, g1, g2) = 2*real(dot(g1, g2))
    function retract(x, g, α)
        Ψ, = x
        Qold, Rsold, vLold, vRold = Ψ
        dQ, dRs, dvL, dvR = g
        Q = Qold + α * dQ
        Rs = Rsold .+ α .* dRs
        if left_boundary == :free
            ndvL = norm(dvL)
            sα, cα = sincos(α*ndvL)
            vL = normalize!(axpby!(sα/ndvL, dvL, cα, copy(vLold)))
            dvL = axpby!(-sα*ndvL, vLold, cα, copy(dvL))
        else
            vL = vLold
        end
        if right_boundary == :free
            ndvR = norm(dvR)
            sα, cα = sincos(α*ndvR)
            vR = normalize!(axpby!(sα/ndvR, dvR, cα, copy(vRold)))
            dvR = axpby!(-sα*ndvR, vRold, cα, copy(dvR))
        else
            vR = vRold
        end
        Ψ = FiniteCMPS(Q, Rs, vL, vR)
        HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR =
            environments!(H, Ψ; Kmax = Kmax, tol = tol)
        infoHL || @warn "HL not converged"
        infoHR || @warn "HR not converged"
        infoρL || @warn "ρL not converged"
        infoρR || @warn "ρR not converged"
        g = FiniteCMPSTangent(Ψ, dQ, dRs, dvL, dvR, optimindices)
        return (Ψ, HL, HR, ρL, ρR, hL, hR, E, e), g
    end
    function transport!(ξ, x, η, α, x′)
        Ψ, = x
        Ψ′, = x′
        Qold, Rolds, vLold, vRold = Ψ
        dQ, dRs, dvL, dvR = ξ
        _, _, dvL′, dvR′ = η
        if left_boundary == :free
            ndvL′ = norm(dvL′)
            sα, cα = sincos(α*ndvL′)
            d1 = dot(vLold, dvL)
            d2 = dot(dvL′, dvL)/ndvL′
            dvL = axpy!((cα-1)*d1 - sα*d2, vLold, dvL)
            dvL = axpy!(((cα-1)*d2 + sα*d1)/ndvL′, dvL′, dvL)
        end
        if right_boundary == :free
            ndvR′ = norm(dvR′)
            sα, cα = sincos(α*ndvR′)
            d1 = dot(vRold, dvR)
            d2 = dot(dvR′, dvR)/ndvR′
            dvR = axpy!((cα-1)*d1 - sα*d2, vRold, dvR)
            dvR = axpy!(((cα-1)*d2 + sα*d1)/ndvR′, dvR′, dvR)
        end
        return FiniteCMPSTangent(Ψ′, dQ, dRs, dvL, dvR, ξ.indices)
    end
    function finalize!(x, f, g, numiter)
        Ψ, HL, HR, ρL, ρR, hL, hR, E, e = x
        normgrad = sqrt(inner(x, g, g))
        callback(numiter, E, normgrad, Ψ, ρL, ρR, e)
        return x, f, g
    end
    function precondition(x, g)
        Ψ, HL, HR, ρL, ρR, hL, hR, E, e = x
        g′, info = linsolve(g, g, CG(; verbosity = 0, maxiter = 5, tol = 1e-9)) do Φ
            dQ, dRs = CMPSKit.metric(Φ, (Ψ, ρL, ρR))
            axpy!(1e-2, Φ, FiniteCMPSTangent(Ψ, dQ, dRs, zero(Ψ.vL), zero(Ψ.vR)))
        end
        return g′
    end

    x = let (Q₀, R₀s, vL₀, vR₀) = Ψ₀
        Q = PiecewiseLinear(collect(Q₀.nodes), copy.(Q₀.values))
        Rs = map(R₀s) do R₀
            PiecewiseLinear(collect(R₀.nodes), copy.(R₀.values))
        end
        vL = normalize(vL₀)
        vR = normalize(vR₀)
        Ψ = FiniteCMPS(Q, Rs, vL, vR)
        HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR =
            environments!(H, Ψ; Kmax = Kmax, tol = tol)
        (Ψ, HL, HR, ρL, ρR, hL, hR, E, e)
    end
    x, _, grad, numfg, history = optimize(fg, x, optalg; retract = retract, inner = inner,
                                            transport! = transport!, finalize! = finalize!,
                                            precondition = precondition)

    return let (Ψ, HL, HR, ρL, ρR, hL, hR, E, e) = x
        normgrad = sqrt(inner(x, grad, grad))
        Ψ, ρL, ρR, E, e, normgrad, numfg, history
    end
end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE, ACE.SphericalHarmonics, StaticArrays, LinearAlgebra
using ACE: _mrange

module RotationCoeffs

   using ACE, ACE.SphericalHarmonics, StaticArrays, LinearAlgebra
   using ACE: _mrange

   struct Rot3DCoeffs{NMAX}
      vals::SVector{NMAX, Dict}
      cg::ClebschGordan{Float64}
   end

   dicttype(N::Integer) = dicttype(Val(N))
   dicttype(::Val{N}) where {N} =
         Dict{Tuple{SVector{N,Int8}, SVector{N,Int8}, SVector{N,Int8}}, Float64}

   function Rot3DCoeffs(Nmax, Lmax)
      vals = SVector{Nmax, Dict}([ dicttype(N)() for N = 1:Nmax ]...)
      cg = ClebschGordan(Lmax)
      return Rot3DCoeffs(vals, cg)
   end

   function get_vals(A::Rot3DCoeffs, ::Val{N}) where {N}
      return A.vals[N]::Dict{Tuple{SVector{N,Int8}, SVector{N,Int8}, SVector{N,Int8}}, Float64}
   end

   _key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
         (SVector{N, Int8}(ll), SVector{N, Int8}(mm), SVector{N, Int8}(kk))

   function (A::Rot3DCoeffs)(ll::StaticVector{N},
                            mm::StaticVector{N},
                            kk::StaticVector{N}) where {N}
      if       sum(mm) != 0 ||
               sum(kk) != 0 ||
               !all(abs.(mm) .<=  ll) ||
               !all(abs.(kk) .<= ll)
         return 0.0
      end
      vals = A.vals[N]::Dict{Tuple{SVector{N,Int8}, SVector{N,Int8}, SVector{N,Int8}}, Float64}
      key = _key(ll, mm, kk)
      if haskey(vals, key)
         val  = vals[key]
      else
         val = _compute_val(A, key...)
         vals[key] = val
      end
      return val
   end

   function _compute_val(A::Rot3DCoeffs, ll::StaticVector{N},
                                        mm::StaticVector{N},
                                        kk::StaticVector{N}) where {N}
      val = 0.0
      llp = ll[1:N-2]
      mmp = mm[1:N-2]
      kkp = kk[1:N-2]
      for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
         if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
            continue
         end
         val += A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N]) *
                A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N]) *
                      A( SVector(llp..., j),
                         SVector(mmp..., mm[N-1]+mm[N]),
                         SVector(kkp..., kk[N-1]+kk[N]) )
      end
      return val
   end

   function _compute_val(A::Rot3DCoeffs, ll::StaticVector{2},
                                        mm::StaticVector{2},
                                        kk::StaticVector{2})
      if ll[1] != ll[2] || sum(mm) != 0 || sum(kk) != 0
         return 0.0
      else
         return 8 * pi^2 / (2*ll[1]+1) * (-1)^(mm[1]-kk[1])
      end
   end

   compute_Al(ll::SVector{N}) where {N} =
      compute_Al(Rot3DCoeffs(N, sum(ll)), ll)

   function compute_Al(A::Rot3DCoeffs, ll::SVector)
      len = 0
      for mm in _mrange(ll)
         len += 1
      end
      CC = zeros(len, len)
      for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
         CC[ik, im] = A(ll, mm, kk)
      end
      return CC
   end

   compute_Al_symm(ll::SVector{N}) where {N} =
      compute_Al_symm(Rot3DCoeffs(N, sum(ll)), ll)

   function compute_Al_symm(A::Rot3DCoeffs, ll::SVector)
      len = 0
      for mm in _mrange(ll)
         len += 1
      end
      CC = zeros(len, len)
      KK = Vector{Any}(undef, len)
      for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
         CC[ik, im] = A(ll, mm, kk)
         KK[ik] = kk
      end

      return CC, KK
   end

end

using Combinatorics
ll = SVector(4,4,4,4)
CC, KK = RotationCoeffs.compute_Al_symm(ll)
DK = Dict([ kk => ik for (ik, kk) in enumerate(KK)]...)
sumCC = zeros(size(CC))
for p in permutations(1:4)
   global sumCC
   # kk = KK[i] then a_kk,j = CC[i,j]
   # kk[p] = KK[ip] then CCp[ip,:] = CC[i,:]
   CCp = zeros(size(CC))
   for i = 1:size(CC,1)
      kk = KK[i]
      # kk[p] is the permutation of p
      ip = DK[SVector(kk[p] ...)]
      CCp[ip, :] = CC[i, :]
   end
   sumCC += CCp
end

rank(CC)
rank(sumCC)
svdvals(sumCC)

function Alkm_old(ll::SVector{2}, mm, kk, cg)
   if ll[1] != ll[2] ||  sum(mm) != 0 || sum(kk) != 0
      return 0.0
   end
   return 8 * pi^2 / (2*ll[1]+1) * (-1)^(mm[1]-kk[1])
end

function Alkm_old(ll::SVector{3}, mm, kk, cg)
   if sum(mm) != 0 || sum(kk) != 0
      return 0.0
   end
   return 8 * pi^2 / (2*ll[3]+1) * (-1)^(mm[3]-kk[3]) *
          cg(ll[1], mm[1], ll[2], mm[2], ll[3], -mm[3]) *
          cg(ll[1], kk[1], ll[2], kk[2], ll[3], -kk[3])
end

function Alkm_old(ll::SVector{4}, mm, kk, cg)
   Alkm = 0.0
   jlo = max(abs(ll[1]-ll[2]), abs(ll[3]-ll[4]))
   jhi = min(ll[1]+ll[2], ll[3]+ll[4])
   for j = jlo:jhi
      if (abs(mm[1]+mm[2]) > j) || (abs(mm[3]+mm[4]) > j) ||
         (abs(kk[1]+kk[2]) > j) || (abs(kk[3]+kk[4]) > j)
         continue
      end
      Alkm += 8 * pi^2 * (-1)^(mm[1]+mm[2]-kk[1]-kk[2]) / (2*j+1) *
                    cg(ll[1], mm[1], ll[2], mm[2], j, mm[1]+mm[2]) *
                    cg(ll[3], mm[3], ll[4], mm[4], j, mm[3]+mm[4]) *
                    cg(ll[1], kk[1], ll[2], kk[2], j, kk[1]+kk[2]) *
                    cg(ll[3], kk[3], ll[4], kk[4], j, kk[3]+kk[4])
   end
   return Alkm
end

##
A = RotationCoeffs.Rot3DCoeffs(5, 12)
cg = ClebschGordan(12)

## len-2
l2a = SVector(2,3)
l2b = SVector(2,2)
m2 = SVector(-1,1)
k2 = SVector(2,-2)
Alkm_old(l2a, m2, k2, cg)
A(l2a, m2, k2)
Alkm_old(l2b, m2, k2, cg)
A(l2b, m2, k2)

## len-3
l3 = SVector(3,2,3)
m3 = SVector(-2,1,1)
k3 = SVector(0,2,-2)
Alkm_old(l3, m3, k3, cg)
A(l3, m3, k3)

##  len-4
l4 = SVector(3,3,3,3)
m4 = SVector(-2,1,3,-2)
k4 = SVector(0,3,-2,-1)
Alkm_old(l4, m4, k4, cg)
A(l4, m4, k4)

A4 = A.vals[4]

ll = SVector(5,3,2,2)
@time RotationCoeffs.compute_Al(A, ll)

ll = SVector(3,2,4)
A = RotationCoeffs.Rot3DCoeffs(5, 12)
ll = SVector(3,2)
@time Cl = RotationCoeffs.compute_Al(A, ll)
ll = SVector(4,3,3,2)
@time Cl = RotationCoeffs.compute_Al(A, ll)
@time Cold = compute_Ckm_old(ll)

Ul = svd(Cl).U[:, 1:5]
Uo = svd(Cold).U[:, 1:5]

rank([Ul Uo])

ll = SVector(2, 2, 2, 2, 2)
Cl = RotationCoeffs.compute_Al(A, ll)
rank(Cl)
Matrix(svd(Cl).U)
A.vals[5]

function compute_Ckm_old(ll::SVector{4})
   cg = ClebschGordan(sum(ll))

   len = 0
   for mm in _mrange(ll)
      len += 1
   end
   @show len

   Ckm = zeros(len, len)

   for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
      jlo = max(abs(ll[1]-ll[2]), abs(ll[3]-ll[4]))
      jhi = min(ll[1]+ll[2], ll[3]+ll[4])
      for j = jlo:jhi
         if (abs(mm[1]+mm[2]) > j) || (abs(mm[3]+mm[4]) > j) ||
            (abs(kk[1]+kk[2]) > j) || (abs(kk[3]+kk[4]) > j)
            continue
         end
         cg1 = cg(ll[1], mm[1], ll[2], mm[2], j, mm[1]+mm[2])
         cg2 = cg(ll[3], mm[3], ll[4], mm[4], j, mm[3]+mm[4])
         cg3 = cg(ll[1], kk[1], ll[2], kk[2], j, kk[1]+kk[2])
         # @show ll[3], kk[3], ll[4], kk[4], j, kk[3]+kk[4]
         cg4 = cg(ll[3], kk[3], ll[4], kk[4], j, kk[3]+kk[4])
         Ckm[ik,im] += 8 * pi^2 * (-1)^(mm[1]+mm[2]-kk[1]-kk[2]) / (2*j+1) *
                       cg(ll[1], mm[1], ll[2], mm[2], j, mm[1]+mm[2]) *
                       cg(ll[3], mm[3], ll[4], mm[4], j, mm[3]+mm[4]) *
                       cg(ll[1], kk[1], ll[2], kk[2], j, kk[1]+kk[2]) *
                       cg(ll[3], kk[3], ll[4], kk[4], j, kk[3]+kk[4])
      end
   end
   return Ckm
end

##

##

# CASE 1
ll1 = SVector(2,1,1,2)
Ckm = compute_Ckm(ll1)
@show rank(Ckm)
svdf = svd(Ckm)
@show svdf.S[1:5]
for i = 1:3
   @info("V[:,$i]")
   @show round.(svdf.Vt[i,:], digits=2)
end

# CASE 2
ll2 = SVector(2,3,4,3)
Ckm = compute_Ckm(ll2)
@show rank(Ckm)
@show svdvals(Ckm)[1:8]

# CASE 3
ll3 = SVector(5,4,4,3)
Ckm = compute_Ckm(ll3)
@show rank(Ckm)
@show svdvals(Ckm)[1:10]

# MAIN EXAMPLE:
ll = SVector(1,1,1,1)
Ckm = compute_Ckm(ll)
Ckm_sym = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 15.7914 0.0 -7.89568 -7.89568 2.63189 5.26379 2.63189 0.0 -7.89568 -7.89568 5.26379 10.5276 5.26379 -7.89568 -7.89568 0.0 2.63189 5.26379 2.63189 -7.89568 -7.89568 0.0 15.7914 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -7.89568 0.0 10.5276 -2.63189 -7.89568 -2.63189 5.26379 0.0 -2.63189 10.5276 -2.63189 -5.26379 -2.63189 10.5276 -2.63189 0.0 5.26379 -2.63189 -7.89568 -2.63189 10.5276 0.0 -7.89568 0.0 0.0; 0.0 0.0 -7.89568 0.0 -2.63189 10.5276 5.26379 -2.63189 -7.89568 0.0 10.5276 -2.63189 -2.63189 -5.26379 -2.63189 -2.63189 10.5276 0.0 -7.89568 -2.63189 5.26379 10.5276 -2.63189 0.0 -7.89568 0.0 0.0; 0.0 0.0 2.63189 0.0 -7.89568 5.26379 15.7914 -7.89568 2.63189 0.0 5.26379 -7.89568 -7.89568 10.5276 -7.89568 -7.89568 5.26379 0.0 2.63189 -7.89568 15.7914 5.26379 -7.89568 0.0 2.63189 0.0 0.0; 0.0 0.0 5.26379 0.0 -2.63189 -2.63189 -7.89568 10.5276 -7.89568 0.0 -2.63189 -2.63189 10.5276 -5.26379 10.5276 -2.63189 -2.63189 0.0 -7.89568 10.5276 -7.89568 -2.63189 -2.63189 0.0 5.26379 0.0 0.0; 0.0 0.0 2.63189 0.0 5.26379 -7.89568 2.63189 -7.89568 15.7914 0.0 -7.89568 5.26379 -7.89568 10.5276 -7.89568 5.26379 -7.89568 0.0 15.7914 -7.89568 2.63189 -7.89568 5.26379 0.0 2.63189 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -7.89568 0.0 -2.63189 10.5276 5.26379 -2.63189 -7.89568 0.0 10.5276 -2.63189 -2.63189 -5.26379 -2.63189 -2.63189 10.5276 0.0 -7.89568 -2.63189 5.26379 10.5276 -2.63189 0.0 -7.89568 0.0 0.0; 0.0 0.0 -7.89568 0.0 10.5276 -2.63189 -7.89568 -2.63189 5.26379 0.0 -2.63189 10.5276 -2.63189 -5.26379 -2.63189 10.5276 -2.63189 0.0 5.26379 -2.63189 -7.89568 -2.63189 10.5276 0.0 -7.89568 0.0 0.0; 0.0 0.0 5.26379 0.0 -2.63189 -2.63189 -7.89568 10.5276 -7.89568 0.0 -2.63189 -2.63189 10.5276 -5.26379 10.5276 -2.63189 -2.63189 0.0 -7.89568 10.5276 -7.89568 -2.63189 -2.63189 0.0 5.26379 0.0 0.0; 0.0 0.0 10.5276 0.0 -5.26379 -5.26379 10.5276 -5.26379 10.5276 0.0 -5.26379 -5.26379 -5.26379 15.7914 -5.26379 -5.26379 -5.26379 0.0 10.5276 -5.26379 10.5276 -5.26379 -5.26379 0.0 10.5276 0.0 0.0; 0.0 0.0 5.26379 0.0 -2.63189 -2.63189 -7.89568 10.5276 -7.89568 0.0 -2.63189 -2.63189 10.5276 -5.26379 10.5276 -2.63189 -2.63189 0.0 -7.89568 10.5276 -7.89568 -2.63189 -2.63189 0.0 5.26379 0.0 0.0; 0.0 0.0 -7.89568 0.0 10.5276 -2.63189 -7.89568 -2.63189 5.26379 0.0 -2.63189 10.5276 -2.63189 -5.26379 -2.63189 10.5276 -2.63189 0.0 5.26379 -2.63189 -7.89568 -2.63189 10.5276 0.0 -7.89568 0.0 0.0; 0.0 0.0 -7.89568 0.0 -2.63189 10.5276 5.26379 -2.63189 -7.89568 0.0 10.5276 -2.63189 -2.63189 -5.26379 -2.63189 -2.63189 10.5276 0.0 -7.89568 -2.63189 5.26379 10.5276 -2.63189 0.0 -7.89568 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 2.63189 0.0 5.26379 -7.89568 2.63189 -7.89568 15.7914 0.0 -7.89568 5.26379 -7.89568 10.5276 -7.89568 5.26379 -7.89568 0.0 15.7914 -7.89568 2.63189 -7.89568 5.26379 0.0 2.63189 0.0 0.0; 0.0 0.0 5.26379 0.0 -2.63189 -2.63189 -7.89568 10.5276 -7.89568 0.0 -2.63189 -2.63189 10.5276 -5.26379 10.5276 -2.63189 -2.63189 0.0 -7.89568 10.5276 -7.89568 -2.63189 -2.63189 0.0 5.26379 0.0 0.0; 0.0 0.0 2.63189 0.0 -7.89568 5.26379 15.7914 -7.89568 2.63189 0.0 5.26379 -7.89568 -7.89568 10.5276 -7.89568 -7.89568 5.26379 0.0 2.63189 -7.89568 15.7914 5.26379 -7.89568 0.0 2.63189 0.0 0.0; 0.0 0.0 -7.89568 0.0 -2.63189 10.5276 5.26379 -2.63189 -7.89568 0.0 10.5276 -2.63189 -2.63189 -5.26379 -2.63189 -2.63189 10.5276 0.0 -7.89568 -2.63189 5.26379 10.5276 -2.63189 0.0 -7.89568 0.0 0.0; 0.0 0.0 -7.89568 0.0 10.5276 -2.63189 -7.89568 -2.63189 5.26379 0.0 -2.63189 10.5276 -2.63189 -5.26379 -2.63189 10.5276 -2.63189 0.0 5.26379 -2.63189 -7.89568 -2.63189 10.5276 0.0 -7.89568 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 15.7914 0.0 -7.89568 -7.89568 2.63189 5.26379 2.63189 0.0 -7.89568 -7.89568 5.26379 10.5276 5.26379 -7.89568 -7.89568 0.0 2.63189 5.26379 2.63189 -7.89568 -7.89568 0.0 15.7914 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

badI = [1,2,4,10, 18, 24, 26, 27]
goodI = setdiff(1:27, badI)
Ckm_sym = Ckm_sym[goodI, goodI]

display(Ckm_sym[:, [1,2,3]])
display(Ckm[:, [1,2,3]])

norm(Ckm_sym - Ckm, Inf)
svdvals(Ckm)
svdvals(Ckm_sym)

ll = SVector(2,1,1,1)
Ckm = compute_Ckm(ll)
rank(Ckm)
c1 = Ckm[:,1]
c2 = Ckm[:,3]
c1 /= norm(c1)
c2 /= norm(c2)
[sort(c1) sort(c2)]

U = svd(Ckm).U[:,1:2]
[ sort(U[:,1]) sort(U[:,2]) ]
U' * U

svd(Ckm)

for l = 1:6
   ll = SVector(l,l,l,l)
   Ckm = compute_Ckm(ll)
   @show l, rank(Ckm)
end

using Profile
A = RotationCoeffs.Rot3DCoeffs(5, 12)
ll = SVector(1,1,1,1,1)
ll = SVector(4,4,2,1,1)
@profile Al = RotationCoeffs.compute_Al(A, ll)

Profile.print()


text:
mutable struct Board
	wonAt::Int32
	rows::Vector{Vector{Int32}}
	columns::Vector{Vector{Int32}}
end

input = readlines("input.txt")

boards = Vector{Board}()
for boardDefinition in Iterators.partition(input[3:end], 6)
	board = Board(0, Vector{Vector{Int32}}(), Vector{Vector{Int32}}())
	for row in boardDefinition[1:end-1]
		push!(board.rows, map(x -> parse(Int32, x), split(strip(row), r"\s+")))
	end
	board.columns = [x[:] for x in eachrow(hcat(board.rows...))]
	push!(boards, board)
end

winners = Vector{Board}()
for number in map(x -> parse(Int32, x), split(input[1], ","))
	for board in boards
		f(x) = x != number
		board.rows = [filter(f, x) for x in board.rows]
		board.columns = [filter(f, x) for x in board.columns]

		e(x) = size(x, 1) == 0
		if findfirst(e, board.rows) != nothing || findfirst(e, board.columns) != nothing
			board.wonAt = number
			push!(winners, board)
		end
	end
	deleteat!(boards, findall(x -> x in winners, boards))
end

function score(board)
	return sum([sum(x) for x in board.rows]) * board.wonAt
end

println("The answer to the first part is: ", score(winners[1]))
println("The answer to the second part is: ", score(winners[end]))

max_stars_repo_name: Kwarf/adventofcode

max_issues_repo_name: Kwarf/adventofcode

max_forks_repo_name: Kwarf/adventofcode


text:
import RigidBodyTools: assign_velocity!

struct LidDrivenCavity <: DirectlySpecifiedMotion
    U :: Float64
end

"""
    assign_velocity!(u,v,x,y,m::DirectlySpecifiedMotion,t)

Given a `DirectlySpecifiedMotion`, vectors of `x` and `y` coordinates,
and time `t`, return the velocities `u` and `v` at those points.
"""
function assign_velocity!(u::AbstractVector{Float64},v::AbstractVector{Float64},
                 b::Body,
                 m::LidDrivenCavity,t::Real)

   fill!(v,0.0)
   fill!(u,0.0)
   # find the elements of the x and y vectors that lie on the
   # top wall, and assign velocity m.U to their u values.

   nx, ny = normalmid(b)
   topindices = findall(x -> (x ≈ 1.0), ny)
   u[topindices] .= m.U

   return u, v

end


text:
# tests sequence choice points arising from reps, mult, and plus

#
# reps(function_as_symbol, [min[, max]])
#

# Can't define generators inside a @testset as it gives rise to error:
#    LoadError: error compiling anonymous: type definition not allowed inside a local scope
# but even though this scope is inside the local scope of a @testset block in runtests.jl, the
# the fact the code is included in the local scope appears to avoid this issue

@generator SCRepsGen begin # prefix SC (for Sequence Choice points) to avoid type name clashes
start() = reps(a,2,4)
a() = gensym()
end

@generator SCRepsNoMaxGen begin
start() = reps(a,2)
a() = 'a'
end

@generator SCRepsNoMinGen begin
start() = reps(a)
a() = 'a'
end

@generator SCRepsNonLiteralMinMaxGen begin
start() = reps(a,2*2,5+1)
a() = 'a'
end

@generator SCRepsRuleNonLiteralMinMaxGen begin
start() = reps(a,x(),y())
    a() = 'a'
x() = 4
y() = 6
end

@generator SCRepsShortFormParenGen begin
start() = reps(a(),2,4)
    a() = gensym()
end

@generator SCRepsRuleParamGen begin
start() = begin
    a = 12
    b = 60
    reps(add(a,b),2,4)
end
    add(x,y) = (gensym(),x+y)
end

@generator SCRepsNonRuleGen begin
start() = reps(gensym(),2,4)
end

@generator SCSubGen() begin
start() = gensym()
end

@generator SCRepsSubGenGen(subgen) begin
start() = reps(choose(subgen),2,4)
end

@generator SCMultGen begin
start() = mult(a)
a() = 'a'
end

@generator SCPlusGen begin
start() = plus(a)
a() = 'a'
end

@testset "sequence choice points" begin

	@testset "reps choice point" begin

		gn = SCRepsGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
	    	td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td) <= 4
		    @test all(map(x->typeof(x),td) .== Symbol)
		    @test td[1] != td[2]
		    @mtest_values_include [2,3,4] length(td)
		end

	end

	@testset "reps choice point with no maximum" begin

		gn = SCRepsNoMaxGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td)
		    @test all(td .== 'a')
		    @mtest_values_include [2,3,4] length(td)
		end

	end

	@testset "reps choice point with no minimum nor maximum" begin

		gn = SCRepsNoMinGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 0 <= length(td)
		    @test all(td .== 'a')
		    @mtest_values_include [0,1,2] length(td)
		end

	end

	@testset "reps choice point with non-literal minimum and maximum" begin

		gn = SCRepsNonLiteralMinMaxGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 4 <= length(td) <= 6
		    @test all(td .== 'a')
		    @mtest_values_are [4,5,6] length(td)
		end

	end

	@testset "reps choice point with non-literal minimum and maximum defined by rules" begin

		gn = SCRepsRuleNonLiteralMinMaxGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 4 <= length(td) <= 6
		    @test all(td .== 'a')
		    @mtest_values_are [4,5,6] length(td)
		end

	end

	@testset "reps choice point where function called uses short form with parentheses" begin

		gn = SCRepsShortFormParenGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td) <= 4
		    @test all(map(x->typeof(x),td) .== Symbol)
		    @test td[1] != td[2]
		    @mtest_values_are [2,3,4] length(td)
		end

	end

	@testset "reps choice point where function called with parameters" begin

		gn = SCRepsRuleParamGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td) <= 4
		    td1 = map(x->x[1],td)
		    td2 = map(x->x[2],td)
		    @test all(map(x->typeof(x),td1) .== Symbol)
		    @test td1[1] != td1[2]
		    @test all(td2 .== 72)
		    @mtest_values_are [2,3,4] length(td)
		end

	end

	@testset "reps choice point where function called is not a rule" begin

		gn = SCRepsNonRuleGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td) <= 4
		    @test all(map(x->typeof(x),td) .== Symbol)
		    @test td[1] != td[2]
		    @mtest_values_are [2,3,4] length(td)
		end

	end

	@testset "reps choice point where function called is not a rule" begin

		sg = SCSubGen()
		gn = SCRepsSubGenGen(sg)

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 2 <= length(td) <= 4
		    @test all(map(x->typeof(x),td) .== Symbol)
		    @test td[1] != td[2]
		    @mtest_values_are [2,3,4] length(td)
		end

	end

	#
	# mult(function_as_symbol)
	#

	@testset "mult choice point" begin

		gn = SCMultGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 0 <= length(td)
		    @test all(td .== 'a')
		    @mtest_values_include [0,1,2] length(td)
		end

	end

	#
	# plus(function_as_symbol)
	#

	@testset "plus choice point" begin

		gn = SCPlusGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) <: Array
		    @test 1 <= length(td)
		    @test all(td .== 'a')
		    @mtest_values_include [1,2,3] length(td)
		end

	end

end


text:
module CommunalHelperMoveBlockRedirect

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/MoveBlockRedirect" MoveBlockRedirect(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    direction::String="Up",
    fastRedirect::Bool=false,
    deleteBlock::Bool=false,
    oneUse::Bool=false,
    modifier::Number=0.0,
    operation::String="Add",
    reskinFolder::String="",
    overrideColor::String="",
    overrideUsedColor::String="",
)

const placements = Ahorn.PlacementDict(
    "Move Block Redirect (Communal Helper)" => Ahorn.EntityPlacement(
        MoveBlockRedirect,
        "rectangle",
    ),
    "Move Block Redirect (One Use) (Communal Helper)" => Ahorn.EntityPlacement(
        MoveBlockRedirect,
        "rectangle",
        Dict{String,Any}(
            "oneUse" => true,
        ),
    ),
    "Move Block Redirect (Delete Block) (Communal Helper)" => Ahorn.EntityPlacement(
        MoveBlockRedirect,
        "rectangle",
        Dict{String,Any}(
            "deleteBlock" => true,
        ),
    ),
)

const opTypes = ["Add", "Subtract", "Multiply"]

Ahorn.editingOptions(entity::MoveBlockRedirect) = Dict{String, Any}(
    "direction" => Maple.move_block_directions,
    "operation" => opTypes,
)

Ahorn.minimumSize(entity::MoveBlockRedirect) = Maple.defaultBlockWidth, Maple.defaultBlockHeight

Ahorn.resizable(entity::MoveBlockRedirect) = true, true

Ahorn.selection(entity::MoveBlockRedirect) = Ahorn.getEntityRectangle(entity)

const defaultColor = (251, 206, 54, 255) ./ 255
const deleteColor = (204, 37, 65, 255) ./ 255
const fastColor = (41, 195, 47, 255) ./ 255
const slowColor = (28, 91, 179, 255) ./ 255

const clockwise = ["Up", "Right", "Down", "Left"]
function Ahorn.rotated(entity::MoveBlockRedirect, steps::Int)
    dir = get(entity.data, "direction", "Up")
    idx = findall(d -> d == dir, clockwise)[1]
    entity.data["direction"] = clockwise[mod1(idx + steps, 4)]

    return entity
end

const rotations = Dict{String,Float64}(
    "Up" => pi * 1.5,
    "Right" => 0,
    "Down" => pi / 2,
    "Left" => pi,
)

function getRotation(dir::String)
    if haskey(rotations, dir)
        return rotations[dir]
    elseif (fAngle = tryparse(Float64, dir)) !== nothing
        return fAngle
    else
        return 0
    end
end

function getTexturesAndColor(entity::MoveBlockRedirect)
    path = "objects/CommunalHelper/moveBlockRedirect/"
    reskinFolder = String(get(entity.data, "reskinFolder", ""))
    if reskinFolder != ""
        if last(reskinFolder, 1) != "/"
            reskinFolder = reskinFolder * "/"
        end
        path = "objects/" * reskinFolder
    end

    block = path * "block"
    overrideColorHex = String(get(entity.data, "overrideColor", ""))
    overrideColor = (overrideColorHex != "" && length(overrideColorHex) == 6) ? hexToRGBA(overrideColorHex) : nothing

    deleteBlock = Bool(get(entity.data, "deleteBlock", false))
    if deleteBlock
        return block, path * "x", overrideColor === nothing ? deleteColor : overrideColor
    end

    operation = get(entity.data, "operation", "Add")
    modifier = abs(get(entity.data, "modifier", 0.0))
    if operation == "Add"
        if modifier != 0.0
            return block, path * "fast", overrideColor === nothing ? fastColor : overrideColor
        end
    elseif operation == "Subtract"
        if modifier != 0
            return block, path * "slow", overrideColor === nothing ? slowColor : overrideColor
        end
    elseif operation == "Multiply"
        if modifier == 0.0
            return block, path * "x", overrideColor === nothing ? deleteColor : overrideColor
        elseif modifier > 1.0
            return block, path * "fast", overrideColor === nothing ? fastColor : overrideColor
        elseif modifier < 1.0
            return block, path * "slow", overrideColor === nothing ? slowColor : overrideColor
        end
    end

    return block, path * "arrow", overrideColor === nothing ? defaultColor : overrideColor
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::MoveBlockRedirect)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    direction = String(get(entity.data, "direction", "Up"))
    block, sprite, color = getTexturesAndColor(entity)

    tileWidth = ceil(Int, width / 8)
    tileHeight = ceil(Int, height / 8)

    for i in -1:tileWidth, j in -1:tileHeight
        tx = (i == -1) ? 0 : ((i == tileWidth) ? 16 : 8)
        ty = (j == -1) ? 0 : ((j == tileHeight) ? 16 : 8)
        if i == -1 || i == tileWidth || j == -1 || j == tileHeight
            Ahorn.drawImage(ctx, block, i * 8, j * 8, tx, ty, 8, 8, tint=color)
        end
    end

    # finicky
    Ahorn.Cairo.save(ctx)
    Ahorn.Cairo.translate(ctx, width / 2, height / 2)
    Ahorn.Cairo.rotate(ctx, getRotation(direction))
    Ahorn.Cairo.translate(ctx, -8, -8) # sprite always has a 16 x 16 size
    Ahorn.drawImage(ctx, sprite, 0, 0, tint=color)
    Ahorn.Cairo.restore(ctx)
end

end


text:
a = 5 - 3im

b = complex(5,3)

c = 5 + (3/4)im      # Parentheses needed since 3/4im == 3/(4im)

a * b

a / b

(1 + im)^2

sqrt(-1.0 + 0im)  # OK, complex input gives complex output

sqrt(-1.0)        # Error - real input but needs complex output

function roots_of_quadratic(a,b,c)
    # Compute the roots of the quadratic ax^2 + bx + c = 0
    d = sqrt(complex(b^2 - 4*a*c))
    r1 = (-b - d) / 2a
    r2 = (-b + d) / 2a
    r1, r2
end

roots_of_quadratic(1,0,1)

roots_of_quadratic(1,4,13)

roots_of_quadratic(-1,5,6)


text:
export NormalMixture, NormalMixtureNode
export GaussianMixture, GaussianMixtureNode

# Normal Mixture Functional Form
struct NormalMixture{N} end

# Special node
# Generic FactorNode implementation does not work with dynamic number of inputs
# We need to reimplement the following set of functions
# functionalform(factornode::FactorNode)
# sdtype(factornode::FactorNode)
# interfaces(factornode::FactorNode)
# factorisation(factornode::FactorNode)
# localmarginals(factornode::FactorNode)
# localmarginalnames(factornode::FactorNode)
# metadata(factornode::FactorNode)
# get_pipeline_stages(factornode::FactorNode)
#
# setmarginal!(factornode::FactorNode, cname::Symbol, marginal)
# getmarginal!(factornode::FactorNode, localmarginal::FactorNodeLocalMarginal)
#
# functional_dependencies(factornode::FactorNode, iindex::Int)
# get_messages_observable(factornode, message_dependencies)
# get_marginals_observable(factornode, marginal_dependencies)
#
# score(::Type{T}, ::FactorBoundFreeEnergy, ::Stochastic, node::AbstractFactorNode, scheduler) where T
#
# Base.show
#
const NormalMixtureNodeFactorisationSupport = Union{MeanField, }

struct NormalMixtureNode{N, F <: NormalMixtureNodeFactorisationSupport, M, P} <: AbstractFactorNode
    factorisation :: F

    # Interfaces
    out    :: NodeInterface
    switch :: NodeInterface
    means  :: NTuple{N, IndexedNodeInterface}
    precs  :: NTuple{N, IndexedNodeInterface}

    meta     :: M
    pipeline :: P
end

const GaussianMixture     = NormalMixture
const GaussianMixtureNode = NormalMixtureNode

functionalform(factornode::NormalMixtureNode{N}) where N = NormalMixture{N}
sdtype(factornode::NormalMixtureNode)                    = Stochastic()
interfaces(factornode::NormalMixtureNode)                = (factornode.out, factornode.switch, factornode.means..., factornode.precs...)
factorisation(factornode::NormalMixtureNode)             = factornode.factorisation
localmarginals(factornode::NormalMixtureNode)            = error("localmarginals() function is not implemented for NormalMixtureNode")
localmarginalnames(factornode::NormalMixtureNode)        = error("localmarginalnames() function is not implemented for NormalMixtureNode")
metadata(factornode::NormalMixtureNode)                  = factornode.meta
getpipeline(factornode::NormalMixtureNode)               = factornode.pipeline

setmarginal!(factornode::NormalMixtureNode, cname::Symbol, marginal)                = error("setmarginal() function is not implemented for NormalMixtureNode")
getmarginal!(factornode::NormalMixtureNode, localmarginal::FactorNodeLocalMarginal) = error("getmarginal() function is not implemented for NormalMixtureNode")

## activate!

struct NormalMixtureNodeFunctionalDependencies <: AbstractNodeFunctionalDependenciesPipeline end

default_functional_dependencies_pipeline(::Type{ <: NormalMixture }) = NormalMixtureNodeFunctionalDependencies()

function functional_dependencies(::NormalMixtureNodeFunctionalDependencies, factornode::NormalMixtureNode{N, F}, iindex::Int) where { N, F <: MeanField }
    message_dependencies = ()

    marginal_dependencies = if iindex === 1
        (factornode.switch, factornode.means, factornode.precs)
    elseif iindex === 2
        (factornode.out, factornode.means, factornode.precs)
    elseif 2 < iindex <= N + 2
        (factornode.out, factornode.switch, factornode.precs[ iindex - 2 ])
    elseif N + 2 < iindex <= 2N + 2
        (factornode.out, factornode.switch, factornode.means[ iindex - N - 2 ])
    else
        error("Bad index in functional_dependencies for NormalMixtureNode")
    end

    return message_dependencies, marginal_dependencies
end

function get_messages_observable(factornode::NormalMixtureNode{N, F}, message_dependencies::Tuple{}) where { N, F <: MeanField }
    return nothing, of(nothing)
end

function get_marginals_observable(
    factornode::NormalMixtureNode{N, F},
    marginal_dependencies::Tuple{ NodeInterface, NTuple{N, IndexedNodeInterface}, NTuple{N, IndexedNodeInterface} }) where { N, F <: MeanField }

    varinterface    = marginal_dependencies[1]
    meansinterfaces = marginal_dependencies[2]
    precsinterfaces = marginal_dependencies[3]

    marginal_names = Val{ (name(varinterface), name(meansinterfaces[1]), name(precsinterfaces[1])) }
    marginals_observable = combineLatest((
        getmarginal(connectedvar(varinterface), IncludeAll()),
        combineLatest(map((prec) -> getmarginal(connectedvar(prec), IncludeAll()), reverse(precsinterfaces)), PushNew()),
        combineLatest(map((mean) -> getmarginal(connectedvar(mean), IncludeAll()), reverse(meansinterfaces)), PushNew()),
    ), PushNew()) |> map_to((
        getmarginal(connectedvar(varinterface), IncludeAll()),
        map((mean) -> getmarginal(connectedvar(mean), IncludeAll()), meansinterfaces),
        map((prec) -> getmarginal(connectedvar(prec), IncludeAll()), precsinterfaces)
    ))

    return marginal_names, marginals_observable
end

function get_marginals_observable(
    factornode::NormalMixtureNode{N, F},
    marginal_dependencies::Tuple{ NodeInterface, NodeInterface, IndexedNodeInterface }) where { N, F <: MeanField }

    outinterface    = marginal_dependencies[1]
    switchinterface = marginal_dependencies[2]
    varinterface    = marginal_dependencies[3]

    marginal_names       = Val{ (name(outinterface), name(switchinterface), name(varinterface)) }
    marginals_observable = combineLatestUpdates((
        getmarginal(connectedvar(outinterface), IncludeAll()),
        getmarginal(connectedvar(switchinterface), IncludeAll()),
        getmarginal(connectedvar(varinterface), IncludeAll()),
    ), PushNew())

    return marginal_names, marginals_observable
end

# FreeEnergy related functions

@average_energy NormalMixture (q_out::Any, q_switch::Any, q_m::NTuple{N, NormalMeanVariance}, q_p::NTuple{N, GammaDistributionsFamily}) where N = begin
    z_bar = probvec(q_switch)
    return mapreduce(+, 1:N, init = 0.0) do i
        return z_bar[i] * score(AverageEnergy(), NormalMeanPrecision, Val{ (:out, :μ, :τ) }, map((q) -> Marginal(q, false, false), (q_out, q_m[i], q_p[i])), nothing)
    end
end

@average_energy NormalMixture (q_out::Any, q_switch::Any, q_m::NTuple{N, MvNormalMeanCovariance}, q_p::NTuple{N, Wishart}) where N = begin
    z_bar = probvec(q_switch)
    return mapreduce(+, 1:N, init = 0.0) do i
        return z_bar[i] * score(AverageEnergy(), MvNormalMeanPrecision, Val{ (:out, :μ, :Λ) }, map((q) -> Marginal(q, false, false), (q_out, q_m[i], q_p[i])), nothing)
    end
end

function score(::Type{T}, objective::BetheFreeEnergy, ::FactorBoundFreeEnergy, ::Stochastic, node::NormalMixtureNode{N, MeanField}, scheduler) where { T <: InfCountingReal, N }

    skip_strategy = marginal_skip_strategy(objective)

    stream = combineLatest((
        getmarginal(connectedvar(node.out), skip_strategy) |> schedule_on(scheduler),
        getmarginal(connectedvar(node.switch), skip_strategy) |> schedule_on(scheduler),
        combineLatest(map((mean) -> getmarginal(connectedvar(mean), skip_strategy) |> schedule_on(scheduler), node.means), PushNew()),
        combineLatest(map((prec) -> getmarginal(connectedvar(prec), skip_strategy) |> schedule_on(scheduler), node.precs), PushNew())
    ), PushNew())

    mapping = let fform = functionalform(node), meta = metadata(node)
        (marginals) -> begin
            average_energy   = score(AverageEnergy(), fform, Val{ (:out, :switch, :m, :p) }, marginals, meta)

            out_entropy     = score(DifferentialEntropy(), marginals[1])
            switch_entropy  = score(DifferentialEntropy(), marginals[2])
            means_entropies = mapreduce((m) -> score(DifferentialEntropy(), m), +, marginals[3])
            precs_entropies = mapreduce((m) -> score(DifferentialEntropy(), m), +, marginals[4])

            return convert(T, average_energy - (out_entropy + switch_entropy + means_entropies + precs_entropies))
        end
    end

    return stream |> map(T, mapping)
end

as_node_functional_form(::Type{ <: NormalMixture }) = ValidNodeFunctionalForm()

# Node creation related functions

sdtype(::Type{ <: NormalMixture }) = Stochastic()

collect_factorisation(::Type{ <: NormalMixture }, factorisation) = factorisation

function ReactiveMP.make_node(::Type{ <: NormalMixture{N} }; factorisation::F = MeanField(), meta::M = nothing, pipeline::P = nothing) where { N, F, M, P }
    @assert N >= 2 "NormalMixtureNode requires at least two mixtures on input"
    @assert typeof(factorisation) <: NormalMixtureNodeFactorisationSupport "NormalMixtureNode supports only following factorisations: [ $(NormalMixtureNodeFactorisationSupport) ]"
    out    = NodeInterface(:out)
    switch = NodeInterface(:switch)
    means  = ntuple((index) -> IndexedNodeInterface(index, NodeInterface(:m)), N)
    precs  = ntuple((index) -> IndexedNodeInterface(index, NodeInterface(:p)), N)
    return NormalMixtureNode{N, F, M, P}(factorisation, out, switch, means, precs, meta, pipeline)
end

function ReactiveMP.make_node(::Type{ <: NormalMixture }, out::AbstractVariable, switch::AbstractVariable, means::NTuple{N, AbstractVariable}, precs::NTuple{N, AbstractVariable}; factorisation = MeanField(), meta = nothing, pipeline = nothing) where { N}
    node = make_node(NormalMixture{N}, factorisation = collect_factorisation(NormalMixture, factorisation), meta = collect_meta(NormalMixture, meta), pipeline = collect_pipeline(NormalMixture, pipeline))

    # out
    out_index = getlastindex(out)
    connectvariable!(node.out, out, out_index)
    setmessagein!(out, out_index, messageout(node.out))

    # switch
    switch_index = getlastindex(switch)
    connectvariable!(node.switch, switch, switch_index)
    setmessagein!(switch, switch_index, messageout(node.switch))

    # means
    foreach(zip(node.means, means)) do (minterface, mvar)
        mean_index = getlastindex(mvar)
        connectvariable!(minterface, mvar, mean_index)
        setmessagein!(mvar, mean_index, messageout(minterface))
    end

    # precs
    foreach(zip(node.precs, precs)) do (pinterface, pvar)
        prec_index = getlastindex(pvar)
        connectvariable!(pinterface, pvar, prec_index)
        setmessagein!(pvar, prec_index, messageout(pinterface))
    end

    return node
end

function ReactiveMP.make_node(fform::Type{ <: NormalMixture }, autovar::AutoVar, args::Vararg{ <: AbstractVariable }; kwargs...)
    var  = randomvar(getname(autovar))
    node = make_node(fform, var, args...; kwargs...)
    return node, var
end


text:
using LightGraphs

struct BudgetedSpanningTreeInstance{T, U}
  graph::AbstractGraph{T}
  rewards::Dict{Edge{T}, Float64}
  weights::Dict{Edge{T}, U}
  budget::U
end

abstract type BudgetedSpanningTreeSolution{T, U}
  # instance::BudgetedSpanningTreeInstance{T, U}
  # tree::Vector{Edge{T}}
end

struct BudgetedSpanningTreeLagrangianSolution{T, U} <: BudgetedSpanningTreeSolution{T, U}
  # Used to store important temporary results from solving the Lagrangian dual.
  instance::BudgetedSpanningTreeInstance{T, U}
  tree::Vector{Edge{T}}
  λ::Float64
  value::Float64
  λmax::Float64 # No dual value higher than this is useful (i.e. they all yield the same solution).
end

struct SimpleBudgetedSpanningTreeSolution{T, U} <: BudgetedSpanningTreeSolution{T, U}
  instance::BudgetedSpanningTreeInstance{T, U}
  tree::Vector{Edge{T}}
end

function _budgeted_spanning_tree_compute_value(i::Union{SpanningTreeInstance{T}, BudgetedSpanningTreeInstance{T}}, solution::Vector{Edge{T}}) where T
  return sum(i.rewards[(e in keys(i.rewards)) ? e : reverse(e)] for e in solution)
end

function _budgeted_spanning_tree_compute_weight(i::BudgetedSpanningTreeInstance{T}, solution::Vector{Edge{T}}) where T
  return sum(i.weights[(e in keys(i.weights)) ? e : reverse(e)] for e in solution)
end

function st_prim_budgeted_lagrangian(i::BudgetedSpanningTreeInstance{T}, λ::Float64) where T
  # Solve the subproblem
  #     l(λ) = \max_{x spanning tree} (rewards + λ weights) x - λ budget.
  sti_rewards = Dict{Edge{T}, Float64}(e => i.rewards[e] + λ * i.weights[e] for e in keys(i.rewards))
  sti = SpanningTreeInstance(i.graph, sti_rewards)
  sti_sol = st_prim(sti)
  sti_value = sum(sti_rewards[(e in keys(i.rewards)) ? e : reverse(e)] for e in sti_sol.tree)
  return sti_value - λ * i.budget, sti_sol.tree
end

function st_prim_budgeted_lagrangian_search(i::BudgetedSpanningTreeInstance{T}, ε::Float64) where T
  # Approximately solve the problem \min_{l ≥ 0} l(λ), where
  #     l(λ) = \max_{x spanning tree} (rewards + λ weights) x - λ budget.
  # This problem is the Lagrangian dual of the budgeted maximum spanning-tree problem:
  #     \max_{x spanning tree} rewards x  s.t.  weights x >= budget.
  # This algorithm provides no guarantee on the optimality of the solution.

  # Initial set of values for λ. The optimum is guaranteed to be contained in this interval.
  weights_norm_inf = maximum(values(i.weights)) # Maximum weight.
  m = nv(i.graph) - 1 # Maximum number of items in a solution. Easy to compute for a spanning tree!
  λmax = Float64(weights_norm_inf * m + 1)

  λlow = 0.0
  λhigh = λmax

  # Perform a golden-ratio search.
  while (λhigh - λlow) > ε
    λmidlow = λhigh - (λhigh - λlow) / MathConstants.φ
    λmidhigh = λlow + (λhigh - λlow) / MathConstants.φ
    vmidlow, _ = st_prim_budgeted_lagrangian(i, λmidlow)
    vmidhigh, _ = st_prim_budgeted_lagrangian(i, λmidhigh)

    if vmidlow < vmidhigh
      λhigh = λmidhigh
      vhigh = vmidhigh
    else
      λlow = λmidlow
      vlow = vmidlow
    end
  end

  vlow, stlow = st_prim_budgeted_lagrangian(i, λlow)
  vhigh, sthigh = st_prim_budgeted_lagrangian(i, λhigh)
  if vlow > vhigh
    return BudgetedSpanningTreeLagrangianSolution(i, stlow, λlow, vlow, λmax)
  else
    return BudgetedSpanningTreeLagrangianSolution(i, sthigh, λhigh, vhigh, λmax)
  end
end

function _solution_symmetric_difference(a::Vector{Edge{T}}, b::Vector{Edge{T}}) where T
  only_in_a = Edge{T}[]
  only_in_b = Edge{T}[]

  for e in a
    if e in b
      continue
    end
    push!(only_in_a, e)
  end

  for e in b
    if e in a
      continue
    end
    push!(only_in_b, e)
  end

  return only_in_a, only_in_b
end

function _solution_symmetric_difference_size(a::Vector{Edge{T}}, b::Vector{Edge{T}}) where T
  only_in_a, only_in_b = _solution_symmetric_difference(a, b)
  return length(only_in_a) + length(only_in_b)
end

function st_prim_budgeted_lagrangian_refinement(i::BudgetedSpanningTreeInstance{T};
                                                ε::Float64=1.0e-3, ζ⁻::Float64=0.2, ζ⁺::Float64=5.0,
                                                stalling⁻::Float64=1.0e-5) where T
  # Check assumptions.
  if ζ⁻ >= 1.0
    error("ζ⁻ must be strictly less than 1.0: the dual multiplier λ is multiplied by ζ⁻ to reach an infeasible solution by less penalising the budget constraint.")
  end
  if ζ⁺ <= 1.0
    error("ζ⁺ must be strictly greater than 1.0: the dual multiplier λ is multiplied by ζ⁺ to reach a feasible solution by penalising more the budget constraint.")
  end

  # Ensure the problem is feasible by only considering the budget constraint.
  feasible_rewards = Dict{Edge{T}, Float64}(e => i.weights[e] for e in keys(i.rewards))
  feasible_instance = SpanningTreeInstance(i.graph, feasible_rewards)
  feasible_solution = st_prim(feasible_instance)
  if _budgeted_spanning_tree_compute_value(feasible_instance, feasible_solution.tree) < i.budget
    # By maximising the left-hand side of the budget constraint, impossible to reach the target budget. No solution!
    return SimpleBudgetedSpanningTreeSolution(i, similar(feasible_solution.tree, 0))
  end

  # Solve the Lagrangian relaxation to optimality.
  lagrangian = st_prim_budgeted_lagrangian_search(i, ε)
  λ0, v0, st0 = lagrangian.λ, lagrangian.value, lagrangian.tree
  λmax = lagrangian.λmax
  b0 = _budgeted_spanning_tree_compute_weight(i, st0) # Budget consumption of this first solution.

  # Find two solutions: one above the budget x⁺ (i.e. respecting the constraint), the other not x⁻.
  x⁺, x⁻ = nothing, nothing

  λi = λ0
  if b0 >= i.budget
    x⁺ = st0
    @assert _budgeted_spanning_tree_compute_weight(i, x⁺) >= i.budget

    stalling = false
    while true
      # Penalise less the constraint: it should no more be satisfied.
      λi *= ζ⁻
      _, sti = st_prim_budgeted_lagrangian(i, λi)
      if _budgeted_spanning_tree_compute_weight(i, sti) < i.budget
        x⁻ = sti
        break
      end

      # Is the process stalling?
      if λi <= stalling⁻
        stalling = true
        break
      end
    end

    # Specific handling of stallings.
    if stalling # First test: don't penalise the constraint at all.
      _, sti = st_prim_budgeted_lagrangian(i, 0.0)
      new_budget = _budgeted_spanning_tree_compute_weight(i, sti)
      if new_budget < i.budget
        x⁻ = sti
        stalling = false
      end
    end

    if stalling # Second test: minimise the left-hand side of the budget constraint, in hope of finding a feasible solution.
      # This process is highly similar to the computation of feasible_solution, but with a reverse objective function.
      infeasible_rewards = Dict{Edge{T}, Float64}(e => - i.weights[e] for e in keys(i.rewards))
      infeasible_solution = st_prim(SpanningTreeInstance(i.graph, infeasible_rewards)).tree

      if _budgeted_spanning_tree_compute_weight(i, infeasible_solution) < i.budget
        x⁻ = infeasible_solution
        stalling = false
      end
    end

    if stalling # Third: decide there is no solution strictly below the budget. No refinement is possible.
      # As x⁺ is feasible, return it.
      return SimpleBudgetedSpanningTreeSolution(i, x⁺)
    end
  else
    x⁻ = st0
    @assert _budgeted_spanning_tree_compute_weight(i, x⁻) < i.budget

    while true
      # Penalise more the constraint: it should become satisfied at some point.
      λi *= ζ⁺
      _, sti = st_prim_budgeted_lagrangian(i, λi)
      if _budgeted_spanning_tree_compute_weight(i, sti) >= i.budget
        x⁺ = sti
        break
      end

      # Is the process stalling? If so, reuse feasible_solution, which is guaranteed to be feasible.
      if λi >= λmax
        x⁺ = feasible_solution.tree
        break
      end
    end
  end

  # Normalise the solutions: the input graph is undirected, the direction of the edges is not important.
  # In case one solution has the edge v -> w and the other one w -> v, make them equal. s
  sort_edge(e::Edge{T}) where T = (src(e) < dst(e)) ? e : reverse(e)
  x⁺ = [sort_edge(e) for e in x⁺]
  x⁻ = [sort_edge(e) for e in x⁻]

  # Iterative refinement. Stop as soon as there is a difference of at most one edge between the two solutions.
  while _solution_symmetric_difference_size(x⁺, x⁻) > 2
    # Enforce the loop invariant.
    @assert x⁺ !== nothing
    @assert x⁻ !== nothing
    @assert _budgeted_spanning_tree_compute_weight(i, x⁺) >= i.budget # Feasible.
    @assert _budgeted_spanning_tree_compute_weight(i, x⁻) < i.budget # Infeasible.

    # Switch elements from one solution to another.
    only_in_x⁺, only_in_x⁻ = _solution_symmetric_difference(x⁺, x⁻)
    e1 = first(only_in_x⁺)
    e2 = first(only_in_x⁻)

    # Create the new solution (don't erase x⁺ nor x⁻: only one of the two will be forgotten, the other will be kept).
    new_x = copy(x⁺)
    filter!(e -> e != e1, new_x)
    push!(new_x, e2)

    # Replace one of the two solutions, depending on whether this solution is feasible (x⁺) or not (x⁻).
    if _budgeted_spanning_tree_compute_weight(i, new_x) >= i.budget
      x⁺ = new_x
    else
      x⁻ = new_x
    end
  end

  # Done!
  return SimpleBudgetedSpanningTreeSolution(i, x⁺)
end


text:
module IDebug

export @debug, @bp
using JSON

macro bp(ex)
    return :($(esc(ex)))
end

function nbinclude_string(path::AbstractString)
    codes = []
    path = abspath(path)

    nb = open(JSON.parse, path, "r")

    # check for an acceptable notebook:
    nb["nbformat"] == 4 || error("unrecognized notebook format ", nb["nbformat"])
    lang = lowercase(nb["metadata"]["language_info"]["name"])
    lang == "julia" || error("notebook is for unregognized language $lang")

    shell_or_help = r"^\s*[;?]" # pattern for shell command or help

    for cell in nb["cells"]
        if cell["cell_type"] == "code" && !isempty(cell["source"])
            s = join(cell["source"])
            isempty(strip(s)) && continue # Jupyter doesn't number empty cellsf
            ismatch(shell_or_help, s) && continue
            ismatch(r"IDebug", s) && continue
            ismatch(r"@debug", s) && continue
            push!(codes,s)
        end
    end
    code = join(codes, '\n')
    lines = []
    for (line, str) in enumerate(split(code,'\n'))
        ismatch(r"@bp", str) && push!(lines, line)
    end
    return code, lines
end

macro debug(ex...)
	println("using Gallium, IDebug")
    ex = quote
    cd($(pwd()))
    code, lines = IDebug.nbinclude_string($(ex[1]))
    write("test.jl", code)

    for line in lines
        Gallium.breakpoint("test.jl", line)
    end
    include("test.jl")
    $(ex[2])
    end
    println(ex)
end

end # end of module


text:
module BenchmarkRANSAC

using Dates
using Logging
using Random

using BenchmarkTools
using Tables
using PrettyTables

using RANSAC

export runbenchmark, savebenchmark
export loadbenchmarks
export printresult, saveresult
export info
export printload

const table_header = ["date" "commit sha" "minimum time" "median time" "mean time" "maximum time" "allocs" "memory" "system"; "" "" "[s]" "[s]" "[s]" "[s]" "" "[MiB]" ""]
const md_table_header = ["date" "commit sha" "minimum time [s]" "median time [s]" "mean time [s]" "maximum time [s]" "allocs" "memory [MiB]" "system"]
const nums = :r
const table_align = [:l, :l, nums, nums, nums, nums, nums, nums, :l]

function setupme(iterations)
    vs, ns, norms4Plot, shape_s = examplepc3()
    pcr = RANSACCloud(vs, ns, 32)
    #p = RANSACParameters{Float64}()
    # plane
    # p = RANSACParameters(p, ϵ_plane=0.3, α_plane=deg2rad(5))
    p_plane = (ϵ=0.3, α=deg2rad(5),)
    # cylidner
    # p = RANSACParameters(p, ϵ_cylinder=0.3, α_cylinder=deg2rad(5))
    p_cyl = (ϵ=0.3, α=deg2rad(5),)
    # sphere
    # p = RANSACParameters(p, ϵ_sphere=0.3, α_sphere=deg2rad(5))
    p_sp = (ϵ=0.3, α=deg2rad(5),)
    # number of minimal subsets drawed in one iteration
	# probability that we found shapes
	# minimum shape size
	# maximum number of iteration
    p_it = (minsubsetN=15, prob_det=0.9, τ=900, itermax=iterations, drawN=3)
    p = ransacparameters(;plane=p_plane, cylinder=p_cyl, sphere=p_sp, iteration=p_it)
    return (pcr, p, true)
end

function makenamedtuple(bmresult)
    comsha = read(`git log -n 1 --pretty=format:"%H"`, String)
    pc = getpc()
    mint = BenchmarkTools.minimum(bmresult).time
    maxt = BenchmarkTools.maximum(bmresult).time
    meant = BenchmarkTools.mean(bmresult).time
    mediant = BenchmarkTools.median(bmresult).time
    alld = allocs(bmresult)
    mmm = memory(bmresult)/1024
    (date=Dates.now(), commitsha=comsha, minimumtime=mint, mediantime=mediant, meantime=meant, maximumtime=maxt, allocated=alld, memory=mmm, system = pc)
end

function getpc()
    if Sys.iswindows()
        username = ENV["UserName"]
        if username == "cstamas"
            return "WorkLaptop"
        elseif username == "Ipse"
            return "HomeLaptop"
		elseif username == "Laci"
            return "HomePC"
        else
            return "unknownWin"
        end
    else
        return "notWin"
    end
end

"""
    runbenchmark(show = true; savetocsv = false, showdebug = false)

Run benchmark.
"""
function runbenchmark(show = true; savetocsv = false, showdebug = false)
    glb = global_logger()
    showdebug ? global_logger(nosource_debuglogger()) : global_logger(nosource_infologger())

    Random.seed!(1234);
    sharp = setupme(100_000)
    @info "Running real benchmark..."
    benched = @benchmark ransac($sharp..., reset_rand=true)
    @info "Benchmark finished."

    benchedtuple = makenamedtuple(benched)
    bmi = Tables.rowtable([benchedtuple])
    show && display(benched)
    savetocsv && savebenchmark(bmi)
    global_logger(glb)
    bmi, benched
end

function nicifyone(onel::Array{T}, divisor) where T<:NamedTuple
    onel = onel[1]
    shas = getproperty(onel, :commitsha)[1:7]
    mint = getproperty(onel, :minimumtime)/divisor
    maxt = getproperty(onel, :maximumtime)/divisor
    meant = getproperty(onel, :meantime)/divisor
    mediant = getproperty(onel, :mediantime)/divisor
    mem = getproperty(onel, :memory)/1024 # KiB to MiB
    sys = getproperty(onel, :system)
    [getproperty(onel, :date) shas mint mediant meant maxt getproperty(onel, :allocated) mem sys]
end

function nicifyone(onel, divisor)
    shas = getproperty(onel, :commitsha)[1][1:7]
    mint = getproperty(onel, :minimumtime)[1]/divisor
    maxt = getproperty(onel, :maximumtime)[1]/divisor
    meant = getproperty(onel, :meantime)[1]/divisor
    mediant = getproperty(onel, :mediantime)[1]/divisor
    mem = getproperty(onel, :memory)[1]/1024 # KiB to MiB
    sys = getproperty(onel, :system)
    [getproperty(onel, :date)[1] shas mint mediant meant maxt getproperty(onel, :allocated)[1] mem sys]
end

function nicify(tb)
    divisor = 1_000_000_000 # nanosec to sec
    if size(tb, 1) < 2
        return nicifyone(tb, divisor)
    end
    shas = [ prop[1:7] for prop in getproperty(tb, :commitsha)]
    mint = [ m/divisor for m in getproperty(tb, :minimumtime)]
    maxt = [ m/divisor for m in getproperty(tb, :maximumtime)]
    meant = [ m/divisor for m in getproperty(tb, :meantime)]
    mediant = [ m/divisor for m in getproperty(tb, :mediantime)]
    mem = [ m/1024 for m in getproperty(tb, :memory)] # KiB to MiB
    sys = getproperty(tb, :system)
    [getproperty(tb, :date) shas mint mediant meant maxt getproperty(tb, :allocated) mem sys]
end

function printresult(tb)
    hmnice = nicify(tb)
    tf = PrettyTableFormat(unicode)
    pretty_table(hmnice, table_header, tf, alignment = table_align, formatter=ft_round(3, [3,4,5,6,8]))
end

function info()
    @info "Runnnig one benchmark takes around 5 minutes."
    @info "`bmark, benched = runbenchmark();` to run the benchmark (and also display it)."
    @info "`printresult(bmark)` to show the result of one or more benchmarks."
end

end # module

using .BenchmarkRANSAC
info()


text:
using SparseArrays
using Random

# split from runtests.jl for docs generating script (gen_imgs.jl)

# sprand or sprandn is not stable across versions (e.g. 1.0 vs 1.6)
function _stable_sprand(r, m::Integer, n::Integer, density::AbstractFloat)
    I = Int[]
    J = Int[]
    for li in randsubseq(r, 1:(m*n), density)
        j, i = divrem(li - 1, m)
        push!(I, i + 1)
        push!(J, j + 1)
    end
    V = rand(r, length(I))
    sparse(I, J, V)
end


text:
using FiSUDiPPuS
using FileIO, JLD2
using Test

include("gendata.jl")

@testset "FiSUDiPPuS.jl" begin
    settings_file = joinpath(@__DIR__, "../data/default.jl")
    r = runfit(settings_file; saveprefix="test")
    @test round(r.minimizer[1], digits=6) ≈ A_ssp[1]
    @test round(r.minimizer[4], digits=6) ≈ A_ppp[1]
    @test round(r.minimizer[end-1], digits=6) ≈ a_pow

    # run the other exported functions on the data to check if they error
    # get latest results
    result_file = readdir(joinpath(@__DIR__, "../data/results/"))[end]
    result_file = joinpath("../data/results", result_file)
    plotresult(result_file)
    printresult(result_file)
    viewsettings(settings_file)
end


text:
f_to_c_offset = 32.0
f_to_c_factor = 0.555555555
temp_c = 0.0

function fahrenheit_to_celsius(temp_f)
    global temp_c
    temp_c = (temp_f - f_to_c_offset) * f_to_c_factor
end

fahrenheit_to_celsius(100.0)
println(temp_c)


text:
Random.seed!(2)

mud = 0.01
mup = 10.0
muneg = 0.1

its = 10^7
nmutsp = zeros(Int64, its)
nmutsd = zeros(Int64, its)
nmutsneg = zeros(Int64, its)
global Rmax = log(2)
sfunc() = 0.1
for i in 1:its
  global c = CancerSeqSim.cancercellM(Int64[], Int64[], Int64[], log(2), 0.0,log(2), 0.0, Float64[], Float64[], mup, mud, muneg, Float64[])
  global mutIDp, mutIDd, mutIDneg = 1, 1, 1
  global c, mutIDp, mutIDd, mutIDneg, Rmax = CancerSeqSim.newmutations(c, mutIDp, mutIDd, mutIDneg, Rmax, 1.0, sfunc, CancerSeqSim.multiplicativefitness)
  nmutsp[i] = length(c.mutationsp)
  nmutsd[i] = length(c.mutationsd)
  nmutsneg[i] = length(c.mutationsneg)
end

@test isapprox(mean(nmutsp), mup, rtol = 1e-2)
@test isapprox(mean(nmutsd), mud, rtol = 1e-2)
@test isapprox(mean(nmutsneg), muneg, rtol = 1e-2)

global mud = 0.1
global mup = 10.0
global muneg = 0.5
global its = 10^5
global Ndivs = 5
nmutsp = zeros(Int64, its)
nmutsd = zeros(Int64, its)
nmutsneg = zeros(Int64, its)
global Rmax = log(2)
for i in 1:its
  global c = CancerSeqSim.cancercellM(Int64[], Int64[], Int64[], log(2), 0.0,log(2), 0.0, Float64[], Float64[], mup, mud, muneg, Float64[])
  global mutIDp, mutIDd, mutIDneg = 1, 1, 1
  for j in 1:Ndivs
      global c, mutIDp, mutIDd, mutIDneg, Rmax = CancerSeqSim.newmutations(c, mutIDp, mutIDd, mutIDneg, Rmax, 1.0, sfunc, CancerSeqSim.multiplicativefitness)
  end
  nmutsp[i] = length(c.mutationsp)
  nmutsd[i] = length(c.mutationsd)
  nmutsneg[i] = length(c.mutationsneg)
end

@test isapprox(mean(nmutsp), mup * Ndivs, rtol = 1e-2)
@test isapprox(mean(nmutsd), mud * Ndivs, rtol = 1e-2)

global mud = 1.0
global mup = 10.0
global muneg = 1.0
global Rmax = log(2)
global c = CancerSeqSim.cancercellM(Int64[], Int64[], Int64[], log(2), 0.0, log(2), 0.0, Float64[], Float64[], mup, mud, muneg, Float64[])
global mutIDp, mutIDd, mutIDneg = 1, 1, 1
global c, mutIDp, mutIDd, mutIDneg, Rmax = CancerSeqSim.newmutations(c, mutIDp, mutIDd, mutIDneg, Rmax, 1.0, sfunc, CancerSeqSim.multiplicativefitness)
@test isapprox(prod(c.fitness .+ 1) .* prod(c.fitnessneg .+ 1) * log(2), c.b)
@test Rmax == c.b .+ c.d

for i in 1:10
  global c, mutIDp, mutIDd, mutIDneg, Rmax = CancerSeqSim.newmutations(c, mutIDp, mutIDd, mutIDneg, Rmax, 1.0, sfunc, CancerSeqSim.multiplicativefitness)
  println(c.b)
end
@test isapprox(prod(c.fitness .+ 1) .* prod(c.fitnessneg .+ 1) .* log(2), c.b)

@time cells, tvec, rm1 = CancerSeqSim.tumourgrow_birthdeath(log(2), 0.0, 10^3, 10.0, 0.0, 0.0)
@test rm1 == log(2)

@time cells, tvec, rm1 = CancerSeqSim.tumourgrow_birthdeath(log(2), 0.0, 10^3, 10.0, 0.1, 0.0)
@test rm1 == maximum(map(x -> x.b, cells))

IP = CancerSeqSim.InputParametersM(
10^2,
0.1,
2,
200.0,
200,
10,
0.001,
0.0,
log(2),
0.0,
0.0,
1.0,
0.1,
CancerSeqSim.exptime,
CancerSeqSim.nonmultiplicativefitness)
#get simulation data
simresult = CancerSeqSim.run1simulation(IP)

# if s = 0.0, mean VAF of passengers and drivers should be equal
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.0
for i in 1:100
  x = simulatedifferentmutations(Nmax = 10^3, μp = 10.0, μd = 10.0, clonalmutations = 0, s = mysfunc)
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@time isapprox(mean(meand), mean(meanp), rtol = 1e-2)

# if s >0.0, mean VAF of drivers > passengers
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.1
for i in 1:100
  x = simulatedifferentmutations(Nmax = 10^3, μp = 0.1, μd = 0.01, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio1 = mean(meand) / mean(meanp)

#if s is large above ratio should be larger
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.25
for i in 1:100
  x = simulatedifferentmutations(Nmax = 10^3, μp = 0.1, μd = 0.01, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio2 = mean(meand) / mean(meanp)

#if s is large above ratio should be larger
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.5
for i in 1:100
  x = simulatedifferentmutations(Nmax = 10^3, μp = 0.1, μd = 0.01, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio3 = mean(meand) / mean(meanp)

#if s is large above ratio should be larger
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.75
for i in 1:100
  x = simulatedifferentmutations(Nmax = 10^3, μp = 0.1, μd = 0.01, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio4 = mean(meand) / mean(meanp)

@test ratio4 > ratio3 > ratio2 > ratio1

###
#Moran model
###
μp = 10.0
μd = 0.01
μneg = 0.0
b = log(2)
d = 0.0
Nmax = 10^3
maxt = 5.0
t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg = CancerSeqSim.initializesim(μp, μd, μneg, b, d, Nmax)

c, t, rm1 = CancerSeqSim.tumourmoran(b, d, Nmax, μp, μd, μneg, maxt)
Mp, Md, Mneg, fitness, tvec, cells, timedrivers = CancerSeqSim.getresults(b, d, μp, μd, μneg, Nmax, maxt)

mysfunc() = 0.1
tmax = 20.0
@time x = simulatedifferentmutationsmoran(tmax; Nmax = 10^3, μp = 0.1, μd = 0.1, clonalmutations = 0, s = mysfunc)

#if s is large above ratio should be larger
tmax = 20.0
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.0
for i in 1:100
  global x = simulatedifferentmutationsmoran(tmax, Nmax = 100, μp = 0.1, μd = 0.05, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@time isapprox(mean(meand), mean(meanp), rtol = 1e-2)
ratio1 = mean(meand) / mean(meanp)

tmax = 20.0
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.1
for i in 1:100
  global x = simulatedifferentmutationsmoran(tmax, Nmax = 100, μp = 0.1, μd = 0.05, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio1 = mean(meand) / mean(meanp)

tmax = 20.0
meanp = Float64[]
meand = Float64[]
mysfunc() = 0.2
for i in 1:100
  global x = simulatedifferentmutationsmoran(tmax, Nmax = 100, μp = 0.1, μd = 0.05, clonalmutations = 0, s = mysfunc);
  push!(meand, mean(x.output.trueVAFd))
  push!(meanp, mean(x.output.trueVAFp))
end
@test mean(meand) > mean(meanp)
ratio2 = mean(meand) / mean(meanp)

@test ratio2 > ratio1


text:
# ============================================================
# Isothermal Wall Functions
# ============================================================

function bc_isothermal!(
    ctr::T,
    ctr1::T,
    γ::RN,
    λw = 1.0::RN,
) where {T<:Union{ControlVolume1D,ControlVolume2D,ControlVolumeUS}}

    ctr.prim[2] = -ctr1.prim[2]
    if length(ctr.prim) == 4
        ctr.prim[3] = -ctr1.prim[3]
    elseif length(ctr.prim) == 5
        ctr.prim[4] = -ctr1.prim[4]
    end

    ctr.prim[end] = 2.0 * λw - ctr1.prim[end]
    fac = (ctr1.prim[end] - λw) / λw
    ctr.prim[1] = (1.0 - fac) / (1.0 + fac) * ctr1.prim[1]

    ctr.w .= prim_conserve(ctr.prim, γ)

    return nothing
end