text:
# NOTE: see also interface/univariate_finite.jl

const UnivariateFiniteUnion =
    Union{UnivariateFinite, UnivariateFiniteArray}

# ABOVE NOT EXPORTED!!!!

"""
    classes(d::UnivariateFinite)
    classes(d::UnivariateFiniteArray)

A list of categorial elements in the common pool of classes used to
construct `d`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    classes(d) # CategoricalArray{String,1,UInt32}["maybe", "no", "yes"]
"""
MMI.classes(d::UnivariateFiniteUnion) = d.decoder.classes

"""
    levels(d::UnivariateFinite)

A list of the raw levels in the common pool of classes used to
construct `d`, equal to
`CategoricalArrays.DataAPI.unwrap.(classes(d))`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    levels(d) # Array{String, 1}["maybe", "no", "yes"]

"""
levels(d::UnivariateFinite)  = CategoricalArrays.DataAPI.unwrap.(classes(d))

function Distributions.params(d::UnivariateFinite)
    raw = raw_support(d) # reflects order of pool at instantiation of d
    pairs = tuple([unwrap.(d.decoder(r))=>d.prob_given_ref[r] for r in raw]...)
    levs = unwrap.(classes(d))
    return (levels=levs, probs=pairs)
end

# get the internal integer representations of the support
raw_support(d::UnivariateFiniteUnion) = collect(keys(d.prob_given_ref))

"""
    Distributions.support(d::UnivariateFinite)
    Distributions.support(d::UnivariateFiniteArray)

Ordered list of classes associated with non-zero probabilities.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    support(d) # CategoricalArray{String,1,UInt32}["maybe", "no"]

"""
Distributions.support(d::UnivariateFiniteUnion) =
    map(d.decoder, raw_support(d))

sample_scitype(d::UnivariateFiniteUnion) = d.scitype

CategoricalArrays.isordered(d::UnivariateFinite) = isordered(classes(d))
CategoricalArrays.isordered(u::UnivariateFiniteArray) = isordered(classes(u))

## DISPLAY

function Base.show(stream::IO, d::UnivariateFinite)
    pairs = Distributions.params(d).probs
    x1, p1 = pairs[1]
    str = "UnivariateFinite{$(d.scitype)}($x1=>$(round(p1, sigdigits=3))"
    for pair in pairs[2:end]
        str *= ", $(pair[1])=>$(round(pair[2], sigdigits=3))"
    end
    str *= ")"
    print(stream, str)
end

show_prefix(u::UnivariateFiniteArray{S,V,R,P,1}) where {S,V,R,P} =
    "$(length(u))-element"
show_prefix(u::UnivariateFiniteArray) = join(size(u),'x')

# function show_compact(io::IO, v)
#     K = 4
#     n = length(v)
#     if !isempty(v)
#         indent = div(length(string(v[1])), 2)
#         L = min(n, K - 1)
#         for i in 1:min(n, L - 1)
#             println(io, " ", v[i])
#         end
#         print(io, " ", v[L])
#         if n > K
#             println(io, "\n", fill(" ", indent)..., " ⋮")
#         else
#             println(io)
#         end
#         if n > K - 1
#             println(io, " ", last(v))
#         end
#     end
# end

# function Base.show(io::IO,
#                    u::UnivariateFiniteArray{S,V,R,P,N}) where {S,V,R,P,N}
#     print(io, show_prefix(u), " UnivariateFiniteArray{$S,$V,$R,$P,$N}(...)")
# end

# function Base.show(io::IO, m::MIME"text/plain",
#                    u::UnivariateFiniteArray{S,V,R,P,1}) where {S,V,R,P}
#     support = unwrap.(Dist.support(u))
#     print(io, show_prefix(u), " UnivariateFiniteArray{$S,$V,$R,$P,1}")
#     if !isempty(u)
#         println(io, ":")
#         show_compact(io, u)
#     end
# end

"""
    isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)

Returns `true` if and only if `d1` and `d2` have the same support and
the corresponding probabilities are approximately equal. The key-word
arguments `kwargs` are passed through to each call of `isapprox` on
probability pairs.  Returns `false` otherwise.

"""
function Base.isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)
    support1 = Dist.support(d1)
    support2 = Dist.support(d2)
    for c in support1
        c in support2 || return false
        isapprox(pdf(d1, c), pdf(d2, c); kwargs...) ||
            return false # pdf defined below
    end
    return true
end
function Base.isapprox(d1::UnivariateFiniteArray,
                       d2::UnivariateFiniteArray; kwargs...)
    support1 = Dist.support(d1)
    support2 = Dist.support(d2)
    for c in support1
        c in support2 || return false
        isapprox(pdf.(d1, c), pdf.(d2, c); kwargs...) ||
            return false
    end
    return true
end

function average(dvec::AbstractVector{UnivariateFinite{S,V,R,P}};
                 weights=nothing) where {S,V,R,P}

    n = length(dvec)

    Dist.@check_args(UnivariateFinite, weights == nothing || n==length(weights))

    # check all distributions have consistent pool:
    first_index = first(dvec).decoder.classes
    for d in dvec
        d.decoder.classes == first_index ||
            error("Averaging UnivariateFinite distributions with incompatible"*
                  " pools. ")
    end

    # get all refs:
    refs = reduce(union, [keys(d.prob_given_ref) for d in dvec]) |> collect

    # initialize the prob dictionary for the distribution sum:
    prob_given_ref = LittleDict{R,P}([refs...], zeros(P, length(refs)))

    # make vector of all the distributions dicts padded to have same common keys:
    prob_given_ref_vec = map(dvec) do d
        merge(prob_given_ref, d.prob_given_ref)
    end

    # sum up:
    if weights == nothing
        scale = 1/n
        for x in refs
            for k in 1:n
                prob_given_ref[x] += scale*prob_given_ref_vec[k][x]
            end
        end
    else
        scale = 1/sum(weights)
        for x in refs
            for k in 1:n
                prob_given_ref[x] +=
                    weights[k]*prob_given_ref_vec[k][x]*scale
            end
        end
    end
    d1 = first(dvec)
    return UnivariateFinite(sample_scitype(d1), d1.decoder, prob_given_ref)
end

"""
    Distributions.pdf(d::UnivariateFinite, x)

Probability of `d` at `x`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    pdf(d, "yes")     # 0.3
    pdf(d, v[1])      # 0.3
    pdf(d, "no")      # 0.0
    pdf(d, "house")   # throws error

Other similar methods are available too:

    mode(d)    # CategoricalValue{String, UInt32} "maybe"
    rand(d, 5) # CategoricalArray{String,1,UInt32}["maybe", "no", "maybe", "maybe", "no"] or similar
    d = fit(UnivariateFinite, v)
    pdf(d, "maybe") # 0.25
    logpdf(d, "maybe") # log(0.25)

One can also do weighted fits:

    w = [1, 4, 5, 1] # some weights
    d = fit(UnivariateFinite, v, w)
    pdf(d, "maybe") ≈ 4/11 # true

See also `classes`, `support`.
"""
function Distributions.pdf(
    d::UnivariateFinite{S,V,R,P},
    cv::CategoricalValue,
) where {S,V,R,P}
    return get(d.prob_given_ref, int(cv), zero(P))
end
Distributions.pdf(d::UnivariateFinite{S,V}, c::V) where {S,V} = _pdf(d, c)

# Avoid method ambiguity errors with Distributions >= 0.24
Distributions.pdf(d::UnivariateFinite{S,V}, c::V) where {S,V<:Real} = _pdf(d, c)

function _pdf(d::UnivariateFinite, c)
    _classes = classes(d)
    c in _classes || throw(DomainError("Value $c not in pool. "))
    pool = CategoricalArrays.pool(_classes)
    class = pool[get(pool, c)]
    return pdf(d, class)
end

Distributions.logpdf(d::UnivariateFinite, cv::CategoricalValue) = log(pdf(d,cv))
Distributions.logpdf(d::UnivariateFinite{S,V}, c::V) where {S,V} = log(pdf(d, c))

# Avoid method ambiguity errors with Distributions >= 0.24
Distributions.logpdf(d::UnivariateFinite{S,V}, c::V) where {S,V<:Real} = log(pdf(d, c))

function Distributions.mode(d::UnivariateFinite)
    dic = d.prob_given_ref
    p = values(dic)
    max_prob = maximum(p)
    m = first(first(dic)) # mode, just some ref for now
    for (x, prob) in dic
        if prob == max_prob
            m = x
            break
        end
    end
    return d.decoder(m)
end

# mode(v::Vector{UnivariateFinite}) = mode.(v)
# mode(u::UnivariateFiniteVector{2}) =
#     [u.support[ifelse(s > 0.5, 2, 1)] for s in u.scores]
# mode(u::UnivariateFiniteVector{C}) where {C} =
#     [u.support[findmax(s)[2]] for s in eachrow(u.scores)]

"""
    _cumulative(d::UnivariateFinite)

Return the cumulative probability vector `[0, ..., 1]` for the
distribution `d`, using only classes in the support of `d`, ordered
according to the categorical elements used at instantiation of
`d`. Used only to implement random sampling from `d`.

"""
function _cumulative(d::UnivariateFinite{S,V,R,P}) where {S,V,R,P<:Real}

    # the keys of `d` are in order; see constructor
    p = collect(values(d.prob_given_ref))
    K = length(p)
    p_cumulative = Array{P}(undef, K + 1)
    p_cumulative[1] = zero(P)
    p_cumulative[K + 1] = one(P)
    for i in 2:K
        p_cumulative[i] = p_cumulative[i-1] + p[i-1]
    end
    return p_cumulative
end

"""
_rand(rng, p_cumulative, R)

Randomly sample the distribution with discrete support `R(1):R(n)`
which has cumulative probability vector `p_cumulative=[0, ..., 1]` (of
length `n+1`). Does not check the first and last elements of
`p_cumulative` but does not use them either.

"""
function _rand(rng, p_cumulative, R)
    real_sample = rand(rng)
    K = R(length(p_cumulative))
    index = K
    for i in R(2):R(K)
        if real_sample < p_cumulative[i]
            index = i - R(1)
            break
        end
    end
    return index
end

function Base.rand(rng::AbstractRNG,
                   d::UnivariateFinite{<:Any,<:Any,R}) where R
    p_cumulative = _cumulative(d)
    return Distributions.support(d)[_rand(rng, p_cumulative, R)]
end

function Base.rand(rng::AbstractRNG,
                   d::UnivariateFinite{<:Any,<:Any,R},
                   dim1::Int, moredims::Int...) where R # ref type
    p_cumulative = _cumulative(d)
    A = Array{R}(undef, dim1, moredims...)
    for i in eachindex(A)
        @inbounds A[i] = _rand(rng, p_cumulative, R)
    end
    support = Distributions.support(d)
    return broadcast(i -> support[i], A)
end

rng(d::UnivariateFinite, args...) = rng(Random.GLOBAL_RNG, d, args...)

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C}) where C
    C <: CategoricalValue ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "`CategoricalValue type`. ")
    y = skipmissing(v) |> collect
    isempty(y) && error("No non-missing data to fit. ")
    N = length(y)
    count_given_class = Dist.countmap(y)
    classes = Tuple(keys(count_given_class))
    probs = values(count_given_class)./N
    prob_given_class = LittleDict(classes, probs)
    return UnivariateFinite(prob_given_class)
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C},
                           weights::Nothing) where C
    return Distributions.fit(d, v)
end

function Distributions.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C},
                           weights::AbstractVector{<:Real}) where C
    C <: CategoricalValue ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "`CategoricalValue` type. ")
    y = broadcast(identity, skipmissing(v))
    isempty(y) && error("No non-missing data to fit. ")
    classes_seen = filter(in(unique(y)), classes(y[1]))

    # instantiate and initialize prob dictionary:
    prob_given_class = LittleDict{C,Float64}()
    for c in classes_seen
        prob_given_class[c] = 0
    end

    # compute unnormalized  probablilities:
    for i in eachindex(y)
        prob_given_class[y[i]] += weights[i]
    end

    # normalize the probabilities:
    S = sum(values(prob_given_class))
    for c in keys(prob_given_class)
        prob_given_class[c] /=S
    end

    return UnivariateFinite(prob_given_class)
end

## NORMAL DISTRIBUTION ARITHMETIC

# *(lambda::Number, d::Dist.Normal) = Dist.Normal(lambda*d.μ, lambda*d.σ)
# function +(ds::Dist.Normal...)
#     μ = sum([d.μ for d in ds])
#     σ = sum([d.σ^2 for d in ds]) |> sqrt
#     return Dist.Normal(μ, σ)
# end


text:
function get_initial_covariance_matrix(obj::CovarianceMatrixProvider, arg0::SpacecraftState)
    return jcall(obj, "getInitialCovarianceMatrix", RealMatrix, (SpacecraftState,), arg0)
end

ext: jl


text:
module FileTransportTests

using Thrift
using Compat
using Base.Test

function testfiletransport()
    println("\nTesting file transport...")

    fname = tempname()
    s1 = "Hello World"
    open(fname, "w") do f
        t = TFileTransport(f)
        p = TCompactProtocol(t)
        for i in 1:5
            write(p, s1)
        end
    end

    open(fname, "r") do f
        t = TFileTransport(f)
        p = TCompactProtocol(t)
        for i in 1:5
            s2 = read(p, typeof(s1))
            @test s2 == s1
        end
        println("passed.")
    end
    rm(fname)
end

testfiletransport()

end


text:
## Main algorithm of AMRVW
## This follows that given in the paper very closely
## should trace algorithm better with `verbose`
function AMRVW_algorithm!(state::AbstractQRRotatorFactorization)
    m = eltype(state.QF) <: Real ? 2 : 1
    AMRVW_algorithm!(state, m)
end

function AMRVW_algorithm!(state::AbstractQRRotatorFactorization, m)

    QF, RF = state.QF, state.RF
    S = eltype(QF)
    T = real(S)

    # set counters
    n = length(QF.Q)
    ctr = AMRVW_Counter(0, 1, n, IT_COUNT, n-1)

    # set up storage
    storage = make_storage(QF, m)
    A = storage.A  # cleanup
    EIGS = zeros(Complex{T}, n+1)

    it_max = 20*n
    kk = 0

    check_deflation(state, ctr)

    while kk <= it_max

        if ctr.stop_index <= 0
            break
        end

        check_deflation(state, ctr)

        ##show_status(state, ctr)

        kk += 1

        ## Delta is number of rotators in piece being considered
        ## m must be greater than delta-1 to run a bulge step
        k = ctr.stop_index
        delta = ctr.stop_index - ctr.zero_index  ## number of rotators in Q factorization

        if delta == 1

            diagonal_block!(A, QF, RF, k, k+1)
            e1::Complex{T}, e2::Complex{T} = eigen_values(A)

            EIGS[k] = e1
            EIGS[k+1] = e2

            # can finish up if near end
            if ctr.stop_index == 2
                diagonal_block!(A, QF, RF, 1, 2)
                e1 = A[1,1]
                EIGS[1] = e1

                break
            end

            ctr.stop_index = ctr.zero_index - 1
            ctr.zero_index = 0
            ctr.start_index = 1
            check_deflation(state, ctr)

            ctr.stop_index <= 0 && break

        elseif delta <= 0

            diagonal_block!(A, QF, RF, k, k+1)
            e1, e2 = A[1,1], A[2,2]  # 0s off diagonal

            if ctr.stop_index == 1
                # finish up
                EIGS[1], EIGS[2] = e1, e2
                ctr.stop_index = 0

                break
            end

            EIGS[k+1] = e2

            ctr.stop_index = ctr.zero_index - 1
            ctr.stop_index <= 0 && break

            ctr.zero_index = 0
            ctr.start_index = 1
            check_deflation(state, ctr)

        elseif m > 2 && m >= delta - 1
            ## When finding shifts we find m eigenvalues from some other means
            ## so here we find the eigvals

            δ, Δ = ctr.start_index, ctr.stop_index
            diagonal_block!(A, QF, RF, δ, Δ+1)

            es = _eigvals(view(A,1:Δ-δ+2,1:Δ-δ+2))
            for (ind, ev) in enumerate(es)
                EIGS[δ+ind-1] = ev
            end

            if  δ == 2
                # finish up
                diagonal_block!(A, QF, RF, 1, 2)

                e1 = A[1,1] # δ == 2 means first rotator a diagonal, so product will be diagonal in 1:2x1:2 block
                EIGS[1] = e1

                break
            end

            ctr.zero_index = 0
            ctr.start_index = 1
            ctr.stop_index = δ - 2

            ctr.stop_index <= 0 && break
            check_deflation(state, ctr)

        else

            bulge_step(QF, RF, storage, ctr, m)
            ctr.it_count -= 1

        end

    end

    LinearAlgebra.sorteig!(EIGS)

    EIGS

end

##################################################
## A container for our counters
mutable struct AMRVW_Counter
    zero_index::Int  # this should be dropped, it is just start_index-1
    start_index::Int
    stop_index::Int
    it_count::Int    # count on use of random rotators
    tr::Int          # count for simple_passthrough
end

function make_counter(state)
    n = length(state.QF.Q)
    AMRVW_Counter(0, 1, n, 0, n-1)
end

##################################################

## Used to  create storage based  on Qfactorization  type
make_storage(state::QRFactorization) = make_storage(state.QF, eltype(state) <: Complex ? 1 : 2)

function make_storage(QF::QFactorization{T,S,V}, m) where {T, S, V}
    A = zeros(S, 4, 4)
    VU = Vector{Rotator{T, S}}(undef, m)   # Ascending Chain
    limb = DescendingChain(copy(QF.Q.x[1:(m-1)]))
    (A=A, VU=VU, limb=limb)
end

function make_storage(QF::QFactorizationTwisted{T,S, Vt,  PVt}, m) where {T, S, Vt, PVt}
    A = zeros(S, max(4, m+2), max(4, m+2))
    VU = Vector{Rotator{T, S}}(undef, m)   # Ascending Chain
    limb = TwistedChain(copy(QF.Q.x[1:(m-1)]), copy(QF.Q.pv[1:(m-2)]))
    (A=A, VU=VU, limb=limb)
end

##################################################
## Deflation
## when a Q[k] matrix becomes a diagonal matrix, we deflate.
## This is checked by the sine term being basically 0.
function check_deflation(state::AbstractQRRotatorFactorization, ctr)

    QF = state.QF
    tol = eps(real(eltype(QF)))

    for k in ctr.stop_index:-1:ctr.start_index

        c, s = vals(QF.Q[k])

        if abs(s) <= tol

            deflate(QF, k, ctr)

            ctr.zero_index = k      # points to a matrix Q[k] either
            ctr.start_index = k + 1
            ctr.it_count = 15        # reset counter

            return
        end
    end

end

## Find eigenvalues from factorization
function eigvals(state::AbstractQRRotatorFactorization)

    es = AMRVW_algorithm!(copy(state))
    LinearAlgebra.sorteig!(es)
    es

end


text:
using OrdinaryDiffEq, Test, DiffEqDevTools, DiffEqOperators
u0 = rand(2)
A = DiffEqArrayOperator([2.0 -1.0; -1.0 2.0])
function (p::typeof(A))(::Type{Val{:analytic}},u0,p,t)
    exp(p.A*t)*u0
end

prob = ODEProblem(A,u0,(0.0,1.0))

x = rand(2)
@test A(0.0,x) == A*x

sol = solve(prob,LinearImplicitEuler())

dts = 1./2.^(8:-1:4) #14->7 good plot
sim  = test_convergence(dts,prob,LinearImplicitEuler())
@test abs(sim.𝒪est[:l2]-1) < 0.2

# using Plots; pyplot(); plot(sim)

B = ones(2)
L = AffineDiffEqOperator{Float64}((A,),(B,),rand(2))
prob = ODEProblem(L,u0,(0.0,1.0))
sol = solve(prob,LinearImplicitEuler())

B = DiffEqArrayOperator(ones(2,2))
L = AffineDiffEqOperator{Float64}((A,B),(),rand(2))
function (p::typeof(L))(::Type{Val{:analytic}},u0,p,t)
    exp((p.As[1].A+p.As[2].A)*t)*u0
end

# Midpoint splitting

prob = ODEProblem(L,u0,(0.0,1.0))
sol = solve(prob,MidpointSplitting(),dt=1/10)
# using Plots; pyplot; plot(sol)

## Midpoint splitting convergence
##
## We use the inhomogeneous Lorentz equation for an electron in a
## time-dependent field. To write this on matrix form and simplify
## comparison with the analytic solution, we introduce two dummy
## variables:
## 1) As the third component, a one is stored to allow the
##    inhomogeneous part to be expressed on matrix form.
## 2) As the fourth component, the initial time t_i is stored,
##    for use by the analytical formula.
## This wastes a lot of space, but simplifies the error analysis.
##
## We can then write the Lorentz equation as q̇ = [A + f(t)B]q.

f = t -> -sin(2pi*t)
F = t -> cos(2pi*t)/2pi # Primitive function of f(t)

A = DiffEqArrayOperator([0 1 0 0
                         0 0 0 0
                         0 0 0 0
                         0 0 0 0])

B = DiffEqArrayOperator([0 0 0 0
                         0 0 1 0
                         0 0 0 0
                         0 0 0 0], f)

H = AffineDiffEqOperator{Float64}((A,B),(),rand(4))
function (p::typeof(H))(::Type{Val{:analytic}},u0,p,t)
    x0,v0 = u0[1:2]
    ti = u0[end]
    x = x0 + (t-ti)*v0 - (f.(t)-f(ti))/(2pi)^2 - (t-ti)*F(ti)
    v = v0 + (F.(t)-F(ti))
    [x, v, 1, ti]
end

x0,v0,ti = rand(3)
prob = ODEProblem(H, [x0, v0, 1, ti], (ti, 5.))
dts = 1./2.^(10:-1:1)
sim  = test_convergence(dts,prob,MidpointSplitting())
@test abs(sim.𝒪est[:l2]-2) < 0.2


text:
### This support BLAS style multiplication
#           α * A * B + β C
# it represents the same thing as Mul(α, A, B) .+ Mul(β, C)
# but avoids the broadcast machinery

# Lazy representation of α*A*B + β*C
struct MulAdd{StyleA, StyleB, StyleC, T, AA, BB, CC}
    style_A::StyleA
    style_B::StyleB
    style_C::StyleC
    α::T
    A::AA
    B::BB
    β::T
    C::CC
end

function MulAdd(styleA::StyleA, styleB::StyleB, styleC::StyleC, α::T, A::AA, B::BB, β::V, C::CC) where {StyleA,StyleB,StyleC,T,V,AA,BB,CC}
    axes(A,2) == axes(B,1) || throw(DimensionMismatch())
    axes(A,1) == axes(C,1) || throw(DimensionMismatch())
    axes(B,2) == axes(C,2) || throw(DimensionMismatch())
    MulAdd{StyleA,StyleB,StyleC,promote_type(T,V),AA,BB,CC}(styleA, styleB, styleC, α, A, B, β, C)
end

@inline MulAdd(α, A, B, β, C) = MulAdd(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), α, A, B, β, C)

eltype(::MulAdd{StyleA,StyleB,StyleC,T,AA,BB,CC}) where {StyleA,StyleB,StyleC,T,AA,BB,CC} =
     promote_type(T, eltype(AA), eltype(BB), eltype(CC))

size(M::MulAdd, p::Int) = size(M)[p]
axes(M::MulAdd, p::Int) = axes(M)[p]
length(M::MulAdd) = prod(size(M))
size(M::MulAdd) = length.(axes(M))
axes(M::MulAdd) = axes(M.C)

const ArrayMulArrayAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray}
const MatMulVecAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractVector,<:AbstractVector}
const MatMulMatAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}

BroadcastStyle(::Type{<:MatMulVecAdd{StyleA,StyleB,StyleC}}) where {StyleA,StyleB,StyleC} =
    ArrayMulArrayStyle{StyleA,StyleB,2,1}()
BroadcastStyle(::Type{<:MatMulMatAdd{StyleA,StyleB,StyleC}}) where {StyleA,StyleB,StyleC} =
    ArrayMulArrayStyle{StyleA,StyleB,2,2}()

broadcastable(M::MatMulMatAdd) = M
broadcastable(M::MatMulVecAdd) = M

const BlasMatMulVec{StyleA,StyleB,StyleC,T} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractVector{T},<:AbstractVector{T}}
const BlasMatMulMat{StyleA,StyleB,StyleC,T} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractMatrix{T},<:AbstractMatrix{T}}

@inline function _copyto!(_, dest::AbstractArray, M::MulAdd)
    M.C ≡ dest || copyto!(dest, M.C)
    materialize!(MulAdd(M.α, M.A, M.B, M.β, dest))
end

@inline copyto!(dest::AbstractArray, M::MulAdd) = _copyto!(MemoryLayout(dest), dest, M)

const BArrayMulArrayAdd{styleA, styleB, styleC, p, q} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q}, <:Any, typeof(identity),
                <:Tuple{<:ArrayMulArrayAdd{styleA,styleB,styleC}}}

@inline function _copyto!(_, dest::AbstractArray, bc::BArrayMulArrayAdd)
    (M,) = bc.args
    copyto!(dest, M)
end

import LinearAlgebra: tilebufsize, Abuf, Bbuf, Cbuf

# Modified from LinearAlgebra._generic_matmatmul!
function tile_size(T, S, R)
    tile_size = 0
    if isbitstype(R) && isbitstype(T) && isbitstype(S)
        tile_size = floor(Int, sqrt(tilebufsize / max(sizeof(R), sizeof(S), sizeof(T))))
    end
    tile_size
end

function tiled_blasmul!(tile_size, α, A::AbstractMatrix{T}, B::AbstractMatrix{S}, β, C::AbstractMatrix{R}) where {S,T,R}
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    @inbounds begin
        sz = (tile_size, tile_size)
        # FIXME: This code is completely invalid!!!
        Atile = unsafe_wrap(Array, convert(Ptr{T}, pointer(Abuf[Threads.threadid()])), sz)
        Btile = unsafe_wrap(Array, convert(Ptr{S}, pointer(Bbuf[Threads.threadid()])), sz)

        z1 = zero(A[1, 1]*B[1, 1] + A[1, 1]*B[1, 1])
        z = convert(promote_type(typeof(z1), R), z1)

        if mA < tile_size && nA < tile_size && nB < tile_size
            copy_transpose!(Atile, 1:nA, 1:mA, 'N', A, 1:mA, 1:nA)
            copyto!(Btile, 1:mB, 1:nB, 'N', B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    C[i,j] = α*s + β*C[i,j]
                end
            end
        else
            # FIXME: This code is completely invalid!!!
            Ctile = unsafe_wrap(Array, convert(Ptr{R}, pointer(Cbuf[Threads.threadid()])), sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    copyto!(Ctile, 1:ilen, 1:jlen, C, ib:ilim, jb:jlim)
                    lmul!(β,Ctile)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        copy_transpose!(Atile, 1:klen, 1:ilen, 'N', A, ib:ilim, kb:klim)
                        copyto!(Btile, 1:klen, 1:jlen, 'N', B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += α*s
                            end
                        end
                    end
                    copyto!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                end
            end
        end
    end

    C
end

function default_blasmul!(α, A::AbstractMatrix, B::AbstractMatrix, β, C::AbstractMatrix)
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    @inbounds for k = 1:mA, j = 1:nB
        z2 = zero(A[k, 1]*B[1, j] + A[k, 1]*B[1, j])
        Ctmp = convert(promote_type(eltype(C), typeof(z2)), z2)
        @simd for ν = 1:size(A,2)
            Ctmp = muladd(A[k, ν],B[ν, j],Ctmp)
        end
        C[k,j] = α*Ctmp + β*C[k,j]
    end
    C
end

function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    mB = length(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    length(C) == mA || throw(DimensionMismatch("Dimensions must match"))

    lmul!(β, C)
    (nA == 0 || mB == 0)  && return C

    z = zero(A[1]*B[1] + A[1]*B[1])
    Astride = size(A, 1) # use size, not stride, since its not pointer arithmetic

    @inbounds for k = 1:mB
        aoffs = (k-1)*Astride
        b = B[k]
        for i = 1:mA
            C[i] += α * A[aoffs + i] * b
        end
    end

    C
end

function materialize!(M::MatMulMatAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    ts = tile_size(eltype(A), eltype(B), eltype(C))
    if iszero(β) # false is a "strong" zero to wipe out NaNs
        ts == 0 ? default_blasmul!(α, A, B, false, C) : tiled_blasmul!(ts, α, A, B, false, C)
    else
        ts == 0 ? default_blasmul!(α, A, B, β, C) : tiled_blasmul!(ts, α, A, B, β, C)
    end
end

function materialize!(M::MatMulVecAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    default_blasmul!(α, A, B, iszero(β) ? false : β, C)
end

# make copy to make sure always works
@inline function _gemv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.gemv!(tA, α, A, copy(x), β, y)
    else
        BLAS.gemv!(tA, α, A, x, β, y)
    end
end

# make copy to make sure always works
@inline function _gemm!(tA, tB, α, A, B, β, C)
    if B ≡ C
        BLAS.gemm!(tA, tB, α, A, copy(B), β, C)
    else
        BLAS.gemm!(tA, tB, α, A, B, β, C)
    end
end

@inline materialize!(M::BlasMatMulVec{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasFloat}) =
    _gemv!('N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVec{<:AbstractRowMajor,<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasFloat}) =
    _gemv!('T', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVec{<:ConjLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _gemv!('C', M.α, M.A', M.B, M.β, M.C)

@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasFloat}) =
    _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasFloat}) =
    _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.A, M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMat{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasFloat}) =
    _gemm!('T', 'N', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMat{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'N', M.α, M.A', M.B, M.β, M.C)

@inline materialize!(M::BlasMatMulMat{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasFloat}) =
    _gemm!('T', 'T', M.α, transpose(M.A), transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMat{<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('T', 'C', M.α, transpose(M.A), M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMat{<:ConjLayout{<:AbstractRowMajor},<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'T', M.α, M.A', M.B', M.β, M.C)
@inline materialize!(M::BlasMatMulMat{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'C', M.α, M.A', M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractRowMajor,<:BlasFloat}) =
    _gemm!('T', 'T', M.α, M.B, M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'C', M.α, M.B, M.A, M.β, M.C')

@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractRowMajor,<:BlasFloat}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, M.C')
@inline materialize!(M::BlasMatMulMat{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.B', M.A, M.β, M.C')

@inline materialize!(M::BlasMatMulMat{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractRowMajor,<:BlasFloat}) =
    _gemm!('T', 'N', M.α, M.B, transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMat{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'N', M.α, M.B, M.A', M.β, M.C')

@inline materialize!(M::BlasMatMulMat{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractRowMajor,<:BlasFloat}) =
    _gemm!('N', 'N', M.α, transpose(M.B), transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMat{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'N', M.α, M.B', M.A', M.β, M.C')

###
# Symmetric
###

# make copy to make sure always works
@inline function _symv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.symv!(tA, α, A, copy(x), β, y)
    else
        BLAS.symv!(tA, α, A, x, β, y)
    end
end

@inline function _hemv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.hemv!(tA, α, A, copy(x), β, y)
    else
        BLAS.hemv!(tA, α, A, x, β, y)
    end
end

materialize!(M::BlasMatMulVec{<:SymmetricLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasFloat}) =
    _symv!(M.style_A.uplo, M.α, symmetricdata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVec{<:SymmetricLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasFloat}) =
    _symv!(M.style_A.uplo == 'L' ? 'U' : 'L', M.α, transpose(symmetricdata(M.A)), M.B, M.β, M.C)

materialize!(M::BlasMatMulVec{<:HermitianLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(M.style_A.uplo, M.α, hermitiandata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVec{<:HermitianLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(M.style_A.uplo == 'L' ? 'U' : 'L', M.α, hermitiandata(M.A)', M.B, M.β, M.C)

###
# Triangular
###

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractColumnMajor},
                                   <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.factors
    x ≡ dest || copyto!(dest, x)
    BLAS.trmv!(UPLO, 'N', UNIT, triangulardata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractRowMajor},
                                   <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.factors
    x ≡ dest || copyto!(dest, x)
    BLAS.trmv!(UPLO, 'T', UNIT, transpose(triangulardata(A)), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,<:ConjLayout{<:AbstractRowMajor}},
                                   <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.factors
    x ≡ dest || copyto!(dest, x)
    BLAS.trmv!(UPLO, 'C', UNIT, triangulardata(A)', dest)
end

# Triangular *\ Matrix

function _copyto!(_, dest::AbstractMatrix, M::MatMulMat{<:TriangularLayout})
    A,X = M.factors
    size(dest,2) == size(X,2) || thow(DimensionMismatch("Dimensions must match"))
    @views for j in axes(dest,2)
        dest[:,j] .= Mul(A, X[:,j])
    end
    dest
end


text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

# checks
function check_sol(sol, D, H)
    sol.prob.D >= D || Base.error("the required Dth-degree derivative is unavailable")
    H in sol.prob.h || Base.error("the required Hth-harmonic phasor is unavailable")
end

# frequency and angular frequency
F(sol)   = sol.prob.F
F(sol,H) = H * F(sol)
ω(sol,H) = 2.0 * pi * F(sol,H)

# amplitude
"""
    TFT.a(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Shorthand function to obtain the amplitude of the Dth-degree derivative of the
Hth-harmonic phasor, dispatching to `amplitude(sol,D,H)`.
"""
a(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) = amplitude(sol,D,H)

"""
    TFT.amplitude(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Function to obtain the amplitude of the Dth-degree derivative of the
Hth-harmonic phasor.

    ∀ h ∈ {0}:
        a₀⁽ᴰ⁾(t) = ξ₀⁽ᴰ⁾(t)
    ∀ h ∉ {0}:
        aₕ⁽⁰⁾(t) = |2 ⋅ ξₕ⁽⁰⁾(t)| ∈ 𝐑⁺
        aₕ⁽¹⁾(t) = ℜ[2 ⋅ ξₕ⁽¹⁾(t) ⋅ exp(-im ⋅ ϕₕ⁽⁰⁾(t))] ∈ 𝐑
        aₕ⁽²⁾(t) = ℜ[2 ⋅ ξₕ⁽²⁾(t) ⋅ exp(-im ⋅ ϕₕ⁽⁰⁾(t))] + aₕ⁽⁰⁾(t) ⋅ ϕₕ⁽¹⁾(t)² ∈ 𝐑

See: [Assessing Synchrophasor Estimates of an Event Captured by a Phasor
Measurement Unit, pg. 3112](https://ieeexplore.ieee.org/document/9239915)

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `D::Int`                      | degree of the derivative [-], default=0
- `H::Int`                      | harmonic number [-], default=1

Output:
- `a::Vector{<:Real}`           | amplitude aₕ⁽ᴰ⁾(t) [?]
"""
function amplitude(sol::AbstractDTFTSolution, D::Int=0, H::Int=1)
    H == 0 && return ξ(sol,D,H)
    D == 0 && return abs.(2.0 .* ξ(sol,0,H))
    D == 1 && return real.(2.0 .* ξ(sol,1,H) .* exp.(-im .* ϕ(sol,0,H)))
    D == 2 && return real.(2.0 .* ξ(sol,2,H) .* exp.(-im .* ϕ(sol,0,H))) .+
                        a(sol,0,H) .* ϕ(sol,1,H).^2
    return nothing
end

# phase
"""
    TFT.ϕ(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Shorthand function to obtain the phase of the Dth-degree derivative
of the Hth-harmonic phasor, dispatching to `phase(sol, D, H)`.
"""
ϕ(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) = phase(sol, D, H)

"""
    TFT.phase(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Function to obtain the alternative angle of the Dth-degree derivative of the
Hth-harmonic phasor.

    ∀ h ∈ {0}:
        ϕ₀⁽ᴰ⁾(t) = 0.0
    ∀ h ∉ {0}:
        ϕₕ⁽⁰⁾(t) = ∠[pₕ⁽⁰⁾(t)] ∈ [-π,π]
        ϕₕ⁽¹⁾(t) = ℑ[2 ⋅ ξₕ⁽¹⁾(t) ⋅ exp(-im ⋅ ϕₕ⁽⁰⁾(t))] / aₕ⁽⁰⁾(t) ∈ 𝐑
        ϕₕ⁽²⁾(t) = {ℑ[2 ⋅ ξₕ⁽²⁾(t) ⋅ exp(-im ⋅ ϕₕ⁽⁰⁾(t))] - 2 ⋅ aₕ⁽¹⁾(t) ⋅ ϕₕ⁽¹⁾(t)} / aₕ⁽⁰⁾(t) ∈ 𝐑

See: [Assessing Synchrophasor Estimates of an Event Captured by a Phasor
Measurement Unit, pg. 3112](https://ieeexplore.ieee.org/document/9239915)

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `D::Int`                      | degree of the derivative [-], default=0
- `H::Int`                      | harmonic number [-], default=1

Output:
- `ϕ::Vector{<:Real}`           | phase ϕₕ⁽ᴰ⁾(t) [(rad)]
"""
function phase(sol::AbstractDTFTSolution, D::Int=0, H::Int=1)
    H == 0 && return zero(ξ(sol,D,H))
    D == 0 && return Base.angle.(ξ(sol,0,H))
    D == 1 && return imag.(2.0 .* ξ(sol,1,H) .* exp.(-im .* ϕ(sol,0,H))) ./
                        a(sol,0,H)
    D == 2 && return (imag.(2.0 .* ξ(sol,2,H) .* exp.(-im .* ϕ(sol,0,H))) .-
                        2.0 .* a(sol,1,H) .* ϕ(sol,1,H)) ./ a(sol,0,H)
    return nothing
end

# anti-rotating phase
"""
    TFT.φ(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Shorthand function to obtain the anti-rotating phase of the Dth-degree
derivative of the Hth-harmonic phasor, dispatching to `angle(sol,D,H)`.
"""
φ(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) = ar_phase(sol,D,H)

"""
    TFT.ar_phase(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Function to obtain the anti-rotating phase of the Dth-degree derivative of the
Hth-harmonic phasor.

    ∀ h ∈ {0}:
        φ₀⁽ᴰ⁾(t) = 0.0
    ∀ h ∉ {0}:
        φₕ⁽⁰⁾(t) = ∠[ψₕ⁽⁰⁾(t)] ∈ [-π,π]
        φₕ⁽¹⁾(t) = ℑ[2 ⋅ ψₕ⁽¹⁾(t) ⋅ exp(-im ⋅ φₕ⁽⁰⁾(t))] / aₕ⁽⁰⁾(t) ∈ 𝐑
        φₕ⁽²⁾(t) = {ℑ[2 ⋅ ψₕ⁽²⁾(t) ⋅ exp(-im ⋅ φₕ⁽⁰⁾(t))] - 2 ⋅ aₕ⁽¹⁾(t) ⋅ φₕ⁽¹⁾(t)} / aₕ⁽⁰⁾(t) ∈ 𝐑

See: [Assessing Synchrophasor Estimates of an Event Captured by a Phasor
Measurement Unit, pg. 3112](https://ieeexplore.ieee.org/document/9239915)

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `D::Int`                      | degree of the derivative [-], default=0
- `H::Int`                      | harmonic number [-], default=1

Output:
- `φ::Vector{<:Real}`           | anti-rotating phase φₕ⁽ᴰ⁾(t) [(rad)]
"""
function ar_phase(sol::AbstractDTFTSolution, D::Int=0, H::Int=1)
    H == 0 && return zero(ψ(sol,D,H))
    D == 0 && return Base.angle.(ψ(sol,0,H))
    D == 1 && return imag.(2.0 .* ψ(sol,1,H) .* exp.(-im .* φ(sol,0,H))) ./
                        a(sol,0,H)
    D == 2 && return (imag.(2.0 .* ψ(sol,2,H) .* exp.(-im .* φ(sol,0,H))) .-
                        2.0 .* a(sol,1,H) .* φ(sol,1,H)) ./ a(sol,0,H)
    return nothing
end

# frequency
"""
    TFT.f(sol::TFT.AbstractDTFTSolution, H::Int=1)

Shorthand function to obtain the frequency of the zeroth-degree derivative of
the Hth-harmonic phasor, dispatching to `frequency(sol,H)`.
"""
f(sol::AbstractDTFTSolution, H::Int=1) = frequency(sol,H)

"""
    TFT.frequency(sol::TFT.AbstractDTFTSolution, H::Int=1)

Function to obtain the frequency of the zeroth-degree derivative of the
Hth-harmonic phasor.

    fₕ(t) = Fₕ + ϕₕ⁽¹⁾(t) / (2 π) ∈ 𝐑⁺

See: [Fast Taylor-Fourier Transform for Monitoring Modern Power Grids with
Real-Time Dynamic Harmonic Estimation](tbp)

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `H::Int`                      | harmonic number [-], default=1

Output:
- `f::Vector{<:Real}`           | frequency fₕ(t) [Hz]
"""
frequency(sol::AbstractDTFTSolution, H::Int=1) =
    F(sol,H) .+ (ϕ(sol,1,H) ./ (2 * pi))

# rocof
"""
    TFT.r(sol::TFT.AbstractDTFTSolution, H::Int=1)

Shorthand function to obtain the rate-of-change-of-frequency of the
zeroth-degree derivative of the Hth-harmonic phasor, dispatching to
`rocof(sol,H)`.
"""
r(sol::AbstractDTFTSolution, H::Int=1) = rocof(sol,H)

"""
    TFT.rocof(sol::TFT.AbstractDTFTSolution, H::Int=1)

Function to obtain the rate-of-change-of-frequency of the zeroth-degree
derivative of the Hth-harmonic phasor.

    rₕ(t) = ϕₕ⁽²⁾(t) / (2 π)² ∈ 𝐑

See: [Fast Taylor-Fourier Transform for Monitoring Modern Power Grids with
Real-Time Dynamic Harmonic Estimation](tbp)

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `H::Int`                      | harmonic number [-], default=1

Output:
- `r::Vector{<:Real}`           | rocof rₕ(t) [Hz²]
"""
rocof(sol::AbstractDTFTSolution, H::Int=1) =
    ϕ(sol,2,H) ./ (2 * pi)^2

# dynamic phasor
"""
    TFT.ξ(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Shorthand function to obtain the Dth-degree derivative of the Hth-harmonic
dynamic phasor, dispatching to `phasor(sol,D,H)`.
"""
ξ(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) = phasor(sol,D,H)

"""
    TFT.phasor(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Function to obtain the D-th-degree derivative of the Hth-harmonic dynamic
phasor. For the zeroth-harmonic dynamic phasor, only the real part of the
dynamic phasor is returned.

    ∀ h ∈ {0}:
        ξ₀⁽ᴰ⁾(t) ∈ ℝ
    ∀ h ∉ {0}:
        ξₕ⁽ᴰ⁾(t) ∈ ℂ

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `D::Int`                      | degree of the derivative [-], default=0
- `H::Int`                      | harmonic number [-], default=1

Output:
- `ξ::Vector{<:Complex}`        | dynamic phasor ξₕ⁽ᴰ⁾(t) [?]
"""
function phasor(sol::AbstractDTFTSolution, D::Int=0, H::Int=1)
    check_sol(sol, D, H)

    H == 0 && return real(sol.X[H][:,D+1])
    return sol.X[H][:,D+1]
end

# anti-rotating dynamic phasor
"""
    TFT.ψ(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Shorthand function to obtain the Dth-degree derivative of the Hth-harmonic anti-
rotating dynamic phasor, dispatching to `ar_phasor(sol,D,H)`.
"""
ψ(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) = ar_phasor(sol,D,H)

"""
    TFT.ar_phasor(sol::TFT.AbstractDTFTSolution, D::Int=0, H::Int=1)

Function to obtain the D-th-degree derivative of the Hth-harmonic anti-rotating
dynamic phasor.

    ∀ h ∈ {0}:
        ψ₀⁽ᴰ⁾(t) = ξ₀⁽ᴰ⁾(t) ∈ ℝ
    ∀ h ∉ {0}:
        ψₕ⁽ᴰ⁾(t) = ξₕ⁽ᴰ⁾(t) exp(-im ωₕ t) ∈ 𝐂

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `D::Int`                      | degree of the derivative [-], default=0
- `H::Int`                      | harmonic number [-], default=1

Output:
- `ψ::Vector{<:Complex}`        | anti-rotating dynamic phasor ψₕ⁽ᴰ⁾(t) [?]
"""
ar_phasor(sol::AbstractDTFTSolution, D::Int=0, H::Int=1) =
    ifelse(H == 0, ξ(sol,D,H), ξ(sol,D,H) .* exp.(-im .* ω(sol,H) .* sol.prob.t))

# signal
"""
    TFT.signal(sol::TFT.AbstractDTFTSolution)

Function to obtain the overall signal.

    s(t) = ∑ₕ ξₕ⁽⁰⁾(t) + conj(ξₕ⁽⁰⁾(t)) ∈ 𝐑

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]

Output:
- `s::Vector{<:Real}`           | signal s(t) [?]
"""
signal(sol::AbstractDTFTSolution) = sum(signal(sol,nh) for nh in sol.prob.h)

"""
    TFT.signal(sol::TFT.AbstractDTFTSolution, H::Int)

Function to obtain the Hth-harmonic signal.

    ∀ h ∈ {0}:
        s₀(t) = ξ₀⁽⁰⁾(t) ∈ 𝐑
    ∀ h ∉ {0}:
        sₕ(t) = ℜ[ξₕ⁽⁰⁾(t) + conj(ξₕ⁽⁰⁾(t))] ∈ 𝐑

Note: In its essence, the real operator `ℜ[]`` is unnecessary, however, it is
used to convert the complex number `x + j0` to a real number `x`.

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]
- `H::Int`                      | harmonic number [-]

Output:
- `s::Vector{<:Real}`           | signal sₕ(t) [?]
"""
signal(sol::AbstractDTFTSolution, H::Int) =
    ifelse(H == 0, ξ(sol,0,H), real(ξ(sol,0,H) + conj(ξ(sol,0,H))))

# error
"""
    TFT.error(sol::TFT.AbstractDTFTSolution)

Function to obtain the error between the input and computed signal.

Input:
- `sol::AbstractDTFTSolution`   | DTFT solution struct [-]

Output:
- `e::Vector{<:Real}`           | error [?]
"""
error(sol::TFT.AbstractDTFTSolution) = sol.prob.s .- signal(sol)


text:
# Motif Search based on Position Weight Matrix
# ============================================
#
# This file defines two types of matrices: PFM and PWM.  PFM is a thin wrapper
# of Matrix and stores non-negative frequency values for each symbol and
# position. PWM is a position weighted matrix and stores score values like PFM.
# PFMs are mutable like usual matrices but PWMs cache partial scores to each
# position and hence must be immutable so that updating operations won't cause
# inconsistency.  PFM is created from a raw frequency matrix or a set of
# sequences. PWM is created from a raw score matrix or a PFM.
#
# The motif search algorithm implemented here uses PWM. It tries to find the
# nearest position from which the total score exceeds a specified threshold. The
# algorithm stops searching immediately when it finds that the total score will
# not exceeds the threshold even if the maximum partial score is achieved from
# that position. This uses the partial score cache stored in a PWM object.

"""
Position frequency matrix.
"""
struct PFM{S,T<:Real} <: AbstractMatrix{T}
    # row: symbol, column: position
    data::Matrix{T}

    function PFM{S,T}(m::AbstractMatrix{T}) where {S<:Union{DNA,RNA},T}
        if size(m, 1) != 4
            throw(ArgumentError("PFM must have four rows"))
        end
        return new{S,T}(m)
    end
end

function PFM{S}(m::AbstractMatrix{T}) where {S, T}
    return PFM{S,T}(m)
end
Base.convert(::Type{PFM{S}}, m::AbstractMatrix{T}) where {S, T} = PFM{S}(m)

function Base.convert(::Type{Matrix{T}}, m::PFM) where T
    return convert(Matrix{T}, m.data)
end

function PFM(set)
    return PFM(collect(set))
end

function PFM(set::Vector)
    if isempty(set)
        throw(ArgumentError("sequence set must be non-empty"))
    end
    S = eltype(eltype(set))
    if S ∉ (DNA, RNA)
        throw(ArgumentError("sequence element must be DNA or RNA"))
    end
    len = length(set[1])
    freq = zeros(Int, (4, len))
    for i in 1:lastindex(set)
        seq = set[i]
        if eltype(seq) != S
            throw(ArgumentError("sequence element must be $(S)"))
        elseif length(seq) != len
            throw(ArgumentError("all sequences must be of the same length"))
        end
        for j in 1:len
            s = seq[j]
            if iscertain(s)  # ignore uncertain symbols
                freq[index_nuc(s, j)] += 1
            end
        end
    end
    return PFM{S}(freq)
end

# Broadcasting
struct PFMBroadcastStyle{S} <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{PFM{S,T}}) where {S,T} = PFMBroadcastStyle{S}()
Base.BroadcastStyle(s1::PFMBroadcastStyle, s2::Base.BroadcastStyle) where {S,T} = s1
function Base.similar(bc::Broadcast.Broadcasted{PFMBroadcastStyle{S}}, elt::Type{T}) where {S, T}
    return PFM{S, T}(similar(Array{T}, axes(bc)))
end

function Base.IndexStyle(::Type{<:PFM})
    return IndexLinear()
end

function Base.size(m::PFM)
    return size(m.data)
end

function Base.getindex(m::PFM, i::Integer)
    return m.data[i]
end

function Base.getindex(m::PFM{S}, s::S, j::Integer) where S<:Union{DNA,RNA}
    return m.data[index_nuc(s, j)]
end

function Base.setindex!(m::PFM, val, i::Integer)
    return setindex!(m.data, val, i)
end

function Base.setindex!(m::PFM, val, i::Integer, j::Integer)
    return setindex!(m.data, val, i, j)
end

function Base.setindex!(m::PFM, val, s::S, j::Integer) where S<:Union{DNA,RNA}
    return setindex!(m.data, val, index_nuc(s, j))
end

function Base.show(io::IO, m::PFM{<:Union{DNA,RNA}})
    show_nuc_matrix(io, m)
end

function Base.show(io::IO, ::MIME"text/plain", m::PFM{<:Union{DNA,RNA}})
    show_nuc_matrix(io, m)
end

"""
Position weight matrix.
"""
struct PWM{S,T<:Real} <: AbstractMatrix{T}
    # symbol type
    symtype::Type{S}

    # score matrix (row: symbols, column: position)
    data::Matrix{T}

    # maxscore[i] == sum(maximum(data, 1)[i:end])
    maxscore::Vector{T}

    function PWM{S,T}(pwm::AbstractMatrix) where {S<:Union{DNA,RNA},T}
        if size(pwm, 1) != 4
            throw(ArgumentError("PWM must have four rows"))
        end
        # make a copy for safety
        return new{S,T}(S, copy(pwm), make_maxscore(pwm))
    end
end

"""
    PWM(pfm::PFM{<:Union{DNA,RNA}}; prior=fill(1/4,4))

Create a position weight matrix from a position frequency matrix `pfm`.

The positive weight matrix will be `log2.((pfm ./ sum(pfm, 1)) ./ prior)`.
"""
function PWM(pfm::PFM{S}; prior=fill(1/4, 4)) where S <: Union{DNA,RNA}
    if !all(x -> x > 0, prior)
        throw(ArgumentError("prior must be positive"))
    elseif sum(prior) ≉ 1
        throw(ArgumentError("prior must be sum to 1"))
    end
    prob = pfm ./ sum(pfm, dims=1)
    return PWM{S,Float64}(log2.(prob ./ prior))
end

"""
    PWM{T}(pwm::AbstractMatrix)

Create a PWM from a raw matrix `pwm` (rows are symbols and columns are
positions).

Examples
--------
```julia-repl
julia> pwm = [
           0.065 -0.007 0.042 0.016 0.003  # A
           0.025 -0.007 0.016 0.058 0.082  # C
           0.066  0.085 0.070 0.054 0.010  # G
           0.016  0.025 0.051 0.058 0.038  # T
       ]
4×5 Array{Float64,2}:
 0.065  -0.007  0.042  0.016  0.003
 0.025  -0.007  0.016  0.058  0.082
 0.066   0.085  0.07   0.054  0.01
 0.016   0.025  0.051  0.058  0.038

julia> PWM{DNA}(pwm)
BioSequences.PWM{BioSymbols.DNA,Float64}:
 A  0.065 -0.007  0.042  0.016  0.003
 C  0.025 -0.007  0.016  0.058  0.082
 G  0.066  0.085  0.07   0.054  0.01
 T  0.016  0.025  0.051  0.058  0.038

```
"""
function PWM{T}(pwm::AbstractMatrix) where T <: Union{DNA,RNA}
    return PWM{T,eltype(pwm)}(pwm)
end

function Base.IndexStyle(::Type{<:PWM})
    return IndexLinear()
end

function Base.size(pwm::PWM)
    return size(pwm.data)
end

function Base.getindex(pwm::PWM, i::Integer)
    return getindex(pwm.data, i)
end

function Base.getindex(m::PWM{S}, s::S, j::Integer) where S<:Union{DNA,RNA}
    return m.data[index_nuc(s, j)]
end

function Base.show(io::IO, pwm::PWM{<:Union{DNA,RNA}})
    show_nuc_matrix(io, pwm)
end

function Base.show(io::IO, ::MIME"text/plain", m::PWM{<:Union{DNA,RNA}})
    show_nuc_matrix(io, m)
end

"""
    maxscore(pwm::PWM)

Return the maximum achievable score of `pwm`.
"""
function maxscore(pwm::PWM)
    if size(pwm, 2) == 0
        return zero(eltype(pwm))
    else
        return pwm.maxscore[1]
    end
end

# Make an accumulated maximum score vector.
function make_maxscore(pwm)
    len = size(pwm, 2)
    maxscore = Vector{eltype(pwm)}(undef, len)
    for j in len:-1:1
        s = pwm[1,j]
        for i in 2:size(pwm, 1)
            s = max(s, pwm[i,j])
        end
        if j == len
            maxscore[j] = s
        else
            maxscore[j] = s + maxscore[j+1]
        end
    end
    return maxscore
end

"""
    scoreat(seq::BioSequence, pwm::PWM, start::Integer)

Calculate the PWM score starting from `seq[start]`.
"""
function scoreat(seq::BioSequence, pwm::PWM, start::Integer)
    check_pwm(seq, pwm)
    pwmlen = size(pwm, 2)
    checkbounds(seq, start:start+pwmlen-1)
    score = zero(eltype(pwm))
    @inbounds for j in 0:pwmlen-1
        x = seq[start+j]
        score += iscertain(x) ? pwm[j<<2+trailing_zeros(x)+1] : zero(score)
    end
    return score
end

function Base.findfirst(pwm::PWM, seq::BioSequence, threshold::Real, start = 1, stop=lastindex(seq))
    if eltype(seq) == DNA || eltype(seq) == RNA
        return search_nuc(seq, start:stop, pwm, convert(eltype(pwm), threshold))
    else
        throw(ArgumentError("no search algorithm for '$(typeof(seq))'"))
    end
end

function check_pwm(seq, pwm::PWM{S}) where S <: Union{DNA,RNA}
    if size(pwm, 1) != 4
        throw(ArgumentError("matrix must have four rows"))
    elseif eltype(seq) != S
        throw(ArgumentError("sequence and PWM must have the same element type"))
    end
end

function search_nuc(seq::BioSequence, range::UnitRange{Int}, pwm::PWM{<:Union{DNA,RNA},S}, threshold::S) where S<:Real
    check_pwm(seq, pwm)
    checkbounds(seq, range)
    pwmlen = size(pwm, 2)
    for p in range.start:range.stop-pwmlen+1
        score = zero(eltype(pwm))
        @inbounds for j in 0:pwmlen-1
            if score + pwm.maxscore[j+1] < threshold
                break
            end
            x = seq[p+j]
            score += iscertain(x) ? pwm[j<<2+trailing_zeros(x)+1] : zero(score)
        end
        if score ≥ threshold
            return p
        end
    end
    return nothing
end

function index_nuc(s::Union{DNA,RNA}, j::Integer)
    if !iscertain(s)
        throw(ArgumentError(string("index symbol must be A, C, G or ", typeof(s) == DNA ? "T" : "U")))
    end
    return (j-1) << 2 + (trailing_zeros(s)) + 1
end

function show_nuc_matrix(io::IO, m::Union{PFM{S},PWM{S}}) where S<:Union{DNA,RNA}
    compact(x) = string(x ≥ 0 ? " " : "", sprint(show, x, context=:compact=>true))
    cells = hcat(['A', 'C', 'G', S == DNA ? 'T' : 'U'], compact.(m.data))
    width = maximum(length.(cells), dims=1)
    print(io, summary(m), ':')
    for i in 1:size(cells, 1)
        print(io, "\n ", rpad(cells[i,1], width[1]+1))
        for j in 2:size(cells, 2)
            print(io, rpad(cells[i,j], width[j]+1))
        end
    end
end


text:
module Dec9Part2

function get_continuous_numbers(x::Array{Int64,1}, ref::Int64)::Union{Nothing,Array{Int64,1}}
    x = collect(filter(y -> y < ref, x))
    n_x = length(x)
    for i in 2:n_x
        for j in 1:n_x - (i + 1)
            if ref == sum(x[j:j + i])
                return x[j:j + i]
            end
        end
    end
    return nothing
end
function main(file_path::String)::Union{Int64,Nothing}
    numbers = parse.(Int64, readlines(file_path))
    ref = 144381670
    numbers = get_continuous_numbers(numbers, ref)
    return max(numbers...) + min(numbers...)
end
end


text:
@testset "expand_poly" begin
    x = 1:6
    @test_throws MethodError expand_poly(x, obsdim = 3)
    @test_throws MethodError expand_poly(x, 5, ObsDim.Constant{3}())

    @testset "ObsDim.Last()" begin
        X = expand_poly(x, degree = 5)
        @test @inferred(expand_poly(x)) ≈ X
        @test @inferred(expand_poly(x, 5)) ≈ X
        @test @inferred(expand_poly(x, 5, ObsDim.Last())) ≈ X
        @test @inferred(expand_poly(x, 5, ObsDim.Constant{2}())) ≈ X

        @test size(X) == (5, 6)
        for i in 1:6, k in 1:5
            @test X[k, i] === Float64(x[i]^k)
        end
    end

    @testset "ObsDim.First()" begin
        X = expand_poly(x, degree = 5, obsdim = 1)
        @test expand_poly(x, obsdim = 1) ≈ X
        @test @inferred(expand_poly(x, 5, ObsDim.First())) ≈ X

        @test size(X) == (6, 5)
        for k in 1:5, i in 1:6
            @test X[i, k] === Float64(x[i]^k)
        end
    end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/concave_increasing.jl
# Contains definitions of functions used to enforce bounds on intermediate
# terms in the computation of relaxations:
# positive, negative, lower_bnd, upper_bnd, bnd
#############################################################################

function cut_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    zL = z.lo
    zU = z.hi
    cvv, ccv, cv_grad, cc_grad = cut(zL, zU, x.cv, x.cc, x.cv_grad, x.cc_grad)
    return MC{N, T}(cvv, ccv, z, cv_grad, cc_grad, x.cnst)
end

"""
positive(x::MC)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `McCormick.MC_DOMAIN_TOL`. (Sub)gradients are adjusted appropriately.
"""
positive(x) = x
function positive(x::Interval)
     x ∩ Interval(MC_DOMAIN_TOL, Inf)
end
function positive_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end
function positive(x::MC{N,T}) where {N, T<:Union{NS,MV}}
    positive_kernel(x, positive(x.Intv))
end

negative(x) = x
function negative(x::Interval)
     x ∩ Interval(-Inf, -MC_DOMAIN_TOL)
end
function negative_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
negative(x::MC)

Sets the upper interval bound and the concave relaxation of `x` to a value of
at most `-McCormick.MC_DOMAIN_TOL`. (Sub)gradients are adjusted appropriately.
"""
function negative(x::MC{N,T}) where {N, T<:Union{NS,MV}}
    negative_kernel(x, negative(x.Intv))
end

lower_bnd(x, lb) = x
function lower_bnd(x::Interval{Float64}, lb::Float64)
     x ∩ Interval(lb, Inf)
end
function lower_bnd_kernel(x::MC{N,T}, lb::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
lower_bnd(x::MC, lb::Float64)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `lb`. (Sub)gradients are adjusted appropriately.
"""
function lower_bnd(x::MC{N,T}, lb::Float64) where {N, T<:Union{NS,MV}}
    lower_bnd_kernel(x, lb, lower_bnd(x.Intv, lb))
end

upper_bnd(x, lb) = x
function upper_bnd(x::Interval{Float64}, ub::Float64)
     x ∩ Interval(-Inf, ub)
end
function upper_bnd_kernel(x::MC{N,T}, ub::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
upper_bnd(x::MC, ub)

Sets the upper interval bound and the concave relaxation of `x` to a value of
at most `ub`. (Sub)gradients are adjusted appropriately.
"""
function upper_bnd(x::MC{N,T}, lb::Float64) where {N, T<:Union{NS,MV}}
    upper_bnd_kernel(x, lb, upper_bnd(x.Intv, lb))
end

bnd(x, lb, ub) = x
function bnd(x::Interval{Float64}, lb::Float64, ub::Float64)
     x ∩ Interval(lb, ub)
end
function bnd_kernel(x::MC{N,T}, lb::Float64, ub::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
bnd(x::MC, lb, ub)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `lb`. Sets the upper interval bound and the concave relaxation of `x`
to a value of at most `ub`. (Sub)gradients are adjusted appropriately.
"""
function bnd(x::MC{N,T}, lb::Float64, ub::Float64) where {N, T<:Union{NS,MV}}
    bnd_kernel(x, lb, ub, bnd(x.Intv, lb, ub))
end

function d_lower_bnd_grad(g, x::T, y::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    return
end
function d_upper_bnd_grad(g, x::T, y::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    return
end
function d_bnd_grad(g, x::T, y::T, z::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    g[3] = zero(T)
    return
end


text:
include("../src/vm_basis.jl")

function test_vm_basis_length()
   print("vm_basis.length...")

   @test length(victor_miller_basis(12)) == dim_Mk(12)
   @test length(victor_miller_basis(12,10,true)) == dim_Sk(12)
   @test length(victor_miller_basis(32)) == dim_Mk(32)
   @test length(victor_miller_basis(32,10,true)) == dim_Sk(32)

   println("PASS")
end

function test_vm_basis_type()
   print("vm_basis.type...")

   @test typeof(victor_miller_basis(12)[1]) == fmpz_rel_series
   @test typeof(victor_miller_basis(32)[2]) == fmpz_rel_series

   println("PASS")
end

function test_vm_basis_cusp()
   print("vm_basis.cusp...")

   vm1 = victor_miller_basis(36,10,true)
   @test coeff(vm1[1],0) == 0
   @test coeff(vm1[1],1) == 1

   vm2 = victor_miller_basis(36,10,false)
   @test coeff(vm2[1],0) == 1
   @test coeff(vm2[1],1) == 0

   println("PASS")
end

function test_vm_basis_correctness()
   print("vm_basis.correctness...")

   vm1 = victor_miller_basis(12,6,true)
   S = parent(vm1[1])
   q = gen(S)
   @test vm1[1] - (q-24q^2+252q^3-1472q^4+4830q^5+O(q^6)) == 0

   vm2 = victor_miller_basis(24,6,false)
   S = parent(vm2[1])
   q = gen(S)
   @test vm2[1] - (1+52416000q^3+39007332000q^4+6609020221440q^5+O(q^6)) == 0
   @test vm2[2] - (q+195660q^3+12080128q^4+44656110q^5+O(q^6)) == 0
   @test vm2[3] - (q^2-48q^3+1080q^4-15040q^5+O(q^6)) == 0

   vm3 = victor_miller_basis(32,6,false)
   S = parent(vm3[1])
   q = gen(S)
   @test vm3[1] - (1+2611200q^3+19524758400q^4+19715347537920q^5+O(q^6)) == 0
   @test vm3[2] - (q+50220q^3+87866368q^4+18647219790q^5+O(q^6)) == 0
   @test vm3[3] - (q^2+432q^3+39960q^4-1418560q^5+O(q^6)) == 0

   println("PASS")
end

function test_vm_basis()
   test_vm_basis_length()
   test_vm_basis_type()
   test_vm_basis_cusp()
   test_vm_basis_correctness()

   println("")
end


text:
using Pkg
Pkg.activate("wo3", shared=true)
using MLJBase
using MLJTuning
include("models.jl")
using .Models
using Random
using Test

# for debugging (rng's have verbose shows):
Base.show(io::IO, ::AbstractRNG)  = print("<RNG>")
Base.show(io::IO, ::MIME"text/plain", ::AbstractRNG)  = print("<RNG>")

N = 30
x1 = rand(N);
x2 = rand(N);
x3 = rand(N);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.4*rand(N);

m(K) = KNNRegressor(K=K)
r = [m(K) for K in 2:13]

@testset "constructor" begin
    @test_throws ErrorException TunedModel(model=first(r), tuning=Explicit(),
                                           measure=rms)
    @test_throws ErrorException TunedModel(tuning=Explicit(),
                                           range=r, measure=rms)
    @test_throws ErrorException TunedModel(model=42, tuning=Explicit(),
                                           range=r, measure=rms)
    TunedModel(model=first(r), tuning=Explicit(), range=r)
end

# old julia_bug.jl:
x1 = rand(100);
x2 = rand(100);
x3 = rand(100);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

atom = FooBarRegressor()
ensemble = EnsembleModel(atom=atom, n=50)
mach = machine(ensemble, X, y)
r_n = range(ensemble, :n, lower=10, upper=100)

curves = learning_curve!(mach; range=r_n, resolution=7,
                         rngs = 3,
                         rng_name=:rng)

true


text:
include("virtual_memory.jl")
using JSON, Printf

IP = 1
Fun = []
Global = GlobalMem(Persistent([], [], []), Persistent([], [], []))
CurrMem = MemoryObj(Persistent([], [], []), Temporary([], [], [], []))
MemoryStack = [CurrMem]
PointerStack = [IP]

#=
    extract():
        Extracts the values inside the obj.json file and returns the constants
        table, function directory and quadruples as a tuple.
=#
function extract()
    try
        file = open("obj.json", "r")
        input = JSON.parse(read(file, String))
        close(file)
        input["constants"], input["modules"], input["code"]
    catch LoadError
        exit()
    end
end

#=
    store_or_fetch(address, do_store, value):
        Given an address, it looks for the relevant memory object in order to
        perform a call to either fetch or store depending on the value of do_store.
=#
function store_or_fetch(address::Int64, do_store::Bool=false, value::Any=false)
    if address ∈ ranges[1]
        do_store && return store(Global.variables, value, convert(UInt16, address - ranges[1][1] + 1), '1')
        return fetch(Global.variables, convert(UInt16, address - ranges[1][1] + 1), '1')
    elseif address ∈ ranges[2]
        do_store && return store(Global.variables, value, convert(UInt16, address - ranges[2][1] + 1), '2')
        return fetch(Global.variables, convert(UInt16, address - ranges[2][1] + 1), '2')
    elseif address ∈ ranges[3]
        do_store && return store(Global.variables, value, convert(UInt16, address - ranges[3][1] + 1), '3')
        return fetch(Global.variables, convert(UInt16, address - ranges[3][1] + 1), '3')
    elseif address ∈ ranges[4]
        do_store && return store(MemoryStack[end].persistent, value, convert(UInt16, address - ranges[4][1] + 1), '1')
        return fetch(MemoryStack[end].persistent, convert(UInt16, address - ranges[4][1] + 1), '1')
    elseif address ∈ ranges[5]
        do_store && return store(MemoryStack[end].persistent, value, convert(UInt16, address - ranges[5][1] + 1), '2')
        return fetch(MemoryStack[end].persistent, convert(UInt16, address - ranges[5][1] + 1), '2')
    elseif address ∈ ranges[6]
        do_store && return store(MemoryStack[end].persistent, value, convert(UInt16, address - ranges[6][1] + 1), '3')
        return fetch(MemoryStack[end].persistent, convert(UInt16, address - ranges[6][1] + 1), '3')
    elseif address ∈ ranges[7]
        do_store && return store(MemoryStack[end].temporary, value, convert(UInt16, address - ranges[7][1] + 1), '1')
        return fetch(MemoryStack[end].temporary, convert(UInt16, address - ranges[7][1] + 1), '1')
    elseif address ∈ ranges[8]
        do_store && return store(MemoryStack[end].temporary, value, convert(UInt16, address - ranges[8][1] + 1), '2')
        return fetch(MemoryStack[end].temporary, convert(UInt16, address - ranges[8][1] + 1), '2')
    elseif address ∈ ranges[9]
        do_store && return store(MemoryStack[end].temporary, value, convert(UInt16, address - ranges[9][1] + 1), '3')
        return fetch(MemoryStack[end].temporary, convert(UInt16, address - ranges[9][1] + 1), '3')
    elseif address ∈ ranges[10]
        do_store && return store(Global.constants, value, convert(UInt16, address - ranges[10][1] + 1), '1')
        return fetch(Global.constants, convert(UInt16, address - ranges[10][1] + 1), '1')
    elseif address ∈ ranges[11]
        do_store && return store(Global.constants, value, convert(UInt16, address - ranges[11][1] + 1), '2')
        return fetch(Global.constants, convert(UInt16, address - ranges[11][1] + 1), '2')
    elseif address ∈ ranges[12]
        do_store && return store(Global.constants, value, convert(UInt16, address - ranges[12][1] + 1), '3')
        fetch(Global.constants, convert(UInt16, address - ranges[12][1] + 1), '3')
    else
        do_store && return store(MemoryStack[end].temporary, value, convert(UInt16, address - ranges[13][1] + 1), '4')
        return fetch(MemoryStack[end].temporary, convert(UInt16, address - ranges[13][1] + 1), '4')
    end
end

function conversion(address::Int64, value::String)
    if address ∈ ranges[1] || address ∈ ranges[4] || address ∈ ranges[7] || address ∈ ranges[10]
        return parse(Int64, value)
    elseif address ∈ ranges[2] || address ∈ ranges[5] || address ∈ ranges[8] || address ∈ ranges[11]
        return parse(Float64, value)
    elseif address ∈ ranges[9]
        return parse(Bool, value)
    elseif address ∈ ranges[end]
        return parse(UInt16, value)
    else
        return value
    end
end

#=
    operations(data, operation):
        Given a quadruple and the operand, it performs the operation with the
        values inside the quadruple and stores the result.
=#
function operations(data::Vector{Any}, operation::String)
    if data[2] ∈ ranges[end]
        new_add = convert(Int64, store_or_fetch(data[2]))
        return operations([data[1], new_add, data[3], data[4]], operation)
    end
    if data[3] ∈ ranges[end]
        new_add = convert(Int64, store_or_fetch(data[3]))
        return operations([data[1], data[2], new_add, data[4]], operation)
    end
    left = data[2] == nothing ? 0 : store_or_fetch(data[2])
    right = store_or_fetch(data[3])
    if left === false || right === false
        println("Semantic error! Variable was never assigned a value!")
        exit()
    end
    # Arithmetic
    operation == "+" && store_or_fetch(data[4], true, left + right)
    operation == "-" && store_or_fetch(data[4], true, left - right)
    operation == "*" && store_or_fetch(data[4], true, left * right)
    try
        operation == "^" && store_or_fetch(data[4], true, left ^ right)
    catch DomainError
        println("Error! Expression returned imaginary result!")
        exit()
    end
    if operation == "/"
        result = left / right
        if result == Inf
            println("Error! Division by zero is undefined!")
            exit()
        end
        result = store_or_fetch(data[4], true, result)
        return
    end
    if operation == "//"
        result = trunc(Int64, left) ÷ trunc(Int64, right)
        if result == Inf
            println("Error! Division by zero is undefined!")
            exit()
        end
        result = store_or_fetch(data[4], true, result)
        return
    end
    operation == "++" && store_or_fetch(data[4], true, right + 1)
    operation == "--" && store_or_fetch(data[4], true, right - 1)
    # Boolean logic
    operation == ">"   && store_or_fetch(data[4], true, left > right)
    operation == ">="  && store_or_fetch(data[4], true, left >= right)
    operation == "<"   && store_or_fetch(data[4], true, left < right)
    operation == "<="  && store_or_fetch(data[4], true, left <= right)
    operation == "=="  && store_or_fetch(data[4], true, left == right)
    operation == "¬="  && store_or_fetch(data[4], true, left != right)
    operation == "and" && store_or_fetch(data[4], true, left && right)
    operation == "or"  && store_or_fetch(data[4], true, left || right)
    # Assignment
    if operation == "<-"
        if data[4] ∈ ranges[end]
            new_add = convert(Int64, store_or_fetch(data[4]))
            return operations([data[1], data[2], data[3], new_add], operation)
        end
        return store_or_fetch(data[4], true, right)
    end
    # Array sum
    operation == "<+>" && store_or_fetch(data[4], true, left + right)
end

#=
    printquad(msg_address, last, jump):
        Given the address of a message, it fetches its value and prints it to the
        terminal with either a newline or a space depending on the value last.
=#
function printquad(msg_address::Union{Int64, Nothing}, last::Bool, jump::Bool=true)
    if msg_address == nothing
        message = ""
    else
        msg_address ∈ ranges[end] && return printquad(convert( Int64, store_or_fetch(msg_address) ), last, jump)
        message = store_or_fetch(msg_address)
        if message === false
            println("Semantic error! Variable was never assigned a value!")
            exit()
        end
    end
    if msg_address ∈ ranges[3] || msg_address ∈ ranges[6] || msg_address ∈ ranges[end-1]
        message = message[2:end-1]
    end
    last || print(message, ' ')
    last && println(message)
    jump && global PointerStack[end] += 1
end

#=
    go_to_eval(data, type):
        Evaluates the truth value of the data received and decides weather to
        adjust the IP to the valua in data or skip to the next quadruple.
=#
function go_to_eval(data::Vector{Any}, type::Bool)
    result = store_or_fetch(data[2])
    if (result && type) || (!result && !type)
        global PointerStack[end] = data[4] + 1
    else
        global PointerStack[end] += 1
    end
end

#=
    statistics(data, operation, storage):
        Given a string dictating the statistic function to perform, it will
        store the value inside the address inside the storage address.
=#
function statistics(data::Vector{Real}, operation::String, storage::Int64)
    operation == "Size"       && store_or_fetch(storage, true, length(data))
    operation == "Mean"       && store_or_fetch(storage, true, mean(data))
    operation == "Median"     && store_or_fetch(storage, true, median(data))
    operation == "Mode"       && store_or_fetch(storage, true, mode(data))
    operation == "Variance"   && store_or_fetch(storage, true, var(data))
    operation == "Std"        && store_or_fetch(storage, true, std(data))
    operation == "Sum"        && store_or_fetch(storage, true, sum(data))
    operation == "Min"        && store_or_fetch(storage, true, minimum(data))
    operation == "Max"        && store_or_fetch(storage, true, maximum(data))
end

#=
    plot_generator(data, <data>, type, filename):
        Given one or two arrays of data, it will generate a chart and export it
        to the received filename. The type of chart to plot is specified by the
        type attribute.
=#
function plot_generator(data::Vector{Real}, type::String, filename::String)
    color = ["rgb(153, 9, 136)", "rgb(131, 76, 167)"]
    if type == "Histogram"
        trace = histogram(x=data, marker=attr(color=color[1], opacity=0.6, line=attr(color=color[2], width=1.5)))
    elseif type == "Violin"
        trace = violin(x=data, marker=attr(color=color[1], opacity=0.6, line=attr(color=color[2], width=1.5)))
    elseif type == "Box"
        trace = box(x=data, marker=attr(color=color[1], opacity=0.6, line=attr(color=color[2], width=1.5)))
    end
    savefig(plot(trace), filename)
end

function plot_generator(xdata::Vector{Real}, ydata::Vector{Real}, type::String, filename::String)
    color = ["rgb(153, 9, 136)", "rgb(131, 76, 167)"]
    if type == "Bar"
        trace = bar(x=xdata, y=ydata, marker=attr(color=color[1], opacity=0.6, line=attr(color=color[2], width=1.5)))
    elseif type == "Scatter"
        trace = scatter(x=xdata, y=ydata, marker=attr(color=color[1], opacity=0.6, line=attr(color=color[2], width=1.5)), mode="lines+markers")
    end
    savefig(plot(trace), filename)
end

# Main
constants, modules, instructions = extract()
for constant ∈ constants
    constant[2] ∈ ranges[end-3] && store(Global.constants, constant[1], convert(UInt16, constant[2] - ranges[end-3][1]+1), '1')
    constant[2] ∈ ranges[end-2] && store(Global.constants, constant[1], convert(UInt16, constant[2] - ranges[end-2][1]+1), '2')
    constant[2] ∈ ranges[end-1] && store(Global.constants, constant[1], convert(UInt16, constant[2] - ranges[end-1][1]+1), '3')
end
constants = nothing

while true
    current = instructions[PointerStack[end]]
    #println(current, '\n', Global, '\n', CurrMem)
    # End runtime
    if current[1] == "EndProgram" break

    # Basic operations
    elseif current[1] ∈ operators
        operations(current, current[1])
        global PointerStack[end] += 1
        continue

    # I/O
    elseif current[1] == "Print"
        printquad(current[4], false)
        continue
    elseif current[1] == "LPrint"
        printquad(current[4], true)
        continue
    elseif current[1] == "Input"
        # Print the message of the input and read the value from the terminal
        printquad(current[3], false, false)
        address = current[4]
        expected = nothing
        input = readline()
        # If the value doesn't match the storage address, raise an error
        try
            if address ∈ ranges[1] || address ∈ ranges[4] || address ∈ ranges[7] || address ∈ ranges[10]
                expected = "int"
                input = parse(Int64, input)
            elseif address ∈ ranges[2] || address ∈ ranges[5] || address ∈ ranges[8] || address ∈ ranges[11]
                expected = "float"
                input = parse(Float64, input)
            end
        catch
            @printf("Error! Type mismatch on input, expected %s.\n", expected)
            exit()
        end
        # Store the value
        store_or_fetch(address, true, input)
        global PointerStack[end] += 1
        continue
    elseif current[1] == "Mirror"
        filename = store_or_fetch(current[2])[2:end-1]
        contents = nothing
        # Attempt to open the file or raise an error if not found
        try
            file = open(filename, "r")
            contents = split( read(file, String)[1:end-1], ';' )
            close(file)
        catch
            println("Error file ", filename, " doesn't exist!")
            exit()
        end
        try
            # Try to save the values inside the given array or raise error
            if current[3] ∈ ranges[1] || current[3] ∈ ranges[4]
                contents = map(s -> parse(Int64, s), contents)
            else
                contents = map(s -> parse(Float64, s), contents)
            end
        catch
            println("Error file ", filename, " contains incorrect values or does not match the type of the storage value!")
            exit()
        end

        for val ∈ current[3]:current[4]
            store_or_fetch(val, true, contents[val - current[3] + 1])
        end
        contents = nothing
        global PointerStack[end] += 1
        continue

    # Jumps
    elseif current[1] == "Goto"
        global PointerStack[end] = current[4] + 1
        continue
    elseif current[1] == "GotoF"
        go_to_eval(current, false)
        continue
    elseif current[1] == "GotoT"
        go_to_eval(current, true)
        continue

    # Modules
    elseif current[1] == "ARE"
        # Find the function in the function directory
        for fun ∈ modules
            fun[1] != current[2] && continue
            push!(Fun, fun)
        end
        # Occupy the memory stablished in the resources in the current memory
        for address ∈ current[3]
            try
                if store_or_fetch(address) !== false
                    continue
                end
            catch
                store_or_fetch(address, true, nothing)
            end
            store_or_fetch(address, true, nothing)
        end
        for address ∈ current[4]
            store_or_fetch(address, true, nothing)
        end
        global PointerStack[end] += 1
        continue
    elseif current[1] == "Parameter"
        # Remember de paramameters given by the function
        push!(Fun, store_or_fetch(current[3]))
        global PointerStack[end] += 1
        continue
    elseif current[1] == "GoSub"
        # Generate a new MemoryObj and IP and put them in their respective stacks
        params = Fun[1][2]
        global Fun = Fun[2:end]
        NewMem = MemoryObj(Persistent([], [], []), Temporary([], [], [], []))
        push!(MemoryStack, NewMem)
        for i = 1:length(params)
            store_or_fetch(params[i][2], true, Fun[i])
        end
        global Fun = []
        global PointerStack[end] += 1
        push!(PointerStack, current[4] + 1)
        continue
    elseif current[1] == "Return"
        # Update the global value of the function and return to the previous memory
        operations(current, "<-")
        pop!(MemoryStack)
        pop!(PointerStack)
        continue
    elseif current[1] == "EndModule"
        # Return to the previous memory
        pop!(MemoryStack)
        pop!(PointerStack)
        global PointerStack[end] += 1
        continue

    # Arrays
    elseif current[1] == "Verify"
        # Evaluate if the index is within the boundaries of the array size
        val = store_or_fetch(current[2])
        inf = store_or_fetch(current[3])
        sup = store_or_fetch(current[4])
        if !(val ∈ inf:sup)
            @printf("Index out of bounds! Expected a value between %d and %d.\n", inf, sup)
            exit()
        end
        global PointerStack[end] += 1
        continue
    elseif current[1] ∈ stats
        using StatsBase

        values = [store_or_fetch(val) for val ∈ current[2]:current[3]]
        values = convert(Vector{Real}, values)
        # If the statistical function is Range, print the values to the terminal
        if current[1] == "Range"
            (current[2] ∈ ranges[1] || current[2] ∈ ranges[4]) && @printf("Range: %d - %d\n", minimum(values), maximum(values))
            (current[2] ∈ ranges[2] || current[2] ∈ ranges[5]) && @printf("Range: %f - %f\n", minimum(values), maximum(values))
        else
            statistics(values, current[1], current[4])
        end
        global PointerStack[end] += 1
        continue

    # Plots
    elseif current[1] ∈ xplots
        println("Generating ", lowercase(current[1]), " plot...")
        using PlotlyJS

        # Get the values for the array to be plotted
        values = [store_or_fetch(val) for val ∈ current[2]:current[3]]
        values = convert(Vector{Real}, values)
        filename = store_or_fetch(current[4])[2:end-1]

        plot_generator(values, current[1], filename)
        println("Plot created as ", filename, '.')
        global PointerStack[end] += 1
        continue
    elseif current[1] ∈ xyplots
        println("Generating ", lowercase(current[1]), " plot...")
        using PlotlyJS

        # Get the x values for the array to be plotted
        xvalues = [store_or_fetch(val) for val ∈ current[2]:current[3]]
        xvalues = convert(Vector{Real}, xvalues)
        global PointerStack[end] += 1
        current = instructions[PointerStack[end]]

        # Get the y values for the array to be plotted after updating the IP
        yvalues = [store_or_fetch(val) for val ∈ current[2]:current[3]]
        yvalues = convert(Vector{Real}, yvalues)
        filename = store_or_fetch(current[4])[2:end-1]

        plot_generator(xvalues, yvalues, current[1], filename)
        println("Plot created as ", filename, '.')
        global PointerStack[end] += 1
        continue
    end
end

max_stars_repo_name: CharlesGovea/Alice

max_issues_repo_name: CharlesGovea/Alice

max_forks_repo_name: CharlesGovea/Alice


text:
# Returns true when currws in a room
function inRoom(currws)
    if !haskey(getRoomfromWS, currws)
        println("Invalid Game Command occured outside Room")
        errormsg(currws, "room-not-in-room")
        return false
    end
    return true
end

# Returns true when JSON has all the fields
function validJSON(JSON, fields::Array{String})
    for i in fields
        if !haskey(JSON, i)
            println("Invalid JSON recieved")
            println("$(String(logdetails)) is missing field $i")
            return false
        end
    end
    return true
end

validJSON(logdetails, field::String) = validJSON(logdetails, [field])

validName(handle) = !(strip(handle) == "")

function getSpectators(roomname)
    (!haskey(SpecDict, roomname)) && return []
    return collect(SpecDict[roomname])
end

function getUsers(roomname)
    return [ getHandlefromWS[i] for i in RoomDict[roomname] ]
end

function getPlayers(roomname)
    (!haskey(SpecDict, roomname)) && return getUsers(roomname)
    return filter(x -> !(x in SpecDict[roomname]), getUsers(roomname))
end


text:
using JuMP
using Cbc
using ChooseOptimizer
using DiscreteFunctions
using LinearAlgebra
import Base.sqrt

"""
`sqrt(g::DiscreteFunction)` returns a function `f` such that
`f*f==g` or throws an error if no such function exists.
This method uses integer linear programming.
"""
function sqrt(g::DiscreteFunction)::DiscreteFunction
    err_msg = "This function does not have a square root."
    n = length(g)
    B = Matrix(g)

    if !quick_sqrt_test(g)
        error(err_msg)
    end

    # options = Dict()
    # options[:logLevel] = 0

    MOD = Model(get_solver())

    @variable(MOD, a[1:n,1:n], Bin)  # entry in A matrix
    @variable(MOD, w[1:n,1:n,1:n], Bin) # w[i,j,k] is a[i,j]*a[j,k]

    # this ensures w[i,j,k] is 0 if either a[i,j] or a[j,k] are
    for i=1:n
        for j=1:n
            for k=1:n
                @constraint(MOD, a[i,j]>=w[i,j,k])
                @constraint(MOD, a[j,k]>=w[i,j,k])
            end
        end
    end

    # this ensures that the A matrix represents a function
    for i=1:n
        @constraint(MOD, sum(a[i,j] for j=1:n) == 1)
    end

    # this ensures that A^2 == B
    for i=1:n
        for k=1:n
            @constraint(MOD, sum(w[i,j,k] for j=1:n) == B[i,k])
        end
    end

    # Turns out, we don't need to maximize the w's to get this to work!
    # @objective(MOD, Max, sum(w[i,j,k] for i=1:n for j=1:n for k=1:n))

    # try to solve
    optimize!(MOD)
    status = Int(termination_status(MOD))
    if status != 1   # if status isn't 1, the IP is not feasible, no sqrt
        error(err_msg)
    end

    # build the sqrt function from the A matrix
    A = Int.(value.(a))
    f = DiscreteFunction(n)
    for i=1:n
        for j=1:n
            if A[i,j]==1
                f[i] = j
            end
        end
    end
    return f
end

"""
`has_sqrt(f::DiscreteFunction)` checks if there is a function `g`
such that `g*g==f`. Returns `true` if so and `false` otherwise.
"""
function has_sqrt(f::DiscreteFunction)::Bool
    try
        g = sqrt(f)
        return true
    catch
        return false
    end
end

"""
`quick_sqrt_test(f::DiscreteFunction)` is a quick check of a necessary
condition for `f` to have a square root. If this returns `false` then `f`
does *not* have a square root. If it returns `true`, it might.
"""
function quick_sqrt_test(f::DiscreteFunction)::Bool
    evals = reduced_eigvals(f)
    d = real(prod(evals))
    return d > 0.5
end

function reduced_eigvals(A::Array{T,2}) where T
    elist = eigvals(A)
    return filter(x -> abs(x)>1e-7, elist)
end

"""
`reduced_eigvals(f::DiscreteFunction)` returns the nonzero eigenvalues of
`f`'s matrix.
"""
function reduced_eigvals(f::DiscreteFunction)
    return reduced_eigvals(Matrix(f))
end


text:
using CSV, DataFrames

# x ← FGMRES(A=Helmholtz, M=V-Cycle, b, x = 0, maxIter = 1)
function generate_vcycle!(n, m, kappa, omega, gamma, b; v2_iter=10, level=3, restrt=1)
    h = r_type(2.0 ./ (n+m))

    sl_matrix_level3, h_matrix_level3 = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))
    kappa_coarse = down(reshape(kappa, n-1, m-1, 1, 1)|>pu)[:,:,1,1]
    gamma_coarse = down(reshape(gamma, n-1, m-1, 1, 1)|>pu)[:,:,1,1]
    sl_matrix_level2, h_matrix_level2 = get_helmholtz_matrices!(kappa_coarse, omega, gamma_coarse; alpha=r_type(0.5))
    kappa_coarse = down(reshape(kappa_coarse, Int64((n/2)-1),  Int64((m/2)-1), 1, 1)|>pu)[:,:,1,1]
    gamma_coarse = down(reshape(gamma_coarse,  Int64((n/2)-1),  Int64((m/2)-1), 1, 1)|>pu)[:,:,1,1]
    sl_matrix_level1, h_matrix_level1 = get_helmholtz_matrices!(kappa_coarse, omega, gamma_coarse; alpha=r_type(0.5))

    A(v::a_type) = vec(helmholtz_chain!(reshape(real(v), n-1, m-1, 1, 1), h_matrix_level3; h=h) + im*helmholtz_chain!(reshape(imag(v), n-1, m-1, 1, 1), h_matrix_level3; h=h))
    function M(v::a_type)
        v = reshape(v, n-1, m-1)
        x = zeros(c_type,n-1,m-1)|>pu
        x, = v_cycle_helmholtz!(n, m, h, x, v, h_matrix_level1, sl_matrix_level1, h_matrix_level2, sl_matrix_level2, h_matrix_level3, sl_matrix_level3;  v2_iter=v2_iter, level=level)
        return vec(x)
    end

    x0 = zeros(c_type,n-1,m-1,1,1)|>pu
    if restrt == -1
        restrt = rand(1:10)
    end
    x_vcycle, = fgmres_func(A, vec(b), restrt, tol=1e-10, maxIter=1, M=M, x=vec(x0), out=-1, flexible=true)
    x_vcycle_channels = complex_grid_to_channels!(x_vcycle)
    return x_vcycle, x_vcycle_channels
end

# x ← FGMRES(A=Helmholtz, M=Jacobi, b, x = 0, maxIter = 1)
function generate_jacobi!(n, m, kappa, omega, gamma, b; v2_iter=10, level=3, restrt=1)
    h = r_type(2.0 ./ (n+m))
    sl_matrix, h_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))

    A(v::a_type) = vec(helmholtz_chain!(reshape(v, n-1, m-1, 1, 1), h_matrix; h=h))
    function M(v::a_type)
        v = reshape(v, n-1, m-1)
        x = zeros(c_type,n-1,m-1)|>pu
        x = jacobi_helmholtz_method!(n, m, h, x, v, sl_matrix)
        return vec(x)
    end

    x0 = zeros(c_type,n-1,m-1,1,1)|>pu
    if restrt == -1
        restrt = rand(1:10)
    end
    x_vcycle, = fgmres_func(A, vec(b), restrt, tol=1e-10, maxIter=1,
                                                    M=M, x=vec(x0), out=-1, flexible=true)
    x_vcycle_channels = complex_grid_to_channels!(x_vcycle)
    return x_vcycle, x_vcycle_channels
end

# r ← Ax - A(FGMRES(A=Helmholtz, M=V-Cycle, b, x = 0, maxIter = 1))
function generate_r_vcycle!(n, m, kappa, omega, gamma, x_true; v2_iter=10, level=3, restrt=1, jac=false)
    h = r_type(2.0 ./ (n+m))
    _, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))
    b_true = helmholtz_chain!(real(x_true), helmholtz_matrix; h=h) + im*helmholtz_chain!(imag(x_true), helmholtz_matrix; h=h)

    if jac == true
        x_vcycle, _ = generate_jacobi!(n, m, kappa, omega, gamma, b_true; v2_iter=v2_iter, level=level, restrt=restrt)
    else
        x_vcycle, _ = generate_vcycle!(n, m, kappa, omega, gamma, b_true; v2_iter=v2_iter, level=level, restrt=restrt)
    end
    x_vcycle = reshape(x_vcycle, n-1, m-1, 1, 1)
    e_true = x_true .- x_vcycle
    r_vcycle = b_true .- (helmholtz_chain!(real(x_vcycle), helmholtz_matrix; h=h) + im*helmholtz_chain!(imag(x_vcycle), helmholtz_matrix; h=h))

    return r_vcycle, e_true
end

function generate_random_data!(data_set_m, n, m, kappa, omega, gamma; e_vcycle_input=true, v2_iter=10, level=3, data_augmentetion=false,
                                                          kappa_type=1, threshold=50, kappa_input=true, kappa_smooth=false, k_kernel=3, axb=false, jac=false, norm_input=false, gmres_restrt=1)

    h = r_type(2.0 ./ (n+m))
    dataset = Tuple[]
    data_set_m = data_augmentetion == true ? floor(Int32,0.75*data_set_m) : data_set_m
    for i = 1:data_set_m
        # Generate Model
        kappa = generate_kappa!(n,m; type=kappa_type, smooth=kappa_smooth, threshold=threshold, kernel=k_kernel)|>pu

        # Generate Random Sample
        x_true = randn(c_type,n-1,m-1, 1, 1)|>pu

        if axb == true
            # Generate b
            _, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))
            r_vcycle = helmholtz_chain!(x_true, helmholtz_matrix; h=h)
            e_true = x_true
        else
            # Generate r,e
            r_vcycle, e_true = generate_r_vcycle!(n, m, kappa, omega, gamma, x_true;restrt=gmres_restrt, jac=jac)
        end

        r_vcycle = r_type(h^2) .* r_vcycle
        if norm_input == true
            norm_r = r_type(norm(r_vcycle))
            r_vcycle = r_vcycle ./ norm_r
            e_true = e_true ./ norm_r
        end

        # Print scale information
        if mod(i,5000) == 0
            if axb == true
                @info "$(Dates.format(now(), "HH:MM:SS")) - i = $(i) norm b = $(norm(r_vcycle)) norm x = $(norm(e_true))"
            else
                @info "$(Dates.format(now(), "HH:MM:SS")) - i = $(i) norm r = $(norm(r_vcycle)) norm e = $(norm(e_true))"
            end
        end

        r_vcycle_channels = complex_grid_to_channels!(r_vcycle)
        e_true_channels = complex_grid_to_channels!(e_true)

        # Generate e-vcycle
        if e_vcycle_input == true
            e_vcycle, e_vcycle_channels = generate_vcycle!(n, m, kappa, omega, gamma, r_vcycle; v2_iter=v2_iter, level=level)
            input = cat(e_vcycle_channels, r_vcycle_channels, dims=3)
        else
            input = r_vcycle_channels
        end

        input = kappa_input == true ? cat(input, reshape(kappa, n-1, m-1, 1, 1), dims=3) : input
        append!(dataset,[(input, e_true_channels)])

        # Data Augmentetion
        if data_augmentetion == true && mod(i,3) == 0
            (input_2,e_2) = dataset[rand(1:size(dataset,1))]
            r_index = e_vcycle_input == true ? 3 : 1
            r_2 = input_2[:,:,r_index:r_index+1,:]

            scalar = abs(rand(r_type))
            r_t = scalar*r_vcycle_channels+(1-scalar)*r_2
            scale = (scalar*norm(r_vcycle_channels) + (1-scalar)*norm(r_2))/norm(r_t)

            input_t = (scalar*input+(1-scalar)*input_2)*scale
            e_t = (scalar*e_true_channels+(1-scalar)*e_2)*scale
            append!(dataset,[(input_t, e_t)])
        end
    end
    return dataset
end

function get_csv_set!(path, data_set_m, n, m)
    h = r_type(2.0 ./ (n+m))
    df_training = CSV.File(path)|> DataFrame
    dataset = Tuple[]
    for i = 1:data_set_m
        input = cat(reshape(df_training.RR[(i-1)*(n-1)*(m-1)+1:i*(n-1)*(m-1)],n-1,m-1,1,1),
                    reshape(df_training.RI[(i-1)*(n-1)*(m-1)+1:i*(n-1)*(m-1)],n-1,m-1,1,1),
                    reshape(df_training.KAPPA[(i-1)*(n-1)*(m-1)+1:i*(n-1)*(m-1)],n-1,n-1,1,1), dims=3)

        output = cat(reshape(df_training.ER[(i-1)*(n-1)*(m-1)+1:i*(n-1)*(m-1)],n-1,n-1,1,1) ./ h^2,
                    reshape(df_training.EI[(i-1)*(n-1)*(m-1)+1:i*(n-1)*(m-1)],n-1,n-1,1,1) ./ h^2, dims=3)

        append!(dataset,[(input, output)])
    end
    return dataset
end

# not in use:
# function generate_gmres_data!(data_set_m, n, m, kappa, omega, gamma; e_vcycle_input=true, v2_iter=10, level=3, data_augmentetion=false,
#                                                           kappa_type=1, threshold=50, kappa_input=true, kappa_smooth=false, axb=false, norm_input=false)
#     h = r_type(2.0 ./ (n+m))
#     coefficient = norm_input == true ? h^2 : 1.0
#     dataset = Tuple[]
#     data_set_m = floor(Int64,data_set_m / 5)
#     for i = 1:data_set_m
#         # Generate Model
#         kappa = generate_kappa!(n, m; type=kappa_type, smooth=kappa_smooth, threshold=threshold)
#
#         # Generate Random Sample
#         x_true = randn(c_type,n-1,m-1, 1, 1)|>pu
#         s_matrix, h_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))
#
#         if axb == true
#             # Generate b
#             b_true = helmholtz_chain!(x_true, h_matrix; h=h)
#         else
#             # Generate r,e
#             b_true, x_true = generate_r_vcycle!(n, m, kappa, omega, gamma, x_true)
#         end
#
#
#         # append!(dataset,[(cat(complex_grid_to_channels!(coefficient .* b_true, n), reshape(kappa, n-1, n-1, 1, 1), dims=3),
#         #                 complex_grid_to_channels!(x_true, n))])
#
#         # if mod(i,500) == 0
#         #     @info "$(Dates.format(now(), "HH:MM:SS")) - i = $(i) norm b = $(norm(b_true)) $(norm((h^2) .* b_true)) norm x = $(norm(x_true))"
#         # end
#
#         index = 1
#         A(r) = vec(helmholtz_chain!(reshape(r, n-1, m-1, 1, 1), h_matrix; h=h))
#         function M(r)
#             r = reshape(r, n-1, m-1)
#             if mod(index,4) == 0 # == 3 || index == 13 || index == 18
#                 e_true = fgmres_v_cycle_helmholtz!(n, m, h, r, kappa, omega, gamma)
#                 append!(dataset,[(cat(complex_grid_to_channels!(r), reshape(kappa, n-1, m-1, 1, 1), dims=3),
#                                 complex_grid_to_channels!(e_true ./ coefficient))])
#
#                 # Print scale information
#                 if mod(i,500) == 0
#                     @info "$(Dates.format(now(), "HH:MM:SS")) - $(i) $(index) norm r = $(norm(r)) norm e = $(norm(e_true)) $(norm(e_true ./ (h^2)))"
#                 end
#             end
#             index = index + 1
#             x = zeros(gmres_type,n-1,m-1)|>pu
#             x, = v_cycle_helmholtz!(n, m, h, x, r, kappa, omega, gamma; u=1,
#                         v1_iter = 1, v2_iter = 20, alpha=r_type(0.5), log = 0, level = 3)
#             return vec(x)
#         end
#         x = zeros(gmres_type,n-1,m-1)|>pu
#         x,flag,err,iter,resvec = gpu_flexible_gmres(A, vec(b_true), 10, tol=1e-20, maxIter=2,
#                                                         M=M, x=vec(x), out=-1, flexible=true)
#
#     end
#     return dataset
# end
#
# function generate_point_source_data!(data_set_m, n, m, kappa, omega, gamma; v2_iter=10, level=3, save=false, path="results")
#     h = r_type(2.0 ./ (n+m))
#     _, helmholtz_matrix = get_helmholtz_matrices!(kappa, omega, gamma; alpha=r_type(0.5))
#
#     dataset = Tuple[]
#     for i = 1:data_set_m
#         src = [rand(1:n-1),rand(1:m-1)]
#         b = zeros(c_type,n-1,m-1, 1, 1)|>pu
#         b[src[1],src[2],1,1] = 1.0 ./mean(h.^2)
#
#         # V-cycle
#         x = zeros(c_type,n-1,m-1,1,1)|>pu
#         A(v) = vec(helmholtz_chain!(reshape(v, n-1, m-1, 1, 1), helmholtz_matrix; h=h))
#         function M(v)
#             v = reshape(v, n-1, m-1)
#             x = zeros(c_type,n-1,m-1)|>pu
#             x, = v_cycle_helmholtz!(n, m, h, x, v, kappa, omega, gamma; v2_iter = v2_iter, level = level)
#             return vec(x)
#         end
#         x,flag,err,iter,residual = flexible_gmres(A, vec(b), 5, tol=1e-10, maxIter=1,
#                                                         M=M, x=vec(x), out=-1, flexible=true)
#
#         x = reshape(x,n-1,m-1,1,1)
#         x_channels = complex_grid_to_channels!(x)
#         b_channels = complex_grid_to_channels!(b)
#         append!(dataset,[(x_channels, b_channels)])
#     end
#     return dataset
# end


text:
function Calibration(filename::AbstractString)
    lines = readlines(filename)
    filter!(e->!startswith(e, "#") && !isempty(strip(e)), lines)

    first_line = lowercase(first(lines))  # version can be v or V, halleluja

    if occursin("v", first_line)
        det_id, version = map(x->parse(Int,x), split(first_line, 'v'))
        n_modules = parse(Int, lines[4])
        idx = 5
    else
        det_id, n_modules = map(x->parse(Int,x), split(first_line))
        version = 1
        idx = 2
    end

    pos = Dict{Int32,Vector{NeRCA.Position}}()
    dir = Dict{Int32,Vector{NeRCA.Direction}}()
    t0s = Dict{Int32,Vector{Float64}}()
    dus = Dict{Int32,UInt8}()
    floors = Dict{Int32,UInt8}()
    omkeys = Dict{OMKey,Int32}()
    pmts = Dict{Int32,Tuple{Int32, UInt8}}()  # DOM ID and channel ID

    max_z = 0.0
    for mod ∈ 1:n_modules
        elements = split(lines[idx])
        module_id, du, floor = map(x->parse(Int, x), elements[1:3])
        if version == 4
            xₘ, yₘ, zₘ, q0, qx, qy, qz, t0ₘ = map(x->parse(Float64, x), elements[4:end-1])
        end
        n_pmts = parse(Int, elements[end])
        pos[module_id] = Vector{NeRCA.Position}()
        dir[module_id] = Vector{NeRCA.Direction}()
        t0s[module_id] = Vector{Float64}()
        dus[module_id] = du
        floors[module_id] = floor

        for pmt in 1:n_pmts
            l = split(lines[idx+pmt])
            pmt_id = parse(Int,first(l))
            x, y, z, dx, dy, dz = map(x->parse(Float64, x), l[2:7])
            max_z = max(max_z, z)
            t0 = parse(Float64,l[8])
            push!(pos[module_id], Position(x, y, z))
            push!(dir[module_id], Direction(dx, dy, dz))
            push!(t0s[module_id], t0)
            omkeys[OMKey(du, floor)] = module_id
            pmts[pmt_id] = (module_id, pmt-1)
        end
        idx += n_pmts + 1
    end
    n_dus = length(unique(values(dus)))
    Calibration(det_id, pos, dir, t0s, dus, floors, omkeys, pmts, max_z, n_dus)
end

@deprecate read_calibration(filename::AbstractString) Calibration(filename::AbstractString)